commit a39295781b7906b2021d2b39a39a4a5747f6c06c Author: Gaspard Jankowiak Date: Fri Feb 23 14:46:25 2024 +0100 import diff --git a/Minimization/FitTools.py b/Minimization/FitTools.py new file mode 100755 index 0000000..04b656b --- /dev/null +++ b/Minimization/FitTools.py @@ -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 ) + + +########################################################################################### +########################################################################################### + diff --git a/Minimization/Fit_TSA.py b/Minimization/Fit_TSA.py new file mode 100755 index 0000000..e9fdadc --- /dev/null +++ b/Minimization/Fit_TSA.py @@ -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) + + + diff --git a/Minimization/Fit_TSA_mp.py b/Minimization/Fit_TSA_mp.py new file mode 100755 index 0000000..5a954ea --- /dev/null +++ b/Minimization/Fit_TSA_mp.py @@ -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) + + + diff --git a/Minimization/PLUV.py b/Minimization/PLUV.py new file mode 100644 index 0000000..d56d7dc --- /dev/null +++ b/Minimization/PLUV.py @@ -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 + +################################################################################################################## +################################################################################################################## + + + + + + + diff --git a/Minimization/Tools.py b/Minimization/Tools.py new file mode 100755 index 0000000..95720cd --- /dev/null +++ b/Minimization/Tools.py @@ -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: \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]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 + + + + + +########################################################################################### +########################################################################################### +########################################################################################### +########################################################################################### + + +