## QobjEvo "advanced" methods¶

Test and present diverse feature of QobjEvo.

Made by Eric Giguere

In [1]:
import qutip as qt
import numpy as np
from qutip import QobjEvo

In [2]:
N = 5
destroy, create, Id = qt.destroy(N), qt.create(N), qt.qeye(N)
def exp_i(t,args):
return np.exp(-1j*t)
def cos_w(t,args):
return np.cos(args["w"]*t)
tlist = np.linspace(0,10,10000)
tlistlog = np.logspace(-3,1,10000)

# state vector as np array
vec = np.arange(N)*.5+.5j
vec_super = np.arange(N**2)*.5+.5j
mat_c = (np.arange(N**2)*.5+.5j).reshape((5,5))
mat_f = np.asfortranarray(mat_c*1.)

# Construct QobjEvo of all type
td_cte1 = QobjEvo(Id)
td_cte2 = QobjEvo([Id])

td_func = QobjEvo([Id,[create,exp_i],[destroy,cos_w]],args={"w":2})
td_str = QobjEvo([Id,[create,"exp(-1j*t)"],[destroy,"cos(w*t)"]],args={'w':2.})
td_array = QobjEvo([Id,[create,np.exp(-1j*tlist)],[destroy,np.cos(2*tlist)]],tlist=tlist)
td_array_log = QobjEvo([Id,[create,np.exp(-1j*tlistlog)],[destroy,np.cos(2*tlistlog)]],tlist=tlistlog)

td_super = qt.liouvillian(td_func, c_ops=td_cte1)


## Product and expectation value¶

QobjEvo.mul_vec(t,state) = spmv(QobjEvo(t), state)
QobjEvo.expect(t, state, real) = cy_expect_psi/cy_expect_rho_vec (QobjEvo(t), state, real)

In [3]:
from qutip.cy.spmatfuncs import spmv, cy_expect_rho_vec, cy_expect_psi

In [4]:
spmv(td_func(2).data, vec) == td_func.mul_vec(2,vec)

Out[4]:
array([ True,  True,  True,  True,  True])
In [5]:
print(td_func(2).data * mat_c == td_func.mul_mat(2,mat_c))
mat_c.flags

[[ True  True  True  True  True]
[ True  True  True  True  True]
[ True  True  True  True  True]
[ True  True  True  True  True]
[ True  True  True  True  True]]

Out[5]:
  C_CONTIGUOUS : True
F_CONTIGUOUS : False
OWNDATA : False
WRITEABLE : True
ALIGNED : True
WRITEBACKIFCOPY : False
UPDATEIFCOPY : False
In [6]:
print(td_func(2).data * mat_f == td_func.mul_mat(2,mat_f))
mat_f.flags

[[ True  True  True  True  True]
[ True  True  True  True  True]
[ True  True  True  True  True]
[ True  True  True  True  True]
[ True  True  True  True  True]]

Out[6]:
  C_CONTIGUOUS : False
F_CONTIGUOUS : True
OWNDATA : True
WRITEABLE : True
ALIGNED : True
WRITEBACKIFCOPY : False
UPDATEIFCOPY : False
In [7]:
cy_expect_psi(td_str(2).data, vec, 0) == td_str.expect(2, vec, 0)

Out[7]:
True
In [8]:
cy_expect_rho_vec(td_super(2).data, vec_super, 0) == td_super.expect(2, vec_super, 0)

Out[8]:
True

## Function with state¶

For the solvers which have the option "rhs_with_state"

the QobjEvo can take coefficient function with the signature (for backward compatibility):

def coeff(t, psi. args)

or use advanced args: args={"psi=vec":psi0}

In [9]:
def coeff_state(t, args):
return np.max(args["psi"])*args["w"]
td_state = QobjEvo([Id, [destroy, coeff_state]],args={"w":1, "psi=vec":qt.basis(N,0)})

In [10]:
td_state(2,state=vec)

Out[10]:
Quantum object: dims = [[5], [5]], shape = (5, 5), type = oper, isherm = False\begin{equation*}\left(\begin{array}{*{11}c}1.0 & (2.0+0.500j) & 0.0 & 0.0 & 0.0\\0.0 & 1.0 & (2.828+0.707j) & 0.0 & 0.0\\0.0 & 0.0 & 1.0 & (3.464+0.866j) & 0.0\\0.0 & 0.0 & 0.0 & 1.0 & (4.0+1.0j)\\0.0 & 0.0 & 0.0 & 0.0 & 1.0\\\end{array}\right)\end{equation*}

## Cython object and "Compiling"¶

There is a cython version of the qobjevo for fast call:

• qutip.cy.cqobjevo.CQobjEvo

The cython is created when the "compile" method is created. The cython object contain fast version of the call, expect and rhs (spmv) methods.

In [11]:
td_str.compiled = False
print("Before compilation")
%timeit td_str(2, data=True)
%timeit td_str.mul_vec(2,vec)
%timeit td_str.mul_mat(2,mat_c)
%timeit td_str.mul_mat(2,mat_f)
%timeit td_str.expect(2,vec,0)
td_str.compile()
print("After compilation")
%timeit td_str(2, data=True)
%timeit td_str.mul_vec(2,vec)
%timeit td_str.mul_mat(2,mat_c)
%timeit td_str.mul_mat(2,mat_f)
%timeit td_str.expect(2,vec,0)

Before compilation
195 µs ± 8.68 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
205 µs ± 1.98 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
219 µs ± 15.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
268 µs ± 41.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
206 µs ± 1.34 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
After compilation
13.2 µs ± 477 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
8.53 µs ± 52.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
8.66 µs ± 48.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
9.18 µs ± 49.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
10.7 µs ± 288 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)


## apply¶

Pass a function ( Qobj, *args, **kwargs) -> Qobj to act on each component of the qobjevo.

Will only be mathematicaly valid if the transformation is linear.

In [12]:
def multiply(qobj, b, factor = 3.):
return qobj*b*factor

print(td_func.apply(multiply,2)(2) == td_func(2)*6)
print(td_func.apply(multiply,2,factor=2)(2) == td_func(2)*4)

True
True


## apply_decorator¶

Transform the functions containing the time dependence using a decorator.

The decorator must return a function of (t, **kwargs).

Do not modify the constant part (the contant part do not have a function f(t) = 1).

In [13]:
def rescale_time_and_scale(f_original, time_scale, factor=2.):
def f(t, *args, **kwargs):
return f_original(time_scale*t, *args, **kwargs)*factor
return f

print(td_func.apply_decorator(rescale_time_and_scale, 2)(2) == td_func(4)*2-Id)
print(td_func.apply_decorator(rescale_time_and_scale, 3, factor=3)(2) ==
td_func(6)*3.0 - 2*Id)

True
True


QobjEvo based on string and np.array are changed to a function then the decorator is applied. There are option so that the type of coefficient stay unchanged:

str_mod : change the string -> str_mod[0] + str + str_mod[1]

inplace_np : modify the array  (array[i] = decorator(lambda v: v)(array[i]))
*any modification that rescale the time will not work properly



Decorator can cause problem when used in parallel. (function cannot be pickled error)

In [14]:
td_func_1 = QobjEvo([[create,exp_i]],args={"w":2})
td_str_1 = QobjEvo([[create,"exp(-1j*t)"]],args={'w':2.})
td_array_1 = QobjEvo([[create,np.exp(-1j*tlist)]],tlist=tlist)

def square_qobj(qobj):
return qobj*qobj

def square_f(f_original):
def f(t, *args, **kwargs):
return f_original(t, *args, **kwargs)**2
return f

t1 = td_func_1.apply(square_qobj).apply_decorator(square_f)
print(t1(2) == td_func_1(2)*td_func_1(2))
print((t1.ops[0][2]))

t1 = td_str_1.apply(square_qobj).apply_decorator(square_f)
print(t1(2) == td_str_1(2)*td_str_1(2))
print("str not updated:", (t1.ops[0][2]))

t1 = td_str_1.apply(square_qobj).apply_decorator(square_f, str_mod=["(",")**2"])
print(t1(2) == td_str_1(2)*td_str_1(2))
print("str  updated:",(t1.ops[0][2]))

t1 = td_array_1.apply(square_qobj).apply_decorator(square_f)
print(t1(2) == td_array_1(2)*td_array_1(2))
print("array not updated:",(t1.ops[0][2]))

t1 = td_array_1.apply(square_qobj).apply_decorator(square_f, inplace_np=1)
print(t1(2) == td_array_1(2)*td_array_1(2))
print("array  updated:",(t1.ops[0][2]))

True
<function exp_i at 0x7f1c1e7bb158>
True
str not updated: <function square_f.<locals>.f at 0x7f1c1dee0c80>
True
str  updated: (exp(-1j*t))**2
True
array not updated: <function square_f.<locals>.f at 0x7f1c1dee09d8>
True
array  updated: [1.        -0.j         0.999998  -0.0020002j  0.999992  -0.00400039j ...
0.41173093-0.91130546j 0.40990732-0.91212718j 0.40808206-0.91294525j]


## Removing redundant Qobj¶

When multiple components of the QobjEvo are made from the same Qobj, you can unite them with the "compress" method. It is only done with the same form of time dependance:

In [15]:
small = qt.destroy(2)
def f1(t,args):
return np.sin(t)
def f2(t,args):
return np.cos(args["w"]*t)
def f3(t,args):
return np.sin(args["w"]*t)
def f4(t,args):
return np.cos(t)

In [16]:
td_redoundance = QobjEvo([qt.qeye(2),[small,"sin(t)"],[small,"cos(w*t)"],[small,"sin(w*t)"],
[small,"cos(t)"]],args={'w':2.})
td_redoundance1 = QobjEvo([qt.qeye(2),[small,"sin(t)"],[small,"cos(w*t)"],[small,"sin(w*t)"],
[small,"cos(t)"]],args={'w':2.})
td_redoundance2 = QobjEvo([qt.qeye(2),[small,f1],[small,f2],[small,f3],[small,f4]],args={'w':2.})
td_redoundance3 = QobjEvo([qt.qeye(2),[small,np.sin(tlist)],[small,np.cos(2*tlist)],
[small,np.sin(2*tlist)],[small,np.cos(tlist)]],tlist=tlist)
td_redoundance4 = QobjEvo([qt.qeye(2),[small,f1],[small,"cos(w*t)"],
[small,np.sin(2*tlist)],[small,"cos(t)"]],args={'w':2.},tlist=tlist)

td_redoundance1.compress()
print(td_redoundance1.to_list())
print(len(td_redoundance1.ops))
print(td_redoundance(1.) == td_redoundance1(1.))
td_redoundance2.compress()
print(len(td_redoundance2.ops))
print(td_redoundance(1.) == td_redoundance2(1.))
td_redoundance3.compress()
print(len(td_redoundance3.ops))
print(td_redoundance(1.) == td_redoundance3(1.))
td_redoundance4.compress()
print(len(td_redoundance4.ops))
print(td_redoundance(1.) == td_redoundance4(1.))

[Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
Qobj data =
[[1. 0.]
[0. 1.]], [Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = False
Qobj data =
[[0. 1.]
[0. 0.]], '(sin(t)) + (cos(w*t)) + (sin(w*t)) + (cos(t))']]
1
True
1
True
1
True
3
True

In [17]:
td_redoundance_v2 = QobjEvo([qt.qeye(2),[qt.qeye(2),"sin(t)"],[small,"sin(t)"],[qt.create(2),"sin(t)"]])
td_redoundance_v2.compress()
td_redoundance_v2.to_list()

Out[17]:
[Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
Qobj data =
[[1. 0.]
[0. 1.]],
[Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
Qobj data =
[[1. 1.]
[1. 1.]], 'sin(t)']]

## cimport cython object¶

The cython object can be 'cimported'.

cdef class CQobjEvo:
cdef void _mul_vec(self, double t, complex* vec, complex* out)
cdef void _mul_matf(self, double t, complex* mat, complex* out,int nrow, int ncols)
cdef void _mul_matc(self, double t, complex* mat, complex* out,int nrow, int ncols)
cdef complex _expect(self, double t, complex* vec, int isherm)
cdef complex _expect_super(self, double t, complex* rho, int isherm)
In [18]:
%%cython
from qutip.cy.cqobjevo cimport CQobjEvo
cimport numpy as np

def rhs_call_from_cy(CQobjEvo qobj, double t, np.ndarray[complex, ndim=1] vec, np.ndarray[complex, ndim=1] out):
qobj._mul_vec(t,&vec[0],&out[0])

def expect_call_from_cy(CQobjEvo qobj, double t, np.ndarray[complex, ndim=1] vec, int isherm):
return qobj._expect(t,&vec[0])

def rhs_cdef_timing(CQobjEvo qobj, double t, np.ndarray[complex, ndim=1] vec, np.ndarray[complex, ndim=1] out):
cdef int i
for i in range(10000):
qobj._mul_vec(t,&vec[0],&out[0])

def expect_cdef_timing(CQobjEvo qobj, double t, np.ndarray[complex, ndim=1] vec, int isherm):
cdef complex aa = 0.
cdef int i
for i in range(10000):
aa = qobj._expect(t, &vec[0])
return aa

