diff --git a/Minimization/FitTools.py b/Minimization/FitTools.py index 04b656b..d34cc5b 100755 --- a/Minimization/FitTools.py +++ b/Minimization/FitTools.py @@ -8,7 +8,7 @@ import math sys.path.append("/mntdirect/_users/semeraro/python_tools") import FormFactor -import PLUV_POPC_RecBuf +#import PLUV_POPC_RecBuf from PLUV import SDP_base_POPC_RecBuf, SDP_POPC_RecBuf ########################################################################################### diff --git a/Minimization/FormFactor.py b/Minimization/FormFactor.py new file mode 100644 index 0000000..c6d9ba8 --- /dev/null +++ b/Minimization/FormFactor.py @@ -0,0 +1,71 @@ +#!/usr/bin/python + +from scipy.special import j1 , hyp1f1 +from scipy.integrate import simps +from scipy.stats import norm , gamma +import numpy as np +import sys +import math +pi=math.pi +e=math.e + +sys.path.append("/mntdirect/_users/semeraro/python_tools") +import FitTools + +N_OrAv = 90*2 +N_CapCyl = 50 + +########################################################################################### +########################################################################################### + +""" SPHERICAL BESSEL FUNCTION j0 & j1 """ + +def spherical_j0_np(x): + return np.where(x==0, 1, np.sin(x)/x ) + +def spherical_j1_np(x): + return np.where(x==0, 1./3, (np.sin(x)-x*np.cos(x))/x**2 ) + +########################################################################################### +########################################################################################### +""" NORMAL DISTRIBUTION """ +def Normal(x,xo,w): + return np.exp(- (x-xo)**2/(2*w**2) ) / (w*np.sqrt(2*pi)) + +########################################################################################### +########################################################################################### +""" FORM FACTORS """ + +################################ SPHERE +def Amp_Sphere(q,R): + if R==0 or q==0: + A=Vol_Sphere(R) + else: + A=Vol_Sphere(R) * 3 * spherical_j1(q*R) / (q*R) + return A + +################################ SPHERES +def Sphere( q, PAR ): + [n,rho,R,C] = PAR + + amp = (4*np.pi*R**3/3) * 3 * spherical_j1_np(q*R) / (q*R) + + return n * (rho * amp)**2 + C + +################################ POLYDISPERSE SPHERES (NORMAl PDF) +def Sphere_Normal (q,PAR): + [n,rho,R,w,C] = PAR + + N = 21 + R_array = np.linspace(R-3*w,R+3*w,N) + Normal = norm.pdf(R_array, R, w) + + I = np.zeros(q.shape[0], dtype=float) + + for r in range(R_array.shape[0]): + if r==0 or r==N-1 : I+= Sphere(q,[1,1,R_array[r],0]) * Normal[r] / 2 + else : I+= Sphere(q,[1,1,R_array[r],0]) * Normal[r] + I*= 6*w/(N-1) + + return n*rho**2*I + C + diff --git a/Minimization/TSA_algorithm.py b/Minimization/TSA_algorithm.py new file mode 100644 index 0000000..d1ba384 --- /dev/null +++ b/Minimization/TSA_algorithm.py @@ -0,0 +1,297 @@ +#!/usr/bin/python + +from __future__ import print_function +import sys +import os.path +import re +import math +import time +import random +import copy +import array +import numpy as np +from numpy.linalg import inv +from math import sqrt +import matplotlib.pyplot as plt +from scipy.optimize import curve_fit + +sys.path.append("/mntdirect/_users/semeraro/python_tools") +import FitTools +import FormFactor +#import Decoupling + +########################################################################################### + +def X2function(Q,IDATA,I,ERR,N): + 'X^2 calculation' + X2 = np.sum( ( ( IDATA-I ) / ERR )**2 ) + return X2/(Q.shape[0]-(N-1)) + +def exponential_decay(x, a, st, tau_iter, Cave) : + return a * np.exp( -(x-st)/tau_iter ) + Cave + +########################################################################################### +########################################################################################### +def Min_likelyhood (prior_ratio) : + min = 1 + if prior_ratio < 1 : + min = prior_ratio + return min + +#------------------------- Partial ratio between priors +def PriorRatio (start, check, rel, prev) : + new = np.exp(-(start-check)**2/(2*(rel*start)**2)) + old = np.exp(-(start-prev)**2/(2*(rel*start)**2)) + return new/old + +########################################################################################### +########################################################################################### +########################################################################################### +########################### SIMULATED ANNEALING ALGORITHM ########################### +########################################################################################### +########################################################################################### +########################################################################################### +#def SimAnnealing ( Q, IDATA, ERR, PAR, FIX, PRIOR, L_LIM, H_LIM, function, temperature, save ): +def SimAnnealing ( input_list ): + 'Simulated Annealing Algorithm' + + Q, IDATA, ERR, PAR, FIX, PRIOR, L_LIM_org, H_LIM_org, function, temperature, save = input_list + + print("START") + ##################################################### + #------------------------ Define Controls + good = 0 + tau = 5 # tau>1 fast / tau<1 slow + maxgood = 12500 + maxcount = 50000 + thermo = 10 # for POPC LUVs + #thermo = 500 + T = temperature + + FREE_MIN = [] + + prBol = 0 + + alpha = 0.15 + alphamax = 1.0 + alphamin = 0.024 + + cumDEnt = 0 + cumDX2 = 0 + dly = 0 + + negWater_check = 0 + + stop = "CONTINUE" + + ##################################################### + #------------------------ Define Free Parameters Array + (FREE,L_LIM,H_LIM) = FitTools.FreePar(PAR,FIX,L_LIM_org,H_LIM_org) + START = FREE.copy() + + ##################################################### + #------------------------ Merge Parameters Arrays + CALC = FitTools.MergePar(PAR,FREE,FIX) + + ##################################################### + #------------------------ Evaluate I(q) & Calculate new CHI squared + + #I = function(Q,CALC) + intensity_init = function(Q,CALC) + I = intensity_init.intensity() + X2 = X2function(Q,IDATA,I,ERR,len(FREE)) + X2_min = X2 + + ##################################################### + #------------------------ Set priors index and sigma + pr_idx = [] + pr_rel = [] + + p = 0 + f = 0 + while p < len(PRIOR) : + if FIX[p] == 'f' : f+= 1 + if PRIOR[p] != '-': + pr_idx.append(p-f) + pr_rel.append(float(PRIOR[p])) + p+=1 + + #print(pr_idx) + #print(pr_rel) + + ##################################################### + #------------------------ Start Loop + loop=0 + t0=time.time() + Nm = 500 + + trace= [[] for i in range(3+len(PAR))] + + while stop == "CONTINUE": + + ##################################################### + #------------------------ Calculate Parameter Perturbations + PERT=[] + for i in range(len(FREE)): + PERT.append( np.random.normal( 0, alpha * ( H_LIM[i] - L_LIM[i] )/6. ) ) + + ##################################################### + #------------------------ Calculate temporary parameters + #------------------------ Check for FIX parameters & boundarys + CHECK=[] + for i in range(len(FREE)): + if (FREE[i]+PERT[i])>H_LIM[i] : + CHECK.append(FREE[i]+(H_LIM[i]-FREE[i])/10) + elif (FREE[i]+PERT[i])= 0 : + T = temperature + else : + T = thermo * cumDX2 / cumDEnt + if T < 1 : + T = 1 + + if (good+1)/(loop+1) < 0.6 : alpha-= alpha*0.1 + elif (good+1)/(loop+1) > 0.7 : alpha+= alpha*0.1 + + if alpha < alphamin : alpha = alphamin + elif alpha > alphamax : alpha = alphamax + + ##################################################### + #------------------------ Update maxgood + if loop % 1000 == 0 and loop < 9001 : + + try: + popt, pcov = curve_fit( exponential_decay, np.array(trace[0]), np.array(trace[2]) ) + a_fit, st_fit, tau_iter_fit, Cave_fit = popt + + except RuntimeError: + #print("\nError - curve_fit failed\n") + tau_iter_fit = 500 + + #fig1, (figA) = plt.subplots(1, 1, figsize=(7.0, 7.0)) + #plt.scatter(trace[0],trace[2], marker='o', facecolors='none', color='b') + #plt.plot(np.array(trace[0]),exponential_decay(np.array(trace[0]),a_fit, st_fit, tau_iter_fit, Cave_fit), '-', color='r', linewidth=2) + #fig1.tight_layout() + #plt.xscale('log') + #plt.yscale('log') + #plt.show() + + if tau_iter_fit > 500 : maxgood = 12500 + elif tau_iter_fit*5 < 250 : maxgood = 10250 + else : maxgood = int(tau_iter_fit*5) + 10000 + + ##################################################### + #------------------------ Time counting + t1=time.time()-t0 + (minutes,seconds,cents)=FitTools.TimeCount(t1) + #success = (good+1)/(loop+1) + + #if loop/1000 == 1: dly = (time.time()-t0)/1001 + #elif loop/3000 == 1: dly = (time.time()-t0)/3001 + #elif loop/6000 == 1: dly = (time.time()-t0)/6001 + #elif loop/9000 == 1: dly = (time.time()-t0)/9001 + #(minutes_ex,seconds_ex,cents_ex)=FitTools.TimeCount( dly * maxgood/success )#*maxgood/((good+1)/(loop+1))-t0) + + ##################################################### + #------------------------ Print Loop Results + #if T>0.9 and loop%10 == 0 : + # out="T = %0.3f \ a = %0.3f \ min.X²= %0.5s \ N. %0.3d/%0.3d (%0.3f) -> %.3d \ %2d:%2d s (exp. %2d:%2d s) \ " %(T,alpha,X2_min,good,loop,success,maxgood,minutes,seconds,minutes_ex,seconds_ex) + #print (out, end='\r') + + ##################################################### + #------------------------ Set the Loop Conditions + if good == maxgood+1 or loop == maxcount+1 : + stop="STOP" + print("STOP") + + ##################################################### + #------------------------ End Loop + #with open("./"+save+"/Trace.dat", 'w') as trc: + # for i in range(len(trace[p])) : + # for p in range(len(trace)) : + # trc.writelines( "%.7f\t" %(trace[p][i]) ) + # trc.writelines( "\n" ) + + + return [ FitTools.MergePar(PAR,FREE_MIN,FIX), X2_min] + + +########################################################################################### +########################################################################################### + + + + +