Cutiepy*

"Cu" stands for Quantum

"tie" stands for Tools

</div> * provisional name

Solving Linear ODEs

You need two tools:

Linear Algebra Libraries
ODE solvers

Numpy for LinAlg, Scipy for ODE solving

zvode solver used internally (1980s-1990s)

small (N<20)	big (N>100)	time dependent
Bad - Severe overhead from python calls	Depends - Оverhead from memory allocation (can be worked around)	Depends - Any python function works, but python functions are very slow



In [ ]:

    
ts = np.linspace(...)
op = -1j*np.array(...)
f_prime = lambda t, state: op*(... t).dot(state)
state0_num = np.array(...)                                             
ode = scipy.integrate.ode(f_prime) ...
res = [state0_num]                                                              
for t in ts[1:]:                                                             
    ode.integrate(t)                                                            
    res.append(ode.y)                                                           
print(res[-1])

Qutip

zvode (from Scipy) is used as well, but linear algebra goes to cython (most of the time)

small (N<20)	big (N>100)	time dependent
Good - Noticeable overhead from cython/python calls	Good - Minor overhead from cython/python calls	Bad - Severely limited in scope and Inconsistent



In [ ]:

    
qutip.sesolve(H,init,ts,[])

if timedependent, H has to be:

a python function for sesolve
a list of lists (not tuples!?) for mesolve and mcsolve
only recent versions try to bring some consistency to the interface (not yet done)

To run $\hat{H} = \hat{H}_0 + \operatorname{sin}(t)\hat{H}_1$ you write



In [ ]:

    
H = [H0, [H1, 'sin(t)']]

Theano

no built-in solver, but you can use (again) zvode (from Scipy)
linear algebra is done internally with various optimizations (including in-place operations)
bad sparse matrix support (getting better)

small (N<20)	big (N>100)	time dependent
Bad - Severe overhead from theano/python calls	Good - Minor overhead from theano/python calls	Great - Symbolic optimizations (and in-place operations!)

Summary

Everybody uses the old zvode through scipy (if it has an ODE solver at all).
- bad because if the slow detour through python
Except Theano time dependence awkward if supported at all.

	small (N<20)	big (N>100)	time dependent
Numpy/Scipy	Bad - Severe overhead from python calls	Depends - Оverhead from memory allocation (can be worked around)	Depends - Any python function works, but python functions are very slow
Qutip	Good - Noticeable overhead from cython/python calls	Good - Minor overhead from cython/python calls	Bad - Severely limited in scope and Inconsistent
Theano	Bad - Severe overhead from theano/python calls	Good - Minor overhead from theano/python calls	Great - Symbolic optimizations (and in-place operations!)
Cutiepy	?	?	?

Cutiepy

Having short simple code was as important as making it fast.

code generator (a simple one, <500 lines of code)
symbolic algebra module (also a simple one, <1000 lines of code)

Symbolic Expressions



In [1]:

    
from cutiepy import *
k = Ket('psi', 4)
k









    Out[1]:





$$\text{Ket }{| {psi}_{} \rangle} \text{ on the space }\mathbb{C}^{4}\text{ without an attached numerical content.}$$



In [2]:

    
H = Operator('H', 4)
H









    Out[2]:





$$\text{Operator }{\hat{H}_{}} \text{ on the space }\mathbb{C}^{4}\text{ without an attached numerical content.}$$



In [3]:

    
tensor(H,H)









    Out[3]:





$${{\hat{H}_{}}\tiny\otimes\normalsize{\hat{H}_{}}}$$



In [4]:

    
k2 = Ket('psi_2', [2,2])
k2









    Out[4]:





$$\text{Ket }{| {psi}_{2} \rangle} \text{ on the space }\mathbb{C}^{2}\otimes\mathbb{C}^{2}\text{ without an attached numerical content.}$$



In [5]:

    
sin(2*t)*H*k









    Out[5]:





$${{\operatorname{sin}\left( {{{2}\tiny\times\normalsize{t}}} \right)}\tiny\times\normalsize{{\hat{H}_{}}{| {psi}_{} \rangle}}}$$

Numerical Content



In [6]:

    
import numpy as np
Ket('phi', 4, np.array([0,0,0,1], dtype=complex))









    Out[6]:





$$\text{Ket }{| {phi}_{} \rangle} \text{ on the space }\mathbb{C}^{4}\text{ with numerical content: }$$
$$\begin{equation*}\left(\begin{array}{*{11}c}0.0\\0.0\\0.0\\1.0\\\end{array}\right)\end{equation*}$$

Built-in Operators (I should work a bit more on this)



In [7]:

    
sigmax()









    Out[7]:





$$\text{Operator }{\hat{σ}_{x}} \text{ on the space }\mathbb{C}^{2}\text{ with numerical content: }$$
$$\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 1.0\\1.0 & 0.0\\\end{array}\right)\end{equation*}$$



In [8]:

    
num(3)









    Out[8]:





$$\text{Operator }{\hat{{n}}_{{3}}} \text{ on the space }\mathbb{C}^{3}\text{ with numerical content: }$$
$$\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0 & 0.0\\0.0 & 1.0 & 0.0\\0.0 & 0.0 & 2.0\\\end{array}\right)\end{equation*}$$

Lazy evaluation of expressions



In [9]:

    
expression = sigmax()*sigmay()
expression









    Out[9]:





$${{\hat{σ}_{x}}{\hat{σ}_{y}}}$$



In [10]:

    
evalf(expression)









    Out[10]:





$$\text{Anonymous }\text{Operator }{\hat{\tiny\boxed{{O}_{abb0e9ed...}}\normalsize}_{}} \text{ on the space }\mathbb{C}^{2}\text{ with numerical content: }$$
$$\begin{equation*}\left(\begin{array}{*{11}c}1.0j & 0.0\\0.0 & -1.0j\\\end{array}\right)\end{equation*}$$

Getting the numerical content



In [11]:

    
numerical(destroy(3))









    Out[11]:





array([[ 0.00000000+0.j,  1.00000000+0.j,  0.00000000+0.j],
       [ 0.00000000+0.j,  0.00000000+0.j,  1.41421356+0.j],
       [ 0.00000000+0.j,  0.00000000+0.j,  0.00000000+0.j]])

Code generation

All these lazy expressions can be compiled to C (with cython) for very fast evaluation.
Coupled to a modern ODE solver (in C, no python detour)
Hidden behind the sesolve, mesolve, etc. intefaces

Low Level Example

No need to worry about this in user code, but interesting if you want to help develop the library



In [12]:

    
y = Ket.anon(2)
op = -1j*num(2)
y_prime = op*y
y_prime









    Out[12]:





$${{-1j}\tiny\times\normalsize{{\hat{{n}}_{{2}}}{| {\tiny\boxed{{K}_{906adb3b...}}\normalsize}_{} \rangle}}}$$



In [13]:

    
f = generate_cython(y_prime, # Expression
                    NDArrayFunction(), # Type of compilation
                                       # (whether to couple to an ODE solver)
                    [y]) # Arguments for the function



In [14]:

    
codegen.DEBUG = True
cf = f.compiled()
codegen.DEBUG = False









    









  
  
  




#cython: boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True
cimport numpy as np
import numpy as np
from libc.math cimport sin, cos, exp, sinh, cosh, tan, tanh, log
from scipy.linalg.cython_blas cimport zgemv, zgemm

cdef int iONE=1, iZERO=0
cdef double complex zONE=1, zZERO=0, zNHALF=-0.5

# Declaration of global variables for
# - intermediate results
# - predefined constants
cdef double complex [:, :] C = np.empty((2, 1), dtype=np.complex)
cdef double complex [:, :] A # initialized in setup

# The generated cython code
cdef generated_function(
    # input arguments
    double complex [:, :] B,
    # output argument
    double complex [:, :] D
    ):
    # indices and intermediate values for various matrix operations
    cdef int i, j, k, ii, jj, kk
    cdef double complex c
    # intermediate results
    global C
    # predefined constants
    global A
    # evaluating the function
    # C = A.B
    ii = A.shape[0]
    jj = A.shape[1]
    kk = B.shape[1]
    zgemm('N', 'N', &kk, &ii, &jj, &zONE, &B[0,0], &kk, &A[0,0], &jj, &zZERO, &C[0,0], &kk)
    # D = -1j*C
    ii = C.shape[0]
    jj = C.shape[1]
    for i in range(ii):
        for j in range(jj):
            D[i,j] = -1j*C[i,j]