def rhs_def_timing(qobj, double t, np.ndarray[complex, ndim=1] vec, complex[::1] out):
cdef int i
for i in range(10000):
out = qobj.mul_vec(t,vec)

def expect_def_timing(qobj, double t, np.ndarray[complex, ndim=1] vec, int isherm):
cdef complex aa = 0.
cdef int i
for i in range(10000):
aa = qobj.expect(t, vec)
return aa

In [19]:
td_str.compile()
print(expect_call_from_cy(td_str.compiled_qobjevo, 2, vec, 0) - td_str.expect(2,vec,0))
%timeit expect_def_timing(td_str.compiled_qobjevo, 2, vec, 0)
%timeit expect_cdef_timing(td_str.compiled_qobjevo, 2, vec, 0)

0j
56.1 ms ± 447 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
43.4 ms ± 365 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [20]:
out = np.zeros(N,dtype=np.complex128)
rhs_call_from_cy(td_str.compiled_qobjevo, 2, vec, out)
print( [a - b for a,b in zip(out, td_str.mul_vec(2,vec))])
%timeit rhs_def_timing(td_str.compiled_qobjevo, 2, vec, out)
%timeit rhs_cdef_timing(td_str.compiled_qobjevo, 2, vec, out)
# Most of the time gained is from allocating the out vector, not the

[0j, 0j, 0j, 0j, 0j]
64.6 ms ± 2.26 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
11 ms ± 24.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [21]:
td_cte = QobjEvo([Id])
td_cte.compile()
out = np.zeros(N,dtype=np.complex128)
rhs_call_from_cy(td_cte.compiled_qobjevo, 2, vec, out)
print( [a - b for a,b in zip(out, td_cte.mul_vec(2,vec))])
%timeit rhs_def_timing(td_cte.compiled_qobjevo, 2, vec, out)
%timeit rhs_cdef_timing(td_cte.compiled_qobjevo, 2, vec, out)

[0j, 0j, 0j, 0j, 0j]
39.7 ms ± 209 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
635 µs ± 1.07 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)


# Compiled string code¶

In [22]:
td_str.compiled = False
print(td_str.compile(code=True))

#!python
#cython: language_level=3
# This file is generated automatically by QuTiP.

import numpy as np
cimport numpy as np
import scipy.special as spe
cimport cython
np.import_array()
cdef extern from "numpy/arrayobject.h" nogil:
void PyDataMem_NEW_ZEROED(size_t size, size_t elsize)
void PyArray_ENABLEFLAGS(np.ndarray arr, int flags)
from qutip.cy.spmatfuncs cimport spmvpy
from qutip.cy.inter cimport _spline_complex_t_second, _spline_complex_cte_second
from qutip.cy.inter cimport _spline_float_t_second, _spline_float_cte_second
from qutip.cy.interpolate cimport (interp, zinterp)
from qutip.cy.cqobjevo_factor cimport StrCoeff
from qutip.cy.cqobjevo cimport CQobjEvo
from qutip.cy.math cimport erf, zerf
from qutip.qobj import Qobj
cdef double pi = 3.14159265358979323

include '/home/eric/anaconda3/lib/python3.7/site-packages/qutip-4.4.0.dev0+8b414dd-py3.7-linux-x86_64.egg/qutip/cy/complex_math.pxi'

cdef class CompiledStrCoeff(StrCoeff):
cdef double w

def set_args(self, args):
self.w=args['w']

cdef void _call_core(self, double t, complex * coeff):
cdef double w = self.w

coeff[0] = exp(-1j*t)
coeff[1] = cos(w*t)


In [23]:
qt.about()

QuTiP: Quantum Toolbox in Python
Copyright (c) 2011 and later.
A. J. Pitchford, P. D. Nation, R. J. Johansson, A. Grimsmo, and C. Granade

QuTiP Version:      4.4.0.dev0+8b414dd
Numpy Version:      1.16.2
Scipy Version:      1.2.1
Cython Version:     0.29.6
Matplotlib Version: 3.0.1
Python Version:     3.7.0
Number of CPUs:     4
BLAS Info:          OPENBLAS
OPENMP Installed:   False
INTEL MKL Ext:      False
Platform Info:      Linux (x86_64)
Installation path:  /home/eric/anaconda3/lib/python3.7/site-packages/qutip-4.4.0.dev0+8b414dd-py3.7-linux-x86_64.egg/qutip
==============================================================================
For your convenience a bibtex reference can be easily generated using qutip.cite()