# 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
# Here are the initial values
# For test use
array12 = np.asarray(np.split(np.random.rand(1,60)[0],12))
# Here is the activation function
def act(x):
return expit(x)
# 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
rho(array12,0,init)
array([ 1., 0., 0., 0.])
# Hamiltonian of the problem, in terms of four real components
hamil = np.array( [ np.cos(2.0),np.sin(2.0) , np.sin(2.0),np.cos(2.0) ] )
#hamil = 1.0/2*np.array( [ -np.cos(2.0),np.sin(2.0) , np.sin(2.0),np.cos(2.0) ] )
print hamil
[-0.41614684 0.90929743 0.90929743 -0.41614684]
# Cost function for each time step
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]* (act(ti*x[i*3+1] + x[i*3+2]) ) * (1.0 - (act(ti*x[i*3+1] + x[i*3+2]) ) )* x[i*3+1] )
return rhoprime
## This is the regularization
regularization = 0.0001
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
# return np.sum(costTemp) + regularization*np.sum(x**2)
costi(array12,0,init)
15.207712841852992
def cost(x,t,initialCondition):
costTotal = map(lambda t: costi(x,t,initialCondition),t)
return 0.5*np.sum(costTotal)
cost(array12,np.array([0,1,2]),init)
#cost(xresult,np.array([0,4,11]),init)
91.171908432261404
# with ramdom initial guess. TO make it more viable, I used (-5,5)
initGuess = np.asarray(np.split( 5.0*(np.random.rand(1,60)[0] - 0.5),12))
#initGuess = np.split(np.zeros(60),12)
endpoint = 4
tlin = np.linspace(0,endpoint,11)
costF = lambda x: cost(x,tlin,init)
start = timeit.default_timer()
costvFResult = minimize(costF,initGuess,method="SLSQP",tol=1e-20)
stop = timeit.default_timer()
print stop - start
print costvFResult
29.0142710209 status: 9 success: False njev: 101 nfev: 6309 fun: 0.54986503002143661 x: array([-8.92058678, -0.13492851, -3.10941404, 2.19527238, -0.15759464, -1.55261971, 1.60513912, -0.31253223, -1.75685996, -4.0134798 , -1.44660797, -1.40727698, 2.34579259, -1.4181895 , -1.20672898, 1.63251151, -2.1915661 , -2.01251504, 1.68690961, 1.27759708, 0.44283011, -1.56722276, 1.23812704, -0.69068384, -1.90243599, -0.53131647, 0.53973357, 1.07152732, 1.03052607, 2.07195604, 1.00096812, -2.20206403, 2.41427559, -0.25938366, 0.59837747, -2.39353717, -0.59077064, -0.41344543, -2.0255061 , -0.87488353, 0.42200697, 1.71708869, -0.13074538, -0.28937053, -0.21290055, -0.82239556, 1.83856369, -1.42455788, 0.25667742, -2.35525549, -1.97730218, -0.11343628, 0.78372443, 0.89597664, 1.06860755, 1.85481695, 1.36804639, -2.16455287, -0.42619305, -1.76387092]) message: 'Iteration limit exceeded' jac: array([ 0.00025084, -0.00299339, -0.00226416, -0.00105266, -0.00668124, -0.0025664 , -0.0012297 , -0.00267243, -0.00213507, 0.00186877, -0.00875118, -0.00752775, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ]) nit: 101
xmid = costvFResult.get("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": 10000000,"maxiter":1000000})
stop = timeit.default_timer()
print stop - start
print costvFResult
xresult = costvFResult.get("x")
#xresult = np.array([-0.01486401, 2.25493868, -1.84543911, -0.07335087, 2.2548026 ,
# -1.84421662, 0.10454078, -1.64040244, 2.99923129, 0.01486399,
# 2.2549433 , -1.84544286, 0.61994826, 0.96756294, 0.60929092,
# 0.23839364, 0.45364599, 0.18426882, 0.30242425, 0.96719724,
# 0.13016584, 0.9192801 , 0.116001 , 0.46777053, 0.17497595,
# 0.96035958, 0.21763616, 0.73997804, 0.88071662, 0.1620245 ,
# 0.66904538, 0.66084959, 0.89772078, 0.49020208, 0.67802378,
# 0.53307714, 0.59867975, 0.16864478, 0.4257949 , 0.5364126 ,
# 0.78476644, 0.4910997 , 0.834945 , 0.45061367, 0.16736545,
# 0.42579168, 0.16877594, 0.98282177, 0.08852038, 0.12633737,
# 0.50922379, 0.93146299, 0.66505978, 0.33157336, 0.05408186,
# 0.04504323, 0.27311737, 0.27651656, 0.47313653, 0.12806564])
#xresult = np.array([-1.37886409, 2.81454922, -0.3571002 , 0.02582831, -1.05414931,
# -1.52308153, -2.24747468, 0.33947049, -0.32310112, -1.43887103,
# 0.81176258, 0.05139705, -1.02669705, -0.97236805, -0.27536667,
# 0.34860447, -1.06962772, 0.89978175, 2.39662887, -1.45165477,
# -1.54636469, -2.79921374, -1.30335793, -0.62844367, -3.04440811,
# -2.74566393, -2.16222918, -1.60535643, -0.77298204, 0.13848754,
# -0.36544212, 1.23901581, -0.80586367, -0.30212561, -1.02818302,
# -2.82928373, -0.80776632, -2.90056107, -2.42432246, -2.87572658,
# -0.8645904 , -0.59526987, -1.87029203, -1.60957508, -1.83106839,
# 1.07020356, -0.84892132, -0.97053555, -0.2005098 , -0.72422578,
# -3.32948549, -4.99349947, -3.46242765, -3.52481528, -3.36820222,
# -4.1848837 , -1.90748847, -2.09206645, -4.10831718, 2.76094325])
print xresult
rho(xresult,2,init)
plttlin=np.linspace(0,endpoint,100)
pltdata11 = np.array([])
for i in plttlin:
pltdata11 = np.append(pltdata11 ,rho(xresult,i,init)[0] )
print pltdata11
#MMA_optmize_Vac_pltdata = np.genfromtxt('./assets/homogen/MMA_optmize_Vac_pltdata.txt', delimiter = ',')
plt.figure(figsize=(16,9.36))
plt.ylabel('P')
plt.xlabel('Time')
#plt.plot(np.linspace(0,15,4501),MMA_optmize_Vac_pltdata,"r-",label="MMAVacrho11")
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")
print scipy.__version__
0.15.1
# with ramdom initial guess
%%snakeviz
devoendpoint = 2
devotlin = np.linspace(0,devoendpoint,11)
devocostF = lambda x: cost(x,devotlin,init)
bounds=np.zeros([3*4*5,2])
for i in range(3*4*5):
bounds[i,0]=-5
bounds[i,1]=5
#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
11998.9440718 nfev: 1701061 success: True fun: 0.50220275532735592 x: array([-0.38533032, -1.64831529, 4.38168257, 0.38533177, -1.64830245, 4.38164573, 0.2341858 , -1.24310539, 2.7120793 , -1.17805166, -2.02416497, 0.58760437, -1.3433545 , -0.83102552, -1.87713645, 1.56420419, 4.00351989, 1.59243772, 0.21001786, -3.2407917 , -0.21969595, -1.73653359, -0.20812335, -2.54570885, -2.22593135, -1.29257834, 3.92793899, -2.36899476, 4.07210649, 2.49986995, 1.79095183, 2.27666122, -4.01403399, 2.6457365 , -2.00015276, -2.41416953, -0.86537594, 2.60821781, -3.30716236, -2.00846314, 0.5856739 , -1.31974757, 4.97549085, 0.85912326, 0.16231823, 4.78895253, -0.894443 , -3.69719569, 0.49323922, 2.90520738, -1.49061519, 0.28778316, 1.9116826 , -4.82550174, -3.78757261, 1.78348816, -2.56329151, -3.4732061 , 2.09031887, 0.17996168]) message: 'Optimization terminated successfully.' nit: 1889
devoxmid = devo.x
devocostFmid = lambda x: cost(devoxmid,devotlin,init)
startdevo = timeit.default_timer()
devo = differential_evolution(devocostFmid,bounds,strategy='best1bin',tol=1e-10,maxiter=10000,polish=True)
stopdevo = timeit.default_timer()
print stopdevo - startdevo
print devo
7.17289209366 nfev: 1861 success: True fun: 0.50220275532735592 x: array([ 1.24025196, -1.88835977, 4.12070508, 3.11592914, 4.42746631, -4.19707414, 3.41163435, 0.39575733, -1.36370922, 0.10272485, -2.50074836, -0.87324897, 2.5963044 , 4.9221439 , 3.10923979, 3.94014244, -4.62658608, -1.89053953, -2.26800682, -3.20219453, 2.64633316, 2.66767086, 4.1280028 , 1.2258348 , 1.24287907, 4.58866753, -0.45173114, -1.03266506, -0.16861406, 1.93959417, -0.96305919, -4.15124674, 3.31256891, -0.26546418, 1.66579996, -4.49053607, -0.82588613, 2.02271267, -4.28116946, 3.91625742, -1.35224539, -3.730362 , 4.52091474, 3.26261582, 3.07828013, -2.58320648, 0.9280175 , 4.66541153, 4.14676708, 0.62891694, -4.81911323, 2.70865828, 1.01728352, -4.96619298, 0.70850557, 2.2863048 , 2.16242189, 2.20548228, 0.70093172, -4.88308632]) message: 'Optimization terminated successfully.' nit: 1
devoxresult=devo.get("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
[ 1. 1.02464006 1.04924836 1.0738231 1.09836243 1.12286438 1.14732688 1.17174779 1.19612485 1.2204557 1.24473788 1.2689688 1.29314577 1.31726596 1.34132643 1.36532411 1.38925577 1.41311807 1.43690751 1.46062044 1.48425305 1.50780138 1.53126131 1.55462852 1.57789856 1.60106676 1.62412828 1.64707808 1.66991095 1.69262145 1.71520394 1.73765259 1.75996132 1.78212387 1.80413371 1.82598413 1.84766816 1.86917859 1.890508 1.91164869 1.93259274 1.953332 1.97385804 1.9941622 2.01423557 2.03406899 2.05365307 2.07297816 2.09203437 2.11081158 2.12929944 2.14748738 2.16536459 2.18292007 2.20014261 2.21702082 2.23354313 2.2496978 2.26547295 2.28085656 2.2958365 2.31040055 2.32453642 2.33823177 2.35147422 2.36425142 2.37655103 2.38836078 2.39966849 2.41046211 2.42072975 2.43045971 2.43964053 2.44826103 2.45631032 2.46377789 2.47065362 2.47692783 2.48259131 2.48763538 2.49205196 2.49583354 2.49897329 2.50146509 2.50330352 2.50448398 2.50500268 2.50485668 2.50404394 2.50256332 2.50041466 2.49759878 2.49411749 2.48997363 2.48517107 2.47971475 2.47361065 2.46686584 2.45948841 2.45148755]
plt.figure(figsize=(16,9.36))
plt.ylabel('MMArho11')
plt.xlabel('Time')
plt.plot(devoplttlin,1-(np.sin(2.0)**2)*(np.sin(0.5*devoplttlin) )**2)
plt.plot(devoplttlin,devopltdata11,"b4-",label="devo_vac_rho11")
plt.show()
#py.iplot_mpl(plt.gcf(),filename="MMA-rho11-Vac-80-60")
print scipy.__version__
print initGuess
[[-0.95067917 -0.09074367 -4.82398457 -0.7615439 -0.27515522] [-3.33365998 -1.4162127 -2.8632105 -2.7298333 -0.62978916] [-2.52407383 -2.57476616 -4.88282221 -4.65025538 -1.49872754] [-3.15147754 -0.06354641 -2.3908632 -0.86447528 -2.33229558] [-0.5739389 -4.61810188 -0.13228954 -2.58397247 -4.33343592] [-4.83035004 -3.06425551 -0.45936255 -2.91194388 -2.87745423] [-3.93824244 -2.08182497 -3.05850129 -1.92361242 -3.81654995] [-4.19771673 -0.53177727 -3.9155934 -0.04361043 -4.73211585] [-3.68483106 -1.71179777 -1.13691215 -4.794143 -0.47377443] [-3.66788546 -2.8080668 -0.33475682 -4.81757391 -1.18803969] [-3.05637085 -3.98693405 -1.00737055 -2.01300012 -4.43150289] [-2.6875922 -4.52969243 -0.0724917 -3.63840275 -0.79722408]]