-
-
Save harshadsurdi/3736067 to your computer and use it in GitHub Desktop.
##==============================================================================## | |
##-----------ALL IMPORTS-----------## | |
from pylab import * | |
import numpy as np | |
import matplotlib.pyplot as plt | |
from mpl_toolkits.mplot3d import Axes3D | |
import math | |
from scipy.optimize import * | |
import cmath as cm | |
##import powell #This is my written powell method | |
##from powell import * | |
from decimal import * | |
##==============================================================================## | |
##-----------ALL INITIALIZATIONS-----------## | |
v1=complex(0,0) | |
v2=complex(0,0) | |
fix_values=-1j*2*math.pi*0.01 | |
L=varr_L=400 | |
c=2.99792458 | |
delta_f= 0.001953125 #---delta_w = 4THz/2048 = 1953125000 | |
constant=fix_values*(delta_f*L/c) | |
constant_varr_L=fix_values*(delta_f*varr_L/c) | |
z=complex(0,1.0) | |
n1=complex(1.0,0) | |
n3=complex(1.0,0) | |
global Pwl_omega | |
##==============================================================================## | |
##-----------ALL DEFINITIONS-----------## | |
def twopower(x): | |
y=int(math.log(x,2)) | |
return 2**y | |
#=========================================== | |
#======T formulated========================= | |
#=========================================== | |
try: | |
def Tformulated(n,k,w): | |
v3=w*constant | |
n2=complex(n,-k) | |
v1=1.00/(1.00-(((n2-n1)*(n2-n3)*cm.exp(2.00*v3*n2))/((n2+n1)*(n2+n3)))) | |
FP=v1 | |
v2=((2.0*n2*(n1+n3))/((n2+n1)*(n2+n3)))*((cm.exp((n2-n1)*v3)))*FP | |
Tform=v2 | |
return(Tform) | |
except OverflowError: | |
pass | |
#=========================================== | |
#======T measured with fixed n and k values===================== | |
def Tmeas_fix(n,k): | |
FP=[] | |
Tform=[] | |
for i in range(2048): | |
v3=i*constant | |
n2=complex(n,-k) | |
v1=1.0/(1.0-(((n2-n1)*(n2-n3)*cm.exp(2*v3*n2))/((n2+n1)*(n2+n3)))) | |
FP.append(v1) | |
v2=((2.0*n2*(n1+n3))/((n2+n1)*(n2+n3)))*((cm.exp((n2-n1)*v3)))*FP[i] | |
Tform.append(v2) | |
return(Tform) | |
#=============================================================================== | |
#======T measured with linear n and k values===================== | |
def Tmeas_lin(): | |
FP=[] | |
Tform=[] | |
for i in range(2048): | |
v3=i*constant | |
n=3.5 + 2.0*i/1000;k=1.0 + 2.0*i/1000 | |
n2=complex(n,-k) | |
v1=1.0/(1.0-(((n2-n1)*(n2-n3)*cm.exp(2.0*v3*n2))/((n2+n1)*(n2+n3)))) | |
FP.append(v1) | |
v2=((2.0*n2*(n1+n3))/((n2+n1)*(n2+n3)))*((cm.exp((n2-n1)*v3)))*FP[i] | |
Tform.append(v2) | |
return(Tform) | |
#=============================================================================== | |
##==============================================================================## | |
##==============================================================================## | |
##-----------MAIN ALGORITHM-----------## | |
#============= Generating T(w) for given n2,k2 and range of w ================== | |
n2=4 #the parameters to | |
k2=3 #generate Tmeas() | |
Tw=Tmeas_lin() #Tmeas() is the generation function | |
Tw1=[abs(i) for i in Tw] | |
plt.figure() | |
plt.plot(range(2048),Tw1) | |
plt.title("Formula Generated T(w)") | |
plt.show() | |
#===========ERROR FUNCTIONS===================================================== | |
error_powell=lambda x: (math.log(abs(Tformulated(x[0],x[1],Pwl_omega))) - math.log(abs(Tw[Pwl_omega])))**2 + (math.atan2(Tformulated(x[0],x[1],Pwl_omega).imag,Tformulated(x[0],x[1],Pwl_omega).real) - math.atan2(Tw[Pwl_omega].imag,Tw[Pwl_omega].real))**2 | |
error=lambda x,w: (math.log(abs(Tformulated(x[0],x[1],w))) - math.log(abs(Tw[w])))**2 + (math.atan2(Tformulated(x[0],x[1],w).imag,Tformulated(x[0],x[1],w).real) - math.atan2(Tw[w].imag,Tw[w].real))**2 | |
#-----------------function trials----------------------------------------------- | |
w_form=1024 | |
error_form=lambda x,w: (math.log(abs(Tformulated(x[0],x[1],w_form))) - math.log(abs(Tw[w])))**2 + (math.atan2(Tformulated(x[0],x[1],w_form).imag/Tformulated(x[0],x[1],w_form).real) - math.atan2(Tw[w].imag/Tw[w].real))**2 | |
function = lambda x : -sin(sqrt(x[0]**2+x[1]**2))/sqrt(x[0]**2+x[1]**2) | |
#------------------------------------------------------------------------------- | |
#========================initial guess vector=================================== | |
v0=[5,3] | |
v0_spow=v0 | |
v0_sbfgs=v0 | |
v0_cg=v0 | |
v0_bfgs=v0 | |
v0_nm=v0 | |
v0_nrpow=v0 | |
v0_ncg=v0 | |
omega=[500] | |
nrpow_k2=[] | |
nrpow_n2=[] | |
spow_n2=[] | |
spow_k2=[] | |
sbfgs_n2=[] | |
sbfgs_k2=[] | |
snm_n2=[] | |
snm_k2=[] | |
sfbfgs_n2=[] | |
sfbfgs_k2=[] | |
sfcg_n2=[] | |
sfcg_k2=[] | |
fminncg_k2=[] | |
fminncg_n2=[] | |
#===========================Scipy POWELL======================================== | |
func_type=raw_input("Minimization Method : ") | |
freq_sel=raw_input("Enter Frequency Range : ") | |
a,b = map(int,freq_sel.split()) | |
for w in range(a,b): | |
omega[0]=w | |
if func_type == "ALL": | |
#print "=======================================================" | |
retval=fmin_powell(error,v0_spow,args=(omega),disp=0) | |
#print "fmin_powell of COUTAZ at %s --> %s "%(retval,error(retval,omega[0])) | |
spow_n2.append(retval[0]) | |
spow_k2.append(retval[1]) | |
v0_spow=[retval[0],retval[1]] | |
#===========================Scipy BFGS========================================== | |
#print "=======================================================" | |
retval2 = minimize(error,v0_sbfgs,args=(omega),method='BFGS') | |
#print "minimize BFGS : ",retval2.x | |
#print "VALUE of COUTAZ at (%f,%f) = %s"%(retval2.x[0],retval2.x[1],error([retval2.x[0],retval2.x[1]],omega[0])) | |
sbfgs_n2.append(retval2.x[0]) | |
sbfgs_k2.append(retval2.x[1]) | |
v0_sbfgs=[retval2.x[0],retval2.x[1]] | |
## #===========================Scipy fmin_BFGS========================================== | |
## #print "=======================================================" | |
retval2 = fmin_bfgs(error,v0_bfgs,args=(omega),disp=0,gtol=1e-5) | |
#print "minimize BFGS : ",retval2 | |
#print "VALUE of COUTAZ at (%f,%f) = %s"%(retval2.x[0],retval2.x[1],error([retval2.x[0],retval2.x[1]],omega[0])) | |
sfbfgs_n2.append(retval2[0]) | |
sfbfgs_k2.append(retval2[1]) | |
v0_bfgs=[retval2[0],retval2[1]] | |
## #===========================Scipy fmin_cg========================================== | |
## print "=======================================================" | |
retval2 = fmin_cg(error,v0,args=(omega),disp=0) | |
#print "minimize BFGS : ",retval2 | |
#print "VALUE of COUTAZ at (%f,%f) = %s"%(retval2.x[0],retval2.x[1],error([retval2.x[0],retval2.x[1]],omega[0])) | |
sfcg_n2.append(retval2[0]) | |
sfcg_k2.append(retval2[1]) | |
v0_cg=[retval2[0],retval2[1]] | |
## #===========================Scipy Nelder-Mead=================================== | |
#print "=======================================================" | |
retval2 = minimize(error,v0_nm,args=(omega),method='Nelder-Mead') | |
#print "minimize Nelder-Mead : ",retval2.x | |
#print "VALUE of COUTAZ at (%f,%f) = %s"%(retval2.x[0],retval2.x[1],error([retval2.x[0],retval2.x[1]],omega[0])) | |
snm_n2.append(retval2.x[0]) | |
snm_k2.append(retval2.x[1]) | |
v0_nm=[retval2.x[0],retval2.x[1]] | |
## #===========================Numerical Reciepes POWELL=========================== | |
## #print "=======================================================" | |
##Pwl_omega = omega[0] | |
##retval,nIter = powell(error_powell,array([v0_nrpow[0],v0_nrpow[1]])) | |
##print "NR_powell of COUTAZ at %s --> %s - "%(xMin,float(error_powell(xMin))) | |
##nrpow_n2.append(retval[0]) | |
## nrpow_k2.append(retval[1]) | |
## v0_nrpow=[retval[0],retval[1]] | |
###===========================fmin Newton CG=========================== | |
retval2 = fmin_ncg(error,v0_ncg,args=(omega)) | |
#print "minimize Nelder-Mead : ",retval2.x | |
#print "VALUE of COUTAZ at (%f,%f) = %s"%(retval2.x[0],retval2.x[1],error([retval2.x[0],retval2.x[1]],omega[0])) | |
fminncg_n2.append(retval2[0]) | |
fminncg_k2.append(retval2[1]) | |
v0_ncg=[retval2[0],retval2[1]] | |
## #print "=======================================================" | |
elif func_type == "Powell": | |
retval=fmin_powell(error,v0_spow,args=(omega),disp=0) | |
#print "fmin_powell of COUTAZ at %s --> %s "%(retval,error(retval,omega[0])) | |
spow_n2.append(retval[0]) | |
spow_k2.append(retval[1]) | |
v0_spow=[retval[0],retval[1]] | |
elif func_type == "mini BFGS": | |
retval2 = minimize(error,v0_sbfgs,args=(omega),method='BFGS') | |
#print "minimize BFGS : ",retval2.x | |
#print "VALUE of COUTAZ at (%f,%f) = %s"%(retval2.x[0],retval2.x[1],error([retval2.x[0],retval2.x[1]],omega[0])) | |
sbfgs_n2.append(retval2.x[0]) | |
sbfgs_k2.append(retval2.x[1]) | |
v0_sbfgs=[retval2.x[0],retval2.x[1]] | |
elif func_type == "fmin BFGS": | |
retval2 = fmin_bfgs(error,v0_bfgs,args=(omega),disp=0,gtol=1e-5) | |
#print "minimize BFGS : ",retval2 | |
#print "VALUE of COUTAZ at (%f,%f) = %s"%(retval2.x[0],retval2.x[1],error([retval2.x[0],retval2.x[1]],omega[0])) | |
sfbfgs_n2.append(retval2[0]) | |
sfbfgs_k2.append(retval2[1]) | |
v0_bfgs=[retval2[0],retval2[1]] | |
elif func_type == "mini Nelder_Mead": | |
retval2 = minimize(error,v0_nm,args=(omega),method='Nelder-Mead') | |
#print "minimize Nelder-Mead : ",retval2.x | |
#print "VALUE of COUTAZ at (%f,%f) = %s"%(retval2.x[0],retval2.x[1],error([retval2.x[0],retval2.x[1]],omega[0])) | |
snm_n2.append(retval2.x[0]) | |
snm_k2.append(retval2.x[1]) | |
v0_nm=[retval2.x[0],retval2.x[1]] | |
##elif func_type == "NR Powell": | |
## Pwl_omega = omega[0] | |
## retval,nIter = powell(error_powell,array([v0_nrpow[0],v0_nrpow[1]])) | |
## #print "NR_powell of COUTAZ at %s --> %s - "%(xMin,float(error_powell(xMin))) | |
## nrpow_n2.append(retval[0]) | |
## nrpow_k2.append(retval[1]) | |
## v0_nrpow=[retval[0],retval[1]] | |
#print "\n=======================================================" #blue - ; green- ; red : ; light blue : spowell | |
n=[] | |
k=[] | |
for i in range(w+1): | |
temp=3.5 + 2.0*i/1000 | |
temp1=1 + 2.0*i/1000 | |
n.append(temp) | |
k.append(temp1) | |
plt.figure(); | |
if func_type == "ALL": | |
p1,=plt.plot(range(len(spow_k2)),spow_k2) | |
p2,=plt.plot(range(len(sbfgs_k2)),sbfgs_k2) | |
p3,=plt.plot(range(len(sfbfgs_k2)),sfbfgs_k2) | |
p4,=plt.plot(range(len(snm_k2)),snm_k2) | |
p5,=plt.plot(range(len(nrpow_k2)),nrpow_k2) | |
p6,=plt.plot(range(w+1),k,'g+') | |
legend([p1,p2,p3,p4,p5,p6], ["spow", "sbfgs","snm","nrpow","sfbfgs","n2"]) | |
plt.title("K2") | |
plt.show() | |
elif func_type=="Powell": | |
p1,=plt.plot(range(len(spow_k2)),spow_k2) | |
p6,=plt.plot(range(w+1),k,'g+') | |
legend([p1,p6], ["Powell_k2","k2"]) | |
plt.title("K2") | |
plt.show() | |
elif func_type=="mini BFGS": | |
p2,=plt.plot(range(len(sbfgs_k2)),sbfgs_k2) | |
p6,=plt.plot(range(w+1),k,'g+') | |
legend([p2,p6], ["mini BFGS_k2","k2"]) | |
plt.title("K2") | |
plt.show() | |
elif func_type=="fmin BFGS": | |
p3,=plt.plot(range(len(sfbfgs_k2)),sfbfgs_k2) | |
p6,=plt.plot(range(w+1),k,'g+') | |
legend([p3,p6], ["fmin BFGS_k2","k2"]) | |
plt.title("K2") | |
plt.show() | |
elif func_type=="mini Nelder_Mead": | |
p4,=plt.plot(range(len(snm_k2)),snm_k2) | |
p6,=plt.plot(range(w+1),k,'g+') | |
legend([p4,p6], ["mini Nelder Mead","k2"]) | |
plt.title("K2") | |
plt.show() | |
##elif func_type=="NR Powell": | |
## p5,=plt.plot(range(len(nrpow_k2)),nrpow_k2) | |
##p6,=plt.plot(range(w+1),k,'g+') | |
##legend([p5,p6], ["NR Powell","k2"]) | |
##plt.title("K2") | |
##plt.show() | |
plt.figure(); | |
if func_type == "ALL": | |
p1,=plt.plot(range(len(spow_k2)),spow_k2) | |
p2,=plt.plot(range(len(sbfgs_k2)),sbfgs_k2) | |
p3,=plt.plot(range(len(sfbfgs_k2)),sfbfgs_k2) | |
p4,=plt.plot(range(len(snm_k2)),snm_k2) | |
p5,=plt.plot(range(len(nrpow_k2)),nrpow_k2) | |
p6,=plt.plot(range(w+1),k,'g+') | |
legend([p1,p2,p3,p4,p5,p6], ["spow", "sbfgs","snm","nrpow","sfbfgs","n2"]) | |
plt.title("K2") | |
plt.show() | |
elif func_type=="Powell": | |
p1,=plt.plot(range(len(spow_n2)),spow_n2) | |
p6,=plt.plot(range(w+1),n,'g+') | |
legend([p1,p6], ["Powell_k2","k2"]) | |
plt.title("K2") | |
plt.show() | |
elif func_type=="mini BFGS": | |
p2,=plt.plot(range(len(sbfgs_n2)),sbfgs_n2) | |
p6,=plt.plot(range(w+1),n,'g+') | |
legend([p2,p6], ["mini BFGS_k2","k2"]) | |
plt.title("K2") | |
plt.show() | |
elif func_type=="fmin BFGS": | |
p3,=plt.plot(range(len(sfbfgs_n2)),sfbfgs_n2) | |
p6,=plt.plot(range(w+1),n,'g+') | |
legend([p3,p6], ["fmin BFGS_k2","k2"]) | |
plt.title("K2") | |
plt.show() | |
elif func_type=="mini Nelder_Mead": | |
p4,=plt.plot(range(len(snm_n2)),snm_n2) | |
p6,=plt.plot(range(w+1),n,'g+') | |
legend([p4,p6], ["mini Nelder Mead","k2"]) | |
plt.title("K2") | |
plt.show() | |
##elif func_type=="NR Powell": | |
## p5,=plt.plot(range(len(nrpow_n2)),nrpow_n2) | |
## p6,=plt.plot(range(w+1),n,'g+') | |
## legend([p5,p6], ["NR Powell","k2"]) | |
## plt.title("K2") | |
## plt.show() |
ok, here's how you do the program :
first it will display a graph, then it will ask for the minimization method , enter either of the following : "ALL" or "Powell" or "mini BFGS" or "fmin BFGS" or " mini Nelder-Mead" (ofcourse w/o the inverted commas). Then it will ask for the frequency range, enter the following : fmin fmax [ fmin can be minimum 0 and fmax can be within 2048]. for example - "Enter Frequency Range : " you enter "0 2048"(ofcourse w/o the inverted commas)
For "ALL" it will perform all the minimization methods and for the rest accordingly.
Where does the 'minimize' function in line 155 come from?
It's a optimization routine in the scipy.optimize module
http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
Sorry, line 155.