提交: d1bd54aad8204ba385a9b1e6d913dbe43b55d86d - myprojects-hg-reborn (hg)

Commit MetaInfo

• 修订版: d1bd54aad8204ba385a9b1e6d913dbe43b55d86d (tree)
• 时间: 2011-03-16 06:50:44
• 作者: lorenzo
• Commiter: lorenzo

Log Message

I added a code to calculate the anisotropy coefficients of a cluster.

更改概述

• add: Python-codes/cluster_anysotropy.py

差异

diff -r 6f131a48a3ab -r d1bd54aad820 Python-codes/cluster_anysotropy.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Python-codes/cluster_anysotropy.py	Tue Mar 15 21:50:44 2011 +0000
@@ -0,0 +1,223 @@
+#! /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"