tsa_saxs/Minimization/Tools.py
2024-03-05 11:13:24 +01:00

587 lines
18 KiB
Python
Executable file

#!/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( r"^DATA_FILE:$", self.OPT[i][0] ):
options['datafile'] = str(self.OPT[i][1])
file_exists(options['datafile'], "Data File")
elif re.match( r"^TITLE:$", self.OPT[i][0] ):
options['title'] = str(self.OPT[i][1])
print("Title:\t\t", options['title'])
elif re.match( r"^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( r"^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( r"^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( r"^MODEL:$", self.OPT[i][0] ):
options['function'] = str(self.OPT[i][1])
print("Function:\t", options['function'])
elif re.match( r"^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( r"^BINNING:$", self.OPT[i][0] ):
options['bing'] = check_int(self.OPT[i][1], "Binning")
print("Binning:\t", options['bing'])
elif re.match( r"^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( r"^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( r"^#" , 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( r"^#" , 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. Format:
1st column: q value
2nd column: I value
3rd column: error
"""
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
###########################################################################################
###########################################################################################
###########################################################################################
###########################################################################################