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

Commit MetaInfo

• 修訂: c6b282d993a0796f7229dede583ae71156ab3418 (tree)
• 時間: 2012-06-18 04:29:59
• 作者: lorenzo
• Commiter: lorenzo

Log Message

I added a code which calculated the Rgeo and Rg of a set of aggregates.

Change Summary

• add: Python-codes/Rg_vs_Rgeo.py

差異

diff -r dec7f903d28c -r c6b282d993a0 Python-codes/Rg_vs_Rgeo.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Python-codes/Rg_vs_Rgeo.py	Sun Jun 17 19:29:59 2012 +0000
@@ -0,0 +1,165 @@
+#! /usr/bin/env python
+
+
+
+from mayavi import mlab
+
+import scipy as s
+import numpy as n
+
+
+# print sys.argv
+
+
+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 compatibility
+
+######################################################
+
+kf=1.3
+df= 1.3   # 1.8
+
+
+
+N_clust_ini=1024 #initial number of cluster that I combine
+N_gen=9 #number of generations of clusters from collisions
+
+gen_array=2**(s.arange(N_gen)+1)
+
+print "gen_array is, ", gen_array
+
+R_geo_arr=s.zeros(0)
+R_g_arr=s.zeros(0)
+N_clust_arr=s.zeros(0).astype("int")
+
+for i in xrange(N_gen):
+    N_clust_in_gen=N_clust_ini/gen_array[i]
+    print "N_clust_in_gen is, ", N_clust_in_gen
+
+    for k in xrange(N_clust_in_gen):
+
+        prefix="../../df-13/kf-13/aggregate_number_%01d"%(k+1)
+        prefix2="_generation_%01d"%(i+1)
+        filename="_.dat"
+        filename=prefix+prefix2+filename
+
+        
+
+
+
+        final_cluster=n.loadtxt(filename)
+
+        dist_mat=euclidean_distances(final_cluster,final_cluster)
+
+
+
+
+        R_geo=max(s.ravel(dist_mat))
+
+
+        x=final_cluster[:,0]
+        y=final_cluster[:,1]
+        z=final_cluster[:,2]
+
+
+        R1_agg_sq=s.var(final_cluster[:,0])+\
+                  s.var(final_cluster[:,1])+\
+                  s.var(final_cluster[:,2]) +1.
+
+        R1_agg=s.sqrt(R1_agg_sq)
+
+        # print "R1agg is, ", R1_agg
+
+
+        N_tot=len(x)
+
+        # print "the total number of monomers is, ", N_tot
+
+        R_g_theory=(N_tot/kf)**(1./df)
+
+        # print "the theoretical R_g is, ", R_g_theory
+
+        # print "R_geo is, ", R_geo
+
+        R_g_arr=s.hstack((R_g_arr,R1_agg))
+
+        R_geo_arr=s.hstack((R_geo_arr,R_geo))
+        N_clust_arr=s.hstack((N_clust_arr,N_tot))
+
+n.savetxt("R_geo_all.dat", R_geo_arr)
+n.savetxt("R_g_all.dat", R_g_arr)
+n.savetxt("N_clust_all.dat", N_clust_arr)
+
+
+
+print "So far so good"