<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Homogeneous Gas
Here is a notebook for homogeneous gas model.
Here we are talking about a homogeneous gas bulk of neutrinos with single energy. The EoM is
$$
i \partial_t \rho_E = \left[ \frac{\delta m^2}{2E}B +\lambda L +\sqrt 2 G_F \int_0^\infty dE' ( \rho_{E'} - \bar \rho_{E'} ) ,\rho_E \right]
$$
while the EoM for antineutrinos is
$$
i \partial_t \bar\rho_E = \left[- \frac{\delta m^2}{2E}B +\lambda L +\sqrt 2 G_F \int_0^\infty dE' ( \rho_{E'} - \bar \rho_{E'} ) ,\bar\rho_E \right]
$$
Initial
Step1:  Using Mathematica, I can find the 4*2 equations
Step2:  I am going to substitute all density matrix elements using their corrosponding network expressions.
So first of all, I need the network expression for the unknown functions.
A function is written as
$$ y_i= 1+t_i v_k f(t_i w_k+u_k) ,$$
while it's derivative is
$$v_k f(t w_k+u_k) + t v_k f(tw_k+u_k) (1-f(tw_k+u_k)) w_k .$$
Now I can write down the equations using these two forms.
Step3:  Minimization
Here is the minimization
Step4:  Functions
Find the solutions to each elements.
Step5:  Practice
<ASSISTANT_TASK:>
Python Code:
# 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.special import expit
import matplotlib.pyplot as plt
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
# hbar=1.054571726*10**(-34)
hbar=1
delm2E=1
lamb=1  ## lambda for neutrinos
lambb=1 ## lambda for anti neutrinos
gF=1
nd=1  ## number density
ndb=1   ## number density
omega=1
omegab=1
## Here are some matrices to be used
elM = np.array([[1,0],[0,0]])
bM = 1/2*np.array( [ [ - 0.38729833462,0.31622776601] , [0.31622776601,0.38729833462] ] )
## sqareroot of 2
sqrt2=np.sqrt(2)
Explanation: Homogeneous Gas
Here is a notebook for homogeneous gas model.
Here we are talking about a homogeneous gas bulk of neutrinos with single energy. The EoM is
$$
i \partial_t \rho_E = \left[ \frac{\delta m^2}{2E}B +\lambda L +\sqrt 2 G_F \int_0^\infty dE' ( \rho_{E'} - \bar \rho_{E'} ) ,\rho_E \right]
$$
while the EoM for antineutrinos is
$$
i \partial_t \bar\rho_E = \left[- \frac{\delta m^2}{2E}B +\lambda L +\sqrt 2 G_F \int_0^\infty dE' ( \rho_{E'} - \bar \rho_{E'} ) ,\bar\rho_E \right]
$$
Initial:
Homogeneous, Isotropic, Monoenergetic $\nu_e$ and $\bar\nu_e$
The equations becomes
$$
i \partial_t \rho_E = \left[ \frac{\delta m^2}{2E} B +\lambda L +\sqrt 2 G_F   ( \rho_{E} - \bar \rho_{E} ) ,\rho_E \right]
$$
$$
i \partial_t \bar\rho_E = \left[- \frac{\delta m^2}{2E}B +\lambda_b L +\sqrt 2 G_F ( \rho_{E} - \bar \rho_{E} ) ,\bar\rho_E \right]
$$
Define $\omega=\frac{\delta m^2}{2E}$,  $\omega = \frac{\delta m^2}{-2E}$, $\mu=\sqrt{2}G_F n_\nu$
$$
i \partial_t \rho_E = \left[ \omega B +\lambda L +\mu   ( \rho_{E} - \bar \rho_{E} ) ,\rho_E \right]
$$
$$
i \partial_t \bar\rho_E = \left[\bar\omega B +\bar\lambda L +\mu ( \rho_{E} - \bar \rho_{E} ) ,\bar\rho_E \right]
$$
where
$$
B = \frac{1}{2} \begin{pmatrix} 
-\cos 2\theta_v & \sin 2\theta_v \
\sin 2\theta_v & \cos 2\theta_v
\end{pmatrix} = 
\begin{pmatrix} 
-0.38729833462 & 0.31622776601\
0.31622776601 & 0.38729833462
\end{pmatrix}
$$
$$
L = \begin{pmatrix}
1 & 0 \
0 & 0
\end{pmatrix}
$$
Initial condition 
$$
\rho(t=0) = \begin{pmatrix}
1 & 0 \
0 & 0
\end{pmatrix}
$$
$$
\bar\rho(t=0) =\begin{pmatrix}
1 & 0 \
0 & 0
\end{pmatrix}
$$
define the following quantities
hbar$=\hbar$
delm2E$= \delta m^2/2E$
lamb $= \lambda$, lambb $= \bar\lambda$
gF $= G_F$
mu $=\mu$
omega $=\omega$, omegab $=\bar\omega$
Numerical
End of explanation
#r11prime(t)
## The matrix eqn for neutrinos. Symplify the equation to the form A.X=0. Here I am only writing down the LHS.
## Eqn for r11'
# 1/2*( r21(t)*( bM12*delm2E - 2*sqrt2*gF*rb12(t) ) + r12(t) * ( -bM21*delm2E + 2*sqrt2*gF*rb21(t) ) - 1j*r11prime(t)  )
## Eqn for r12'
# 1/2*( r22(t)* ( bM12 ) )
### wait a minute I don't actually need to write down this. I can just do this part in numpy.
Explanation: Using Mathematica, I can find the 4*2 equations
End of explanation
def trigf(x):
    #return 1/(1+np.exp(-x)) # It's not bad to define this function here for people could use other functions other than expit(x).
    return expit(x)
## The time derivative part
### Here are the initial conditions
init = np.array( [[1,0],[0,0]] )
initb = np.array([[0,0],[0,0]])
### For neutrinos
def rho(x,ti,initialCondition): # x is the input structure arrays, ti is a time point
    v11,w11,u11,v12,w12,u12,v21,w21,u21,v22,w22,u22 = np.split(x,12)[:12]
        
    elem11= np.sum(ti * v11 * trigf( ti*w11 +u11 ) )
    elem12= np.sum(ti * v12 * trigf( ti*w12 +u12 ) )
    elem21= np.sum(ti * v21 * trigf( ti*w21 +u21 ) )
    elem22= np.sum(ti * v22 * trigf( ti*w22 +u22 ) )
    
    return initialCondition + np.array([[ elem11 , elem12 ],[elem21, elem22]])
def rhob(xb,ti,initialConditionb): # x is the input structure arrays, ti is a time point
    vb11,wb11,ub11,vb12,wb12,ub12,vb21,wb21,ub21,vb22,wb22,ub22 = np.split(xb,12)[:12]
    elem11= np.sum(ti * vb11 * trigf( ti*wb11 +ub11 ) )
    elem12= np.sum(ti * vb12 * trigf( ti*wb12 +ub12 ) )
    elem21= np.sum(ti * vb21 * trigf( ti*wb21 +ub21 ) )
    elem22= np.sum(ti * vb22 * trigf( ti*wb22 +ub22 ) )
    
    return initialConditionb + np.array([[ elem11 , elem12 ],[elem21, elem22]])
## Test
xtemp=np.ones(120)
rho(xtemp,1,init)
## Define Hamiltonians for both
def hamil(x,xb,ti,initialCondition,initialConditionb):
    
    return delm2E*bM + lamb*elM + sqrt2*gF*( rho(x,ti,initialCondition) - rhob(xb,ti,initialConditionb) )
def hamilb(x,xb,ti,initialCondition,initialConditionb):
    
    return -delm2E*bM + lambb*elM + sqrt2*gF*( rho(x,ti,initialCondition) - rhob(xb,ti,initialConditionb) )
## The commutator
def comm(x,xb,ti,initialCondition,initialConditionb):
    
    return np.dot(hamil(x,xb,ti,initialCondition,initialConditionb), rho(x,ti,initialCondition) ) - np.dot(rho(x,ti,initialCondition), hamil(x,xb,ti,initialCondition,initialConditionb) )
def commb(x,xb,ti,initialCondition,initialConditionb):
    
    return np.dot(hamilb(x,xb,ti,initialCondition,initialConditionb), rhob(xb,ti,initialConditionb) ) - np.dot(rhob(xb,ti,initialConditionb), hamilb(x,xb,ti,initialCondition,initialConditionb) )
## Test
print "neutrino\n",comm(xtemp,xtemp,1,init,initb)
print "antineutrino\n",commb(xtemp,xtemp,0.5,init,initb)
## The COST of the eqn set
def costTi(x,xb,ti,initialCondition,initialConditionb):
    
    v11,w11,u11,v12,w12,u12,v21,w21,u21,v22,w22,u22 = np.split(x,12)[:12]
    vb11,wb11,ub11,vb12,wb12,ub12,vb21,wb21,ub21,vb22,wb22,ub22 = np.split(xb,12)[:12]
    
    fvec11 = np.array(trigf(ti*w11 + u11) )  # This is a vector!!!
    fvec12 = np.array(trigf(ti*w12 + u12) )
    fvec21 = np.array(trigf(ti*w21 + u21) )
    fvec22 = np.array(trigf(ti*w22 + u22) )
    
    fvecb11 = np.array(trigf(ti*wb11 + ub11) )  # This is a vector!!!
    fvecb12 = np.array(trigf(ti*wb12 + ub12) )
    fvecb21 = np.array(trigf(ti*wb21 + ub21) )
    fvecb22 = np.array(trigf(ti*wb22 + ub22) )
    
    
    costi11= ( np.sum (v11*fvec11 + ti * v11* fvec11 * ( 1 -  fvec11  ) * w11 ) + 1j*  ( comm(x,xb,ti,initialCondition,initialConditionb)[0,0] )  )  
    costi12= ( np.sum (v12*fvec12 + ti * v12* fvec12 * ( 1 -  fvec12  ) * w12 ) + 1j*  ( comm(x,xb,ti,initialCondition,initialConditionb)[0,1] )  )  
    costi21= ( np.sum (v21*fvec21 + ti * v21* fvec21 * ( 1 -  fvec21  ) * w21 ) + 1j*  ( comm(x,xb,ti,initialCondition,initialConditionb)[1,0] )  )  
    costi22= ( np.sum (v22*fvec22 + ti * v22* fvec22 * ( 1 -  fvec22  ) * w22 ) + 1j*  ( comm(x,xb,ti,initialCondition,initialConditionb)[1,1] )  )  
    costbi11= ( np.sum (vb11*fvecb11 + ti * vb11* fvecb11 * ( 1 -  fvecb11  ) * wb11 ) + 1j*  ( commb(x,xb,ti,initialCondition,initialConditionb)[0,0] )  )  
    costbi12= ( np.sum (vb12*fvecb12 + ti * vb12* fvecb12 * ( 1 -  fvecb12  ) * wb12 ) + 1j*  ( commb(x,xb,ti,initialCondition,initialConditionb)[0,1] )  )  
    costbi21= ( np.sum (vb21*fvecb21 + ti * vb21* fvecb21 * ( 1 -  fvecb21  ) * wb21 ) + 1j*  ( commb(x,xb,ti,initialCondition,initialConditionb)[1,0] )  )  
    costbi22= ( np.sum (vb22*fvecb22 + ti * vb22* fvecb22 * ( 1 -  fvecb22  ) * wb22 ) + 1j*  ( commb(x,xb,ti,initialCondition,initialConditionb)[1,1] )  )  
    
    return (np.real(costi11))**2 + (np.real(costi12))**2+ (np.real(costi21))**2 +  (np.real(costi22))**2 + (np.real(costbi11))**2 + (np.real(costbi12))**2 +(np.real(costbi21))**2 + (np.real(costbi22))**2 + (np.imag(costi11))**2 + (np.imag(costi12))**2+ (np.imag(costi21))**2 +  (np.imag(costi22))**2 + (np.imag(costbi11))**2 + (np.imag(costbi12))**2 +(np.imag(costbi21))**2 + (np.imag(costbi22))**2    
costTi(xtemp,xtemp,0,init,initb)
## Calculate the total cost
def cost(xtot,t,initialCondition,initialConditionb):
    
    x,xb = np.split(xtot,2)[:2]
    t = np.array(t)
    
    costTotal = np.sum( costTList(x,xb,t,initialCondition,initialConditionb)  )
        
    return costTotal
    
def costTList(x,xb,t,initialCondition,initialConditionb):  ## This is the function WITHOUT the square!!! 
        
    t = np.array(t)
    
    costList = np.asarray([])
    
    for temp in t:
        tempElement = costTi(x,xb,temp,initialCondition,initialConditionb)
        costList = np.append(costList, tempElement)
        
    return np.array(costList)
    
ttemp = np.linspace(0,10)
print ttemp
ttemp = np.linspace(0,10)
print costTList(xtemp,xtemp,ttemp,init,initb)
print cost(xtemp,ttemp,init,initb)
Explanation: I am going to substitute all density matrix elements using their corrosponding network expressions.
So first of all, I need the network expression for the unknown functions.
A function is written as
$$ y_i= 1+t_i v_k f(t_i w_k+u_k) ,$$
while it's derivative is
$$v_k f(t w_k+u_k) + t v_k f(tw_k+u_k) (1-f(tw_k+u_k)) w_k .$$
Now I can write down the equations using these two forms.
End of explanation
tlin = np.linspace(0,0.5,3)
initGuess = np.ones(120)
# initGuess = np.random.rand(1,30)+2
costF = lambda x: cost(x,tlin,init,initb)
cost(initGuess,tlin,init,initb)
## %%snakeviz
# startCG = timeit.default_timer()
#costFResultCG = minimize(costF,initGuess,method="CG")
#stopCG = timeit.default_timer()
#print stopCG - startCG
#print costFResultCG
%%snakeviz
startSLSQP = timeit.default_timer()
costFResultSLSQP = minimize(costF,initGuess,method="SLSQP")
stopSLSQP = timeit.default_timer()
print stopSLSQP - startSLSQP
print costFResultSLSQP
costFResultSLSQP.get('x')
np.savetxt('./assets/homogen/optimize_ResultSLSQPT2120.txt', costFResultSLSQP.get('x'), delimiter = ',')
Explanation: Minimization
Here is the minimization
End of explanation
# costFResultSLSQPx = np.genfromtxt('./assets/homogen/optimize_ResultSLSQP.txt', delimiter = ',')
## The first element of neutrino density matrix
xresult = np.split(costFResultSLSQP.get('x'),2)[0]
xresultb = np.split(costFResultSLSQP.get('x'),2)[1]
#xresult = np.split(costFResultSLSQPx,2)[0]
#xresultb = np.split(costFResultSLSQPx,2)[1]
## print xresult11
plttlin=np.linspace(0,5,100)
pltdata11 = np.array([])
pltdata22 = np.array([])
for i in plttlin:
    pltdata11 = np.append(pltdata11 ,rho(xresult,i,init)[0,0] )
    
print pltdata11
for i in plttlin:
    pltdata22 = np.append(pltdata22 ,rho(xresult,i,init)[1,1] )
    
print pltdata22
print "----------------------------------------"
pltdatab11 = np.array([])
pltdatab22 = np.array([])
for i in plttlin:
    pltdatab11 = np.append(pltdatab11 ,rho(xresultb,i,init)[0,0] )
    
print pltdatab11
for i in plttlin:
    pltdatab22 = np.append(pltdatab22 ,rho(xresultb,i,init)[1,1] )
    
print pltdatab22
#np.savetxt('./assets/homogen/optimize_pltdatar11.txt', pltdata11, delimiter = ',')
#np.savetxt('./assets/homogen/optimize_pltdatar22.txt', pltdata22, delimiter = ',')
plt.figure(figsize=(20,12.36))
plt.ylabel('Time')
plt.xlabel('rho11')
plt.plot(plttlin,pltdata11,"b4-",label="rho11")
py.iplot_mpl(plt.gcf(),filename="HG-rho11")
# tls.embed("https://plot.ly/~emptymalei/73/")
plt.figure(figsize=(20,12.36))
plt.ylabel('Time')
plt.xlabel('rho22')
plt.plot(plttlin,pltdata22,"r4-",label="rho22")
py.iplot_mpl(plt.gcf(),filename="HG-rho22")
plt.figure(figsize=(20,12.36))
plt.ylabel('Time')
plt.xlabel('rhob11')
plt.plot(plttlin,pltdatab11,"b*-",label="rhob11")
py.iplot_mpl(plt.gcf(),filename="HG-rhob11")
plt.figure(figsize=(20,12.36))
plt.ylabel('Time')
plt.xlabel('rhob22')
plt.plot(plttlin,pltdatab11,"b*-",label="rhob22")
py.iplot_mpl(plt.gcf(),filename="HG-rhob22")
MMA_optmize_pltdata = np.genfromtxt('./assets/homogen/MMA_optmize_pltdata.txt', delimiter = ',')
plt.figure(figsize=(20,12.36))
plt.ylabel('MMArho11')
plt.xlabel('Time')
plt.plot(np.linspace(0,5,501),MMA_optmize_pltdata,"r-",label="MMArho11")
plt.plot(plttlin,pltdata11,"b4-",label="rho11")
py.iplot_mpl(plt.gcf(),filename="MMA-rho11")
Explanation: Functions
Find the solutions to each elements.
End of explanation
xtemp1 = np.arange(4)
xtemp1.shape = (2,2)
print xtemp1
xtemp1[0,1]
np.dot(xtemp1,xtemp1)
xtemp1[0,1]
Explanation: Practice
End of explanation
<END_TASK>