#!/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)