import
This commit is contained in:
commit
a39295781b
5 changed files with 1535 additions and 0 deletions
201
Minimization/FitTools.py
Executable file
201
Minimization/FitTools.py
Executable file
|
|
@ -0,0 +1,201 @@
|
||||||
|
#!/usr/bin/python
|
||||||
|
|
||||||
|
import sys
|
||||||
|
import os.path
|
||||||
|
import re
|
||||||
|
import time
|
||||||
|
import math
|
||||||
|
|
||||||
|
sys.path.append("/mntdirect/_users/semeraro/python_tools")
|
||||||
|
import FormFactor
|
||||||
|
import PLUV_POPC_RecBuf
|
||||||
|
from PLUV import SDP_base_POPC_RecBuf, SDP_POPC_RecBuf
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
def TimeCount(t):
|
||||||
|
'Time counting'
|
||||||
|
minutes=int(t/60)
|
||||||
|
seconds=int(t-minutes*60)
|
||||||
|
cents=(t-seconds-minutes*60)*100
|
||||||
|
return (minutes,seconds,cents)
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
#def FreePar (PAR,FIX,L_LIM,H_LIM):
|
||||||
|
# 'Redefine all the sub-array related to free parameters'
|
||||||
|
# FREE = PAR.copy() #list(PAR)
|
||||||
|
# shift=0
|
||||||
|
#
|
||||||
|
# for i in range (len(PAR)):
|
||||||
|
# if FIX[i]=="f":
|
||||||
|
# del FREE[i-shift]
|
||||||
|
# del L_LIM[i-shift]
|
||||||
|
# del H_LIM[i-shift]
|
||||||
|
# shift=shift+1
|
||||||
|
#
|
||||||
|
# return ( FREE , L_LIM , H_LIM )
|
||||||
|
|
||||||
|
def FreePar (PAR,FIX,L_LIM_org,H_LIM_org):
|
||||||
|
'Redefine all the sub-array related to free parameters'
|
||||||
|
FREE = []
|
||||||
|
L_LIM = []
|
||||||
|
H_LIM = []
|
||||||
|
|
||||||
|
for i in range (len(PAR)):
|
||||||
|
if FIX[i]!="f":
|
||||||
|
FREE.append(PAR[i])
|
||||||
|
L_LIM.append(L_LIM_org[i])
|
||||||
|
H_LIM.append(H_LIM_org[i])
|
||||||
|
|
||||||
|
return ( FREE , L_LIM , H_LIM )
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
def MergePar (PAR,FREE,FIX):
|
||||||
|
'Define the merged parameter list to evaluate the function'
|
||||||
|
CALC=[]
|
||||||
|
TEMP_PAR=list(PAR)
|
||||||
|
TEMP_FREE=list(FREE)
|
||||||
|
for el in (FIX):
|
||||||
|
if el=="f":
|
||||||
|
CALC.append(TEMP_PAR.pop(0))
|
||||||
|
else:
|
||||||
|
CALC.append(TEMP_FREE.pop(0))
|
||||||
|
TEMP_PAR.pop(0)
|
||||||
|
return CALC
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
def PrintResults(NAME,PAR,FIX,FREE,X2):
|
||||||
|
'Print Results'
|
||||||
|
|
||||||
|
if X2>1e+7 or X2<1e-4:
|
||||||
|
print("Reduced X^2 ...... %.1e\n" %X2 )
|
||||||
|
else:
|
||||||
|
print("Reduced X^2 ...... %.7s\n" %X2 )
|
||||||
|
|
||||||
|
shift=0
|
||||||
|
TEMP_NAME=list(NAME)
|
||||||
|
TEMP_PAR=list(PAR)
|
||||||
|
for i in range (len(TEMP_PAR)):
|
||||||
|
if FIX[i]=="f":
|
||||||
|
del TEMP_NAME[i-shift]
|
||||||
|
del TEMP_PAR[i-shift]
|
||||||
|
shift=shift+1
|
||||||
|
for i in range(len(FREE)):
|
||||||
|
print(" %d.\t%4s ...... %.4e" %(i,TEMP_NAME[i],FREE[i]) )
|
||||||
|
print ()
|
||||||
|
return
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
def PrintResultsFile(title,folder,datafile,function,NAME,PAR_RES,FIX,X2):
|
||||||
|
'Print Results on a File'
|
||||||
|
localtime = time.asctime( time.localtime(time.time()) )
|
||||||
|
|
||||||
|
|
||||||
|
# res_file = open(resfile,"a")
|
||||||
|
# print(resfile)
|
||||||
|
# res_file.write(localtime+"\n")
|
||||||
|
# res_file.write("Title: %s\n" %title)
|
||||||
|
# res_file.write("Data File: %s\n" %datafile)
|
||||||
|
# res_file.write("Model: %s\n" %function)
|
||||||
|
|
||||||
|
resfile="./"+folder+"/Results.dat"
|
||||||
|
with open(resfile, 'a') as fl:
|
||||||
|
fl.writelines(localtime+"\n")
|
||||||
|
fl.writelines("Title: %s\n" %title)
|
||||||
|
fl.writelines("Data File: %s\n" %datafile)
|
||||||
|
fl.writelines("Model: %s\n\n" %function)
|
||||||
|
for i in range(len(NAME)) :
|
||||||
|
fl.writelines( "%d\t%4s\t%.7f\n" %(i,NAME[i],PAR_RES[i]) )
|
||||||
|
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
def Convolution(q0,function,PAR,Dq):
|
||||||
|
" Convolution "
|
||||||
|
|
||||||
|
N= 10
|
||||||
|
step= 6 * Dq / N
|
||||||
|
Qv=[ (q0-3*Dq + q*step) for q in range(N) ]
|
||||||
|
PAR[-1]=0.0
|
||||||
|
|
||||||
|
CI = 0
|
||||||
|
CI = CI + FormFactor.Normal(q0,Qv[0],Dq)*function(Qv[0],PAR)/2.0
|
||||||
|
for i in range(1,len(Qv)-2):
|
||||||
|
CI = CI + FormFactor.Normal(q0,Qv[i],Dq)*function(Qv[i],PAR)
|
||||||
|
CI = CI + FormFactor.Normal(q0,Qv[len(Qv)-1],Dq)*function(Qv[len(Qv)-1],PAR)/2.0
|
||||||
|
return CI*step
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
def Convolution2(q0,function,PAR,Dql,Dqt):
|
||||||
|
" Convolution "
|
||||||
|
lam=0.0955
|
||||||
|
Dq = lambda Dql,Dqt,q: math.sqrt( ( q*Dql )**2 + ( 2*math.pi*Dqt )**2 )/lam
|
||||||
|
N= 10
|
||||||
|
step= 6 * Dq(Dql,Dqt,q0) / N
|
||||||
|
Qv=list()
|
||||||
|
for i in range(0,N):
|
||||||
|
Qv.append( q0-3*Dq(Dql,Dqt,q0) + i*step )
|
||||||
|
TEMP_PAR=list(PAR)
|
||||||
|
TEMP_PAR.pop(-1)
|
||||||
|
TEMP_PAR.append(0.0)
|
||||||
|
CI = 0
|
||||||
|
CI = CI + FormFactor.Normal(q0,Qv[0],Dq(Dql,Dqt,Qv[0]))*function(Qv[0],TEMP_PAR)/2.0
|
||||||
|
for i in range(1,len(Qv)-2):
|
||||||
|
CI = CI + FormFactor.Normal(q0,Qv[i],Dq(Dql,Dqt,Qv[i]))*function(Qv[i],TEMP_PAR)
|
||||||
|
CI = CI + FormFactor.Normal(q0,Qv[len(Qv)-1],Dq(Dql,Dqt,Qv[len(Qv)-1]))*function(Qv[len(Qv)-1],TEMP_PAR)/2.0
|
||||||
|
return CI*step
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
def ChooseFunction (function ):
|
||||||
|
" Choose Function "
|
||||||
|
|
||||||
|
pltOptions = { 'Form': 0, 'strucCtr': 0, 'parDiff': 0, 'porodCtr': 0, 'debyeCtr': 0 }
|
||||||
|
|
||||||
|
########### Sphere
|
||||||
|
|
||||||
|
if function=="Sphere":
|
||||||
|
intensity = FormFactor.Sphere
|
||||||
|
|
||||||
|
elif function=="Sphere_Normal":
|
||||||
|
intensity = FormFactor.Sphere_Normal
|
||||||
|
|
||||||
|
elif function=="SDP_POPC_RecBuf":
|
||||||
|
#intensity = PLUV.SDP_POPC_RecBuf
|
||||||
|
intensity = SDP_POPC_RecBuf
|
||||||
|
|
||||||
|
elif function=="SDP_base_POPC_RecBuf":
|
||||||
|
intensity = SDP_base_POPC_RecBuf
|
||||||
|
|
||||||
|
elif function=="SDP_POPC_RecBuf_conv":
|
||||||
|
intensity = PLUV_POPC_RecBuf.SDP_POPC_RecBuf_conv
|
||||||
|
|
||||||
|
########### Sticky Hard Sphere Potential
|
||||||
|
|
||||||
|
# elif function=="SHS_MSA_a0_Schulz":
|
||||||
|
# pltOptions['strucCtr']=1
|
||||||
|
# pltOptions['parDiff'] = 3
|
||||||
|
# if conv == 0:
|
||||||
|
# intensity = StickyHardSphere.SHS_MSA_a0_Schulz
|
||||||
|
# pltOptions['Form'] = FormFactor.Sphere_Schulz
|
||||||
|
# else:
|
||||||
|
# intensity = StickyHardSphere.SHS_MSA_a0_Schulz_CONV
|
||||||
|
# pltOptions['Form'] = FormFactor.Sphere_Schulz_CONV
|
||||||
|
# pltOptions['parDiff']+= 1
|
||||||
|
|
||||||
|
###########
|
||||||
|
|
||||||
|
else:
|
||||||
|
sys.exit("--- This function name does not exist")
|
||||||
|
|
||||||
|
return ( intensity , pltOptions )
|
||||||
|
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
|
||||||
156
Minimization/Fit_TSA.py
Executable file
156
Minimization/Fit_TSA.py
Executable file
|
|
@ -0,0 +1,156 @@
|
||||||
|
#!/usr/bin/python
|
||||||
|
|
||||||
|
import sys
|
||||||
|
import os.path
|
||||||
|
import time
|
||||||
|
import numpy as np
|
||||||
|
import os
|
||||||
|
|
||||||
|
localtime = time.asctime( time.localtime(time.time()) )
|
||||||
|
print(localtime)
|
||||||
|
|
||||||
|
sys.path.append("/home/semeraro/Programming/PythonStuff/python_tools")
|
||||||
|
from Tools import Reading
|
||||||
|
from Tools import ReadData
|
||||||
|
from Tools import PlotData
|
||||||
|
import FitTools
|
||||||
|
import TSA_algorithm
|
||||||
|
|
||||||
|
|
||||||
|
############################################################
|
||||||
|
############################################################
|
||||||
|
############################################################
|
||||||
|
|
||||||
|
def X2function(Q,IDATA,I,ERR,N):
|
||||||
|
'X^2 calculation'
|
||||||
|
X2 = np.sum( ( ( IDATA-I ) / ERR )**2 )
|
||||||
|
return X2/(Q.shape[0]-(N-1))
|
||||||
|
|
||||||
|
|
||||||
|
############################################################
|
||||||
|
|
||||||
|
t0=time.time()
|
||||||
|
|
||||||
|
############################################################ Read Options & Parameters
|
||||||
|
|
||||||
|
############### Read File & Options
|
||||||
|
print("\n-- Reading File:")
|
||||||
|
read_par = Reading(sys.argv)
|
||||||
|
iters = int(sys.argv[2])
|
||||||
|
print("\n-- Reading Options")
|
||||||
|
options = read_par.Options()
|
||||||
|
|
||||||
|
############### Read Function
|
||||||
|
"-- Reading Function"
|
||||||
|
(function, pltOptions) = FitTools.ChooseFunction( options['function'] )
|
||||||
|
|
||||||
|
############### Read Parameters
|
||||||
|
print("\n-- Reading Parameters\n")
|
||||||
|
(NAME, PAR, FIX, PRIOR, L_LIM, H_LIM) = read_par.Parameters()
|
||||||
|
|
||||||
|
|
||||||
|
############################################################ Read & Select Data
|
||||||
|
(Q, IDATA ,ERR) = ReadData( options['datafile'] ).SelectColumns( options['qmin'],
|
||||||
|
options['qmax'],
|
||||||
|
options['bing'],
|
||||||
|
options['err_mul'] )
|
||||||
|
|
||||||
|
############################################################ Fit Routine
|
||||||
|
|
||||||
|
t0=time.time()
|
||||||
|
|
||||||
|
collection = []
|
||||||
|
collection_X2 = []
|
||||||
|
|
||||||
|
if options['plot']==0:
|
||||||
|
for it in range(0,iters) :
|
||||||
|
print("\n---- Iteration N.", it+1," -----\n")
|
||||||
|
|
||||||
|
# Define the inputs for the function
|
||||||
|
inputs = [Q, IDATA, ERR, PAR, FIX, PRIOR, L_LIM, H_LIM, function, options['temp'], options['folder']]
|
||||||
|
|
||||||
|
(PAR_RES, X2) = TSA_algorithm.SimAnnealing( inputs )
|
||||||
|
print("----- Print Results from Iteration N. ",it+1," -----\n")
|
||||||
|
|
||||||
|
print("X^2\t%.7f\n" %X2 )
|
||||||
|
for i in range(len(NAME)) :
|
||||||
|
print("%d\t%4s\t%.7f" %(i,NAME[i],PAR_RES[i]) )
|
||||||
|
|
||||||
|
for it in range(0,iters) :
|
||||||
|
collection.append(PAR_RES)
|
||||||
|
collection_X2.append(X2)
|
||||||
|
|
||||||
|
#os.rename("./"+options['folder']+"/Trace.dat", "./"+options['folder']+"/Trace_"+str(it)+".dat")
|
||||||
|
|
||||||
|
t1=time.time()-t0
|
||||||
|
(minutes,seconds,cents) = FitTools.TimeCount(t1)
|
||||||
|
print ( "%2d:%2d.%2d" %(minutes,seconds,cents) )
|
||||||
|
(minutes,seconds,cents) = FitTools.TimeCount(t1/iters)
|
||||||
|
print ( "%2d:%2d.%2d per iteration" %(minutes,seconds,cents) )
|
||||||
|
|
||||||
|
else:
|
||||||
|
intensity_init = function(np.array(Q),PAR)
|
||||||
|
I = intensity_init.intensity()
|
||||||
|
############################################################ Save Results & Plots
|
||||||
|
if options['plot']==0:
|
||||||
|
collection = np.array(collection).transpose()
|
||||||
|
collection_X2 = np.array(collection_X2)
|
||||||
|
|
||||||
|
mean = np.mean(collection,axis=1)
|
||||||
|
stde = np.std(collection,axis=1)
|
||||||
|
median = np.median(collection,axis=1)
|
||||||
|
Q1 = np.quantile(collection,0.25,axis=1)
|
||||||
|
Q3 = np.quantile(collection,0.75,axis=1)
|
||||||
|
|
||||||
|
intensity_init = function(np.array(Q),mean)
|
||||||
|
I = intensity_init.intensity()
|
||||||
|
N_free=collection.shape[0]
|
||||||
|
for p in range(collection.shape[0]) :
|
||||||
|
if FIX[p] =="f" : N_free-= 1
|
||||||
|
X2_mean = X2function(Q,IDATA,I,ERR,N_free)
|
||||||
|
|
||||||
|
print("\n----- Print statitics -----\n")
|
||||||
|
|
||||||
|
print("\n# Mean\tst.dev.\tX^2 from mean\t",X2_mean,"\n")
|
||||||
|
for p in range(collection.shape[0]) :
|
||||||
|
if FIX[p] !="f": print(NAME[p],"\t",mean[p],"\t",stde[p])
|
||||||
|
else : print(NAME[p],"\t",mean[p],"\t-")
|
||||||
|
print("\n# Median\tQ1\tQ3\n")
|
||||||
|
for p in range(collection.shape[0]) :
|
||||||
|
if FIX[p] !="f": print(NAME[p],"\t",median[p],"\t",Q1[p],"\t",Q3[p])
|
||||||
|
else : print(NAME[p],"\t",median[p],"\t-")
|
||||||
|
|
||||||
|
np.savetxt("./"+options['folder']+"/Results_collection-X2.dat", collection_X2)
|
||||||
|
np.savetxt("./"+options['folder']+"/Results_collection.dat", collection)
|
||||||
|
|
||||||
|
X2_mean_to_print = np.empty(1,dtype=float)
|
||||||
|
X2_mean_to_print[0] = X2_mean
|
||||||
|
|
||||||
|
np.savetxt("./"+options['folder']+"/Results_mean-X2.dat", X2_mean_to_print)
|
||||||
|
np.savetxt("./"+options['folder']+"/Results_mean-values.dat", np.column_stack((mean,stde)) )
|
||||||
|
np.savetxt("./"+options['folder']+"/Results_IQR-values.dat", np.column_stack((median,Q1,Q3)) )
|
||||||
|
|
||||||
|
print("-- Save Results in ...")
|
||||||
|
if options['plot'] == 0:
|
||||||
|
FitTools.PrintResultsFile( options['title'], options['folder'],
|
||||||
|
options['datafile'], options['function'],
|
||||||
|
NAME, mean, FIX, X2_mean )
|
||||||
|
|
||||||
|
if options['plot'] == 1:
|
||||||
|
PAR_RES = PAR.copy()
|
||||||
|
|
||||||
|
plotObj = PlotData( options['title'],
|
||||||
|
options['folder'],
|
||||||
|
options['qmin'],
|
||||||
|
options['qmax'] )
|
||||||
|
|
||||||
|
if pltOptions['strucCtr'] == 0:
|
||||||
|
if pltOptions['porodCtr'] == 0:
|
||||||
|
plotObj.Form(Q, IDATA, ERR, I, mean, pltOptions)
|
||||||
|
if pltOptions['porodCtr'] == 1:
|
||||||
|
plotObj.Form_Porod(Q, IDATA, ERR, I, mean)
|
||||||
|
else:
|
||||||
|
plotObj.Structure(Q, IDATA, ERR, I, mean, pltOptions)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
208
Minimization/Fit_TSA_mp.py
Executable file
208
Minimization/Fit_TSA_mp.py
Executable file
|
|
@ -0,0 +1,208 @@
|
||||||
|
#!/usr/bin/python
|
||||||
|
|
||||||
|
import sys
|
||||||
|
import os.path
|
||||||
|
import time
|
||||||
|
import numpy as np
|
||||||
|
import os
|
||||||
|
import glob
|
||||||
|
import multiprocessing
|
||||||
|
|
||||||
|
localtime = time.asctime( time.localtime(time.time()) )
|
||||||
|
print(localtime)
|
||||||
|
|
||||||
|
sys.path.append("/home/semeraro/Programming/PythonStuff/python_tools")
|
||||||
|
from Tools import Reading
|
||||||
|
from Tools import ReadData
|
||||||
|
from Tools import PlotData
|
||||||
|
import FitTools
|
||||||
|
import TSA_algorithm
|
||||||
|
|
||||||
|
|
||||||
|
############################################################
|
||||||
|
############################################################
|
||||||
|
############################################################
|
||||||
|
|
||||||
|
def X2function(Q,IDATA,I,ERR,N):
|
||||||
|
'X^2 calculation'
|
||||||
|
X2 = np.sum( ( ( IDATA-I ) / ERR )**2 )
|
||||||
|
return X2/(Q.shape[0]-(N-1))
|
||||||
|
|
||||||
|
|
||||||
|
############################################################
|
||||||
|
|
||||||
|
t0=time.time()
|
||||||
|
|
||||||
|
############################################################ Read Options & Parameters
|
||||||
|
|
||||||
|
############### Read File & Options
|
||||||
|
print("\n-- Reading File:")
|
||||||
|
read_par = Reading(sys.argv)
|
||||||
|
iters = int(sys.argv[2])
|
||||||
|
print("\n-- Reading Options")
|
||||||
|
options = read_par.Options()
|
||||||
|
|
||||||
|
############### Read Function
|
||||||
|
"-- Reading Function"
|
||||||
|
(function, pltOptions) = FitTools.ChooseFunction( options['function'] )
|
||||||
|
|
||||||
|
############### Read Parameters
|
||||||
|
print("\n-- Reading Parameters\n")
|
||||||
|
(NAME, PAR, FIX, PRIOR, L_LIM, H_LIM) = read_par.Parameters()
|
||||||
|
|
||||||
|
|
||||||
|
############################################################ Read & Select Data
|
||||||
|
(Q, IDATA ,ERR) = ReadData( options['datafile'] ).SelectColumns( options['qmin'],
|
||||||
|
options['qmax'],
|
||||||
|
options['bing'],
|
||||||
|
options['err_mul'] )
|
||||||
|
|
||||||
|
############################################################ Fit Routine
|
||||||
|
|
||||||
|
t0=time.time()
|
||||||
|
|
||||||
|
collection = []
|
||||||
|
collection_X2 = []
|
||||||
|
|
||||||
|
if options['plot']==0:
|
||||||
|
|
||||||
|
if __name__ == "__main__":
|
||||||
|
# Define the inputs for the function
|
||||||
|
inputs = [(Q, IDATA, ERR, PAR, FIX, PRIOR, L_LIM, H_LIM, function, options['temp'], options['folder'])]
|
||||||
|
|
||||||
|
# Create a pool of processes with the number of processors
|
||||||
|
num_processors = int(multiprocessing.cpu_count()*0.75)
|
||||||
|
pool = multiprocessing.Pool(processes=num_processors)
|
||||||
|
|
||||||
|
# Map the calculate function to the inputs
|
||||||
|
results = pool.map(TSA_algorithm.SimAnnealing, inputs * iters)
|
||||||
|
|
||||||
|
# Close the pool to free up resources
|
||||||
|
pool.close()
|
||||||
|
pool.join()
|
||||||
|
|
||||||
|
# Print the results
|
||||||
|
for result in results:
|
||||||
|
collection.append(result[0])
|
||||||
|
collection_X2.append(result[1])
|
||||||
|
|
||||||
|
#os.rename("./"+options['folder']+"/Trace.dat", "./"+options['folder']+"/Trace_"+str(it)+".dat")
|
||||||
|
|
||||||
|
t1=time.time()-t0
|
||||||
|
(minutes,seconds,cents) = FitTools.TimeCount(t1)
|
||||||
|
print ( "%2d:%2d.%2d" %(minutes,seconds,cents) )
|
||||||
|
(minutes,seconds,cents) = FitTools.TimeCount(t1/iters)
|
||||||
|
print ( "%2d:%2d.%2d per iteration" %(minutes,seconds,cents) )
|
||||||
|
|
||||||
|
else:
|
||||||
|
intensity_init = function(np.array(Q),PAR)
|
||||||
|
I = intensity_init.intensity()
|
||||||
|
############################################################ Save Results & Plots
|
||||||
|
|
||||||
|
if options['plot']==0:
|
||||||
|
collection = np.array(collection).transpose()
|
||||||
|
collection_X2 = np.array(collection_X2)
|
||||||
|
|
||||||
|
mean = np.mean(collection,axis=1)
|
||||||
|
stde = np.std(collection,axis=1)
|
||||||
|
median = np.median(collection,axis=1)
|
||||||
|
Q1 = np.quantile(collection,0.25,axis=1)
|
||||||
|
Q3 = np.quantile(collection,0.75,axis=1)
|
||||||
|
outliers = []
|
||||||
|
|
||||||
|
intensity_init = function(np.array(Q),mean)
|
||||||
|
I = intensity_init.intensity()
|
||||||
|
N_free=collection.shape[0]
|
||||||
|
for p in range(collection.shape[0]) :
|
||||||
|
if FIX[p] =="f" : N_free-= 1
|
||||||
|
X2_mean = X2function(Q,IDATA,I,ERR,N_free)
|
||||||
|
|
||||||
|
print("\n----- Print statitics -----\n")
|
||||||
|
|
||||||
|
print("\n# Mean\t\tst.dev.\t\tX^2 from mean\t",X2_mean,"\n")
|
||||||
|
for p in range(collection.shape[0]) :
|
||||||
|
if FIX[p] !="f": print(NAME[p],"\t",mean[p],"\t",stde[p])
|
||||||
|
else : print(NAME[p],"\t",mean[p],"\t-")
|
||||||
|
|
||||||
|
print("\n# Median\t\tQ3\t\tQ1\n")
|
||||||
|
for p in range(collection.shape[0]) :
|
||||||
|
if FIX[p] !="f": print(NAME[p],"\t",median[p],"\t",Q3[p],"\t",Q1[p])
|
||||||
|
else : print(NAME[p],"\t",median[p],"\t-")
|
||||||
|
|
||||||
|
for p in range(collection.shape[0]) :
|
||||||
|
out = 0
|
||||||
|
for item in (collection[p]):
|
||||||
|
if item < (Q1[p]-1.5*(Q3[p]-Q1[p])) or item >( Q3[p]+1.5*(Q3[p]-Q1[p])) :
|
||||||
|
out = 1
|
||||||
|
outliers.append(out)
|
||||||
|
|
||||||
|
print("\n# Check\t\tz+/z'\t\tz-/z'\t\toutliers\n")
|
||||||
|
for p in range(collection.shape[0]) :
|
||||||
|
if FIX[p] !="f": print(NAME[p],"\t",(Q3[p]-mean[p])/stde[p]/0.67,"\t",(Q1[p]-mean[p])/stde[p]/0.67,"\t", outliers[p])
|
||||||
|
else : print(NAME[p],"\t-")
|
||||||
|
|
||||||
|
if iters < 10 :
|
||||||
|
bin_his = iters
|
||||||
|
elif iters >= 10 and iters < 60 :
|
||||||
|
bin_his = 10
|
||||||
|
elif iters >= 60 :
|
||||||
|
bin_his = int(iters/5)
|
||||||
|
|
||||||
|
full_pattern = os.path.join("./"+options['folder'], "Histogram_")
|
||||||
|
|
||||||
|
# Check if files exist
|
||||||
|
files_to_delete = glob.glob(full_pattern)
|
||||||
|
if files_to_delete:
|
||||||
|
# Files exist, delete them
|
||||||
|
for file in files_to_delete:
|
||||||
|
os.remove(file)
|
||||||
|
#print("Histrogram_* files updated")
|
||||||
|
|
||||||
|
for p in range(collection.shape[0]) :
|
||||||
|
if FIX[p] !="f":
|
||||||
|
hist = np.empty(bin_his,dtype=float)
|
||||||
|
bin_edges = np.empty(bin_his,dtype=float)
|
||||||
|
hist, bin_edges = np.histogram(collection[p], bins=bin_his, density=False)
|
||||||
|
bin_center = bin_edges[:hist.shape[0]]+(bin_edges[1]-bin_edges[0])/2
|
||||||
|
np.savetxt("./"+options['folder']+"/Histogram_"+NAME[p]+".dat", np.column_stack((bin_center,hist)))
|
||||||
|
|
||||||
|
np.savetxt("./"+options['folder']+"/Results_collection-X2.dat", collection_X2)
|
||||||
|
np.savetxt("./"+options['folder']+"/Results_collection.dat", collection.transpose())
|
||||||
|
|
||||||
|
X2_mean_to_print = np.empty(1,dtype=float)
|
||||||
|
X2_mean_to_print[0] = X2_mean
|
||||||
|
|
||||||
|
np.savetxt("./"+options['folder']+"/Results_mean-X2.dat", X2_mean_to_print)
|
||||||
|
np.savetxt("./"+options['folder']+"/Results_mean-values.dat", np.column_stack((mean,stde)) )
|
||||||
|
np.savetxt("./"+options['folder']+"/Results_IQR-values.dat", np.column_stack((median,Q1,Q3)) )
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
print("\n-- Save Results in ...")
|
||||||
|
if options['plot'] == 0:
|
||||||
|
FitTools.PrintResultsFile( options['title'], options['folder'],
|
||||||
|
options['datafile'], options['function'],
|
||||||
|
NAME, mean, FIX, X2_mean )
|
||||||
|
|
||||||
|
if options['plot'] == 1:
|
||||||
|
PAR_RES = PAR.copy()
|
||||||
|
|
||||||
|
plotObj = PlotData( options['title'],
|
||||||
|
options['folder'],
|
||||||
|
options['qmin'],
|
||||||
|
options['qmax'] )
|
||||||
|
|
||||||
|
if pltOptions['strucCtr'] == 0:
|
||||||
|
if pltOptions['porodCtr'] == 0:
|
||||||
|
plotObj.Form(Q, IDATA, ERR, I, mean, pltOptions)
|
||||||
|
if pltOptions['porodCtr'] == 1:
|
||||||
|
plotObj.Form_Porod(Q, IDATA, ERR, I, mean)
|
||||||
|
else:
|
||||||
|
plotObj.Structure(Q, IDATA, ERR, I, mean, pltOptions)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
388
Minimization/PLUV.py
Normal file
388
Minimization/PLUV.py
Normal file
|
|
@ -0,0 +1,388 @@
|
||||||
|
#!/usr/bin/python
|
||||||
|
|
||||||
|
import sys
|
||||||
|
import os.path
|
||||||
|
import re
|
||||||
|
import time
|
||||||
|
import numpy as np
|
||||||
|
from scipy.stats import norm , gamma
|
||||||
|
from scipy import special
|
||||||
|
from scipy import signal
|
||||||
|
from scipy import ndimage
|
||||||
|
#from numba import njit, prange
|
||||||
|
|
||||||
|
######################################################################################################
|
||||||
|
######################################################################################################
|
||||||
|
######################################################################################################
|
||||||
|
############### POPC/POPG 95/5 mol/mol LUVs ###############
|
||||||
|
############### Reconstitution buffer: TRIS 20 mM, EDTA 2 mM ###############
|
||||||
|
############### Buffer used for PLUV exp. in DESY in 2017: ###############
|
||||||
|
######################################################################################################
|
||||||
|
######################################################################################################
|
||||||
|
######################################################################################################
|
||||||
|
|
||||||
|
|
||||||
|
############### GLOBAL VARIABLES ###############
|
||||||
|
|
||||||
|
# POPG molar ratio
|
||||||
|
x_PG = 0.05
|
||||||
|
|
||||||
|
# ############## POPC and POPG ###############
|
||||||
|
# 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine
|
||||||
|
# 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoglycerol
|
||||||
|
|
||||||
|
# number of chain groups
|
||||||
|
n_CH = 2
|
||||||
|
n_CH2 = 28
|
||||||
|
n_CH3 = 2
|
||||||
|
|
||||||
|
#X-ray scattering length chain groups (nm)
|
||||||
|
b_CH = 1.97256E-05 ;
|
||||||
|
b_CH2 = 2.25435E-05 ;
|
||||||
|
b_CH3 = 2.53615E-05 ;
|
||||||
|
|
||||||
|
### POPC
|
||||||
|
# Lipid-head volume
|
||||||
|
V_HL_PC = 0.331 # 0.320
|
||||||
|
# X-ray scattering length of head groups (nm)
|
||||||
|
b_PC = 2.73340E-04
|
||||||
|
b_CG = 1.88802E-04
|
||||||
|
b_PCN = 1.97256E-04
|
||||||
|
b_Chol = 7.60844E-05
|
||||||
|
# Lipid-volume temperature-dependencies a0 + a1*T (nm^3)
|
||||||
|
a0_V_POPC = 1.22810311835285
|
||||||
|
a1_V_POPC = 0.000934915795086395
|
||||||
|
|
||||||
|
### POPG
|
||||||
|
# Lipid-head volume
|
||||||
|
V_HL_PG = 0.289 ;
|
||||||
|
# X-ray Scattering length of head groups (nm)
|
||||||
|
b_PG = 2.47979E-04
|
||||||
|
b_PG1 = 1.32443E-04
|
||||||
|
b_PG2 = 1.15536E-04
|
||||||
|
|
||||||
|
# Lipid-volume temperature-dependencies a0 + a1*T (nm^3)
|
||||||
|
a0_V_POPG = 1.17881068602663
|
||||||
|
a1_V_POPG = 0.00108364914520327
|
||||||
|
|
||||||
|
|
||||||
|
############### Other variables ###############
|
||||||
|
### Water
|
||||||
|
# V_PW = 24.5e-3 #(nm^3) Perkins 2001 (Hydration shell in proteins)
|
||||||
|
# polynome coefficient for T-dependency of bulk-water-molecule volume (V_HW)
|
||||||
|
# Units in degree Celsius
|
||||||
|
p0_VW = 0.0299218
|
||||||
|
p1_VW = -2.25941e-06
|
||||||
|
p2_VW = 2.5675e-07
|
||||||
|
p3_VW = -1.69661e-09
|
||||||
|
p4_VW = 6.52029e-12
|
||||||
|
# polynome coefficient for T-dependency of bulk-water molar concentration (Cw)
|
||||||
|
#Units in degree Celsius
|
||||||
|
p0_Cw = 55.5052
|
||||||
|
p1_Cw = 0.00131894
|
||||||
|
p2_Cw = -0.000334396
|
||||||
|
p3_Cw = 9.10861e-07
|
||||||
|
|
||||||
|
b_HW = 2.8179E-05
|
||||||
|
d_shl = 0.31 # (nm) Perkins 2001 (Hydration shell in proteins)
|
||||||
|
|
||||||
|
# Composition of the Reconstitution buffer (M)
|
||||||
|
ctris = 0.02
|
||||||
|
cEDTA = 0.002
|
||||||
|
|
||||||
|
### Extra molecules
|
||||||
|
|
||||||
|
# TRIS buffer
|
||||||
|
b_tris = 1.860E-04
|
||||||
|
V_tris = 0.15147 # (nm^3)
|
||||||
|
# EDTA
|
||||||
|
b_EDTA = 4.340E-04 ;
|
||||||
|
V_EDTA = 0.56430 # (nm^3)
|
||||||
|
|
||||||
|
|
||||||
|
##################################################################################################################
|
||||||
|
##################################################################################################################
|
||||||
|
|
||||||
|
#########################################################
|
||||||
|
#@njit(parallel=True)
|
||||||
|
def PDF_normal(x, mu, sig) :
|
||||||
|
return np.exp(-(x-mu)**2 / (2*sig**2) ) /( sig*np.sqrt(2*np.pi))
|
||||||
|
|
||||||
|
#########################################################
|
||||||
|
def lipid_volume(T) :
|
||||||
|
return (1-x_PG) * (a0_V_POPC +T * a1_V_POPC) + x_PG * (a0_V_POPG + T * a1_V_POPG)
|
||||||
|
|
||||||
|
#########################################################
|
||||||
|
def water_volume(T) :
|
||||||
|
return p0_VW + p1_VW*T + p2_VW*T**2 + p3_VW*T**3 + p4_VW*T**4
|
||||||
|
|
||||||
|
#########################################################
|
||||||
|
#@njit(parallel=True)
|
||||||
|
def FTreal_erf(q, mu, d, sig) :
|
||||||
|
return np.where(q==0, 1, np.sin(q*d/2.)/(q*d/2.) * np.exp(-(q*sig)**2/2.) * np.cos(q*mu) )
|
||||||
|
|
||||||
|
#########################################################
|
||||||
|
#@njit(parallel=True)
|
||||||
|
def FTreal_gauss(q, mu, sig) :
|
||||||
|
return np.exp(-(q*sig)**2/2.) * np.cos(q*mu)
|
||||||
|
|
||||||
|
#########################################################
|
||||||
|
def Slab(x, mu, L, sig) :
|
||||||
|
return 0.5 * ( special.erf( (x - (mu-L/2.))/(np.sqrt(2)*sig) ) - special.erf( (x - (mu+L/2))/(np.sqrt(2)*sig) ) )
|
||||||
|
|
||||||
|
#########################################################
|
||||||
|
def Gauss( x, V, mu, sig, A_L ) :
|
||||||
|
return V * PDF_normal(x, mu, sig) / A_L
|
||||||
|
|
||||||
|
#################################### j0^2 MOMENTS (SCHULZ PDF) #########################################
|
||||||
|
#################################### (checked!) #########################################
|
||||||
|
|
||||||
|
#@njit(parallel=True)
|
||||||
|
def mu4(q, Z, a) :
|
||||||
|
return np.where(q==0, a**4*(Z+1)*(Z+2)*(Z+3)*(Z+3), a**2 * (Z+2)*(Z+1) * ( 1 - (1+4*q*q*a*a)**(-(Z+3)/2.) * np.cos((Z+3)*np.arctan(2*q*a)) ) / ( 2*q**2 ) )
|
||||||
|
|
||||||
|
##################################################################################################################
|
||||||
|
##################################################################################################################
|
||||||
|
|
||||||
|
################################ SYMMETRIC VESICLE FOR X-RAY SLDs #########################################
|
||||||
|
|
||||||
|
######################################### SDP MODELLING #########################################
|
||||||
|
#################################### SEPARATED FORM FACTOR #########################################
|
||||||
|
|
||||||
|
##################################################################################################################
|
||||||
|
##################################################################################################################
|
||||||
|
|
||||||
|
#########################################################
|
||||||
|
######### Symmetric POPC bilayer ########################
|
||||||
|
######### liposomes and proteoliposomes #################
|
||||||
|
#########################################################
|
||||||
|
|
||||||
|
class SDP_base_POPC_RecBuf:
|
||||||
|
|
||||||
|
##################
|
||||||
|
def __init__(self, q, PAR) :
|
||||||
|
self.q = q
|
||||||
|
[self.Norm, self.nv,
|
||||||
|
self.Rm, self.Z,
|
||||||
|
self.n_TR, self.d_TR, self.s_TR,
|
||||||
|
self.d_Chol, self.s_Chol, self.d_PCN, self.s_PCN, self.d_CG, self.s_CG,
|
||||||
|
self.A_L,
|
||||||
|
self.s_CH2, self.d_CH, self.s_CH, self.s_CH3,
|
||||||
|
self.r_PCN, self.r_CG, self.r12, self.r32,
|
||||||
|
self.T, self.V_BW,
|
||||||
|
self.Con] = PAR
|
||||||
|
|
||||||
|
Cw = p0_Cw + p1_Cw*self.T + p2_Cw*self.T**2 + p3_Cw*self.T**3
|
||||||
|
xtris = ctris / Cw # mole fraction of free TRIS in bulk
|
||||||
|
xEDTA = cEDTA / Cw # mole fraction of free EDTA in bulk
|
||||||
|
|
||||||
|
# Volumes
|
||||||
|
self.V_L = lipid_volume(self.T)
|
||||||
|
V_HW = water_volume(self.T)
|
||||||
|
V_HC = self.V_L - ( (1-x_PG) * V_HL_PC + x_PG * V_HL_PG )
|
||||||
|
|
||||||
|
# Quasi-molecular volumes
|
||||||
|
V_CH2 = V_HC / ( n_CH2 + n_CH*self.r12 + n_CH3*self.r32 ) # Volume of CH2 groups
|
||||||
|
V_CH = V_CH2 * self.r12 # Volume of CH groups
|
||||||
|
V_CH3 = V_CH2 * self.r32 # Volume of CH3 groups
|
||||||
|
|
||||||
|
self.V_CG = V_HL_PC * self.r_CG # Volume of CG group
|
||||||
|
self.V_PCN = V_HL_PC * self.r_PCN # Volume of PCN group
|
||||||
|
self.V_Chol = V_HL_PC * (1-self.r_PCN-self.r_CG) # Volume of CholCH3 group
|
||||||
|
|
||||||
|
V_PG1 = V_HL_PG * 0.16 # Kucerka 2012
|
||||||
|
V_PG2 = V_HL_PG * ( 1 - 0.51 - 0.16) # Kucerka 2012
|
||||||
|
|
||||||
|
# Calculation of mean D_C
|
||||||
|
self.D_C = V_HC / self.A_L
|
||||||
|
|
||||||
|
# X-ray scattering lengths (nm)
|
||||||
|
rho_sol = ( b_HW + xtris*b_tris + xEDTA*b_EDTA ) / V_HW
|
||||||
|
drho_Chol = ( (1-x_PG)*b_Chol/self.V_Chol + x_PG*b_PG2/V_PG2 ) - rho_sol
|
||||||
|
drho_PCN = ( (1-x_PG)*b_PCN/self.V_PCN + x_PG*b_PG1/V_PG1 ) - rho_sol
|
||||||
|
drho_CG = b_CG / self.V_CG - rho_sol
|
||||||
|
drho_TR = b_tris/ V_tris - rho_sol
|
||||||
|
drho_CH = b_CH / V_CH - rho_sol
|
||||||
|
drho_CH2 = b_CH2 / V_CH2 - rho_sol
|
||||||
|
drho_CH3 = b_CH3 / V_CH3 - rho_sol
|
||||||
|
drho_HW = b_HW / self.V_BW - rho_sol
|
||||||
|
|
||||||
|
# c-prefactors
|
||||||
|
c_Chol = ( (1-x_PG)*self.V_Chol + x_PG*V_PG2 ) / self.A_L
|
||||||
|
c_PCN = ( (1-x_PG)*self.V_PCN + x_PG*V_PG1 ) / self.A_L
|
||||||
|
c_CG = self.V_CG / self.A_L
|
||||||
|
c_TR = V_tris*self.n_TR / self.A_L
|
||||||
|
c_CH = V_CH * n_CH / self.A_L
|
||||||
|
c_CH3 = V_CH3 * n_CH3 / self.A_L
|
||||||
|
|
||||||
|
# calculating scattering amplitude
|
||||||
|
self.Am=np.zeros(q.shape[0],dtype=float)
|
||||||
|
|
||||||
|
# Adding hydrocarbon-chain envelope
|
||||||
|
self.Am += 2 * drho_CH2 *self.D_C * FTreal_erf(self.q, 0, 2*self.D_C, self.s_CH2)
|
||||||
|
# Adding CH and CH3 groups
|
||||||
|
self.Am += 2 * (drho_CH - drho_CH2) * c_CH * FTreal_gauss(self.q, self.d_CH, self.s_CH)
|
||||||
|
self.Am += 2 * (drho_CH3 - drho_CH2) * c_CH3 * FTreal_gauss(self.q, 0, self.s_CH3)
|
||||||
|
|
||||||
|
# Adding hydration-water envelope
|
||||||
|
self.Am += 4 * drho_HW * ( self.d_CG + self.d_PCN + self.d_Chol + d_shl) * FTreal_erf(self.q, (self.D_C+(self.d_CG+self.d_PCN+self.d_Chol+d_shl)/2.), (self.d_CG+self.d_PCN+self.d_Chol+d_shl), self.s_CH2)
|
||||||
|
# Adding CG, PCN and CholCH3 groups
|
||||||
|
self.Am += 2 * (drho_TR - drho_HW) * c_TR * FTreal_gauss(self.q, (self.D_C+self.d_TR/2.), self.s_TR)
|
||||||
|
self.Am += 2 * (drho_CG - drho_HW) * c_CG * FTreal_gauss(self.q, (self.D_C+self.d_CG/2.), self.s_CG)
|
||||||
|
self.Am += 2 * (drho_PCN - drho_HW) * c_PCN * FTreal_gauss(self.q, (self.D_C+self.d_CG+self.d_PCN/2.), self.s_PCN)
|
||||||
|
self.Am += 2 * (drho_Chol - drho_HW) * c_Chol * FTreal_gauss(self.q, (self.D_C+self.d_CG+self.d_PCN+self.d_Chol/2.), self.s_Chol)
|
||||||
|
|
||||||
|
##################
|
||||||
|
def amplitude(self):
|
||||||
|
return self.Am
|
||||||
|
|
||||||
|
##################
|
||||||
|
def intensity(self):
|
||||||
|
alp = self.Rm/(self.Z+1)
|
||||||
|
return ( self.Norm * self.nv*1e-6 ) * self.Am**2 * ( 16*np.pi**2*mu4(self.q,self.Z,alp) ) + self.Con*( 0.99*(1./(1+np.exp(-8*(self.q-1.)))) + 0.01 )
|
||||||
|
|
||||||
|
##################
|
||||||
|
def negative_water(self):
|
||||||
|
|
||||||
|
self.check = 0
|
||||||
|
z_array = np.linspace(0.,4.,81)
|
||||||
|
|
||||||
|
CG = Gauss(z_array, self.V_CG, self.D_C+self.d_CG/2., self.s_CG, self.A_L)
|
||||||
|
PCN = Gauss(z_array, self.V_PCN, self.D_C+self.d_CG+self.d_PCN/2., self.s_PCN, self.A_L)
|
||||||
|
Chol = Gauss(z_array, self.V_Chol, self.D_C+self.d_CG+self.d_PCN+self.d_Chol/2., self.s_Chol, self.A_L)
|
||||||
|
TRIS = Gauss(z_array, self.n_TR*V_tris, self.D_C+self.d_TR/2., self.s_TR, self.A_L)
|
||||||
|
BW = Slab(z_array, self.D_C+(self.d_CG+self.d_PCN+self.d_Chol+d_shl)/2., self.d_CG+self.d_PCN+self.d_Chol+d_shl, self.s_CH2) - CG - PCN - Chol - TRIS
|
||||||
|
|
||||||
|
for i in(BW) :
|
||||||
|
if i <-0.001 : self.check+= 1
|
||||||
|
|
||||||
|
return self.check
|
||||||
|
##################################################################################################################
|
||||||
|
##################################################################################################################
|
||||||
|
|
||||||
|
class SDP_POPC_RecBuf:
|
||||||
|
|
||||||
|
##################
|
||||||
|
def __init__(self, q, PAR) :
|
||||||
|
self.q = q
|
||||||
|
[self.Norm, self.nv,
|
||||||
|
self.Rm, self.Z,
|
||||||
|
self.n_TR, self.d_TR, self.s_TR,
|
||||||
|
self.d_Chol, self.s_Chol, self.d_PCN, self.s_PCN, self.d_CG, self.s_CG,
|
||||||
|
self.A_L, self.s_D_C,
|
||||||
|
self.s_CH2, self.d_CH, self.s_CH, self.s_CH3,
|
||||||
|
self.r_PCN, self.r_CG, self.r12, self.r32,
|
||||||
|
self.T, self.V_BW,
|
||||||
|
self.Con] = PAR
|
||||||
|
|
||||||
|
Cw = p0_Cw + p1_Cw*self.T + p2_Cw*self.T**2 + p3_Cw*self.T**3
|
||||||
|
xtris = ctris / Cw # mole fraction of free TRIS in bulk
|
||||||
|
xEDTA = cEDTA / Cw # mole fraction of free EDTA in bulk
|
||||||
|
|
||||||
|
# Volumes
|
||||||
|
self.V_L = lipid_volume(self.T)
|
||||||
|
V_HW = water_volume(self.T)
|
||||||
|
V_HC = self.V_L - ( (1-x_PG) * V_HL_PC + x_PG * V_HL_PG )
|
||||||
|
|
||||||
|
# Calculation of mean D_C
|
||||||
|
self.D_C = V_HC / self.A_L
|
||||||
|
|
||||||
|
# Quasi-molecular volumes
|
||||||
|
V_CH2 = V_HC / ( n_CH2 + n_CH*self.r12 + n_CH3*self.r32 ) # Volume of CH2 groups
|
||||||
|
V_CH = V_CH2 * self.r12 # Volume of CH groups
|
||||||
|
V_CH3 = V_CH2 * self.r32 # Volume of CH3 groups
|
||||||
|
|
||||||
|
self.V_CG = V_HL_PC * self.r_CG # Volume of CG group
|
||||||
|
self.V_PCN = V_HL_PC * self.r_PCN # Volume of PCN group
|
||||||
|
self.V_Chol = V_HL_PC * (1-self.r_PCN-self.r_CG) # Volume of CholCH3 group
|
||||||
|
|
||||||
|
V_PG1 = V_HL_PG * 0.16 # Kucerka 2012
|
||||||
|
V_PG2 = V_HL_PG * ( 1 - 0.51 - 0.16) # Kucerka 2012
|
||||||
|
|
||||||
|
############### X-ray scattering lengths (nm)
|
||||||
|
rho_sol = ( b_HW + xtris*b_tris + xEDTA*b_EDTA ) / V_HW
|
||||||
|
drho_Chol = ( (1-x_PG)*b_Chol/self.V_Chol + x_PG*b_PG2/V_PG2 ) - rho_sol
|
||||||
|
drho_PCN = ( (1-x_PG)*b_PCN/self.V_PCN + x_PG*b_PG1/V_PG1 ) - rho_sol
|
||||||
|
drho_CG = b_CG / self.V_CG - rho_sol
|
||||||
|
drho_TR = b_tris/ V_tris - rho_sol
|
||||||
|
drho_CH = b_CH / V_CH - rho_sol
|
||||||
|
drho_CH2 = b_CH2 / V_CH2 - rho_sol
|
||||||
|
drho_CH3 = b_CH3 / V_CH3 - rho_sol
|
||||||
|
drho_HW = b_HW / self.V_BW - rho_sol
|
||||||
|
|
||||||
|
############### D_C polydispersity
|
||||||
|
N = 21
|
||||||
|
HC_array = np.linspace(self.D_C-3*self.s_D_C, self.D_C+3*self.s_D_C, N)
|
||||||
|
Normal = PDF_normal(HC_array, self.D_C, self.s_D_C)
|
||||||
|
|
||||||
|
############### calculating scattering amplitude -----------------------------------------------
|
||||||
|
self.Am = np.zeros([HC_array.shape[0],self.q.shape[0]],dtype=float)
|
||||||
|
c_CH = np.zeros(HC_array.shape[0],dtype=float)
|
||||||
|
c_CH3 = np.zeros(HC_array.shape[0],dtype=float)
|
||||||
|
|
||||||
|
############### c-prefactors
|
||||||
|
c_Chol = ( (1-x_PG)*self.V_Chol + x_PG*V_PG2 ) / self.A_L
|
||||||
|
c_PCN = ( (1-x_PG)*self.V_PCN + x_PG*V_PG1 ) / self.A_L
|
||||||
|
c_CG = self.V_CG / self.A_L
|
||||||
|
c_TR = V_tris*self.n_TR / self.A_L
|
||||||
|
|
||||||
|
for hc in range(HC_array.shape[0]):
|
||||||
|
c_CH[hc] = V_CH * n_CH / (V_HC / HC_array[hc] )
|
||||||
|
c_CH3[hc] = V_CH3 * n_CH3 / (V_HC / HC_array[hc] )
|
||||||
|
|
||||||
|
for hc in range(HC_array.shape[0]):
|
||||||
|
|
||||||
|
# Adding hydrocarbon-chain envelope
|
||||||
|
self.Am[hc] += 2 * drho_CH2 *HC_array[hc] * FTreal_erf(self.q, 0, 2*HC_array[hc], self.s_CH2)
|
||||||
|
# Adding CH and CH3 groups
|
||||||
|
self.Am[hc] += 2 * (drho_CH - drho_CH2) * c_CH[hc] * FTreal_gauss(self.q, self.d_CH, self.s_CH)
|
||||||
|
self.Am[hc] += 2 * (drho_CH3 - drho_CH2) * c_CH3[hc] * FTreal_gauss(self.q, 0, self.s_CH3)
|
||||||
|
|
||||||
|
# Adding hydration-water envelope
|
||||||
|
self.Am[hc] += 4 * drho_HW * ( self.d_CG + self.d_PCN + self.d_Chol + d_shl) * FTreal_erf(self.q, (HC_array[hc]+(self.d_CG+self.d_PCN+self.d_Chol+d_shl)/2.), (self.d_CG+self.d_PCN+self.d_Chol+d_shl), self.s_CH2)
|
||||||
|
# Adding CG, PCN and CholCH3 groups
|
||||||
|
self.Am[hc] += 2 * (drho_TR - drho_HW) * c_TR * FTreal_gauss(self.q, (HC_array[hc]+self.d_TR/2.), self.s_TR)
|
||||||
|
self.Am[hc] += 2 * (drho_CG - drho_HW) * c_CG * FTreal_gauss(self.q, (HC_array[hc]+self.d_CG/2.), self.s_CG)
|
||||||
|
self.Am[hc] += 2 * (drho_PCN - drho_HW) * c_PCN * FTreal_gauss(self.q, (HC_array[hc]+self.d_CG+self.d_PCN/2.), self.s_PCN)
|
||||||
|
self.Am[hc] += 2 * (drho_Chol - drho_HW) * c_Chol * FTreal_gauss(self.q, (HC_array[hc]+self.d_CG+self.d_PCN+self.d_Chol/2.), self.s_Chol)
|
||||||
|
|
||||||
|
############### Ensemble average
|
||||||
|
self.I = np.zeros(self.q.shape[0], dtype=float)
|
||||||
|
|
||||||
|
for hc in range(HC_array.shape[0]):
|
||||||
|
if hc==0 or hc==N-1 : self.I+= self.Am[hc]**2 * Normal[hc] / 2
|
||||||
|
else : self.I+= self.Am[hc]**2 * Normal[hc]
|
||||||
|
|
||||||
|
self.I*= 6*self.s_D_C/(N-1)
|
||||||
|
|
||||||
|
##################
|
||||||
|
def intensity(self):
|
||||||
|
alp = self.Rm/(self.Z+1)
|
||||||
|
return ( self.Norm * self.nv*1e-6 ) * self.I * ( 16*np.pi**2*mu4(self.q,self.Z,alp) ) + self.Con*( 0.99*(1./(1+np.exp(-8*(self.q-1.)))) + 0.01 )
|
||||||
|
|
||||||
|
##################
|
||||||
|
def negative_water(self):
|
||||||
|
|
||||||
|
self.check = 0
|
||||||
|
z_array = np.linspace(0.,4.,81)
|
||||||
|
|
||||||
|
CG = Gauss(z_array, self.V_CG, self.D_C+self.d_CG/2., self.s_CG, self.A_L)
|
||||||
|
PCN = Gauss(z_array, self.V_PCN, self.D_C+self.d_CG+self.d_PCN/2., self.s_PCN, self.A_L)
|
||||||
|
Chol = Gauss(z_array, self.V_Chol, self.D_C+self.d_CG+self.d_PCN+self.d_Chol/2., self.s_Chol, self.A_L)
|
||||||
|
TRIS = Gauss(z_array, self.n_TR*V_tris, self.D_C+self.d_TR/2., self.s_TR, self.A_L)
|
||||||
|
BW = Slab(z_array, self.D_C+(self.d_CG+self.d_PCN+self.d_Chol+d_shl)/2., self.d_CG+self.d_PCN+self.d_Chol+d_shl, self.s_CH2) - CG - PCN - Chol - TRIS
|
||||||
|
|
||||||
|
for i in(BW) :
|
||||||
|
if i <-0.001 : self.check+= 1
|
||||||
|
|
||||||
|
return self.check
|
||||||
|
|
||||||
|
##################################################################################################################
|
||||||
|
##################################################################################################################
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
582
Minimization/Tools.py
Executable file
582
Minimization/Tools.py
Executable file
|
|
@ -0,0 +1,582 @@
|
||||||
|
#!/usr/bin/python
|
||||||
|
|
||||||
|
import sys
|
||||||
|
import os.path
|
||||||
|
import re
|
||||||
|
import time
|
||||||
|
import numpy as np
|
||||||
|
|
||||||
|
import FormFactor
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
|
||||||
|
###############
|
||||||
|
############### Miscellaneous
|
||||||
|
|
||||||
|
def file_exists(arg, name):
|
||||||
|
if os.path.exists(arg):
|
||||||
|
print(name,":\t",arg)
|
||||||
|
else:
|
||||||
|
raise IOError(arg,"The file does not exist")
|
||||||
|
|
||||||
|
###############
|
||||||
|
|
||||||
|
def check_int(strg, name):
|
||||||
|
try:
|
||||||
|
return int(strg)
|
||||||
|
except ValueError:
|
||||||
|
raise ValueError(name," must be an int number")
|
||||||
|
|
||||||
|
###############
|
||||||
|
|
||||||
|
def check_binary(arg, name):
|
||||||
|
if arg ==0 or arg == 1:
|
||||||
|
pass
|
||||||
|
else:
|
||||||
|
raise Exception(name," must be 0 or 1")
|
||||||
|
|
||||||
|
###############
|
||||||
|
|
||||||
|
def check_float_positive(strg, name):
|
||||||
|
try:
|
||||||
|
float(strg)
|
||||||
|
except ValueError:
|
||||||
|
raise ValueError(name," must be a float number")
|
||||||
|
if float(strg) >= 0:
|
||||||
|
return float(strg)
|
||||||
|
else:
|
||||||
|
raise Exception(name," must be a positive number")
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
|
||||||
|
class Reading:
|
||||||
|
'Read the Options & Parameters file'
|
||||||
|
EMPTY=[]
|
||||||
|
|
||||||
|
def __init__ ( self , ARGV ):
|
||||||
|
|
||||||
|
if len(ARGV)>3:
|
||||||
|
sys.exit("--- Too many arguments")
|
||||||
|
elif len(ARGV)==1:
|
||||||
|
sys.exit("--- Specify the parameters file")
|
||||||
|
|
||||||
|
file_exists(ARGV[1], "Parameters File")
|
||||||
|
|
||||||
|
with open(ARGV[1], 'r') as f:
|
||||||
|
LINES = f.read().splitlines()
|
||||||
|
|
||||||
|
TEMP=[]
|
||||||
|
OPT=[]
|
||||||
|
PARA=[]
|
||||||
|
for i in range(len(LINES)):
|
||||||
|
if not LINES[i]=='':
|
||||||
|
TEMP.append(LINES[i])
|
||||||
|
else:
|
||||||
|
TEMP.append("#")
|
||||||
|
Reading.EMPTY.append(i)
|
||||||
|
|
||||||
|
line = 0
|
||||||
|
while line < Reading.EMPTY[0]:
|
||||||
|
OPT.append(TEMP[line].split())
|
||||||
|
line += 1
|
||||||
|
line = Reading.EMPTY[1]+1
|
||||||
|
while line <= len(TEMP)-1:
|
||||||
|
PARA.append(TEMP[line].split())
|
||||||
|
line += 1
|
||||||
|
|
||||||
|
self.OPT=OPT
|
||||||
|
self.PARA=PARA
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
|
||||||
|
def Options ( self ):
|
||||||
|
'Read Options'
|
||||||
|
|
||||||
|
err_str = "\n--- Something missing or wrong.\n--- You need:\n"
|
||||||
|
err_str += "DATA_FILE: <filename>\n"
|
||||||
|
err_str += "TITLE: string\n"
|
||||||
|
err_str += "SAVE_FOLDER: string\n"
|
||||||
|
err_str += "q_MIN(nm^-1): float\n"
|
||||||
|
err_str += "q_MAX(nm^-1): float\n"
|
||||||
|
err_str += "MODEL: string\n"
|
||||||
|
err_str += "TEMPERATURE(none=0): float or 0\n"
|
||||||
|
err_str += "BINNING: 4\n"
|
||||||
|
err_str += "ERROR(0..1): float\n"
|
||||||
|
err_str += "PLOT(0-1): int 0/1"
|
||||||
|
|
||||||
|
options={}
|
||||||
|
|
||||||
|
for i in range(len(self.OPT)):
|
||||||
|
if re.match( "^DATA_FILE:$", self.OPT[i][0] ):
|
||||||
|
options['datafile'] = str(self.OPT[i][1])
|
||||||
|
file_exists(options['datafile'], "Data File")
|
||||||
|
|
||||||
|
elif re.match( "^TITLE:$", self.OPT[i][0] ):
|
||||||
|
options['title'] = str(self.OPT[i][1])
|
||||||
|
print("Title:\t\t", options['title'])
|
||||||
|
|
||||||
|
elif re.match( "^SAVE_FOLDER:$", self.OPT[i][0] ):
|
||||||
|
options['folder'] = str(self.OPT[i][1])
|
||||||
|
if os.path.exists(options['folder'])==0:
|
||||||
|
os.makedirs(options['folder'])
|
||||||
|
print("--- Creating new folder: ", options['folder'])
|
||||||
|
print("Saving Folder:\t", options['folder'])
|
||||||
|
elif os.path.exists(options['folder'])==1:
|
||||||
|
print("Saving Folder:\t", options['folder'])
|
||||||
|
|
||||||
|
elif re.match( "^q_MIN\(nm\^-1\):$", self.OPT[i][0] ):
|
||||||
|
options['qmin'] = check_float_positive(self.OPT[i][1], "qmin")
|
||||||
|
print("q range:\t(%s .." %options['qmin'], )
|
||||||
|
|
||||||
|
elif re.match( "^q_MAX\(nm\^-1\):$", self.OPT[i][0] ):
|
||||||
|
options['qmax'] = check_float_positive(self.OPT[i][1], "qmax")
|
||||||
|
print("%s) nm^-1\n" %options['qmax'], )
|
||||||
|
|
||||||
|
elif re.match( "^MODEL:$", self.OPT[i][0] ):
|
||||||
|
options['function'] = str(self.OPT[i][1])
|
||||||
|
print("Function:\t", options['function'])
|
||||||
|
|
||||||
|
elif re.match( "^TEMPERATURE\(none=0\):$", self.OPT[i][0] ):
|
||||||
|
options['temp'] = check_float_positive(self.OPT[i][1], "Temperarure")
|
||||||
|
if options['temp'] != 0:
|
||||||
|
print("Temperarure:\t", options['temp'])
|
||||||
|
else:
|
||||||
|
print("Temperarure:\tNo Simulated Annealing")
|
||||||
|
|
||||||
|
elif re.match( "^BINNING:$", self.OPT[i][0] ):
|
||||||
|
options['bing'] = check_int(self.OPT[i][1], "Binning")
|
||||||
|
print("Binning:\t", options['bing'])
|
||||||
|
|
||||||
|
elif re.match( "^ERROR\(0\.\.1\):$", self.OPT[i][0] ):
|
||||||
|
options['err_mul'] = check_float_positive(self.OPT[i][1], "Error Mult.")
|
||||||
|
print("Error Mult.:\t", options['err_mul'])
|
||||||
|
|
||||||
|
elif re.match( "^PLOT\(0-1\):$", self.OPT[i][0] ):
|
||||||
|
options['plot'] = check_int(self.OPT[i][1], "Plot")
|
||||||
|
check_binary(options['plot'], "Plot")
|
||||||
|
print("Plot:\t\t", options['plot'])
|
||||||
|
|
||||||
|
else:
|
||||||
|
raise Exception(err_str)
|
||||||
|
|
||||||
|
if len(options)!=10:
|
||||||
|
raise Exception(err_str)
|
||||||
|
|
||||||
|
return options
|
||||||
|
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
|
||||||
|
def Parameters ( self ):
|
||||||
|
'Read Parameters'
|
||||||
|
|
||||||
|
############### Define parameters arrays
|
||||||
|
NAME=[]
|
||||||
|
PAR=[]
|
||||||
|
FIX=[]
|
||||||
|
PRIOR=[]
|
||||||
|
L_LIM=[]
|
||||||
|
H_LIM=[]
|
||||||
|
for i in range(len(self.PARA)):
|
||||||
|
if not re.match( "^#" , self.PARA[i][0] ):
|
||||||
|
NAME.append(self.PARA[i][0])
|
||||||
|
PAR.append(float(self.PARA[i][1]))
|
||||||
|
FIX.append(self.PARA[i][2])
|
||||||
|
PRIOR.append(self.PARA[i][3])
|
||||||
|
L_LIM.append(float(self.PARA[i][4]))
|
||||||
|
H_LIM.append(float(self.PARA[i][5]))
|
||||||
|
|
||||||
|
############### Check if there are free parameters
|
||||||
|
n=0
|
||||||
|
for i in range(len(FIX)):
|
||||||
|
if FIX[i]=="f":
|
||||||
|
n += 1
|
||||||
|
if n==len(PAR):
|
||||||
|
sys.exit("--- All the parameters are fixed; Choose the \"Plot\" option")
|
||||||
|
|
||||||
|
############### Check & Print parameters features
|
||||||
|
for i in range(len(NAME)):
|
||||||
|
if FIX[i]=="f":
|
||||||
|
f="Fixed"
|
||||||
|
print(i,"\t%4s\t%.3e\t%s" %(NAME[i],PAR[i],f))
|
||||||
|
else:
|
||||||
|
if PAR[i]==0:
|
||||||
|
sys.exit("--- 0 value as free parameter is not allowed")
|
||||||
|
if FIX[i]=="-" and PRIOR[i]=="-":
|
||||||
|
f="Free"
|
||||||
|
print(i,"\t%4s\t%.3e\t%s\t\thard boundaries (%.3e, %.3e)" %(NAME[i],PAR[i],f,L_LIM[i],H_LIM[i]))
|
||||||
|
elif FIX[i]=="-" and PRIOR[i]!="-":
|
||||||
|
f="Prior"
|
||||||
|
print(i,"\t%4s\t%.3e\t%s %.5s\thard boundaries (%.3e, %.3e)" %(NAME[i],PAR[i],f,PRIOR[i],L_LIM[i],H_LIM[i]))
|
||||||
|
else:
|
||||||
|
sys.exit("--- Free/Fixed -> -/f")
|
||||||
|
|
||||||
|
if H_LIM[i]-L_LIM[i]<=0:
|
||||||
|
sys.exit("--- Check low & high parameter boundary")
|
||||||
|
if PAR[i]<L_LIM[i]:
|
||||||
|
sys.exit("--- Check low & high parameter boundary")
|
||||||
|
if PAR[i]>H_LIM[i]:
|
||||||
|
sys.exit("--- Check low & high parameter boundary")
|
||||||
|
|
||||||
|
return ( NAME , PAR , FIX, PRIOR, L_LIM , H_LIM )
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
|
||||||
|
def Parameters_refresh ( self ):
|
||||||
|
'Read Parameters'
|
||||||
|
|
||||||
|
############### Define parameters arrays
|
||||||
|
NAME=[]
|
||||||
|
PAR=[]
|
||||||
|
FIX=[]
|
||||||
|
PRIOR=[]
|
||||||
|
L_LIM=[]
|
||||||
|
H_LIM=[]
|
||||||
|
for i in range(len(self.PARA)):
|
||||||
|
if not re.match( "^#" , self.PARA[i][0] ):
|
||||||
|
NAME.append(self.PARA[i][0])
|
||||||
|
PAR.append(float(self.PARA[i][1]))
|
||||||
|
FIX.append(self.PARA[i][2])
|
||||||
|
PRIOR.append(self.PARA[i][3])
|
||||||
|
L_LIM.append(float(self.PARA[i][4]))
|
||||||
|
H_LIM.append(float(self.PARA[i][5]))
|
||||||
|
|
||||||
|
return ( NAME , PAR , FIX, PRIOR, L_LIM , H_LIM )
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
|
||||||
|
class ReadData:
|
||||||
|
'Read raw data from file'
|
||||||
|
|
||||||
|
def __init__ ( self , datafile ):
|
||||||
|
|
||||||
|
DATA = []
|
||||||
|
|
||||||
|
with open(datafile, 'r') as d:
|
||||||
|
DLINES = d.read().splitlines()
|
||||||
|
for row in DLINES:
|
||||||
|
DATA.append(row.split())
|
||||||
|
|
||||||
|
DATA.pop(0)
|
||||||
|
|
||||||
|
for i in range(len(DATA)):
|
||||||
|
for j in range(len(DATA[i])):
|
||||||
|
DATA[i][j] = float(DATA[i][j])
|
||||||
|
|
||||||
|
self.DATA = DATA
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
|
||||||
|
def SelectColumns ( self , qmin , qmax , bing , err_mul ):
|
||||||
|
'Select q / I(q) / err(q) data - Change I(q) value from mm^-1 to nm^-1'
|
||||||
|
|
||||||
|
Q = []
|
||||||
|
IDATA = []
|
||||||
|
ERR = []
|
||||||
|
|
||||||
|
for i in range(len(self.DATA)):
|
||||||
|
if ( self.DATA[i][0] >= qmin and self.DATA[i][0] <= qmax ) :
|
||||||
|
if i%bing == 0 :
|
||||||
|
Q.append(self.DATA[i][0])
|
||||||
|
if self.DATA[i][1] > 0 :
|
||||||
|
IDATA.append(1e-6*self.DATA[i][1])
|
||||||
|
if err_mul == 0 :
|
||||||
|
ERR.append(0.01*1e-6*self.DATA[i][1])
|
||||||
|
elif self.DATA[i][2]*err_mul > 0.01*self.DATA[i][1] :
|
||||||
|
ERR.append(1e-6*self.DATA[i][2]*err_mul)
|
||||||
|
else:
|
||||||
|
ERR.append(0.01*1e-6*self.DATA[i][1])
|
||||||
|
else:
|
||||||
|
Q.pop()
|
||||||
|
|
||||||
|
return ( np.array(Q), np.array(IDATA), np.array(ERR) )
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
|
||||||
|
class PlotData:
|
||||||
|
'Plot Results'
|
||||||
|
|
||||||
|
def __init__ ( self , title , folder , qmin , qmax ):
|
||||||
|
|
||||||
|
self.fitFile="./"+folder+"/"+title+".fit"
|
||||||
|
self.gnuFile="./"+folder+"/"+title+".gp"
|
||||||
|
self.intensityFile="./"+folder+"/"+title+"_I.pdf"
|
||||||
|
self.structureFile="./"+folder+"/"+title+"_Sq.pdf"
|
||||||
|
self.errorFile="./"+folder+"/"+title+"_err.pdf"
|
||||||
|
self.qmin=str(qmin)
|
||||||
|
self.qmax=str(qmax)
|
||||||
|
self.gnuStr0= str( "set term pdfcairo enhanced solid color font \"Helvetica,14\" \
|
||||||
|
lw 2 rounded size 5 inch, 3.5 inch\n \
|
||||||
|
set encoding iso_8859_1\n \
|
||||||
|
set grid\n \
|
||||||
|
set key box\n \
|
||||||
|
set bars small\n \
|
||||||
|
set xtics add ("+str(qmin)+","+str(qmax)+")\n \
|
||||||
|
set ytics format \"10^\{%T\}\"\n \
|
||||||
|
set mytics 10\n \
|
||||||
|
set title \""+title+"\"\n" )
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
|
||||||
|
def Form ( self, Q, IDATA, ERR, I, PAR, pltOptions ):
|
||||||
|
'Plot the best fit of a Form Factor'
|
||||||
|
|
||||||
|
############### Check Cluster Structure Factor
|
||||||
|
if pltOptions['debyeCtr'] != 0:
|
||||||
|
|
||||||
|
############### Save constant value
|
||||||
|
if conv == 0:
|
||||||
|
C = str(1e6*PAR[-1])
|
||||||
|
else:
|
||||||
|
C = str(1e6*PAR[-2])
|
||||||
|
|
||||||
|
############### Prepare Form Factor
|
||||||
|
PAR_FORM=list(PAR)
|
||||||
|
for i in range(pltOptions['parDiff']):
|
||||||
|
PAR_FORM.pop(-1)
|
||||||
|
PAR_FORM.append(0.0)
|
||||||
|
if conv == 1:
|
||||||
|
PAR_FORM.append(PAR[-1])
|
||||||
|
FORM = [ pltOptions['Form'](q, PAR_FORM) for q in Q ]
|
||||||
|
|
||||||
|
############### Prepare Cluster Structure Factor
|
||||||
|
if conv == 0:
|
||||||
|
Io = PAR[-4]
|
||||||
|
Rg = PAR[-3]
|
||||||
|
D = PAR[-2]
|
||||||
|
else:
|
||||||
|
Io = PAR[-5]
|
||||||
|
Rg = PAR[-4]
|
||||||
|
D = PAR[-3]
|
||||||
|
DEBYE = [ pltOptions['debyeCtr'](q,Io,Rg,D) for q in Q ]
|
||||||
|
|
||||||
|
############### Write raw data file
|
||||||
|
fit = open ( self.fitFile, "w" )
|
||||||
|
if pltOptions['debyeCtr'] == 0:
|
||||||
|
for i in range(len(Q)):
|
||||||
|
fit.write( "%f\t%f\t%f\t%f\n" %( Q[i], 1e6*IDATA[i], 1e6*ERR[i], 1e6*I[i] ) )
|
||||||
|
else:
|
||||||
|
for i in range(len(Q)):
|
||||||
|
fit.write( "%f\t%f\t%f\t%f\t%f\t%f\n" %( Q[i], 1e6*IDATA[i], 1e6*ERR[i], 1e6*I[i], 1e6*FORM[i], DEBYE[i] ) )
|
||||||
|
fit.close()
|
||||||
|
|
||||||
|
############### Write GNUPLOT script
|
||||||
|
gnu = open ( self.gnuFile, "w" )
|
||||||
|
gnu.write( self.gnuStr0 )
|
||||||
|
|
||||||
|
gnu.write( "set logscale\n" )
|
||||||
|
gnu.write( "set xlabel \"q (nm^{-1})\"\n" )
|
||||||
|
gnu.write( "set ylabel \"I (mm^{-1})\"\n" )
|
||||||
|
gnu.write( "set out '"+self.intensityFile+"'\n" )
|
||||||
|
gnu.write( "p ["+self.qmin+":"+self.qmax+"][] '"+self.fitFile+"' \
|
||||||
|
u 1:2 w p pt 6 ps 0.5 lt -1 t '', '' u 1:2:3 w errorbars pt 7 ps 0.3 lt 1 t 'data', \
|
||||||
|
'' u 1:4 w l lt 3 lw 2.5 t 'Fit'\n" )
|
||||||
|
|
||||||
|
gnu.write( "set ylabel \"| I^{data}-I^{fit} |/I^{data} (mm^{-1})\"\n" )
|
||||||
|
gnu.write( "set out '"+self.errorFile+"'\n" )
|
||||||
|
gnu.write( "p ["+self.qmin+":"+self.qmax+"][] '"+self.fitFile+"' \
|
||||||
|
u 1:($3/$2) w l lt rgb \"web-green\" lw 2.5 t 'err.', \
|
||||||
|
'' u 1:(abs($2-$4)/$2) w p pt 7 ps 0.6 lt 3 t '', \
|
||||||
|
'' u 1:(abs($2-$4)/$2) w p pt 6 ps 0.6 lt -1 t ''\n" )
|
||||||
|
|
||||||
|
if pltOptions['debyeCtr'] != 0:
|
||||||
|
gnu.write( "unset logscale y\n" )
|
||||||
|
gnu.write( "set xtics auto\n" )
|
||||||
|
gnu.write( "set xtics add ("+self.qmax+")\n" )
|
||||||
|
gnu.write( "set ylabel \"S(q)\"\n" )
|
||||||
|
gnu.write( "set out '"+self.structureFile+"'\n" )
|
||||||
|
gnu.write( "p ["+self.qmin+":"+self.qmax+"][] '"+self.fitFile+"' \
|
||||||
|
u 1:(($2-"+C+")/$5) w p pt 6 ps 0.5 lt -1 t '', \
|
||||||
|
'' u 1:(($2-"+C+")/$5) w p pt 7 ps 0.3 lt 1 t 'data', \
|
||||||
|
'' u 1:(1+$6) w l lt rgb \"purple\" lw 2.5 t '1+D(q)'\n" )
|
||||||
|
|
||||||
|
|
||||||
|
gnu.close()
|
||||||
|
|
||||||
|
############### Run GNUPLOT script
|
||||||
|
os.system("gnuplot "+self.gnuFile)
|
||||||
|
############### Remove GNUPLOT script
|
||||||
|
os.remove(self.gnuFile)
|
||||||
|
############### Print saved files
|
||||||
|
print(self.fitFile)
|
||||||
|
print(self.intensityFile)
|
||||||
|
print(self.errorFile)
|
||||||
|
|
||||||
|
return
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
|
||||||
|
def Form_Porod ( self , Q , IDATA , ERR , I , PAR ):
|
||||||
|
'Plot the best fit of a Form Factor'
|
||||||
|
|
||||||
|
############### Porod function values
|
||||||
|
if conv==0:
|
||||||
|
Pval = PAR[-3]
|
||||||
|
Pexp = -PAR[-2]
|
||||||
|
else:
|
||||||
|
Pval = PAR[-4]
|
||||||
|
Pexp = -PAR[-3]
|
||||||
|
|
||||||
|
############### Write raw data file
|
||||||
|
fit = open ( self.fitFile, "w" )
|
||||||
|
for i in range(len(Q)):
|
||||||
|
fit.write( "%f\t%f\t%f\t%f\t%f\n" %( Q[i], 1e6*IDATA[i], 1e6*ERR[i], 1e6*I[i], 1e6*Pval*Q[i]**Pexp ) )
|
||||||
|
fit.close()
|
||||||
|
|
||||||
|
############### Write GNUPLOT script
|
||||||
|
gnu = open ( self.gnuFile, "w" )
|
||||||
|
gnu.write( self.gnuStr0 )
|
||||||
|
|
||||||
|
gnu.write( "set logscale\n" )
|
||||||
|
gnu.write( "set xlabel \"q (nm^{-1})\"\n" )
|
||||||
|
gnu.write( "set ylabel \"I (mm^{-1})\"\n" )
|
||||||
|
gnu.write( "set out '"+self.intensityFile+"'\n" )
|
||||||
|
gnu.write( "p ["+self.qmin+":"+self.qmax+"][] '"+self.fitFile+"' \
|
||||||
|
u 1:2 w p pt 7 ps 0.6 lt -1 t '', '' u 1:2:3 w errorbars pt 7 ps 0.5 lt 1 t 'data', \
|
||||||
|
'' u 1:4 w l lt 3 lw 2.5 t 'Fit', \
|
||||||
|
'' u 1:5 w l lt rgb \"web-green\" lw 2 t 'Porod', \
|
||||||
|
'' u 1:($4-$5) w l lt 0 lw 3 t ''\n" )
|
||||||
|
|
||||||
|
gnu.write( "set ylabel \"| I^{data}-I^{fit} |/I^{data} (mm^{-1})\"\n" )
|
||||||
|
gnu.write( "set out '"+self.errorFile+"'\n" )
|
||||||
|
gnu.write( "p ["+self.qmin+":"+self.qmax+"][] '"+self.fitFile+"' \
|
||||||
|
u 1:($3/$2) w l lt -1 lw 2.5 t 'err.', \
|
||||||
|
'' u 1:(abs($2-$4)/$2) w p pt 6 ps 0.6 lt 1 t ''\n" )
|
||||||
|
|
||||||
|
gnu.close()
|
||||||
|
|
||||||
|
############### Run GNUPLOT script
|
||||||
|
os.system("gnuplot "+self.gnuFile)
|
||||||
|
############### Remove GNUPLOT script
|
||||||
|
os.remove(self.gnuFile)
|
||||||
|
############### Print saved files
|
||||||
|
print(self.fitFile)
|
||||||
|
print(self.intensityFile)
|
||||||
|
print(self.errorFile)
|
||||||
|
|
||||||
|
return
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
|
||||||
|
def Structure (self, Q, IDATA, ERR, I, PAR, pltOptions):
|
||||||
|
'Plot the best fit of a Form Factor'
|
||||||
|
|
||||||
|
############### Save constant value
|
||||||
|
if conv == 0:
|
||||||
|
C = str(1e6*PAR[-1])
|
||||||
|
else:
|
||||||
|
C = str(1e6*PAR[-2])
|
||||||
|
|
||||||
|
############### Prepare Form Factor
|
||||||
|
PAR_FORM=list(PAR)
|
||||||
|
for i in range(pltOptions['parDiff']):
|
||||||
|
PAR_FORM.pop(-1)
|
||||||
|
PAR_FORM.append(0.0)
|
||||||
|
if conv == 1:
|
||||||
|
PAR_FORM.append(PAR[-1])
|
||||||
|
FORM = [ pltOptions['Form'](q, PAR_FORM) for q in Q ]
|
||||||
|
|
||||||
|
############### Check Cluster Structure Factor
|
||||||
|
if pltOptions['debyeCtr'] != 0:
|
||||||
|
############### Prepare Cluster Structure Factor
|
||||||
|
if conv == 0:
|
||||||
|
Io = PAR[-4]
|
||||||
|
Rg = PAR[-3]
|
||||||
|
D = PAR[-2]
|
||||||
|
else:
|
||||||
|
Io = PAR[-5]
|
||||||
|
Rg = PAR[-4]
|
||||||
|
D = PAR[-3]
|
||||||
|
DEBYE = [ pltOptions['debyeCtr'](q,Io,Rg,D) for q in Q ]
|
||||||
|
|
||||||
|
############### Write raw data file
|
||||||
|
fit = open ( self.fitFile, "w" )
|
||||||
|
if pltOptions['debyeCtr'] == 0:
|
||||||
|
for i in range(len(Q)):
|
||||||
|
fit.write( "%f\t%f\t%f\t%f\t%f\n" %( Q[i], 1e6*IDATA[i], 1e6*ERR[i], 1e6*I[i], 1e6*FORM[i] ) )
|
||||||
|
else:
|
||||||
|
for i in range(len(Q)):
|
||||||
|
fit.write( "%f\t%f\t%f\t%f\t%f\t%f\n" %( Q[i], 1e6*IDATA[i], 1e6*ERR[i], 1e6*I[i], 1e6*FORM[i], DEBYE[i] ) )
|
||||||
|
fit.close()
|
||||||
|
|
||||||
|
############### Write GNUPLOT script
|
||||||
|
gnu = open ( self.gnuFile, "w" )
|
||||||
|
gnu.write( self.gnuStr0 )
|
||||||
|
|
||||||
|
gnu.write( "set logscale\n" )
|
||||||
|
gnu.write( "set xlabel \"q (nm^{-1})\"\n" )
|
||||||
|
gnu.write( "set ylabel \"I (mm^{-1})\"\n" )
|
||||||
|
gnu.write( "set out '"+self.intensityFile+"'\n" )
|
||||||
|
gnu.write( "p ["+self.qmin+":"+self.qmax+"][] '"+self.fitFile+"' \
|
||||||
|
u 1:2 w p pt 6 ps 0.5 lt -1 t '', '' u 1:2:3 w errorbars pt 7 ps 0.3 lt 1 t 'data', \
|
||||||
|
'' u 1:5 w l lt rgb \"web-green\" lw 2.5 t 'FormFactor', \
|
||||||
|
'' u 1:4 w l lt 3 lw 2.5 t 'Fit'\n" )
|
||||||
|
|
||||||
|
gnu.write( "set ylabel \"| I^{data}-I^{fit} |/I^{data} (mm^{-1})\"\n" )
|
||||||
|
gnu.write( "set out '"+self.errorFile+"'\n" )
|
||||||
|
gnu.write( "p ["+self.qmin+":"+self.qmax+"][] '"+self.fitFile+"' \
|
||||||
|
u 1:($3/$2) w l lt rgb \"web-green\" lw 2.5 t 'err.', \
|
||||||
|
'' u 1:(abs($2-$4)/$2) w p pt 7 ps 0.6 lt 3 t '', \
|
||||||
|
'' u 1:(abs($2-$4)/$2) w p pt 6 ps 0.6 lt -1 t ''\n" )
|
||||||
|
|
||||||
|
gnu.write( "unset logscale\n" )
|
||||||
|
gnu.write( "set xtics auto\n" )
|
||||||
|
gnu.write( "set xtics add ("+self.qmax+")\n" )
|
||||||
|
gnu.write( "set ylabel \"S(q)\"\n" )
|
||||||
|
gnu.write( "set out '"+self.structureFile+"'\n" )
|
||||||
|
if pltOptions['debyeCtr'] == 0:
|
||||||
|
gnu.write( "p ["+self.qmin+":"+self.qmax+"][] '"+self.fitFile+"' \
|
||||||
|
u 1:(($2-"+C+")/$5) w p pt 6 ps 0.5 lt -1 t '', \
|
||||||
|
'' u 1:(($2-"+C+")/$5) w p pt 7 ps 0.3 lt 1 t 'data', \
|
||||||
|
'' u 1:(($4-"+C+")/$5) w l lt 3 lw 2.5 t 'S(q)'\n" )
|
||||||
|
else:
|
||||||
|
gnu.write( "p ["+self.qmin+":"+self.qmax+"][] '"+self.fitFile+"' \
|
||||||
|
u 1:(($2-"+C+")/$5) w p pt 6 ps 0.5 lt -1 t '', \
|
||||||
|
'' u 1:(($2-"+C+")/$5) w p pt 7 ps 0.3 lt 1 t 'data', \
|
||||||
|
'' u 1:(1+$6) w l lt rgb \"purple\" lw 2.5 t '1+D(q)', \
|
||||||
|
'' u 1:(($4-"+C+")/$5-$6) w l lt rgb \"dark-green\" lw 2.5 t 'S(q)', \
|
||||||
|
'' u 1:(($4-"+C+")/$5) w l lt 3 lw 2.5 t 'S(q)+D(q)'\n" )
|
||||||
|
|
||||||
|
gnu.close()
|
||||||
|
|
||||||
|
############### Run GNUPLOT script
|
||||||
|
os.system("gnuplot "+self.gnuFile)
|
||||||
|
############### Remove GNUPLOT script
|
||||||
|
os.remove(self.gnuFile)
|
||||||
|
############### Print saved files
|
||||||
|
print(self.fitFile)
|
||||||
|
print(self.intensityFile)
|
||||||
|
print(self.errorFile)
|
||||||
|
print(self.structureFile)
|
||||||
|
|
||||||
|
return
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
###########################################################################################
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
Loading…
Reference in a new issue