Skip to content

Commit

Permalink
Merge pull request #2207 from Ericgig/misc.cy_warnings
Browse files Browse the repository at this point in the history
Fix import with cython 3
  • Loading branch information
Ericgig committed Jul 31, 2023
2 parents a2a289e + e53e5e1 commit c2d7eaa
Show file tree
Hide file tree
Showing 10 changed files with 28 additions and 238 deletions.
2 changes: 1 addition & 1 deletion doc/apidoc/functions.rst
Original file line number Diff line number Diff line change
Expand Up @@ -210,7 +210,7 @@ Time-dependent problems
-----------------------

.. automodule:: qutip.rhs_generate
:members: rhs_generate, rhs_clear
:members: rhs_clear

Scattering in Quantum Optical Systems
-------------------------------------
Expand Down
1 change: 1 addition & 0 deletions doc/changes/2207.misc
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
Run in no cython mode with cython >=3.0.0
64 changes: 0 additions & 64 deletions doc/guide/dynamics/dynamics-time.rst
Original file line number Diff line number Diff line change
Expand Up @@ -385,70 +385,6 @@ When repeatedly simulating a system where only the time-dependent variables, or
The second call to :func:`qutip.mcsolve` does not reorganize the data, and in the case of the string format, does not recompile the Cython code. For the small system here, the savings in computation time is quite small, however, if you need to call the solvers many times for different parameters, this savings will obviously start to add up.


.. _time-parallel:

Running String-Based Time-Dependent Problems using Parfor
==========================================================

.. note:: This section covers a specialized topic and may be skipped if you are new to QuTiP.

In this section we discuss running string-based time-dependent problems using the :func:`qutip.parallel.parfor` function. As the :func:`qutip.mcsolve` function is already parallelized, running string-based time-dependent problems inside of parfor loops should be restricted to the :func:`qutip.mesolve` function only. When using the string-based format, the system Hamiltonian and collapse operators are converted into C code with a specific file name that is automatically genrated, or supplied by the user via the ``rhs_filename`` property of the :class:`qutip.solver.Options` class. Because the :func:`qutip.parallel.parfor` function uses the built-in Python multiprocessing functionality, in calling the solver inside a parfor loop, each thread will try to generate compiled code with the same file name, leading to a crash. To get around this problem you can call the :func:`qutip.rhs_generate` function to compile simulation into C code before calling parfor. You **must** then set the :class:`qutip.solver.Options` object ``rhs_reuse=True`` for all solver calls inside the parfor loop that indicates that a valid C code file already exists and a new one should not be generated. As an example, we will look at the Landau-Zener-Stuckelberg interferometry example that can be found in the notebook "Time-dependent master equation: Landau-Zener-Stuckelberg inteferometry" in the tutorials section of the QuTiP web site.

To set up the problem, we run the following code:

.. plot::
:context:
:nofigs:

delta = 0.1 * 2 * np.pi # qubit sigma_x coefficient
w = 2.0 * 2 * np.pi # driving frequency
T = 2 * np.pi / w # driving period
gamma1 = 0.00001 # relaxation rate
gamma2 = 0.005 # dephasing rate

eps_list = np.linspace(-10.0, 10.0, 51) * 2 * np.pi # epsilon
A_list = np.linspace(0.0, 20.0, 51) * 2 * np.pi # Amplitude

sx = sigmax(); sz = sigmaz(); sm = destroy(2); sn = num(2)

c_ops = [np.sqrt(gamma1) * sm, np.sqrt(gamma2) * sz] # relaxation and dephasing
H0 = -delta / 2.0 * sx
H1 = [sz, '-eps / 2.0 + A / 2.0 * sin(w * t)']
H_td = [H0, H1]
Hargs = {'w': w, 'eps': eps_list[0], 'A': A_list[0]}


where the last code block sets up the problem using a string-based Hamiltonian, and ``Hargs`` is a dictionary of arguments to be passed into the Hamiltonian. In this example, we are going to use the :func:`qutip.propagator` and :func:`qutip.propagator.propagator_steadystate` to find expectation
values for different values of :math:`\epsilon` and :math:`A` in the
Hamiltonian :math:`H = -\frac{1}{2}\Delta\sigma_x -\frac{1}{2}\epsilon\sigma_z- \frac{1}{2}A\sin(\omega t)`.

We must now tell the :func:`qutip.mesolve` function, that is called by :func:`qutip.propagator` to reuse a
pre-generated Hamiltonian constructed using the :func:`qutip.rhs_generate` command:

.. plot::
:context:
:nofigs:

opts = Options(rhs_reuse=True)
rhs_generate(H_td, c_ops, Hargs, name='lz_func')

Here, we have given the generated file a custom name ``lz_func``, however this is not necessary as a generic name will automatically be given. Now we define the function ``task`` that is called by :func:`qutip.parallel.parfor` with the m-index parallelized in loop over the elements of ``p_mat[m,n]``:

.. code-block:: python
def task(args):
m, eps = args
p_mat_m = np.zeros(len(A_list))
for n, A in enumerate(A_list):
# change args sent to solver, w is really a constant though.
Hargs = {'w': w, 'eps': eps,'A': A}
U = propagator(H_td, T, c_ops, Hargs, opts) #<- IMPORTANT LINE
rho_ss = propagator_steadystate(U)
p_mat_m[n] = expect(sn, rho_ss)
return [m, p_mat_m]
Notice the Options ``opts`` in the call to the :func:`qutip.propagator` function. This tells the :func:`qutip.mesolve` function used in the propagator to call the pre-generated file ``lz_func``. If this was missing then the routine would fail.

.. plot::
:context: reset
:include-source: false
Expand Down
16 changes: 12 additions & 4 deletions qutip/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -55,15 +55,23 @@
else:
from qutip.utilities import _version2int
_cy_require = "0.29.20"
_cy_unsupported = "3.0.0"
if _version2int(_Cython.__version__) < _version2int(_cy_require):
warnings.warn(
"Old version of Cython detected: needed {}, got {}."
.format(_cy_require, _Cython.__version__)
)
# Setup pyximport
import qutip.cy.pyxbuilder as _pyxbuilder
_pyxbuilder.install()
del _pyxbuilder, _Cython, _version2int
if _version2int(_Cython.__version__) >= _version2int(_cy_unsupported):
warnings.warn(
"The new version of Cython, (>= 3.0.0) is not supported."
.format(_Cython.__version__)
)
else:
# Setup pyximport
import qutip.cy.pyxbuilder as _pyxbuilder
_pyxbuilder.install()
del _pyxbuilder, _Cython, _version2int
qutip.settings.has_cython = True


# -----------------------------------------------------------------------------
Expand Down
5 changes: 2 additions & 3 deletions qutip/cy/br_codegen.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,11 +27,10 @@ def __init__(self, h_terms=None, h_td_terms=None, h_obj=None,
omp_thresh=None,
omp_threads=None,
atol=None):
try:
import cython
except (ImportError, ModuleNotFoundError):
if not qset.has_cython:
raise ModuleNotFoundError("Cython is needed for "
"time-depdendent brmesolve")
import cython
import sys
import os
sys.path.append(os.getcwd())
Expand Down
2 changes: 1 addition & 1 deletion qutip/propagator.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
from qutip.qobj import Qobj
from qutip.tensor import tensor
from qutip.operators import qeye
from qutip.rhs_generate import (rhs_generate, rhs_clear, _td_format_check)
from qutip.rhs_generate import rhs_clear, _td_format_check
from qutip.superoperator import (vec2mat, mat2vec,
vector_to_operator, operator_to_vector)
from qutip.sparse import sp_reshape
Expand Down
5 changes: 3 additions & 2 deletions qutip/qobjevo.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,10 +32,11 @@
if sys.platform == 'win32':
safePickle[0] = True

try:

if qset.has_cython:
import cython
use_cython = [True]
except:
else:
use_cython = [False]


Expand Down
163 changes: 1 addition & 162 deletions qutip/rhs_generate.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
__all__ = ['rhs_generate', 'rhs_clear']
__all__ = ['rhs_clear']

import os
import numpy as np
Expand Down Expand Up @@ -37,167 +37,6 @@ def rhs_clear():
if "mcsolve" in solver_safe:
del solver_safe["mcsolve"]

def rhs_generate(H, c_ops, args={}, options=Options(), name=None,
cleanup=True):
"""
Generates the Cython functions needed for solving the dynamics of a
given system using the mesolve function inside a parfor loop.
Parameters
----------
H : qobj
System Hamiltonian.
c_ops : list
``list`` of collapse operators.
args : dict
Arguments for time-dependent Hamiltonian and collapse operator terms.
options : Options
Instance of ODE solver options.
name: str
Name of generated RHS
cleanup: bool
Whether the generated cython file should be automatically removed or
not.
Notes
-----
Using this function with any solver other than the mesolve function
will result in an error.
"""
config.reset()
config.options = options

