Created
September 17, 2012 07:46
-
-
Save harshadsurdi/3736067 to your computer and use it in GitHub Desktop.
The math range annoying error
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
##==============================================================================## | |
##-----------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 a optimization routine in the scipy.optimize module
http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Where does the 'minimize' function in line 155 come from?