myprojects-hg-reborn

Fork

#! /usr/bin/env python

import scipy as s
import numpy as n
import pylab as p
import rpy as r
#from rpy import r
#import distance_calc as d_calc

# now I try performing a least-square fitting
import scipy.optimize as sopt


my_config=12 #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

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>1.) #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


# Here I define my own error function i.e. the function I want to minimize

def residuals(p, y, x):
   a,b,c = p
   err = abs(y-(a+b*x**c))
   return err

def residuals2(p, y, x):
   b,c = p
   err = abs(y-(b*x**c))
   return err


def residuals_lin(p, y, x):
   a,b = p
   err = abs(y-(a+b*x))
   return err



r_gyr_dist=p.load(cluster_name)

r_gyr_dist=r_gyr_dist[my_choice]
print "n_comp is,", n_comp
print "r_gyr_dist is, ", r_gyr_dist



#p.loglog(r_gyr_dist,n_comp, "bo")
#p.show()

p.loglog(r_gyr_dist,n_comp,"ro")
p.xlabel('R [radius of gyration]')
p.ylabel('N(R)')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Cluster-size vs radius of gyration')
p.grid(True)
p.savefig("N_vs_radius_gyration_unprocessed.pdf")
p.hold(False)
p.clf()


print "the initial guesses are"
p0 = [1., 1.8]
print s.array(p0)



# now I try actually solving the problem
#print "ok up to here"
plsq = sopt.leastsq(residuals2, p0, args=(n_comp, r_gyr_dist))

#print "ok up to here2"

print "the optimized values of the parameters are:"

print plsq[0]

coeff=plsq[0]

print "the fitted value of the fractal dimension (no averaging the R_g) is, ", coeff[1]

my_r=s.linspace(r_gyr_dist.min(),r_gyr_dist.max(), 100)

my_fit=coeff[0]*my_r**coeff[1]

p.loglog(r_gyr_dist,n_comp,"ro",my_r,my_fit,linewidth=2.)
p.xlabel('R [radius of gyration]')
p.ylabel('N(R)')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Cluster-size vs radius of gyration')
p.grid(True)
p.savefig("N_vs_radius_gyration_unprocessed_fitting.pdf")
p.hold(False)
p.clf()


p.plot(s.log10(r_gyr_dist),s.log10(n_comp),"ro")
p.xlabel('R [radius of gyration]')
p.ylabel('N(R)')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Cluster-size vs radius of gyration')
p.grid(True)
p.savefig("N_vs_radius_gyration_unprocessed_loglog_explicit.pdf")
p.hold(False)
p.clf()





print "the initial guesses are"
p0 = [3., 2.8]
print s.array(p0)



# now I try actually solving the problem
#print "ok up to here"
plsq = sopt.leastsq(residuals_lin, p0, args=(s.log10(n_comp), s.log10(r_gyr_dist)))

#print "ok up to here2"

print "the optimized values of the parameters are:"

print plsq[0]

coeff=plsq[0]


print "the fitted value of the fractal dimension (no averaging the R_g, from log10) is, ", coeff[1]


p.plot(s.log10(r_gyr_dist),s.log10(n_comp),"ro")
p.xlabel('R [radius of gyration]')
p.ylabel('N(R)')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Cluster-size vs radius of gyration')
p.grid(True)
p.savefig("N_vs_radius_gyration_unprocessed_loglog_explicit.pdf")
p.hold(False)
p.clf()




my_fit=coeff[0]+coeff[1]*s.log10(my_r)

p.plot(s.log10(r_gyr_dist),s.log10(n_comp),"ro",s.log10(my_r),my_fit,linewidth=2.)
p.xlabel('R [radius of gyration]')
p.ylabel('N(R)')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Cluster-size vs radius of gyration')
p.grid(True)
p.savefig("N_vs_radius_gyration_unprocessed_fitting_loglog_explicit.pdf")
p.hold(False)
p.clf()



#now I try working with the cluster population

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

print "len(time) now is, ", len(time)

p.save("time_red.dat",time)

#I initialize the coefficients for the other fitting

p0 = [1.0, -1.0]
print s.array(p0)



print "the new initial ansatz is, ", p0

plsq2 = sopt.leastsq(residuals2, p0, args=(n_clus, time))

coeff=plsq2[0]

print "plsq2 is, ", plsq2

my_t=s.linspace(time.min(),time.max(), 100)
#print 'my_r is, ', my_r

#my_conc=coeff[0]+coeff[1]*my_t**coeff[2]


my_conc=coeff[0]*my_t**coeff[1]

# print "my_conc is, ", my_conc
# print "my_t is, ", my_t


#p.plot(time,n_clus, "bo",my_t,my_conc,linewidth=2.)
#p.show()


p.loglog(time,n_clus, "bo",my_t,my_conc,linewidth=2.)
p.xlabel('Time')
p.ylabel('Number of clusters')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Number of clusters vs time')
p.grid(True)
p.savefig("number_cluster_vs_time2.pdf")
p.hold(False)
p.clf()

p.plot(time,n_clus, "bo",my_t,my_conc,linewidth=2.)
p.xlabel('Time')
p.ylabel('Number of clusters')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Number of clusters vs time')
p.grid(True)
p.savefig("number_cluster_vs_time-linear.pdf")
p.hold(False)
p.clf()




#print "ok up to here2"

print "the optimized values of the for fitting the decay in the cluster number are:"



#print plsq[0]
#coeff[1]=coeff[1]**2.
#coeff[2]=coeff[2]**2.


print " b is", coeff[0]


#print "b is", coeff[1]

print "the exponent  of the power-law decay is", coeff[1]


#now another calculation

n_clus2=n_clus[0]/n_clus-1.
my_choose=s.where(n_clus2>0.)
n_clus2=n_clus2[my_choose]
time=time[my_choose]

p.loglog(time,n_clus2, "bo")
p.xlabel('Time')
p.ylabel('N(0)/N(t)-1')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Number of clusters vs time')
p.grid(True)
p.savefig("coefficient_coagulation.pdf")
p.hold(False)
p.clf()

#print "n_clus2 is, ", n_clus2






my_ord=s.argsort(n_comp)
# print "my_ord is, ", my_ord
# print "n_comp[my_ord] is, ", n_comp[my_ord]

n_comp_ord=n_comp[my_ord]
r_gyr_dist_ord=r_gyr_dist[my_ord]
#Now I get ready for the averaging the radius of gyration for each repeated cluster size





#print "len(n_comp) is, ", len(n_comp)

# my_min=n.int(n_comp.min())
# my_max=n.int(n_comp.max())
# n_iter=my_max-my_min+1
# print "max is, ", my_max
# print "min is, ", my_min
# print "n_iter is, ", n_iter

n_iter=n.int(s.sum(n_comp))

R_gyr_aver=s.zeros(n_iter)  
n_clu_rep=s.zeros(n_iter) #the length of these two arrays is for the worst case
#in which n_clusters=n_particles (just to be sure they are long enough for what I
# want to calculate.

R_gyr_aver[:]=2.*n_iter  
n_clu_rep[:]=2.*n_iter  #I am initializing these arrays with "absurd" values
#and at the end I will remove them in order to get the mean radius of gyration
# for a given number of particle in a cluster.
#Right now, for several clusters e.g. all with 10 particles there are different radia
# of gyration; I will end up with a single cluster of 10 particles corresponding
# to different mean radius of gyration.



my_count=0
for i in xrange(0,n_iter):
   #print "i+1 is, ", i+1
   group=s.where(n_comp_ord==((i+1)*1.))
   group=group[0]
   if (len(group)>0):
      #print "my_count is, ", my_count
      #print "group is, ", group
      R_gyr_aver[my_count]=s.mean(r_gyr_dist_ord[group])
      n_clu_rep[my_count]=i+1
      #print "group[0] is, ", group[0]
      # print "R_gyr_aver[my_count],n_clu_rep[my_count] are, ",R_gyr_aver[my_count]\
#             ,n_clu_rep[my_count]
      my_count=my_count+1


#print "the length of n_clu_rep is, ", len(n_clu_rep),len(R_gyr_aver)

#print "n_clu_rep is, ", n_clu_rep
#print "R_gyr_aver is, ", R_gyr_aver


#Now I get back what I wanted by removing the bits with the too high value I
#initialized the arrays to.

R_gyr_aver2=R_gyr_aver[s.where(R_gyr_aver!=(2.*n_iter))]
n_clu_rep2=n_clu_rep[s.where(n_clu_rep!=(2.*n_iter))]

#print "n_clu_rep2 is, ", n_clu_rep2
#print "R_gyr_aver2 is, ", R_gyr_aver2

# print "the length of n_clu_rep2 is, ", len(n_clu_rep2),len(R_gyr_aver2)






p.loglog(R_gyr_aver2,n_clu_rep2,"bo")
p.xlabel('mean radius of gyration')
p.ylabel('Number of particles in cluster')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Cluster-size vs average radius of gyration')
p.grid(True)
p.savefig("N_vs_mean_radius_gyration_processed.pdf")
p.hold(False)
p.clf()


print "the initial guesses are"
p0 = [3., 2.8]
print s.array(p0)



# now I try actually solving the problem
#print "ok up to here"
plsq = sopt.leastsq(residuals2, p0, args=(n_clu_rep2, R_gyr_aver2))

#print "ok up to here2"

print "the optimized values of the parameters are:"

