<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Test for different solvers for stochastic equation
Based on development-smesolver-new-methods by Manuel Grimm, Niels Lörch, and Denis V. Vasilyev.
Eric Giguere, March 2018
Step1:  Just check that analytical solution coincides with the solution of ODE for the variance
Step2:  Test of different SME solvers
Step3:  Plotting the figure - Constant case
Step4:  Deterministic part time dependent
Step5:  Both d1 and d2 time-dependent
Step6:  Multiple sc_ops, time-dependent
Step7:  Versions
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
%config InlineBackend.figure_formats = ['svg']
from qutip import *
from qutip.ui.progressbar import BaseProgressBar
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
y_sse = None
import time
Explanation: Test for different solvers for stochastic equation
Based on development-smesolver-new-methods by Manuel Grimm, Niels Lörch, and Denis V. Vasilyev.
Eric Giguere, March 2018
End of explanation
def arccoth(x):
    return 0.5*np.log((1.+x)/(x-1.))
############ parameters #############
th = 0.1 # Interaction parameter
alpha = np.cos(th)
beta = np.sin(th)
gamma = 1.
def gammaf(t):
    return 0.25+t/12+t*t/6
def f_gamma(t,*args):
    return (0.25+t/12+t*t/6)**(0.5)
################# Solution of the differential equation for the variance Vc ####################
T = 6.
N_store = 200
tlist = np.linspace(0,T,N_store)
y0 = 0.5
def func(y, t):
    return -(gammaf(t) - alpha*beta)*y - 2*alpha*alpha*y*y + 0.5*gammaf(t)
y_td = odeint(func, y0, tlist)
def func(y, t):
    return -(gamma - alpha*beta)*y - 2*alpha*alpha*y*y + 0.5*gamma
y = odeint(func, y0, tlist)
############ Exact steady state solution for Vc #########################
Vc = (alpha*beta - gamma + np.sqrt((gamma-alpha*beta)**2 + 4*gamma*alpha**2))/(4*alpha**2)
#### Analytic solution
A = (gamma**2 + alpha**2 * (beta**2 + 4*gamma) - 2*alpha*beta*gamma)**0.5
B = arccoth((-4*alpha**2*y0 + alpha*beta - gamma)/A)
y_an = (alpha*beta - gamma + A / np.tanh(0.5*A*tlist - B))/(4*alpha**2)
f, (ax, ax2) = plt.subplots(2, 1, sharex=True)
ax.set_title('Variance as a function of time')
ax.plot(tlist,y)
ax.plot(tlist,Vc*np.ones_like(tlist))
ax.plot(tlist,y_an)
ax.set_ylim(0,0.5)
ax2.set_title('Deviation of odeint from analytic solution')
ax2.set_xlabel('t')
ax2.set_ylabel(r'$\epsilon$')
ax2.plot(tlist,y_an - y.T[0]);
Explanation: Just check that analytical solution coincides with the solution of ODE for the variance
End of explanation
####################### Model ###########################
N = 30                   # number of Fock states
Id = qeye(N)
a = destroy(N)
s = 0.5*((alpha+beta)*a + (alpha-beta)*a.dag())
x = (a + a.dag())/np.sqrt(2)
H = Id
c_op = [np.sqrt(gamma)*a]
c_op_td = [[a,f_gamma]]
sc_op = [s]
e_op = [x, x*x]
rho0 = fock_dm(N,0)       # initial vacuum state
#sc_len=1                 # one stochastic operator
############## time steps and trajectories ###################
ntraj = 1 #100           # number of trajectories
T = 6.                   # final time          
N_store = 200            # number of time steps for which we save the expectation values/density matrix
tlist = np.linspace(0,T,N_store)
ddt = (tlist[1]-tlist[0])
Nsubs = list((13*np.logspace(0,1,10)).astype(np.int))
stepsizes = [ddt/j for j in Nsubs] # step size is doubled after each evaluation 
Nt = len(Nsubs)                    # number of step sizes that we compare
Nsubmax = Nsubs[-1]                # Number of intervals for the smallest step size; 
dtmin = (tlist[1]-tlist[0])/(Nsubmax)
Explanation: Test of different SME solvers
End of explanation
# Analetical solution not available: 
# Compute the evolution with the best solver and very small step size and use it as the reference
sol = ssesolve(H, fock(N), tlist, [sc_op[0]+c_op[0]], e_op, nsubsteps=2000, method="homodyne",solver="taylor2.0")
y_sse = sol.expect[1]-sol.expect[0]*sol.expect[0].conj()
ntraj = 1
def run_sse(**kwargs):
    epsilon = np.zeros(Nt)
    std = np.zeros(Nt)
    print(kwargs)
    for jj in range(0,Nt):
        for j in range(0,ntraj):
            Nsub = Nsubs[jj]#int(Nsubmax/(2**jj))
            sol = ssesolve(H, fock(N), tlist, [sc_op[0]+c_op[0]], e_op, nsubsteps=Nsub, **kwargs)
            epsilon_j = 1/T * np.sum(np.abs(y_sse - (sol.expect[1]-sol.expect[0]*sol.expect[0].conj())))*ddt
            epsilon[jj] += epsilon_j
            std[jj] += epsilon_j
    epsilon/= ntraj
    std = np.sqrt(1/ntraj * (1/ntraj * std - epsilon**2))
    return epsilon
def get_stats(**kw):
    start = time.time()
    y = run_sse(**kw)
    tag = str(kw["solver"])
    x = np.log(stepsizes)
    ly = np.log(y)
    fit = np.polyfit(x, ly, 1)[0]
    return y,tag,fit,time.time()-start
stats_cte = []
stats_cte.append(get_stats(solver='euler-maruyama'))
stats_cte.append(get_stats(solver='platen'))
stats_cte.append(get_stats(solver='pred-corr'))
stats_cte.append(get_stats(solver='milstein'))
stats_cte.append(get_stats(solver='milstein-imp', tol=1e-9))
stats_cte.append(get_stats(solver='pred-corr-2'))
stats_cte.append(get_stats(solver='explicit1.5'))
stats_cte.append(get_stats(solver="taylor1.5"))
stats_cte.append(get_stats(solver="taylor1.5-imp", tol=1e-9))
stats_cte.append(get_stats(solver="taylor2.0"))
stats_cte.append(get_stats(solver="taylor2.0", noiseDepth=500))
fig = plt.figure()
ax = plt.subplot(111)
mark = "o*vspx+^<>1hdD"
for i,run in enumerate(stats_cte):
    ax.loglog(stepsizes, run[0], mark[i], label=run[1]+": " + str(run[2]))
ax.loglog(stepsizes, 0.003*np.array(stepsizes)**0.5, label="$\propto\Delta t^{1/2}$")
ax.loglog(stepsizes, 0.01*np.array(stepsizes)**1, label="$\propto\Delta t$")
ax.loglog(stepsizes, 0.001*np.array(stepsizes)**1, label="$\propto\Delta t$")
ax.loglog(stepsizes, 0.01*np.array(stepsizes)**1.5, label="$\propto\Delta t^{3/2}$")
ax.loglog(stepsizes, 0.05*np.array(stepsizes)**2.0, label="$\propto\Delta t^{2}$")
ax.set_xlabel(r'$\Delta t$  $\left[\gamma^{-1}\right]$')
ax.set_ylabel('deviation')
lgd=ax.legend(loc='center left', bbox_to_anchor=(1, 0.64), prop={'size':12})
Explanation: Plotting the figure - Constant case
End of explanation
def H_f(t,args):
    return 0.125+t/12+t*t/72
sol = ssesolve([H,[c_op[0].dag()*c_op[0]/2,H_f]], fock(N), tlist, sc_op, e_op, 
                       nsubsteps=2500, method="homodyne",solver="taylor2.0")
y_sse_td = sol.expect[1]-sol.expect[0]*sol.expect[0].conj()
plt.plot(y_sse_td)
ntraj = 1
def run_sse_td(**kwargs):
    epsilon = np.zeros(Nt)
    std = np.zeros(Nt)
    print(kwargs)
    for jj in range(0,Nt):
        for j in range(0,ntraj):
            Nsub = Nsubs[jj]#int(Nsubmax/(2**jj))
            sol = ssesolve([H,[c_op[0].dag()*c_op[0]/2,H_f]], fock(N), tlist, sc_op, e_op, nsubsteps=Nsub, **kwargs)
            epsilon_j = 1/T * np.sum(np.abs(y_sse_td - (sol.expect[1]-sol.expect[0]*sol.expect[0].conj())))*ddt
            epsilon[jj] += epsilon_j
            std[jj] += epsilon_j
    epsilon/= ntraj
    std = np.sqrt(1/ntraj * (1/ntraj * std - epsilon**2))
    return epsilon
