-
-
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() |
It's showing a plot and then asking for minimization method and frequency range. The 'minimize' method in line 255 is not defined.
Sorry, line 155.
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
The powell module is my written module; but now I have commented it and the parts related to that module