# Function to deliver the constants from python to cython
cpdef setup_generated_function(
    double complex [:, :] _A
    ):
    global A
    A = _A

cpdef np.ndarray[np.complex_t, ndim=2] pythoncall(_0):
    cdef np.ndarray[np.complex_t, ndim=2] result = np.empty((2, 1), complex)
    generated_function(_0, result)
    return result



In [15]:

    
arg = np.ones((2,1), dtype=complex)
arg









    Out[15]:





array([[ 1.+0.j],
       [ 1.+0.j]])



In [16]:

    
cf.pythoncall(arg)









    Out[16]:





array([[ 0.-0.j],
       [ 0.-1.j]])

High Level Example



In [17]:

    
initial_state = basis(2, 0)
initial_state









    Out[17]:





$$\text{Ket }{| {{0}}_{{\tiny N\normalsize 2}} \rangle} \text{ on the space }\mathbb{C}^{2}\text{ with numerical content: }$$
$$\begin{equation*}\left(\begin{array}{*{11}c}1.0\\0.0\\\end{array}\right)\end{equation*}$$



In [18]:

    
ω0 = 1
Δ = 0.002
Ω = 0.005
H = ω0/2 * sigmaz() + Ω * sigmax() * sin((ω0+Δ)*t)
H









    Out[18]:





$${\left( {{0.5}\tiny\times\normalsize{\hat{σ}_{z}}}+{{0.005}\tiny\times\normalsize{\operatorname{sin}\left( {{{1.002}\tiny\times\normalsize{t}}} \right)}\tiny\times\normalsize{\hat{σ}_{x}}} \right)}$$



In [19]:

    
ts = 2*np.pi/Ω*np.linspace(0,1,40)
res = sesolve(H, initial_state, ts)









    



Generating cython code...
Compiling cython code...
Running cython code...
Starting at 10/29 15:35:05.
Finishing at 10/29 15:35:05.
Total time: 0 seconds.
Formatting the output...



In [20]:

    
%matplotlib inline
import matplotlib.pyplot as plt
σz_expect = expect(sigmaz(), res)
plt.plot(ts*Ω/np.pi, σz_expect, 'r.', label='numerical result')
Ωp = (Ω**2+Δ**2)**0.5
plt.plot(ts*Ω/np.pi, 1-(Ω/Ωp)**2*2*np.sin(Ωp*ts/2)**2, 'b-',
         label=r'$1-2(\Omega^\prime/\Omega)^2\sin^2(\Omega^\prime t/2)$')
plt.title(r'$\langle\sigma_z\rangle$-vs-$t\Omega/\pi$ at '
          r'$\Delta/\Omega=%.2f$, $\omega_0/\Omega=%.2f$'%(Δ/Ω, ω0/Ω))
plt.ylim(-1,1)
plt.legend(loc=3);

Benchmarks against Qutip

Dense

Sparse

TODO (i.e. HELP!?)

Sparse matrices not fully supported
Simplectic/Geometric ODE solvers
- qutip uses a "handwavy" approach to preserving the norm. It works, but ...
The "details"
- built-in operators and states
- visualization
- others
Make compilation faster
- Compiling takes us ~6s. qutip does it in <1s.

	small (N<20)	big (N>100)	time dependent
Numpy/Scipy	Bad - Severe overhead from python calls	Depends - Оverhead from memory allocation (can be worked around)	Depends - Any python function works, but python functions are very slow
Qutip	Good - Noticeable overhead from cython/python calls	Good - Minor overhead from cython/python calls	Bad - Severely limited in scope and Inconsistent
Theano	Bad - Severe overhead from theano/python calls	Good - Minor overhead from theano/python calls	Great - Symbolic optimizations (and in-place operations!)
Cutiepy	?	?	?