Here is a test of another parameterization, a simple piecewise function.
In [1]:
# Routine
# This line configures matplotlib to show figures embedded in the notebook,
# instead of opening a new window for each figure. More about that later.
# If you are using an old version of IPython, try using '%pylab inline' instead.
%matplotlib inline
%load_ext snakeviz
import numpy as np
from scipy.optimize import minimize
from scipy.optimize import rosen, differential_evolution
from scipy.special import expit
import matplotlib.pyplot as plt
import scipy
from matplotlib.lines import Line2D
import timeit
import pandas as pd
import plotly.plotly as py
from plotly.graph_objs import *
import plotly.tools as tls
In [2]:
# For test use. This is some random structure.
array8 = np.asarray(np.split(np.random.rand(1,40)[0],8))
array8zero = np.asarray(np.split(np.zeros(40),8))
In [20]:
# Density matrix in the forms that I wrote down on my Neutrino Physics notebook
# x is a real array of 8 lenghth mmax arrays. x[i] is for A_m\sin(m x) while x[i+4] is for B_m \cos(m x)
init = np.array([1.0,0.0,0.0,0.0])
def rho(x,ti):
elem = np.zeros(4)
for i in np.linspace(0,3,4):
elem[i] = np.sum( x[i]*np.sin(x[i]*ti) + x[i+4]*np.cos(x[i+4]*ti) )
return elem
In [4]:
rho(array8,0.3)
Out[4]:
In [39]:
def rhop(x,ti,length):
rhoprime = np.zeros(4)
# mlen = len(x[0])
mlen = length
marray = np.linspace(0,mlen-1,mlen)
for i in np.linspace(0,3,4):
rhoprime[i] = np.sum( ( x[i]*np.cos(x[i]*ti) - x[i+4]*np.sin(x[i+4]*ti) )*marray )
return rhoprime
In [40]:
rhop(array8,0.3,5)
Out[40]:
In [17]:
tri = np.array([np.cos(2.0),np.sin(2.0)])
In [63]:
def costi0(x,ti,initialCondition):
costTemp = np.zeros(4)
for i in np.linspace(0,3,4):
costTemp[i] = (np.sum( x[i+4] ) - initialCondition[i])**2
return np.sum(costTemp)
def costi(x,ti,initialCondition,length):
rhoi = rho(x,ti)
rhopi = rhop(x,ti,length)
costTemp = np.zeros(4)
costTemp[0] = ( rhopi[0] - 2.0*rhoi[3]*tri[1] )**2
costTemp[1] = ( rhopi[1] + 2.0*rhoi[3]*tri[1] )**2
costTemp[2] = ( rhopi[2] - 2.0*rhoi[3]*tri[0] )**2
costTemp[3] = ( rhopi[3] + 2.0*rhoi[2]*tri[0] - tri[1] * (rhoi[1] - rhoi[0] ) )**2
return np.sum(costTemp) + costi0(x,ti,initialCondition)
In [64]:
costi(array8,0,init,5)
Out[64]:
In [65]:
def cost(x,t,initialCondition,length):
costTotal = map(lambda t: costi(x,t,initialCondition,length),t)
return 0.5*np.sum(costTotal)
In [66]:
cost(array8,np.array([0,1,2]),init,5)
Out[66]:
In [67]:
# initial guess
initGuess = np.array(np.split(np.random.rand(1,40)[0],8))
endpoint = 2
tlin = np.linspace(0,endpoint,11)
lengthx = len(initGuess[0])
costF = lambda x: cost(x,tlin,init,lengthx)
start = timeit.default_timer()
costvFResult = minimize(costF,initGuess,method="SLSQP",tol=1e-10)
stop = timeit.default_timer()
print stop - start
print costvFResult
In [68]:
xmid = costvFResult.x
start = timeit.default_timer()
#costvFResult = minimize(costF,xmid,method="SLSQP",tol=1e-30,options={"ftol":1e-30,"maxiter":100000})
costvFResult = minimize(costF,xmid,method='Nelder-Mead',tol=1e-15,options={"ftol":1e-15, "maxfev": 1000000,"maxiter":1000000})
stop = timeit.default_timer()
print stop - start
print costvFResult
In [69]:
plttlin=np.linspace(0,endpoint,100)
xresult = costvFResult.x
pltdata11 = np.array([])
for i in plttlin:
pltdata11 = np.append(pltdata11,rho(xresult,i)[0] )
In [70]:
plt.figure(figsize = (16,9.36))
plt.plot(np.linspace(0,2,10),1-(np.sin(2.0)**2)*(np.sin(0.5*np.linspace(0,2,10)) )**2,"r-")
plt.plot(plttlin,pltdata11)
plt.show()
In [ ]: