<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Runtime Cythonize in QuTiP
Prepared for EuroSciPy 2019
Alex Pitchford (alex.pitchford@gmail.com)
Two types of time-dependent dynamics solving
'function' type and 'string' type
pros and cons to both
string type is a key differentiating feature of QuTiP
recently much improved
scope for further enhancement
Ordinary differential equation
of the form
$$ \frac{\textrm{d} x}{\textrm{d} t} = L(t) \; x $$
$x$ -- vector
Step1:  Time dependent control functions
Step2:  Hamiltonians, initial state and measurements
Step3:  Solving the dynamics
Step4:  Function type
Step5:  String type
Step6:  Comparing execution times
<ASSISTANT_TASK:>
Python Code:
# Imports and utility functions
import time
import numpy as np
import matplotlib.pyplot as plt
from qutip.sesolve import sesolve
from qutip.solver import Options, solver_safe
from qutip import sigmax, sigmay, sigmaz, identity, tensor, basis, Bloch
def timing_val(func):
    def wrapper(*arg, **kw):
        '''source: http://www.daniweb.com/code/snippet368.html'''
        t1 = time.time()
        res = func(*arg, **kw)
        t2 = time.time()
        return res, (t2 - t1), func.__name__
    return wrapper
def plot_exp(tlist, expects, lbls, title):
    fig = plt.figure()
    ax = fig.add_subplot(111)
    for i, e in enumerate(expects):#
        ax.plot(tlist, e, label=lbls[i])
    ax.set_xlabel(r"$t$")
    ax.set_title(title)
    ax.legend()
Explanation: Runtime Cythonize in QuTiP
Prepared for EuroSciPy 2019
Alex Pitchford (alex.pitchford@gmail.com)
Two types of time-dependent dynamics solving
'function' type and 'string' type
pros and cons to both
string type is a key differentiating feature of QuTiP
recently much improved
scope for further enhancement
Ordinary differential equation
of the form
$$ \frac{\textrm{d} x}{\textrm{d} t} = L(t) \; x $$
$x$ -- vector:
state $\newcommand{\ket}[1]{\left|{#1}\right\rangle} \ket{\psi}$ 
or vectorised density matrix $\vec{\rho}$
or $x$ -- matrix:
unitary transformation $U$ 
or dynamical map $M$
and $L$ -- sparse matrix that drives the dynamics.
In Python these are Qobj
wrapper to fast_csr_matrix -- custom to QuTiP.
Time dependent dynamics generator
$$ L(t) = L_0 + g_1(t) \, L_1 + g_2(t) \, L_2 + \ldots $$
$L(t)$ defined in Python as list of matrices 
and (optionally) time-dependent functions
L = [L0, [L1, g1], [L2, g2] ....]
L -- Qobj
g -- Python function, e.g. def g(t, args): return np.cos(args['w']*t)
or string e.g. "cos(w*t)"
Function vs string type time-dependence
| Python function | string type |
| :-: | :-: |
| Python interpreter | Compiled C++ |
| slower execution | Compile time overhead |
Note: C++ compiler required at runtime for string type. Non-standard in MS Windows
String type compilation
Dynamic generation of Cython code .pyx file
right hand side of ODE cqobjevo_compiled_coeff_&lt;hash&gt;.pyx
Dynamic Python to import time-dependent QobjEvo
exec(compile('from ' + filename + ' import CompiledStrCoeff, locals()))
Triggers compilation
Temporary file .pyx deleted
Compiled C++ time-dependent QobjEvo can be reused
Example - control coupled qubits
Unitary dynamics - Schrödinger's equation
$$ \frac{\textrm{d} \ket{\psi}}{\textrm{d} t} = -\textrm{i} \, H(t) \, \ket{\psi} $$
$H(t)$ -- combined Hamiltonian of the system
$$ H(t) = \omega_1 \, H_{\textrm{d}}^{(1)} + \omega_2 \, H_{\textrm{d}}^{(2)} + \omega_{\textrm{i}} \, H_{\textrm{i}} + g_1(t) \, H_{\textrm{c}}^{(1)} + g_2(t) \, H_{\textrm{c}}^{(2)} $$
End of explanation
def g_sine(t, args):
    return args['A_s']*np.sin(args['w_s']*t)
def g_decay(t, args):
    return args['A_d']*np.exp(-t/args['t_d'])
g_sine_str = "A_s*sin(w_s*t)"
g_decay_str = "A_d*exp(-t/t_d)"
t_tot = 10.0
w_1 = 0.3
w_2 = 0.3
w_i = 0.02
A_s = 0.3
A_d = 0.3
w_s = np.pi/t_tot
t_d = 5.0
tlist = np.linspace(0.0, t_tot, 200)
args = {'A_s': A_s, 'A_d': A_s, 'w_s': w_s, 't_d': t_d}
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(tlist, g_sine(tlist, args), label=r"$g_{\rm{sin}}$") #$")
ax.plot(tlist, g_decay(tlist, args), label=r"$g_{\rm{decay}}$")
ax.set_xlabel(r"$t$")
ax.set_title("Control functions")
ax.legend()
Explanation: Time dependent control functions
End of explanation
Id2 = identity(2)
Sz1 = tensor(sigmaz(), Id2)
Sz2 = tensor(Id2, sigmaz())
Sx1 = tensor(sigmax(), Id2)
Sy1 = tensor(sigmay(), Id2)
Sx2 = tensor(Id2, sigmax())
Sy2 = tensor(Id2, sigmay())
H_d1 = w_1*Sz1
H_d2 = w_2*Sz2
H_c1 = Sx1
H_c2 = Sy2
H_i = w_i*tensor(sigmaz(), sigmaz())
H_func_type = [H_d1, H_d2, [H_c1, g_sine], [H_c2, g_decay], H_i]
H_str_type = [H_d1, H_d2, [H_c1, g_sine_str], [H_c2, g_decay_str], H_i]
up_state = basis(2, 0)
b = Bloch()
b.add_states(up_state)
b.show()
init_state = tensor(up_state, up_state)
meas = [Sz1, Sx1, Sz2, Sx2]
Explanation: Hamiltonians, initial state and measurements
End of explanation
@timing_val
def repeat_solve(H, init_state, tlist, num_reps=1,
                 e_ops=None, args=None, options=None):
    if options is None:
        options = Options()
    out = sesolve(H, init_state, tlist,
              e_ops=meas, args=args, options=options)
    if num_reps > 1:
        options.rhs_reuse = True
        tl = np.array([0, tlist[-1]])
        for i in range(num_reps - 1):
            sesolve(H, init_state, tl,
                          e_ops=meas, args=args, options=options)
    return out
Explanation: Solving the dynamics
End of explanation
n_reps = 1
out, t_func, fname = repeat_solve(H_func_type, init_state, tlist, num_reps=n_reps,
                                  e_ops=meas, args=args)
print("{} execution of func type took {} seconds.".format(n_reps, t_func))
# Plot qubit 1 expectations
plot_exp(tlist, out.expect[:2], lbls=["E[Z]", "E[X]"],
         title="Qubit 1 - func type")
Explanation: Function type
End of explanation
n_reps = 1
out, t_func, fname = repeat_solve(H_str_type, init_state, tlist, num_reps=n_reps,
                                  e_ops=meas, args=args)
print("{} execution of string type took {} seconds.".format(n_reps, t_func))
# Plot qubit 1 expectations
plot_exp(tlist, out.expect[:2], lbls=["E[Z]", "E[X]"],
         title="Qubit 1 - string type")
Explanation: String type
End of explanation
n_rep_list = [1, 2, 5, 10, 20] #, 100] #, 1000, 20000, 100000]
# n_rep_list = [100, 1000, 20000, 100000, 500000]
t_per_exec_f = []
t_per_exec_s = []
H_zero = H_i*0.0
for i, n_reps in enumerate(n_rep_list):
    out, t_func, fname = repeat_solve(H_func_type, init_state, tlist,
                                      num_reps=n_reps,
                                      e_ops=meas, args=args)
    t_per_exec_f.append(t_func / n_reps)
    #print("{} execution of func type took {} seconds.".format(n_reps, t_func))
    # twisted method of making the code change to force new hash and
    # hence recompile
    key = 'nreps{}'.format(i)
    args[key] = n_reps
    H = list(H_str_type)
    H.append([H_zero, key])
    out, t_func, fname = repeat_solve(H, init_state, tlist,
                                      num_reps=n_reps,
                                      e_ops=meas, args=args)
    #print("{} execution of string type took {} seconds.".format(n_reps, t_func))
    t_per_exec_s.append(t_func / n_reps)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(n_rep_list, t_per_exec_f, 'o', label="func type")
ax.plot(n_rep_list, t_per_exec_s, 'x', label="string type")
ax.set_xlabel(r"$N_{\rm{reps}}$")
ax.set_ylabel("exec time per rep")
ax.set_title("Comparing Method Exec Time")
ax.legend()
n_rep_list = [1000, 5000, 10000, 15000, 20000, 50000, 100000]
t_per_exec_f = []
t_per_exec_s = []
H_zero = H_i*0.0
for i, n_reps in enumerate(n_rep_list):
    out, t_func, fname = repeat_solve(H_func_type, init_state, tlist,
                                      num_reps=n_reps,
                                      e_ops=meas, args=args)
    t_per_exec_f.append(t_func / n_reps)
    #print("{} execution of func type took {} seconds.".format(n_reps, t_func))
    # twisted method of making the code change to force new hash and
    # hence recompile
    key = 'nreps{}'.format(i)
    args[key] = n_reps
    H = list(H_str_type)
    H.append([H_zero, key])
    out, t_func, fname = repeat_solve(H, init_state, tlist,
                                      num_reps=n_reps,
                                      e_ops=meas, args=args)
    #print("{} execution of string type took {} seconds.".format(n_reps, t_func))
    t_per_exec_s.append(t_func / n_reps)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(n_rep_list, t_per_exec_f, 'o', label="func type")
ax.plot(n_rep_list, t_per_exec_s, 'x', label="string type")
ax.set_xlabel(r"$N_{\rm{reps}}$")
ax.set_ylabel("exec time per rep")
ax.set_title("Comparing Method Exec Time")
ax.legend()
Explanation: Comparing execution times
End of explanation
<END_TASK>