1. qutip-notebooks
  2. development

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
%load_ext cython
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
==============================================================================
Please cite QuTiP in your publication.
==============================================================================
For your convenience a bibtex reference can be easily generated using `qutip.cite()`
In [ ]: