The math range annoying error

Coutaz.py

##==============================================================================##

##-----------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