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.
array12 = np.asarray(np.split(np.random.rand(1,60)[0],12))
In [3]:
# Here is the activation function
# Using sigmoid
def act(x):
x = np.asarray(x, dtype=float)
# np.piecewise complexity
if not x.shape:
x = np.asarray([x])
return np.piecewise(x, [x<=-1.0, (x>-1.0)&(x<1.0), x>=1.0], [0.0, lambda x: x+1.0, 2.0])
def actp(x):
x = np.asarray(x, dtype=float)
# np.piecewise complexity
if not x.shape:
x = np.asarray([x])
return np.piecewise(x, [x<=-1.0,(x>-1.0)&(x<1.0), x>=1.0], [0.0, 1.0, 0.0])
Here comes the density matrix
In [4]:
# Density matrix in the forms that I wrote down on my Neutrino Physics notebook
# x is a real array of 12 arrays.
init = np.array([1.0,0.0,0.0,0.0])
def rho(x,ti,initialCondition):
elem = np.ones(4)
for i in np.linspace(0,3,4):
elem[i] = np.sum(ti*x[i*3]*act(ti*x[i*3+1] + x[i*3+2]) )
return init + elem
In [5]:
xtest = np.linspace(-5, 5, 10)
#np.piecewise(xtest, [xtest < 0, 0 <= xtest , xtest< 10], [lambda x: -x, 0,lambda x: x])
#[act(i) for i in xtest]
print act(xtest), actp(xtest),act(np.linspace(-5, 5, 10))
In [6]:
# Hamiltonian of the problem, in terms of four real components
hamil = np.array( [ np.cos(2.0),np.sin(2.0) ] )
print hamil
In [7]:
def rhop(x,ti,initialCondition):
rhoprime = np.zeros(4)
for i in np.linspace(0,3,4):
rhoprime[i] = np.sum(x[i*3] * (act(ti*x[i*3+1] + x[i*3+2]) ) ) + np.sum( ti*x[i*3]* (actp(ti*x[i*3+1] + x[i*3+2]) ) * x[i*3+1] )
return rhoprime
def costi(x,ti,initialCondition):
rhoi = rho(x,ti,initialCondition)
rhopi = rhop(x,ti,initialCondition)
costTemp = np.zeros(4)
costTemp[0] = ( rhopi[0] - 2.0*rhoi[3]*hamil[1] )**2
costTemp[1] = ( rhopi[1] + 2.0*rhoi[3]*hamil[1] )**2
costTemp[2] = ( rhopi[2] - 2.0*rhoi[3]*hamil[0] )**2
costTemp[3] = ( rhopi[3] + 2.0*rhoi[2]*hamil[0] - hamil[1] * (rhoi[1] - rhoi[0] ) )**2
return np.sum(costTemp)# + 2.0*(rhoi[0]+rhoi[1]-1.0)**2
In [8]:
costi(array12,0,init)
Out[8]:
In [9]:
def cost(x,t,initialCondition):
costTotal = map(lambda t: costi(x,t,initialCondition),t)
return 0.5*np.sum(costTotal)
In [10]:
cost(array12,np.array([0,1,2]),init)
Out[10]:
In [ ]:
In [ ]:
# with ramdom initial guess. TO make it more viable, I used (-5,5)
#initGuess = np.asarray(np.split(10.0*(np.random.rand(1,60)[0] - 0.5),12))
initGuess = np.split(np.zeros(60),12)
endpoint = 2
tlin = np.linspace(0,endpoint,11)
def costF(x):
return cost(x,tlin,init)
start = timeit.default_timer()
costvFResult = minimize(costF,initGuess,method="Nelder-Mead",tol=1e-20,options={"maxiter":100000,"maxfev":1000000})
stop = timeit.default_timer()
print stop - start
print costvFResult
In [ ]:
xresult = costvFResult.x
In [ ]:
plttlin=np.linspace(0,endpoint,100)
pltdata11 = np.array([])
for i in plttlin:
pltdata11 = np.append(pltdata11 ,rho(xresult,i,init)[0] )
print pltdata11
plt.figure(figsize=(16,9.36))
plt.ylabel('P')
plt.xlabel('Time')
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,"b4-",label="vac_rho11")
plt.show()
#py.iplot_mpl(plt.gcf(),filename="MMA-rho11-Vac-80-60")
In [ ]:
In [19]:
# with ramdom initial guess
# %snakeviz
devoendpoint = 3
devotlin = np.linspace(0,devoendpoint,11)
devocostF = lambda x: cost(x,devotlin,init)
bounds=np.zeros([60,2])
for i in range(60):
bounds[i,0]=-1.0
bounds[i,1]=1.0
#print bounds
startdevo = timeit.default_timer()
devo = differential_evolution(devocostF,bounds,strategy='best1bin',tol=1e-10,maxiter=10000,polish=True)
stopdevo = timeit.default_timer()
print stopdevo - startdevo
print devo
In [ ]:
devoxresult = devo.x
devoplttlin=np.linspace(0,devoendpoint,100)
devopltdata11 = np.array([])
for i in devoplttlin:
devopltdata11 = np.append(devopltdata11 ,rho(devoxresult,i,init)[0] )
print devopltdata11
plt.figure(figsize=(16,9.36))
plt.ylabel('P')
plt.xlabel('Time')
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(devoplttlin,devopltdata11,"b4-",label="vac_rho11")
plt.show()
#py.iplot_mpl(plt.gcf(),filename="MMA-rho11-Vac-80-60")
In [ ]: