#!/usr/bin/env python
# coding: utf-8

# <a id='toc'></a>
# -----
# ### TOC
# 
# 
# * [NumpyOPT](#numpyopt)
# * [Differential Evolution](#evo)
# 

# In[1]:


# 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.
get_ipython().run_line_magic('matplotlib', 'inline')
get_ipython().run_line_magic('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]:


# Here are the initial values


# For test use
array12 = np.asarray(np.split(np.random.rand(1,60)[0],12))


# In[3]:


# Here is the activation function

def act(x):
    return expit(x)


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


rho(array12,0,init)


# In[6]:


# 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


# In[8]:


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


# In[9]:


costi(array12,0,init)


# In[10]:


def cost(x,t,initialCondition):
    
    costTotal = map(lambda t: costi(x,t,initialCondition),t)
    
    return 0.5*np.sum(costTotal)


# In[11]:


cost(array12,np.array([0,1,2]),init)
#cost(xresult,np.array([0,4,11]),init)


# <a id='numpyopt'></a>
# -----
# ## NUMPY OPT

# In[12]:


# 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


# - **Should think about the eps(stepsize)**

# In[ ]:


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


# In[ ]:


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


# In[ ]:


rho(xresult,2,init)


# 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


# In[ ]:


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


# In[ ]:





# In[18]:


print scipy.__version__


# <a id='evo'></a>
# -----
# ## Test Differential Evolution

# In[27]:


# with ramdom initial guess
get_ipython().run_line_magic('%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


# In[28]:


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


# In[29]:


devoxresult=devo.get("x")


# In[30]:


devoplttlin=np.linspace(0,devoendpoint,100)
devopltdata11 = np.array([])
for i in devoplttlin:
    devopltdata11 = np.append(devopltdata11 ,rho(devoxresult,i,init)[0] )
    
print devopltdata11


# In[31]:


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")


# In[ ]:





# In[ ]:


print scipy.__version__


# In[80]:


print initGuess


# In[ ]: