of the form
$$ \frac{\textrm{d} x}{\textrm{d} t} = L(t) \; x $$$x$ -- vector:
or $x$ -- matrix:
and $L$ -- sparse matrix that drives the dynamics.
Qobj
fast_csr_matrix
-- custom to QuTiP.L = [L0, [L1, g1], [L2, g2] ....]
L
-- Qobj
g
-- Python function, e.g. def g(t, args): return np.cos(args['w']*t)
"cos(w*t)"
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
.pyx
filecqobjevo_compiled_coeff_<hash>.pyx
QobjEvo
exec(compile('from ' + filename + ' import CompiledStrCoeff, locals()))
.pyx
deletedQobjEvo
can be reused$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)} $$# 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()
You are in the installation directory. Change directories before running QuTiP.
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()
<matplotlib.legend.Legend at 0x7fa8896cc9d0>
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]
@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
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")
1 execution of func type took 0.06014871597290039 seconds.
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")
1 execution of string type took 4.107965469360352 seconds.
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()
<matplotlib.legend.Legend at 0x7fa881966bd0>
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()
<matplotlib.legend.Legend at 0x7fa881e85510>
.pyx
file?