<jupyter_start><jupyter_text>Benchmarks<jupyter_code># NBVAL_IGNORE_OUTPUT
%load_ext watermark
import os
import qutip
import numpy as np
from numpy import pi
import matplotlib
import matplotlib.pylab as plt
import newtonprop
import snakeviz
import numba
%watermark -v --iversions<jupyter_output>qutip       4.3.1
numpy       1.15.4
matplotlib  3.0.2
matplotlib.pylab  1.15.4
newtonprop  0.1.0
snakeviz    1.0.0
numba       0.41.0
CPython 3.6.7
IPython 7.2.0<jupyter_text>While only available interactively, `snakeviz` is a very good way to look at profiling data:<jupyter_code>%load_ext snakeviz<jupyter_output><empty_output><jupyter_text>The benchmarks are based on the Examples, so we rerun the entire Example notebook in this context:<jupyter_code>%%capture
%run example.ipynb<jupyter_output><empty_output><jupyter_text>Propagation runtime<jupyter_code>tlist = np.linspace(0, 10, 100)<jupyter_output><empty_output><jupyter_text>First, we measure low long the propagation with QuTiP takes. This happens mostly in compiled code, so it is pretty fast.<jupyter_code>%%timeit
qutip.mesolve(L, rho0, tlist)
%%timeit
propagate_expm(L, rho0, tlist)<jupyter_output>39.1 ms ± 1.52 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)<jupyter_text>The Newton propagator, being a reference implementation, is implemented in pure Python, and is several orders of magnitude slower. To make the comparison fair, we limit the precision to $10^{-8}$, which is roughly the precision of `mesolve`.<jupyter_code>%%timeit
propagate(L, rho0, tlist, zero_qutip, norm_qutip, inner_qutip, tol=1e-8)<jupyter_output>3.85 s ± 282 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)<jupyter_text>When using lower-level data types, things get considerably faster:<jupyter_code>%%timeit
propagate(apply_cythonized_L, rho0_data, tlist, zero_vectorized, norm_vectorized, inner_vectorized, tol=1e-8)
%%timeit
propagate(L_vectorized, rho0_vectorized, tlist, zero_vectorized, norm_vectorized, inner_vectorized, tol=1e-8)<jupyter_output>210 ms ± 14.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)<jupyter_text>Profiling We can profile the how much time is spent in the various routines, comparing `mesolve` and different variations of the Newton propagator. See https://docs.python.org/3/library/profile.htmlinstant-user-s-manual for the meaning of the colums.  mesolve First, we look the QuTiP's `mesolve`:<jupyter_code>stats = %prun -q -r qutip.mesolve(L, rho0, tlist);<jupyter_output><jupyter_text>We can look at which top-level routines we spent the most time in *cumulativly*, that is including sub-calls:<jupyter_code>stats.sort_stats('cumtime').print_stats(10);<jupyter_output>8378 function calls (8174 primitive calls) in 0.023 seconds

   Ordered by: cumulative time
   List reduced from 105 to 10 due to restriction <10>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.023    0.023 {built-in method builtins.exec}
        1    0.000    0.000    0.023    0.023 <string>:1(<module>)
        1    0.000    0.000    0.023    0.023 mesolve.py:80(mesolve)
        1    0.000    0.000    0.023    0.023 mesolve.py:775(_mesolve_const)
        1    0.001    0.001    0.023    0.023 mesolve.py:985(_generic_ode_solve)
       99    0.000    0.000    0.009    0.000 _ode.py:396(integrate)
       99    0.008    0.000    0.008    0.000 _ode.py:989(run)
      102    0.001    0.000    0.008    0.000 qobj.py:211(__init__)
      112    0.000    0.000    0.004    0.000 qobj.py:1927(type)
      310    0.001    0.000    0.004    0.000 fromnumeric.py:64(_wrapreduction)<jupyter_text>Or, the bottom-level routines where we *actually* spent time:<jupyter_code>stats.sort_stats('tottime').print_stats(10);<jupyter_output>8378 function calls (8174 primitive calls) in 0.023 seconds

   Ordered by: internal time
   List reduced from 105 to 10 due to restriction <10>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
       99    0.008    0.000    0.008    0.000 _ode.py:989(run)
      310    0.002    0.000    0.002    0.000 {method 'reduce' of 'numpy.ufunc' objects}
      100    0.002    0.000    0.003    0.000 {qutip.cy.spconvert.dense2D_to_fastcsr_fmode}
        1    0.001    0.001    0.023    0.023 mesolve.py:985(_generic_ode_solve)
      611    0.001    0.000    0.001    0.000 {built-in method numpy.core.multiarray.array}
      203    0.001    0.000    0.003    0.000 fastsparse.py:47(__init__)
      310    0.001    0.000    0.004    0.000 fromnumeric.py:64(_wrapreduction)
      100    0.001    0.000    0.001    0.000 superoperator.py:227(vec2mat)
      102    0.001    0.000    0.008    0.000 qobj.py:211(__init__)
      301    0.001    0.000    0.001    0.000 qobj.py:1934([...]<jupyter_text>This is dominated by the ODE solver and sparse matrix operations If we're working interactively, we could use `snakeviz` to analyze further details:<jupyter_code>#%snakeviz qutip.mesolve(L, rho0, tlist)<jupyter_output><empty_output><jupyter_text>qutip-propagation Next, we look at the Newton propagator operating on high-level qutip objects:<jupyter_code>#%snakeviz propagate(L, rho0, tlist, zero_qutip, norm_qutip, inner_qutip, tol=1e-8)
stats = %prun -q -r propagate(L, rho0, tlist, zero_qutip, norm_qutip, inner_qutip, tol=1e-8);
stats.sort_stats('cumtime').print_stats(10);
stats.sort_stats('tottime').print_stats(10);<jupyter_output>3902775 function calls (3810705 primitive calls) in 4.882 seconds

   Ordered by: internal time
   List reduced from 182 to 10 due to restriction <10>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
   110385    0.552    0.000    0.552    0.000 {method 'reduce' of 'numpy.ufunc' objects}
    74547    0.357    0.000    0.737    0.000 fastsparse.py:47(__init__)
   272547    0.289    0.000    0.289    0.000 {built-in method numpy.core.multiarray.array}
   109395    0.233    0.000    0.860    0.000 fromnumeric.py:64(_wrapreduction)
     7623    0.157    0.000    0.509    0.000 fastsparse.py:74(_binopt)
    35937    0.141    0.000    1.640    0.000 qobj.py:211(__init__)
     6435    0.139    0.000    0.217    0.000 {qutip.cy.spmath.zcsr_mult}
       99    0.137    0.001    3.991    0.040 propagator.py:26(_arnoldi)
    91476    0.122    0.000    1.106    0.000 dimensions.py:44(is_scalar)
   548262    0.117    0.000    0.162    0.000 {built-in method builtins.[...]<jupyter_text>cythonized qutip-propagation Lastly, we can look at the more efficient propagation using a cythonized application of the QuTiP objects:<jupyter_code>#%snakeviz propagate(apply_cythonized_L, rho0_data, tlist, zero_vectorized, norm_vectorized, inner_vectorized, tol=1e-8)
stats = %prun -q -r propagate(apply_cythonized_L, rho0_data, tlist, zero_vectorized, norm_vectorized, inner_vectorized, tol=1e-8);
stats.sort_stats('cumtime').print_stats(10);
stats.sort_stats('tottime').print_stats(10);<jupyter_output>95441 function calls in 0.289 seconds

   Ordered by: internal time
   List reduced from 72 to 10 due to restriction <10>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
       99    0.067    0.001    0.201    0.002 propagator.py:26(_arnoldi)
       99    0.041    0.000    0.284    0.003 propagator.py:314(_step)
      990    0.033    0.000    0.058    0.000 linalg.py:953(eigvals)
     3762    0.028    0.000    0.029    0.000 {built-in method numpy.core.multiarray.dot}
      990    0.016    0.000    0.016    0.000 {built-in method qutip.cy.spmatfuncs.cy_ode_rhs}
     1386    0.013    0.000    0.038    0.000 linalg.py:2203(norm)
     8613    0.011    0.000    0.014    0.000 defmatrix.py:186(__getitem__)
     5445    0.008    0.000    0.008    0.000 {built-in method numpy.core.multiarray.vdot}
       99    0.008    0.000    0.008    0.000 propagator.py:142(_extend_leja)
     1089    0.006    0.000    0.006    0.000 {method 'reduce' of 'numpy.ufunc' object[...]<jupyter_text>We now see that the (Python-level) arnoldi and step routines of Newton come out as the bottle neck, as the application of the Liouvillian has become efficient level (running at C speed). Note that the Newton implementation, in particular the `_extend_leja` function has been sped up through the use of numba. Without that, `_extend_leja` would dominate.  Number of function evaluations Since the Newton propagator is ultimately still limited by being implemented in Python, it is more fair to measure the runtime in terms of the number of average number of applications of the Liouvillian per time step. This is under the assumption that in an efficient implementation (and for large Hilbert spaces), this is the dominating factor.The number of applications depends on the chosen precision and on the length of the time step: longer time steps tend do be more efficient, as only then we're in a regime where the fast convergence of the Newton series kicks in. We construct a dummy `Qobj` that counts its own applications to a state:<jupyter_code>class CountingQobj(qutip.Qobj):
    def __init__(self, *args, **kwargs):
        self.counter = 0
        super().__init__(*args, **kwargs)
    def __mul__(self, other):
        if isinstance(other, qutip.Qobj):
            self.counter += 1
        return super().__mul__(other)
def count_applications_newton(nt, tol=1e-8):
    tlist = np.linspace(0, 10, nt)
    L_count = CountingQobj(L)
    propagate(L_count, rho0, tlist, zero_qutip, norm_qutip, inner_qutip, tol=1e-8)
    return L_count.counter / len(tlist)
count_applications_newton(nt=10)
count_applications_newton(nt=100)<jupyter_output><empty_output><jupyter_text>To compare this to the average number of applications in `mesolve`, we use a trimmed-down version of the `mesolve` routine:<jupyter_code>import scipy
from qutip.solver import Options
from qutip.superoperator import mat2vec
from qutip.cy.spmatfuncs import cy_ode_rhs
from qutip.mesolve import _generic_ode_solve
from qutip.ui.progressbar import BaseProgressBar


def mesolve(L, rho0, tlist):
    opt = Options()

    def func(t, rho, data, ind, ptr):
        func.counter += 1
        return cy_ode_rhs(t, rho, data, ind, ptr)

    func.counter = 0

    r = scipy.integrate.ode(func)
    r.set_f_params(L.data.data, L.data.indices, L.data.indptr)
    r.set_integrator('zvode', method=opt.method, order=opt.order,
                     atol=opt.atol, rtol=opt.rtol, nsteps=opt.nsteps,
                     first_step=opt.first_step, min_step=opt.min_step,
                     max_step=opt.max_step)
    initial_vector = mat2vec(rho0.full()).ravel('F')
    r.set_initial_value(initial_vector, tlist[0])
    dt = tlist[1] - tlist[0]
    for step in range(len(tlist)-1):
        r.integrate(r.t + dt)
    return func.counter


def count_applications_mesolve(nt):
    tlist = np.linspace(0, 10, nt)
    return mesolve(L, rho0, tlist) / len(tlist)
count_applications_mesolve(10)
count_applications_mesolve(100)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09575247162211645, 'fraction_numerical': 0.09675942878066643, 'mean_word_length': 3.3910771704180065, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 58, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '644608', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4501, 'n_tokens_neox': 3895, 'n_words': 1078, 'path': 'docs/benchmarks.ipynb', 'repo_name': 'qucontrol/newtonprop', 'size': '10924'}