if name:
config.tdname = name
else:
config.tdname = "rhs" + str(os.getpid()) + str(config.cgen_num)

Lconst = 0

Ldata = []
Linds = []
Lptrs = []
Lcoeff = []
Lobj = []

# loop over all hamiltonian terms, convert to superoperator form and
# add the data of sparse matrix represenation to

msg = "Incorrect specification of time-dependence: "

for h_spec in H:
if isinstance(h_spec, Qobj):
h = h_spec

if not isinstance(h, Qobj):
raise TypeError(msg + "expected Qobj")

if h.isoper:
Lconst += -1j * (spre(h) - spost(h))
elif h.issuper:
Lconst += h
else:
raise TypeError(msg + "expected operator or superoperator")

elif isinstance(h_spec, list):
h = h_spec[0]
h_coeff = h_spec[1]

if not isinstance(h, Qobj):
raise TypeError(msg + "expected Qobj")

if h.isoper:
L = -1j * (spre(h) - spost(h))
elif h.issuper:
L = h
else:
raise TypeError(msg + "expected operator or superoperator")

Ldata.append(L.data.data)
Linds.append(L.data.indices)
Lptrs.append(L.data.indptr)
if isinstance(h_coeff, Cubic_Spline):
Lobj.append(h_coeff.coeffs)
Lcoeff.append(h_coeff)

else:
raise TypeError(msg + "expected string format")

# loop over all collapse operators
for c_spec in c_ops:
if isinstance(c_spec, Qobj):
c = c_spec

if not isinstance(c, Qobj):
raise TypeError(msg + "expected Qobj")

if c.isoper:
cdc = c.dag() * c
Lconst += spre(c) * spost(c.dag()) - 0.5 * spre(cdc) \
- 0.5 * spost(cdc)
elif c.issuper:
Lconst += c
else:
raise TypeError(msg + "expected operator or superoperator")

elif isinstance(c_spec, list):
c = c_spec[0]
c_coeff = c_spec[1]

if not isinstance(c, Qobj):
raise TypeError(msg + "expected Qobj")

if c.isoper:
cdc = c.dag() * c
L = spre(c) * spost(c.dag()) - 0.5 * spre(cdc) \
- 0.5 * spost(cdc)
c_coeff = "(" + c_coeff + ")**2"
elif c.issuper:
L = c
else:
raise TypeError(msg + "expected operator or superoperator")

Ldata.append(L.data.data)
Linds.append(L.data.indices)
Lptrs.append(L.data.indptr)
Lcoeff.append(c_coeff)

else:
raise TypeError(msg + "expected string format")

# add the constant part of the lagrangian
if Lconst != 0:
Ldata.append(Lconst.data.data)
Linds.append(Lconst.data.indices)
Lptrs.append(Lconst.data.indptr)
Lcoeff.append("1.0")

# the total number of liouvillian terms (hamiltonian terms + collapse
# operators)
n_L_terms = len(Ldata)

cgen = Codegen(h_terms=n_L_terms, h_tdterms=Lcoeff, args=args,
config=config)
cgen.generate(config.tdname + ".pyx")

code = compile('from ' + config.tdname +
' import cy_td_ode_rhs', '<string>', 'exec')
exec(code, globals())

config.tdfunc = cy_td_ode_rhs

if cleanup:
try:
os.remove(config.tdname + ".pyx")
except:
pass


def _td_format_check(H, c_ops, solver='me'):
"""
Expand Down
2 changes: 2 additions & 0 deletions qutip/settings.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@
has_mkl = False
# Has OPENMP
has_openmp = False
# Has usable cython
has_cython = False
# debug mode for development
debug = False
# Running on mac with openblas make eigh unsafe
Expand Down
6 changes: 5 additions & 1 deletion qutip/tests/conftest.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
import os
import tempfile
import numpy as np
from qutip.utilities import _version2int


def _add_repeats_if_marked(metafunc):
Expand All @@ -27,7 +28,10 @@ def _skip_cython_tests_if_unavailable(item):
if item.get_closest_marker("requires_cython"):
# importorskip rather than mark.skipif because this way we get pytest's
# version-handling semantics.
pytest.importorskip('Cython', minversion='0.14')
_Cython = pytest.importorskip('Cython', minversion='0.14')
# importorskip does not have maxversion
if _version2int(_Cython.__version__) >= _version2int("3.0.0"):
pytest.skip("cython 3.0.0 not supported")


@pytest.hookimpl(trylast=True)
Expand Down

0 comments on commit c2d7eaa

Please sign in to comment.