提交: 3a9fb5caa578dbd583d5c63ef9c2ec4f347dd921 - myprojects-hg-reborn (hg)

Commit MetaInfo

• 修訂: 3a9fb5caa578dbd583d5c63ef9c2ec4f347dd921 (tree)
• 時間: 2011-03-04 09:23:51
• 作者: lorenzo
• Commiter: lorenzo

Log Message

This code combines two clusters into a new one.

Change Summary

• add: Python-codes/compose_two_clusters.py

差異

diff -r 84d1c044288a -r 3a9fb5caa578 Python-codes/compose_two_clusters.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Python-codes/compose_two_clusters.py	Fri Mar 04 00:23:51 2011 +0000
@@ -0,0 +1,302 @@
+#! /usr/bin/env python
+
+from enthought.mayavi import mlab
+
+import scipy as s
+import numpy as n
+import scipy.spatial as sp
+
+import scipy.linalg as sl
+
+
+import sys
+
+
+def random_rot():
+    theta=s.arccos(1.-2.*s.random.uniform(0.,1.,1)[0])-s.pi/2.
+    phi=s.random.uniform(-s.pi,s.pi,1)[0]
+    psi=s.random.uniform(-s.pi,s.pi,1)[0]
+
+
+    oneone=s.cos(theta)*s.cos(psi)
+    onetwo=-s.cos(phi)*s.sin(psi)+s.sin(phi)*s.sin(theta)*s.cos(psi)
+    onethree=s.sin(phi)*s.sin(psi)+s.cos(phi)*s.sin(theta)*s.cos(psi)
+
+    twoone= s.cos(theta)*s.sin(psi)
+    twotwo=s.cos(phi)*s.cos(psi)+s.sin(phi)*s.sin(theta)*s.sin(psi)
+    twothree=-s.sin(phi)*s.cos(psi)+s.cos(phi)*s.sin(theta)*s.sin(psi)
+
+    threeone=-s.sin(theta)
+    threetwo=s.sin(phi)*s.cos(theta)
+    threethree=s.cos(phi)*s.cos(theta)
+
+
+    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 accept_reject_rotation(cluster_1, cluster_2, epsi):
+
+        while(True):
+            random_rot_mat= random_rot()
+            n_row_col=s.shape(cluster_1)
+
+            cluster_rot=s.zeros(s.prod(n_row_col)).reshape((n_row_col[0],\
+                                                            n_row_col[1]))
+
+            for i in s.arange(n_row_col[0]):
+                
+                cluster_rot[i,:]=s.dot(random_rot_mat, cluster_1[i,:])
+
+        
+            dist_list = euclidean_distances(cluster_rot, cluster_2)
+
+        
+
+            # print "dist_list is, ", dist_list
+        
+        
+            # if (not (dist_list < (2.-epsi)).any()) and \
+            if (not (dist_list < 2.).any()) and \
+                   (dist_list<=(2.+epsi)).any():
+                cluster_new=s.vstack((cluster_rot, cluster_2))
+    
+                return cluster_new
+        
+        
+
+
+
+    
+
+
+
+
+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)
+
+euclidian_distances = euclidean_distances # both spelling for backward compat
+
+
+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 dist_gamma_clusters(cluster_1, cluster_2, kf, df):
+    N1=s.shape(cluster_1)[0]*1.
+    N2=s.shape(cluster_2)[0]*1.
+
+    print "N1 and N2 are, ", N1, N2
+    
+    R1sq=s.var(cluster_1[:,0])+s.var(cluster_1[:,1])+\
+          s.var(cluster_1[:,2]) +1.
+    R2sq=s.var(cluster_2[:,0])+s.var(cluster_2[:,1])+\
+          s.var(cluster_2[:,2]) +1.
+    R1=s.sqrt(R1sq)
+    R2=s.sqrt(R2sq)
+
+    print "R1 is, ", R1
+    print "R2 is, ", R2
+
+    gamma_sq=((N1+N2)**2.)/(N1*N2)*((N1+N2)/kf)**(2./df)\
+              -(N1+N2)/N2*(R1**2.)-(N1+N2)/N1*(R2**2.)
+
+    a=((N1+N2)**2.)/(N1*N2)*((N1+N2)/kf)**(2./df)
+    b=(N1+N2)/N2*R1**2.
+    c=(N1+N2)/N1*R2**2.
+
+    print "a,b,c are, ", a,b,c
+
+    print "gamma_sq is, ", gamma_sq
+
+    my_gamma=s.sqrt(gamma_sq)
+    return my_gamma
+
+
+
+kf=1.3
+df= 1.8   # 1.8
+epsi=0.01
+
+
+cluster_1=n.loadtxt("aggregate_number_21_.dat")
+
+cluster_1=relocate_cluster(cluster_1)
+
+cm1=find_CM(cluster_1)
+
+print "cm1 is, ", cm1
+
+print "cluster_1 is, ", cluster_1
+
+cluster_2=n.loadtxt("aggregate_number_87_.dat")
+
+
+cluster_2=relocate_cluster(cluster_2)
+
+
+print "cluster_2 is, ", cluster_2
+
+
+gamma=dist_gamma_clusters(cluster_1, cluster_2, kf, df)
+
+print "gamma is, ", gamma
+
+
+cluster_2[:,0]=cluster_2[:,0]+gamma
+
+
+cm2=find_CM(cluster_2)
+
+print "cm2 is, ", cm2
+
+
+list_dist=euclidean_distances(cluster_1, cluster_2)
+
+print "len(list_dist) is, ", s.prod(s.shape(list_dist))
+
+print "list_dist is, ", list_dist
+
+
+# x=cluster_2[:,0]
+# y=cluster_2[:,1]
+# z=cluster_2[:,2]
+
+
+# mlab.clf()
+# pts = mlab.points3d(x, y, z, scale_mode='none', resolution=20,\
+#                         color=(0,0,1),scale_factor=2.)
+# #mlab.axes(pts)
+
+# mlab.show()
+
+
+my_rot=random_rot()
+
+mat_calc=s.dot( my_rot,s.transpose(my_rot))
+
+my_det=sl.det(my_rot)
+
+print "mat_calc is, ", mat_calc
+
+print "my_det is, ", my_det
+
+cluster_agglomerate=accept_reject_rotation(cluster_1, cluster_2, epsi)
+
+n.savetxt("aggregate_agglomerate.dat", cluster_agglomerate)
+
+
+x=cluster_agglomerate[:,0]
+y=cluster_agglomerate[:,1]
+z=cluster_agglomerate[:,2]
+
+
+mlab.clf()
+pts = mlab.points3d(x, y, z, scale_mode='none', resolution=20,\
+                        color=(0,0,1),scale_factor=2.)
+#mlab.axes(pts)
+
+mlab.show()
+
+
+R1_agg_sq=s.var(cluster_agglomerate[:,0])+s.var(cluster_agglomerate[:,1])+\
+          s.var(cluster_agglomerate[:,2]) +1.
+
+R1_agg=s.sqrt(R1_agg_sq)
+
+print "R1agg is, ", R1_agg
+
+
+
+print "So far so good"