Skip to content

Instantly share code, notes, and snippets.

@harshadsurdi
Created September 17, 2012 07:46
Show Gist options
  • Save harshadsurdi/3736067 to your computer and use it in GitHub Desktop.
Save harshadsurdi/3736067 to your computer and use it in GitHub Desktop.
The math range annoying error
##==============================================================================##
##-----------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()
@jaidevd
Copy link

jaidevd commented Sep 17, 2012

Sorry, line 155.

@harshadsurdi
Copy link
Author

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.

@jaidevd
Copy link

jaidevd commented Sep 23, 2012

Where does the 'minimize' function in line 155 come from?

@harshadsurdi
Copy link
Author

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