myprojects-hg-reborn

Fork

#! /usr/bin/env python

import scipy as s
from scipy import stats #I need this module for the linear fit
import numpy as n
import pylab as p
#import rpy as r


ini_conf=50  
fin_conf=1000   #configurations I read initially in read_test.tcl








#my_config=19 #here labels the saved configurations but it is not directly the number
# of the file with the saved configuration

by=50 #I need it to select the elements from time.dat

ini_conf=ini_conf/by
fin_conf=fin_conf/by

# time=p.load("time.dat")
# #print "time (initially read) is, ", time
# print "by is, ", by
# print "len(time) is, ", len(time)
# pick_time=s.arange(0,len(time),by)
# print "pick_time is, ", pick_time
# time=time[pick_time]



cluster_long=s.zeros(1)
r_gyr_long=s.zeros(1)

for my_config in xrange(ini_conf,fin_conf):


    my_config=my_config*by

    cluster_name="cluster_dist%05d"%my_config

    n_comp=p.load(cluster_name)

    #print "n_comp is,", n_comp

    my_choice=s.where(n_comp>6.) #I do not calculate the radius of gyration of monomers!
    n_comp=n_comp[my_choice] # select from dimer on

    #print 'my_choice is, ', my_choice

    cluster_name="R_gyr_dist%05d"%my_config

    cluster_long=s.concatenate((cluster_long,n_comp))
    r_gyr_dist=p.load(cluster_name)

    r_gyr_dist=r_gyr_dist[my_choice]
    r_gyr_long=s.concatenate((r_gyr_long,r_gyr_dist))
    

cluster_long=cluster_long[1:]
r_gyr_long=r_gyr_long[1:]

lin_fit=stats.linregress(s.log10(cluster_long),s.log10(r_gyr_long))

print "the fractal dimension is, ", 1/lin_fit[0]
    

p.loglog(cluster_long,r_gyr_long, "ro")
p.grid(True)
p.savefig("all_together.pdf")
p.clf()



cluster_ord=s.sort(cluster_long)
r_gyr_ord=r_gyr_long[s.argsort(cluster_long)]

cluster_uni=n.unique1d(cluster_ord)
r_bound=cluster_ord.searchsorted(cluster_uni,side='right')
l_bound=cluster_ord.searchsorted(cluster_uni,side='left')

print "the length of r_bound is, ", len(r_bound)

r_gyr_ari_mean=s.zeros(len(cluster_uni))
r_gyr_ari_mean_log=s.zeros(len(cluster_uni))

for i in xrange(0,len(cluster_uni)):
    r_gyr_ari_mean[i]=r_gyr_ord[l_bound[i]:r_bound[i]].mean()
    r_gyr_ari_mean_log[i]=s.exp(s.sum(s.log(r_gyr_ord[l_bound[i]:r_bound[i]]))\
                                    /len(r_gyr_ord[l_bound[i]:r_bound[i]]))


lin_fit2=stats.linregress(s.log10(cluster_uni),s.log10(r_gyr_ari_mean))

my_fit2=lin_fit2[1]+lin_fit2[0]*s.log10(cluster_uni)

print "the fractal dimension [arithmetic mean] is, ", 1/lin_fit2[0]


lin_fit3=stats.linregress(s.log10(cluster_uni),s.log10(r_gyr_ari_mean_log))

my_fit3=lin_fit3[1]+lin_fit3[0]*s.log10(cluster_uni)

print "the fractal dimension [log mean] is, ", 1/lin_fit3[0]



p.plot(s.log10(cluster_uni),s.log10(r_gyr_ari_mean), "bo",\
       s.log10(cluster_uni), my_fit2,linewidth=2.)
p.xlabel('N')
p.ylabel('R(N)')
#p.legend(('from power-law fitting','from linear fitting on log-log'))
p.title('Evolution Fractal Dimension from R_g')
p.grid(True)
    #cluster_name="number_cluster_vs_time2%05d"%my_config
cluster_name="Fractal_dimension_from_mean_and_time.pdf"
p.savefig(cluster_name)
p.hold(False)
p.clf()    






p.plot(s.log10(cluster_uni),s.log10(r_gyr_ari_mean_log), "bo",\
       s.log10(cluster_uni), my_fit3,linewidth=2.)
p.xlabel('N')
p.ylabel('R(N)')
#p.legend(('from power-law fitting','from linear fitting on log-log'))
p.title('Evolution Fractal Dimension from R_g')
p.grid(True)
    #cluster_name="number_cluster_vs_time2%05d"%my_config
cluster_name="Fractal_dimension_from_mean_log_and_time.pdf"
p.savefig(cluster_name)
p.hold(False)
p.clf()    




count_number=s.zeros(len(r_bound))

for i in xrange(0,len(r_bound)):
    count_number[i]=r_bound[i]-l_bound[i]  #Now I counted how many clusters I have
                                          #for each existing size
survival=s.where(count_number>10.)

cluster_surv=cluster_uni[survival]
r_log_surv=r_gyr_ari_mean_log[survival]



lin_fit4=stats.linregress(s.log10(cluster_surv),s.log10(r_log_surv))

my_fit4=lin_fit4[1]+lin_fit4[0]*s.log10(cluster_surv)

print "the fractal dimension [selective log mean] is, ", 1./lin_fit4[0]


p.plot(s.log10(cluster_surv),s.log10(r_log_surv), "bo",\
       s.log10(cluster_surv), my_fit4,linewidth=2.)
p.xlabel('N')
p.ylabel('R(N)')
#p.legend(('from power-law fitting','from linear fitting on log-log'))
p.title('Evolution Fractal Dimension from R_g')
p.grid(True)
    #cluster_name="number_cluster_vs_time2%05d"%my_config
cluster_name="Fractal_dimension_from_mean_log_and_time_selection.pdf"
p.savefig(cluster_name)
p.hold(False)
p.clf()    



print "So far so good"