209 lines
6.1 KiB
Python
209 lines
6.1 KiB
Python
|
#!/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)
|
||
|
|
|
||
|
|
|
||
|
|
|