Source code for qutip.rhs_generate

# This file is part of QuTiP: Quantum Toolbox in Python.
#
#    Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson.
#    All rights reserved.
#
#    Redistribution and use in source and binary forms, with or without
#    modification, are permitted provided that the following conditions are
#    met:
#
#    1. Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#
#    2. Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#
#    3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names
#       of its contributors may be used to endorse or promote products derived
#       from this software without specific prior written permission.
#
#    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
#    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
#    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
#    PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
#    HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
#    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
#    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
#    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
#    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
#    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
#    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
###############################################################################

__all__ = ['rhs_generate', 'rhs_clear']

import os
import numpy as np
from types import FunctionType, BuiltinFunctionType
from functools import partial

from qutip.cy.codegen import Codegen
from qutip.solver import Options, config, solver_safe
from qutip.qobj import Qobj
from qutip.superoperator import spre, spost
from qutip.interpolate import Cubic_Spline


[docs]def rhs_clear(): """ Resets the string-format time-dependent Hamiltonian parameters. Parameters ---------- Returns ------- Nothing, just clears data from internal config module. """ # time-dependent (TD) function stuff config.tdfunc = None # Placeholder for TD RHS function. config.colspmv = None # Placeholder for TD col-spmv function. config.colexpect = None # Placeholder for TD col_expect function. config.string = None # Holds string of variables to be passed to solver config.tdname = None # Name of td .pyx file (used in parallel mc code) if "sesolve" in solver_safe: del solver_safe["sesolve"] if "mesolve" in solver_safe: del solver_safe["mesolve"] if "mcsolve" in solver_safe: del solver_safe["mcsolve"]
[docs]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'): """ Checks on time-dependent format. """ h_const = [] h_func = [] h_str = [] h_obj = [] # check H for incorrect format if isinstance(H, Qobj): pass elif isinstance(H, (FunctionType, BuiltinFunctionType, partial)): pass # n_func += 1 elif isinstance(H, list): for k, H_k in enumerate(H): if isinstance(H_k, Qobj): h_const.append(k) elif isinstance(H_k, list): if len(H_k) != 2 or not isinstance(H_k[0], Qobj): raise TypeError("Incorrect hamiltonian specification") else: if isinstance(H_k[1], (FunctionType, BuiltinFunctionType, partial)): h_func.append(k) elif isinstance(H_k[1], str): h_str.append(k) elif isinstance(H_k[1], Cubic_Spline): h_obj.append(k) elif hasattr(H_k[1], '__call__'): h_obj.append(k) elif isinstance(H_k[1], np.ndarray): h_str.append(k) else: raise TypeError("Incorrect hamiltonian specification") else: raise TypeError("Incorrect hamiltonian specification") # the the whole thing again for c_ops c_const = [] c_func = [] c_str = [] c_obj = [] if isinstance(c_ops, list): for k in range(len(c_ops)): if isinstance(c_ops[k], Qobj): c_const.append(k) elif isinstance(c_ops[k], list): if len(c_ops[k]) != 2: raise TypeError( "Incorrect collapse operator specification.") else: if isinstance(c_ops[k][1], (FunctionType, BuiltinFunctionType, partial)): c_func.append(k) elif isinstance(c_ops[k][1], str): c_str.append(k) elif isinstance(c_ops[k][1], Cubic_Spline): c_obj.append(k) elif hasattr(c_ops[k][1], '__call__'): c_func.append(k) elif isinstance(c_ops[k][1], np.ndarray): c_str.append(k) elif isinstance(c_ops[k][1], tuple): c_str.append(k) else: raise TypeError( "Incorrect collapse operator specification") else: raise TypeError("Incorrect collapse operator specification") # # if n_str == 0 and n_func == 0: # # no time-dependence at all # if ((len(h_str) > 0 and len(h_func) > 0) or (len(c_str) > 0 and len(c_func) > 0)): raise TypeError( "Cannot mix string and function type time-dependence formats") # check to see if Cython is installed and version is high enough. if len(h_str) > 0 or len(c_str) > 0: pass #try: # import Cython #except: # raise Exception( # "Unable to load Cython. Use Python function format.") #else: # if Cython.__version__ < '0.21': # raise Exception("Cython version (%s) is too old. Upgrade to" + # " version 0.21+" % Cython.__version__) # If only time-dependence is in Objects, then prefer string based format if (len(h_func) + len(c_func) + len(h_str) + len(c_str)) == 0: h_str += h_obj #Does nothing if not objects c_str += c_obj else: # Combine Hamiltonian objects if len(h_func) > 0: h_func += h_obj elif len(h_str) > 0: h_str += h_obj #Combine collapse objects if len(c_func) > 0: c_func += c_obj elif len(c_str) > 0: c_str += c_obj if solver == 'me': return (len(h_const + c_const), len(h_func) + len(c_func), len(h_str) + len(c_str)) elif solver == 'mc': # H C_ops # # -- ----- -- # NO NO 00 # NO STR 01 # NO FUNC 02 # # STR NO 10 # STR STR 11 # # FUNC NO 20 # # FUNC FUNC 22 if isinstance(H, FunctionType): time_type = 3 # Time-indepdent problems elif ((len(h_func) == len(h_str) == 0) and (len(c_func) == len(c_str) == 0)): time_type = 0 # constant Hamiltonian, time-dependent collapse operators elif len(h_func) == len(h_str) == 0: if len(c_str) > 0: time_type = 1 elif len(c_func) > 0: time_type = 2 else: raise Exception("Error determining time-dependence.") # list style Hamiltonian elif len(h_str) > 0: if len(c_func) == len(c_str) == 0: time_type = 10 elif len(c_str) > 0: time_type = 11 else: raise Exception("Error determining time-dependence.") # Python function style Hamiltonian elif len(h_func) > 0: if len(c_func) == len(c_str) == 0: time_type = 20 elif len(c_func) > 0: time_type = 22 else: raise Exception("Error determining time-dependence.") return time_type, [h_const, h_func, h_str], [c_const, c_func, c_str] def _td_wrap_array_str(H, c_ops, args, times): """ Wrap numpy-array based time-dependence in the string-based time dependence format """ n = 0 H_new = [] c_ops_new = [] args_new = {} if not isinstance(H, list): H_new = H else: for Hk in H: if isinstance(Hk, list) and isinstance(Hk[1], np.ndarray): H_op, H_td = Hk td_array_name = "_td_array_%d" % n H_td_str = '(0 if (t > %f) else %s[int(round(%d * (t/%f)))])' %\ (times[-1], td_array_name, len(times) - 1, times[-1]) args_new[td_array_name] = H_td H_new.append([H_op, H_td_str]) n += 1 else: H_new.append(Hk) if not isinstance(c_ops, list): c_ops_new = c_ops else: for ck in c_ops: if isinstance(ck, list) and isinstance(ck[1], np.ndarray): c_op, c_td = ck td_array_name = "_td_array_%d" % n c_td_str = '(0 if (t > %f) else %s[int(round(%d * (t/%f)))])' %\ (times[-1], td_array_name, len(times) - 1, times[-1]) args_new[td_array_name] = c_td c_ops_new.append([c_op, c_td_str]) n += 1 else: c_ops_new.append(ck) if not args_new: args_new = args elif isinstance(args, dict): args_new.update(args) else: raise ValueError("Time-dependent array format requires args to " + "be a dictionary") return H_new, c_ops_new, args_new