myprojects-hg-reborn

Fork

#! /usr/bin/env python

# from enthought.mayavi import mlab

import scipy as s
import numpy as n
import scipy.linalg as sl

import sys


def inertia_tensor(cluster):

    x=cluster[:,0]
    y=cluster[:,1]
    z=cluster[:,2]



    oneone=s.sum(y**2.+z**2.)
    onetwo=-s.sum(x*y)
    onethree=-s.sum(x*z)

    twoone= -s.sum(x*y)
    twotwo=s.sum(x**2.+z**2.)
    twothree=-s.sum(y*z)

    threeone=-s.sum(x*z)
    threetwo=-s.sum(y*z)
    threethree=s.sum(x**2.+y**2.)


    my_mat=s.zeros(9).reshape((3,3))

    my_mat[0,0]=oneone
    my_mat[0,1]=onetwo
    my_mat[0,2]=onethree

    my_mat[1,0]=twoone
    my_mat[1,1]=twotwo
    my_mat[1,2]=twothree

    my_mat[2,0]=threeone
    my_mat[2,1]=threetwo
    my_mat[2,2]=threethree

    return my_mat




def find_CM(cluster):
    CM=s.mean(cluster, axis=0)
    return CM


def relocate_cluster(cluster):
    cluster_shift=find_CM(cluster)
    cluster[:,0]=cluster[:,0]-cluster_shift[0]
    cluster[:,1]=cluster[:,1]-cluster_shift[1]
    cluster[:,2]=cluster[:,2]-cluster_shift[2]

    return cluster



def move_cluster(cluster, vector):
    cluster_new=s.copy(cluster)
    cluster_new[:,0]=cluster_new[:,0]+vector[0]
    cluster_new[:,1]=cluster_new[:,1]+vector[1]
    cluster_new[:,2]=cluster_new[:,2]+vector[2]

    return cluster_new



def calc_rg(cluster):
    x=cluster[:,0]
    y=cluster[:,1]
    z=cluster[:,2]

    rg=s.var(x)+s.var(y)+s.var(z)
    rg=s.sqrt(rg)

    return rg


def euclidean_distances(X, Y, Y_norm_squared=None, squared=False):
    """
Considering the rows of X (and Y=X) as vectors, compute the
distance matrix between each pair of vectors.

Parameters
----------
X: array of shape (n_samples_1, n_features)

Y: array of shape (n_samples_2, n_features)

Y_norm_squared: array [n_samples_2], optional
pre-computed (Y**2).sum(axis=1)

squared: boolean, optional
This routine will return squared Euclidean distances instead.

Returns
-------
distances: array of shape (n_samples_1, n_samples_2)

Examples
--------
>>> from scikits.learn.metrics.pairwise import euclidean_distances
>>> X = [[0, 1], [1, 1]]
>>> # distrance between rows of X
>>> euclidean_distances(X, X)
array([[ 0., 1.],
[ 1., 0.]])
>>> # get distance to origin
>>> euclidean_distances(X, [[0, 0]])
array([[ 1. ],
[ 1.41421356]])
"""
    # should not need X_norm_squared because if you could precompute that as
    # well as Y, then you should just pre-compute the output and not even
    # call this function.
    if X is Y:
        X = Y = n.asanyarray(X)
    else:
        X = n.asanyarray(X)
        Y = n.asanyarray(Y)

    if X.shape[1] != Y.shape[1]:
        raise ValueError("Incompatible dimension for X and Y matrices")

    XX = n.sum(X * X, axis=1)[:, n.newaxis]
    if X is Y: # shortcut in the common case euclidean_distances(X, X)
        YY = XX.T
    elif Y_norm_squared is None:
        YY = Y.copy()
        YY **= 2
        YY = n.sum(YY, axis=1)[n.newaxis, :]
    else:
        YY = n.asanyarray(Y_norm_squared)
        if YY.shape != (Y.shape[0],):
            raise ValueError("Incompatible dimension for Y and Y_norm_squared")
        YY = YY[n.newaxis, :]

    # TODO:
    # a faster cython implementation would do the dot product first,
    # and then add XX, add YY, and do the clipping of negative values in
    # a single pass over the output matrix.
    distances = XX + YY # Using broadcasting
    distances -= 2 * n.dot(X, Y.T)
    distances = n.maximum(distances, 0)
    if squared:
        return distances
    else:
        return n.sqrt(distances)

######################################################################

kf=1.
df= 1.78   # 1.8

print sys.argv



if (len(sys.argv)==2):

    prefix="aggregate_number_"
    middle=sys.argv[-1]
    filename="_.dat"
    filename=prefix+middle+filename
elif (len(sys.argv)==3):
    prefix="aggregate_number_"
    middle=sys.argv[-2]
    prefix2="_generation_"
    middle2=sys.argv[-1]
    filename="_.dat"
    filename=prefix+middle+prefix2+middle2+filename
    





final_cluster=n.loadtxt(filename)
N=len(final_cluster[:,0])
print "N is, ", N

#ensure its CM is in (0,0,0)
final_cluster=relocate_cluster(final_cluster)

inertia_mat=inertia_tensor(final_cluster)

print "inertia_mat is, ", inertia_mat

eigenvalues=s.real(sl.eigvals(inertia_mat)/N)

eigenvalues=s.sort(eigenvalues)[::-1]

print "eigenvalues are, ", eigenvalues

r123=s.sqrt(eigenvalues)

print "r123, ", r123

Rg=calc_rg(final_cluster)

print "Rg is, ", Rg
print "Rg**2. is, ", Rg**2.


print "0.5*(sum(r123**2.)) is, ", 0.5*(sum(r123**2.))

anisotropy=s.zeros(2)
anisotropy[0]=eigenvalues[0]/eigenvalues[2]
anisotropy[1]=eigenvalues[1]/eigenvalues[2]


print "anisotropy is, ", anisotropy

print "So far so good"