Revisão | a0ba34107497f072988cf7a7c9cf7cd73550b97d (tree) |
---|---|
Hora | 2008-01-30 07:46:00 |
Autor | iselllo |
Commiter | iselllo |
I added r_g_time_aver.py which calculates the scaling law for N vs R_g
taking also a time average and selecting only the clusters which appear
at least a certain (to be specified by the user) number of times.
@@ -0,0 +1,182 @@ | ||
1 | +#! /usr/bin/env python | |
2 | + | |
3 | +import scipy as s | |
4 | +from scipy import stats #I need this module for the linear fit | |
5 | +import numpy as n | |
6 | +import pylab as p | |
7 | +#import rpy as r | |
8 | + | |
9 | + | |
10 | +ini_conf=50 | |
11 | +fin_conf=1000 #configurations I read initially in read_test.tcl | |
12 | + | |
13 | + | |
14 | + | |
15 | + | |
16 | + | |
17 | + | |
18 | + | |
19 | + | |
20 | +#my_config=19 #here labels the saved configurations but it is not directly the number | |
21 | +# of the file with the saved configuration | |
22 | + | |
23 | +by=50 #I need it to select the elements from time.dat | |
24 | + | |
25 | +ini_conf=ini_conf/by | |
26 | +fin_conf=fin_conf/by | |
27 | + | |
28 | +# time=p.load("time.dat") | |
29 | +# #print "time (initially read) is, ", time | |
30 | +# print "by is, ", by | |
31 | +# print "len(time) is, ", len(time) | |
32 | +# pick_time=s.arange(0,len(time),by) | |
33 | +# print "pick_time is, ", pick_time | |
34 | +# time=time[pick_time] | |
35 | + | |
36 | + | |
37 | + | |
38 | +cluster_long=s.zeros(1) | |
39 | +r_gyr_long=s.zeros(1) | |
40 | + | |
41 | +for my_config in xrange(ini_conf,fin_conf): | |
42 | + | |
43 | + | |
44 | + my_config=my_config*by | |
45 | + | |
46 | + cluster_name="cluster_dist%05d"%my_config | |
47 | + | |
48 | + n_comp=p.load(cluster_name) | |
49 | + | |
50 | + #print "n_comp is,", n_comp | |
51 | + | |
52 | + my_choice=s.where(n_comp>6.) #I do not calculate the radius of gyration of monomers! | |
53 | + n_comp=n_comp[my_choice] # select from dimer on | |
54 | + | |
55 | + #print 'my_choice is, ', my_choice | |
56 | + | |
57 | + cluster_name="R_gyr_dist%05d"%my_config | |
58 | + | |
59 | + cluster_long=s.concatenate((cluster_long,n_comp)) | |
60 | + r_gyr_dist=p.load(cluster_name) | |
61 | + | |
62 | + r_gyr_dist=r_gyr_dist[my_choice] | |
63 | + r_gyr_long=s.concatenate((r_gyr_long,r_gyr_dist)) | |
64 | + | |
65 | + | |
66 | +cluster_long=cluster_long[1:] | |
67 | +r_gyr_long=r_gyr_long[1:] | |
68 | + | |
69 | +lin_fit=stats.linregress(s.log10(cluster_long),s.log10(r_gyr_long)) | |
70 | + | |
71 | +print "the fractal dimension is, ", 1/lin_fit[0] | |
72 | + | |
73 | + | |
74 | +p.loglog(cluster_long,r_gyr_long, "ro") | |
75 | +p.grid(True) | |
76 | +p.savefig("all_together.pdf") | |
77 | +p.clf() | |
78 | + | |
79 | + | |
80 | + | |
81 | +cluster_ord=s.sort(cluster_long) | |
82 | +r_gyr_ord=r_gyr_long[s.argsort(cluster_long)] | |
83 | + | |
84 | +cluster_uni=n.unique1d(cluster_ord) | |
85 | +r_bound=cluster_ord.searchsorted(cluster_uni,side='right') | |
86 | +l_bound=cluster_ord.searchsorted(cluster_uni,side='left') | |
87 | + | |
88 | +print "the length of r_bound is, ", len(r_bound) | |
89 | + | |
90 | +r_gyr_ari_mean=s.zeros(len(cluster_uni)) | |
91 | +r_gyr_ari_mean_log=s.zeros(len(cluster_uni)) | |
92 | + | |
93 | +for i in xrange(0,len(cluster_uni)): | |
94 | + r_gyr_ari_mean[i]=r_gyr_ord[l_bound[i]:r_bound[i]].mean() | |
95 | + r_gyr_ari_mean_log[i]=s.exp(s.sum(s.log(r_gyr_ord[l_bound[i]:r_bound[i]]))\ | |
96 | + /len(r_gyr_ord[l_bound[i]:r_bound[i]])) | |
97 | + | |
98 | + | |
99 | +lin_fit2=stats.linregress(s.log10(cluster_uni),s.log10(r_gyr_ari_mean)) | |
100 | + | |
101 | +my_fit2=lin_fit2[1]+lin_fit2[0]*s.log10(cluster_uni) | |
102 | + | |
103 | +print "the fractal dimension [arithmetic mean] is, ", 1/lin_fit2[0] | |
104 | + | |
105 | + | |
106 | +lin_fit3=stats.linregress(s.log10(cluster_uni),s.log10(r_gyr_ari_mean_log)) | |
107 | + | |
108 | +my_fit3=lin_fit3[1]+lin_fit3[0]*s.log10(cluster_uni) | |
109 | + | |
110 | +print "the fractal dimension [log mean] is, ", 1/lin_fit3[0] | |
111 | + | |
112 | + | |
113 | + | |
114 | +p.plot(s.log10(cluster_uni),s.log10(r_gyr_ari_mean), "bo",\ | |
115 | + s.log10(cluster_uni), my_fit2,linewidth=2.) | |
116 | +p.xlabel('N') | |
117 | +p.ylabel('R(N)') | |
118 | +#p.legend(('from power-law fitting','from linear fitting on log-log')) | |
119 | +p.title('Evolution Fractal Dimension from R_g') | |
120 | +p.grid(True) | |
121 | + #cluster_name="number_cluster_vs_time2%05d"%my_config | |
122 | +cluster_name="Fractal_dimension_from_mean_and_time.pdf" | |
123 | +p.savefig(cluster_name) | |
124 | +p.hold(False) | |
125 | +p.clf() | |
126 | + | |
127 | + | |
128 | + | |
129 | + | |
130 | + | |
131 | + | |
132 | +p.plot(s.log10(cluster_uni),s.log10(r_gyr_ari_mean_log), "bo",\ | |
133 | + s.log10(cluster_uni), my_fit3,linewidth=2.) | |
134 | +p.xlabel('N') | |
135 | +p.ylabel('R(N)') | |
136 | +#p.legend(('from power-law fitting','from linear fitting on log-log')) | |
137 | +p.title('Evolution Fractal Dimension from R_g') | |
138 | +p.grid(True) | |
139 | + #cluster_name="number_cluster_vs_time2%05d"%my_config | |
140 | +cluster_name="Fractal_dimension_from_mean_log_and_time.pdf" | |
141 | +p.savefig(cluster_name) | |
142 | +p.hold(False) | |
143 | +p.clf() | |
144 | + | |
145 | + | |
146 | + | |
147 | + | |
148 | +count_number=s.zeros(len(r_bound)) | |
149 | + | |
150 | +for i in xrange(0,len(r_bound)): | |
151 | + count_number[i]=r_bound[i]-l_bound[i] #Now I counted how many clusters I have | |
152 | + #for each existing size | |
153 | +survival=s.where(count_number>10.) | |
154 | + | |
155 | +cluster_surv=cluster_uni[survival] | |
156 | +r_log_surv=r_gyr_ari_mean_log[survival] | |
157 | + | |
158 | + | |
159 | + | |
160 | +lin_fit4=stats.linregress(s.log10(cluster_surv),s.log10(r_log_surv)) | |
161 | + | |
162 | +my_fit4=lin_fit4[1]+lin_fit4[0]*s.log10(cluster_surv) | |
163 | + | |
164 | +print "the fractal dimension [selective log mean] is, ", 1./lin_fit4[0] | |
165 | + | |
166 | + | |
167 | +p.plot(s.log10(cluster_surv),s.log10(r_log_surv), "bo",\ | |
168 | + s.log10(cluster_surv), my_fit4,linewidth=2.) | |
169 | +p.xlabel('N') | |
170 | +p.ylabel('R(N)') | |
171 | +#p.legend(('from power-law fitting','from linear fitting on log-log')) | |
172 | +p.title('Evolution Fractal Dimension from R_g') | |
173 | +p.grid(True) | |
174 | + #cluster_name="number_cluster_vs_time2%05d"%my_config | |
175 | +cluster_name="Fractal_dimension_from_mean_log_and_time_selection.pdf" | |
176 | +p.savefig(cluster_name) | |
177 | +p.hold(False) | |
178 | +p.clf() | |
179 | + | |
180 | + | |
181 | + | |
182 | +print "So far so good" |