def get_stats(**kw):
    y = run_sse_td(**kw)
    tag = str(kw["solver"])
    x = np.log(stepsizes)
    ly = np.log(y)
    fit = np.polyfit(x, ly, 1)[0]
    return y,tag,fit
stats_td = []
stats_td.append(get_stats(solver='euler-maruyama'))
stats_td.append(get_stats(solver='platen'))
stats_td.append(get_stats(solver='pred-corr'))
stats_td.append(get_stats(solver='milstein'))
stats_td.append(get_stats(solver='milstein-imp'))
stats_td.append(get_stats(solver='pred-corr-2'))
stats_td.append(get_stats(solver='explicit1.5'))
stats_td.append(get_stats(solver="taylor1.5"))
stats_td.append(get_stats(solver="taylor1.5-imp", tol=1e-9))
stats_td.append(get_stats(solver="taylor2.0"))
stats_td.append(get_stats(solver="taylor2.0", noiseDepth=500))
fig = plt.figure()
ax = plt.subplot(111)
mark = "o*vspx+^<>1hdD"
for i,run in enumerate(stats_td):
    ax.loglog(stepsizes, run[0], mark[i], label=run[1]+": " + str(run[2]))
    
ax.loglog(stepsizes, 0.1*np.array(stepsizes)**0.5, label="$\propto\Delta t^{1/2}$")
ax.loglog(stepsizes, 0.1*np.array(stepsizes)**1, label="$\propto\Delta t$")
ax.loglog(stepsizes, 0.1*np.array(stepsizes)**1.5, label="$\propto\Delta t^{3/2}$")
ax.loglog(stepsizes, 0.5*np.array(stepsizes)**2.0, label="$\propto\Delta t^{2}$")
ax.set_xlabel(r'$\Delta t$  $\left[\gamma^{-1}\right]$')
ax.set_ylabel('deviation')
lgd=ax.legend(loc='center left', bbox_to_anchor=(1, 0.64), prop={'size':12})
Explanation: Deterministic part time dependent
End of explanation
def H_f(t,args):
    return 0.125+t/12+t*t/72
def H_bf(t,args):
    return 0.125+t/10+t*t/108
sc_op_td = [[sc_op[0],H_bf]]
sol = ssesolve([H,[c_op[0].dag()*c_op[0]/2,H_f]], fock(N), tlist, sc_op_td, e_op, 
                       nsubsteps=2000, method="homodyne",solver="taylor15")
y_sse_btd = sol.expect[1]-sol.expect[0]*sol.expect[0].conj()
plt.plot(y_sse_btd)
ntraj = 1
def run_sse_td(**kwargs):
    epsilon = np.zeros(Nt)
    std = np.zeros(Nt)
    print(kwargs)
    for jj in range(0,Nt):
        for j in range(0,ntraj):
            Nsub = Nsubs[jj]#int(Nsubmax/(2**jj))
            sol = ssesolve([H,[c_op[0].dag()*c_op[0]/2,H_f]], fock(N), tlist, sc_op_td, e_op, nsubsteps=Nsub, **kwargs)
            epsilon_j = 1/T * np.sum(np.abs(y_sse_btd - (sol.expect[1]-sol.expect[0]*sol.expect[0].conj())))*ddt
            epsilon[jj] += epsilon_j
            std[jj] += epsilon_j
    epsilon/= ntraj
    std = np.sqrt(1/ntraj * (1/ntraj * std - epsilon**2))
    return epsilon
def get_stats_b(**kw):
    y = run_sse_td(**kw)
    tag = str(kw["solver"])
    x = np.log(stepsizes)
    ly = np.log(y)
    fit = np.polyfit(x, ly, 1)[0]
    return y,tag,fit
stats_d2_td = []
stats_d2_td.append(get_stats_b(solver='euler-maruyama'))
stats_d2_td.append(get_stats_b(solver='platen'))
stats_d2_td.append(get_stats_b(solver='pred-corr'))
stats_d2_td.append(get_stats_b(solver='milstein'))
stats_d2_td.append(get_stats_b(solver='milstein-imp'))
stats_d2_td.append(get_stats_b(solver='pred-corr-2'))
stats_d2_td.append(get_stats_b(solver='explicit1.5'))
stats_d2_td.append(get_stats_b(solver="taylor1.5"))
stats_d2_td.append(get_stats_b(solver="taylor1.5-imp", tol=1e-9))
fig = plt.figure()
ax = plt.subplot(111)
mark = "o*vspx+^<>1hdD"
for i,run in enumerate(stats_d2_td):
    ax.loglog(stepsizes, run[0], mark[i], label=run[1]+": " + str(run[2]))