print plsq[0]

coeff=plsq[0]

print "the fitted value of the fractal dimension (averaging on R_g) is, ", coeff[1]


my_R_gyr_aver_fit=s.logspace(s.log10(R_gyr_aver2.min()),s.log10(R_gyr_aver2.max()),100)

#print "my_R_gyr_aver_fit is, ", my_R_gyr_aver_fit

my_n_clus_fit=coeff[0]*my_R_gyr_aver_fit**coeff[1]


p.loglog(R_gyr_aver2,n_clu_rep2,"bo",my_R_gyr_aver_fit,my_n_clus_fit,linewidth=2.)
p.xlabel('mean radius of gyration')
p.ylabel('Number of particles in cluster')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Cluster-size vs average radius of gyration')
p.grid(True)
p.savefig("N_vs_mean_radius_gyration_processed-fitting.pdf")
p.hold(False)
p.clf()

# p.loglog(R_gyr_aver2,n_clu_rep2,"bo",my_R_gyr_aver_fit,my_n_clus_fit,linewidth=2.)
# p.xlabel('mean radius of gyration')
# p.ylabel('Number of particles in cluster')
# #p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
# p.title('Cluster-size vs average radius of gyration')
# p.grid(True)
# p.show()
# p.hold(False)
# p.clf()


print "the initial guesses are"
p0 = [3., 2.8]
print s.array(p0)



# now I try actually solving the problem
#print "ok up to here"
plsq = sopt.leastsq(residuals_lin, p0, args=(s.log10(n_clu_rep2), s.log10(R_gyr_aver2)))

#print "ok up to here2"

print "the optimized values of the parameters are:"

print plsq[0]

coeff=plsq[0]


print "the fitted value of the fractal dimension (after averaging the R_g, from log10) is, ", coeff[1]


my_n_clus_fit=coeff[0]+coeff[1]*s.log10(R_gyr_aver2)



p.plot(s.log10(R_gyr_aver2),s.log10(n_clu_rep2),"ro",s.log10(R_gyr_aver2),my_n_clus_fit,linewidth=2.)
p.xlabel('R [radius of gyration]')
p.ylabel('N(R)')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Cluster-size vs radius of gyration')
p.grid(True)
p.savefig("N_vs_mean_radius_gyration_processed_fitting_loglog_explicit.pdf")
p.hold(False)
p.clf()






R_sq=R_gyr_aver2**2.


p.loglog(n_clu_rep2,R_sq,"bo")
p.xlabel('particles in cluster')
p.ylabel('R square')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Mean radius of gyration squared vs cluster size')
p.grid(True)
p.savefig("Mean_R_gyr_squared_vs_cluster_size.pdf")
p.hold(False)
p.clf()

   

print "the initial guesses are"
p0 = [1.0,1.0]
print s.array(p0)



# now I try actually solving the problem
#print "ok up to here"
plsq = sopt.leastsq(residuals2, p0, args=(R_gyr_aver2**2.,n_clu_rep2))

#print "ok up to here2"

print "the optimized values of the parameters are:"

print plsq[0]

coeff=plsq[0]
#coeff[1]=coeff[1]**2.
#coeff[2]=coeff[2]**2.


print "a is", coeff[0]

print "and the exponent of the power law is", coeff[1]


ini=s.log10(n_clu_rep2.min())
print "ini is, ", ini
fin=s.log10(n_clu_rep2.max())
print "fin is, ", fin

my_n_clus_fit=s.logspace(ini,fin, 100, base=10.0)

my_r_sq_fit=coeff[0]*my_n_clus_fit**coeff[1]



# print "my_n_clus_fit is, ", my_n_clus_fit
# print "my_r_sq_fit",my_r_sq_fit
# print "n_clu_rep2 is, ", n_clu_rep2
# print "R_sq is", R_sq


p.plot(n_clu_rep2,R_sq,"bo",my_n_clus_fit,my_r_sq_fit)
p.xlabel('particles in cluster')
p.ylabel('R square')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Cluster-size vs average radius of gyration')
p.grid(True)
p.savefig("N_R_gyr_average-linear.pdf")
p.hold(False)
p.clf()


my_r_sq_fit=my_r_sq_fit[s.where(my_r_sq_fit>0.)]
my_n_clus_fit=my_n_clus_fit[s.where(my_r_sq_fit>0.)]


p.loglog(n_clu_rep2,R_sq,"bo",my_n_clus_fit,my_r_sq_fit)
p.xlabel('particles in cluster')
p.ylabel('R square')
#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
p.title('Cluster-size vs average radius of gyration')
p.grid(True)
p.savefig("Mean_R_gyr_squared_vs_cluster_size-fitting.pdf")
p.hold(False)
p.clf()




print "So far so good"