Commit: 6b09f38516dcef8ae9c06ffd7d90a8aa67e0a102 - myprojects-hg-reborn (hg)

Commit MetaInfo

• Revision: 6b09f38516dcef8ae9c06ffd7d90a8aa67e0a102 (tree)
• Zeit: 2011-03-16 06:51:29
• Autor: lorenzo
• Commiter: lorenzo

Log Message

This code now iterates the calculations on a set of clusters.

Ändern Zusammenfassung

• add: Python-codes/iterate_cluster_anisotropy.py

Diff

diff -r d1bd54aad820 -r 6b09f38516dc Python-codes/iterate_cluster_anisotropy.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Python-codes/iterate_cluster_anisotropy.py	Tue Mar 15 21:51:29 2011 +0000
@@ -0,0 +1,256 @@
+#! /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)
+
+
+
+def calc_anisotropy(final_cluster):
+    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(3)
+    anisotropy[0]=eigenvalues[0]/eigenvalues[2]
+    anisotropy[1]=eigenvalues[1]/eigenvalues[2]
+    anisotropy[2]=N
+
+
+    #print "anisotropy is, ", anisotropy
+
+    return anisotropy
+
+
+
+
+######################################################################
+
+kf=1.
+df= 1.78   # 1.8
+
+print sys.argv
+
+N_iter=sys.argv[-1]
+gen_number=sys.argv[-2]
+
+
+# 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
+    
+
+
+prefix2="_generation_"
+middle2=sys.argv[-2]
+filename_end="_.dat"
+
+# print "middle2+filename is, ",middle2+filename_end
+
+overall=s.zeros(3)
+
+for i in xrange(int(N_iter)):
+    middle=i+1
+    prefix="aggregate_number_%01d"%(middle)
+    # print "prefix is,", prefix
+    filename=prefix+prefix2+middle2+filename_end
+
+    final_cluster=n.loadtxt(filename)
+    anis=calc_anisotropy(final_cluster)
+
+    print "filename is, ", filename
+    print "anisotropy is, ", anis
+    overall=s.vstack((overall, anis))
+
+
+overall=overall[1:, :]
+
+n.savetxt("anisotropy_data.dat", overall)
+
+print "So far so good"