ax.loglog(stepsizes, 0.03*np.array(stepsizes)**0.5, label="$\propto\Delta t^{1/2}$")
ax.loglog(stepsizes, 0.03*np.array(stepsizes)**1, label="$\propto\Delta t$")
ax.loglog(stepsizes, 0.03*np.array(stepsizes)**1.5, label="$\propto\Delta t^{3/2}$")
ax.set_xlabel(r'$\Delta t$  $\left[\gamma^{-1}\right]$')
ax.set_ylabel('deviation')
lgd=ax.legend(loc='center left', bbox_to_anchor=(1, 0.64), prop={'size':12})
Explanation: Both d1 and d2 time-dependent
End of explanation
def H_f(t,args):
    return 0.125+t/12+t*t/36
def H_bf(t,args):
    return 0.125+t/10+t*t/108
sc_op_td = [[sc_op[0]],[sc_op[0],H_bf],[sc_op[0],H_f]]
sol = ssesolve([H,[c_op[0].dag()*c_op[0]/2,H_f]], fock(N), tlist/3, sc_op_td, e_op, 
                       nsubsteps=2000, method="homodyne",solver="taylor15")
y_sse_multi = sol.expect[1]-sol.expect[0]*sol.expect[0].conj()
plt.plot(y_sse_multi)
ntraj = 1
def run_sss_multi(**kwargs):
    epsilon = np.zeros(Nt)
    std = np.zeros(Nt)
    print(kwargs)
    for jj in range(0,Nt):
        for j in range(0,ntraj):
            Nsub = Nsubs[jj]#int(Nsubmax/(2**jj))
            sol = ssesolve([H,[c_op[0].dag()*c_op[0]/2,H_f]], fock(N), tlist/3, sc_op_td, e_op, nsubsteps=Nsub, **kwargs)
            epsilon_j = 1/T * np.sum(np.abs(y_sse_multi - (sol.expect[1]-sol.expect[0]*sol.expect[0].conj())))*ddt
            epsilon[jj] += epsilon_j
            std[jj] += epsilon_j
    epsilon/= ntraj
    std = np.sqrt(1/ntraj * (1/ntraj * std - epsilon**2))
    
    return epsilon
def get_stats_multi(**kw):
    y = run_sss_multi(**kw)
    tag = str(kw["solver"])
    x = np.log(stepsizes)
    ly = np.log(y)
    fit = np.polyfit(x, ly, 1)[0]
    return (y,tag,fit)
stats_multi = []
stats_multi.append(get_stats_multi(solver='euler-maruyama'))
stats_multi.append(get_stats_multi(solver="platen"))
stats_multi.append(get_stats_multi(solver='pred-corr'))
stats_multi.append(get_stats_multi(solver='milstein'))
stats_multi.append(get_stats_multi(solver='milstein-imp'))
stats_multi.append(get_stats_multi(solver='pred-corr-2'))
stats_multi.append(get_stats_multi(solver='explicit1.5'))
stats_multi.append(get_stats_multi(solver="taylor1.5"))
stats_multi.append(get_stats_multi(solver="taylor1.5-imp", tol=1e-9))
fig = plt.figure()
ax = plt.subplot(111)
mark = "o*vspx+^<>Dd"
for run in stats_multi:
    ax.loglog(stepsizes, run[0], 'o', label=run[1]+": " + str(run[2]))
ax.loglog(stepsizes, 0.05*np.array(stepsizes)**0.5, label="$\propto\Delta t^{1/2}$")
ax.loglog(stepsizes, 0.05*np.array(stepsizes)**1, label="$\propto\Delta t$")
ax.loglog(stepsizes, 0.05*np.array(stepsizes)**1.5, label="$\propto\Delta t^{3/2}$")
ax.set_xlabel(r'$\Delta t$  $\left[\gamma^{-1}\right]$')
ax.set_ylabel('deviation')
lgd=ax.legend(loc='center left', bbox_to_anchor=(1, 0.64), prop={'size':12})
Explanation: Multiple sc_ops, time-dependent
End of explanation
from qutip.ipynbtools import version_table
version_table()
Explanation: Versions
End of explanation
<END_TASK>