# 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.
###############################################################################
from __future__ import print_function
__all__ = ['Options', 'Odeoptions', 'Odedata']
import sys
import datetime
from collections import OrderedDict
import os
import warnings
from qutip import __version__
from qutip.qobj import Qobj
from types import FunctionType, BuiltinFunctionType
[docs]class Options():
"""
Class of options for evolution solvers such as :func:`qutip.mesolve` and
:func:`qutip.mcsolve`. Options can be specified either as arguments to the
constructor::
opts = Options(order=10, ...)
or by changing the class attributes after creation::
opts = Options()
opts.order = 10
Returns options class to be used as options in evolution solvers.
Attributes
----------
atol : float {1e-8}
Absolute tolerance.
rtol : float {1e-6}
Relative tolerance.
method : str {'adams','bdf'}
Integration method.
order : int {12}
Order of integrator (<=12 'adams', <=5 'bdf')
nsteps : int {2500}
Max. number of internal steps/call.
first_step : float {0}
Size of initial step (0 = automatic).
min_step : float {0}
Minimum step size (0 = automatic).
max_step : float {0}
Maximum step size (0 = automatic)
tidy : bool {True,False}
Tidyup Hamiltonian and initial state by removing small terms.
num_cpus : int
Number of cpus used by mcsolver (default = # of cpus).
norm_tol : float
Tolerance used when finding wavefunction norm in mcsolve.
norm_steps : int
Max. number of steps used to find wavefunction norm to within norm_tol
in mcsolve.
average_states : bool {False}
Average states values over trajectories in stochastic solvers.
average_expect : bool {True}
Average expectation values over trajectories for stochastic solvers.
mc_corr_eps : float {1e-10}
Arbitrarily small value for eliminating any divide-by-zero errors in
correlation calculations when using mcsolve.
ntraj : int {500}
Number of trajectories in stochastic solvers.
rhs_reuse : bool {False,True}
Reuse Hamiltonian data.
rhs_with_state : bool {False,True}
Whether or not to include the state in the Hamiltonian function
callback signature.
rhs_filename : str
Name for compiled Cython file.
seeds : ndarray
Array containing random number seeds for mcsolver.
store_final_state : bool {False, True}
Whether or not to store the final state of the evolution in the
result class.
store_states : bool {False, True}
Whether or not to store the state vectors or density matrices in the
result class, even if expectation values operators are given. If no
expectation are provided, then states are stored by default and this
option has no effect.
"""
def __init__(self, atol=1e-8, rtol=1e-6, method='adams', order=12,
nsteps=1000, first_step=0, max_step=0, min_step=0,
average_expect=True, average_states=False, tidy=True,
num_cpus=0, norm_tol=1e-3, norm_steps=5, rhs_reuse=False,
rhs_filename=None, ntraj=500, gui=False, rhs_with_state=False,
store_final_state=False, store_states=False, seeds=None,
steady_state_average=False, normalize_output=True):
# Absolute tolerance (default = 1e-8)
self.atol = atol
# Relative tolerance (default = 1e-6)
self.rtol = rtol
# Integration method (default = 'adams', for stiff 'bdf')
self.method = method
# Max. number of internal steps/call
self.nsteps = nsteps
# Size of initial step (0 = determined by solver)
self.first_step = first_step
# Minimal step size (0 = determined by solver)
self.min_step = min_step
# Max step size (0 = determined by solver)
self.max_step = max_step
# Maximum order used by integrator (<=12 for 'adams', <=5 for 'bdf')
self.order = order
# Average expectation values over trajectories (default = True)
self.average_states = average_states
# average expectation values
self.average_expect = average_expect
# Number of trajectories (default = 500)
self.ntraj = ntraj
# Holds seeds for rand num gen
self.seeds = seeds
# tidyup Hamiltonian before calculation (default = True)
self.tidy = tidy
# include the state in the function callback signature
self.rhs_with_state = rhs_with_state
# Use preexisting RHS function for time-dependent solvers
self.rhs_reuse = rhs_reuse
# Use filename for preexisting RHS function (will default to last
# compiled function if None & rhs_exists=True)
self.rhs_filename = rhs_filename
# small value in mc solver for computing correlations
self.mc_corr_eps = 1e-10
# Number of processors to use (mcsolve only)
if num_cpus:
self.num_cpus = num_cpus
else:
self.num_cpus = 0
# Tolerance for wavefunction norm (mcsolve only)
self.norm_tol = norm_tol
# Max. number of steps taken to find wavefunction norm to within
# norm_tol (mcsolve only)
self.norm_steps = norm_steps
# store final state?
self.store_final_state = store_final_state
# store states even if expectation operators are given?
self.store_states = store_states
# average mcsolver density matricies assuming steady state evolution
self.steady_state_average = steady_state_average
# Normalize output of solvers (turned off for batch unitary propagator mode)
self.normalize_output = normalize_output
def __str__(self):
if self.seeds is None:
seed_length = 0
else:
seed_length = len(self.seeds)
s = ""
s += "Options:\n"
s += "-----------\n"
s += "atol: " + str(self.atol) + "\n"
s += "rtol: " + str(self.rtol) + "\n"
s += "method: " + str(self.method) + "\n"
s += "order: " + str(self.order) + "\n"
s += "nsteps: " + str(self.nsteps) + "\n"
s += "first_step: " + str(self.first_step) + "\n"
s += "min_step: " + str(self.min_step) + "\n"
s += "max_step: " + str(self.max_step) + "\n"
s += "tidy: " + str(self.tidy) + "\n"
s += "num_cpus: " + str(self.num_cpus) + "\n"
s += "norm_tol: " + str(self.norm_tol) + "\n"
s += "norm_steps: " + str(self.norm_steps) + "\n"
s += "rhs_filename: " + str(self.rhs_filename) + "\n"
s += "rhs_reuse: " + str(self.rhs_reuse) + "\n"
s += "seeds: " + str(seed_length) + "\n"
s += "rhs_with_state: " + str(self.rhs_with_state) + "\n"
s += "average_expect: " + str(self.average_expect) + "\n"
s += "average_states: " + str(self.average_states) + "\n"
s += "ntraj: " + str(self.ntraj) + "\n"
s += "store_states: " + str(self.store_states) + "\n"
s += "store_final_state: " + str(self.store_final_state) + "\n"
return s
[docs]class Result():
"""Class for storing simulation results from any of the dynamics solvers.
Attributes
----------
solver : str
Which solver was used [e.g., 'mesolve', 'mcsolve', 'brmesolve', ...]
times : list/array
Times at which simulation data was collected.
expect : list/array
Expectation values (if requested) for simulation.
states : array
State of the simulation (density matrix or ket) evaluated at ``times``.
num_expect : int
Number of expectation value operators in simulation.
num_collapse : int
Number of collapse operators in simualation.
ntraj : int/list
Number of trajectories (for stochastic solvers). A list indicates
that averaging of expectation values was done over a subset of total
number of trajectories.
col_times : list
Times at which state collpase occurred. Only for Monte Carlo solver.
col_which : list
Which collapse operator was responsible for each collapse in
``col_times``. Only for Monte Carlo solver.
"""
def __init__(self):
self.solver = None
self.times = None
self.states = []
self.expect = []
self.num_expect = 0
self.num_collapse = 0
self.ntraj = None
self.seeds = None
self.col_times = None
self.col_which = None
def __str__(self):
s = "Result object "
if self.solver:
s += "with " + self.solver + " data.\n"
else:
s += "missing solver information.\n"
s += "-" * (len(s) - 1) + "\n"
if self.states is not None and len(self.states) > 0:
s += "states = True\n"
elif self.expect is not None and len(self.expect) > 0:
s += "expect = True\nnum_expect = " + str(self.num_expect) + ", "
else:
s += "states = True, expect = True\n" + \
"num_expect = " + str(self.num_expect) + ", "
s += "num_collapse = " + str(self.num_collapse)
if self.solver == 'mcsolve':
s += ", ntraj = " + str(self.ntraj)
return s
def __repr__(self):
return self.__str__()
def __getstate__(self):
# defines what happens when Qobj object gets pickled
self.__dict__.update({'qutip_version': __version__[:5]})
return self.__dict__
def __setstate__(self, state):
# defines what happens when loading a pickled Qobj
if 'qutip_version' in state.keys():
del state['qutip_version']
(self.__dict__).update(state)
class SolverConfiguration():
def __init__(self):
self.cgen_num = 0
self.reset()
def reset(self):
# General stuff
self.tlist = None # evaluations times
self.ntraj = None # number / list of trajectories
self.options = None # options for solvers
self.norm_tol = None # tolerance for wavefunction norm
self.norm_steps = None # max. number of steps to take in finding
# Initial state stuff
self.psi0 = None # initial state
self.psi0_dims = None # initial state dims
self.psi0_shape = None # initial state shape
# flags for setting time-dependence, collapse ops, and number of times
# codegen has been run
self.cflag = 0 # Flag signaling collapse operators
self.tflag = 0 # Flag signaling time-dependent problem
# time-dependent (TD) function stuff
self.tdfunc = None # Placeholder for TD RHS function.
self.tdname = None # Name of td .pyx file
self.colspmv = None # Placeholder for TD col-spmv function.
self.colexpect = None # Placeholder for TD col_expect function.
self.string = None # Holds string of variables passed to td solver
self.soft_reset()
def soft_reset(self):
# Hamiltonian stuff
self.h_td_inds = [] # indicies of time-dependent Hamiltonian operators
self.h_tdterms = [] # List of td strs and funcs
self.h_data = None # List of sparse matrix data
self.h_ind = None # List of sparse matrix indices
self.h_ptr = None # List of sparse matrix ptrs
# Expectation operator stuff
self.e_num = 0 # number of expect ops
self.e_ops_data = [] # expect op data
self.e_ops_ind = [] # expect op indices
self.e_ops_ptr = [] # expect op indptrs
self.e_ops_isherm = [] # expect op isherm
# Collapse operator stuff
self.c_num = 0 # number of collapse ops
self.c_const_inds = [] # indicies of constant collapse operators
self.c_td_inds = [] # indicies of time-dependent collapse operators
self.c_ops_data = [] # collapse op data
self.c_ops_ind = [] # collapse op indices
self.c_ops_ptr = [] # collapse op indptrs
self.c_args = [] # store args for time-dependent collapse func.
# Norm collapse operator stuff
self.n_ops_data = [] # norm collapse op data
self.n_ops_ind = [] # norm collapse op indices
self.n_ops_ptr = [] # norm collapse op indptrs
# holds executable strings for time-dependent collapse evaluation
self.col_expect_code = None
self.col_spmv_code = None
# hold stuff for function list based time dependence
self.h_td_inds = []
self.h_td_data = []
self.h_td_ind = []
self.h_td_ptr = []
self.h_funcs = None
self.h_func_args = None
self.c_funcs = None
self.c_func_args = None
def _format_time(t, tt=None, ttt=None):
time_str = str(datetime.timedelta(seconds=t))
if tt is not None and ttt is not None:
sect_percent = 100*t/tt
solve_percent = 100*t/ttt
time_str += " ({:03.2f}% section, {:03.2f}% total)".format(
sect_percent, solve_percent)
elif tt is not None:
sect_percent = 100*t/tt
time_str += " ({:03.2f}% section)".format(sect_percent)
elif ttt is not None:
solve_percent = 100*t/ttt
time_str += " ({:03.2f}% total)".format(solve_percent)
return time_str
[docs]class Stats(object):
"""
Statistical information on the solver performance
Statistics can be grouped into sections.
If no section names are given in the the contructor, then all statistics
will be added to one section 'main'
Parameters
----------
section_names : list
list of keys that will be used as keys for the sections
These keys will also be used as names for the sections
The text in the output can be overidden by setting the header property
of the section
If no names are given then one section called 'main' is created
Attributes
----------
sections : OrderedDict of _StatsSection
These are the sections that are created automatically on instantiation
or added using add_section
header : string
Some text that will be used as the heading in the report
By default there is None
total_time : float
Time in seconds for the solver to complete processing
Can be None, meaning that total timing percentages will be reported
Methods
-------
add_section
Add another section
add_count
Add some stat that is an integer count
add_timing
Add some timing statistics
add_message
Add some text type for output in the report
report:
Output the statistics report to console or file.
"""
def __init__(self, section_names=None):
self._def_section_name = 'main'
self.sections = OrderedDict()
self.total_time = None
self.header = None
if isinstance(section_names, list):
c = 0
for name in section_names:
self.sections[name] = _StatsSection(name, self)
if c == 0:
self._def_section_name = name
c += 1
else:
self.sections[self._def_section_name] = \
_StatsSection(self._def_section_name)
def _get_section(self, section):
if section is None:
return self.sections[self._def_section_name]
elif isinstance(section, _StatsSection):
return section
else:
sect = self.sections.get(section, None)
if sect is None:
raise ValueError("Unknown section {}".format(section))
else:
return sect
[docs] def add_section(self, name):
"""
Add another section with the given name
Parameters
----------
name : string
will be used as key for sections dict
will also be the header for the section
Returns
-------
section : `class` : _StatsSection
The new section
"""
sect = _StatsSection(name, self)
self.sections[name] = sect
return sect
[docs] def add_count(self, key, value, section=None):
"""
Add value to count. If key does not already exist in section then
it is created with this value.
If key already exists it is increased by the give value
value is expected to be an integer
Parameters
----------
key : string
key for the section.counts dictionary
reusing a key will result in numerical addition of value
value : int
Initial value of the count, or added to an existing count
section: string or `class` : _StatsSection
Section which to add the count to.
If None given, the default (first) section will be used
"""
self._get_section(section).add_count(key, value)
[docs] def add_timing(self, key, value, section=None):
"""
Add value to timing. If key does not already exist in section then
it is created with this value.
If key already exists it is increased by the give value
value is expected to be a float, and given in seconds.
Parameters
----------
key : string
key for the section.timings dictionary
reusing a key will result in numerical addition of value
value : int
Initial value of the timing, or added to an existing timing
section: string or `class` : _StatsSection
Section which to add the timing to.
If None given, the default (first) section will be used
"""
self._get_section(section).add_timing(key, value)
[docs] def add_message(self, key, value, section=None, sep=";"):
"""
Add value to message. If key does not already exist in section then
it is created with this value.
If key already exists the value is added to the message
The value will be converted to a string
Parameters
----------
key : string
key for the section.messages dictionary
reusing a key will result in concatenation of value
value : int
Initial value of the message, or added to an existing message
sep : string
Message will be prefixed with this string when concatenating
section: string or `class` : _StatsSection
Section which to add the message to.
If None given, the default (first) section will be used
"""
self._get_section(section).add_message(key, value, sep=sep)
[docs] def set_total_time(self, value, section=None):
"""
Sets the total time for the complete solve or for a specific section
value is expected to be a float, and given in seconds
Parameters
----------
value : float
Time in seconds to complete the solver section
section : string or `class` : _StatsSection
Section which to set the total_time for
If None given, the total_time for complete solve is set
"""
if not isinstance(value, float):
try:
value = float(value)
except:
raise TypeError("value is expected to be a float")
if section is None:
self.total_time = value
else:
sect = self._get_section(section)
sect.total_time = value
[docs] def report(self, output=sys.stdout):
"""
Report the counts, timings and messages from the sections.
Sections are reported in the order that the names were supplied
in the constructor.
The counts, timings and messages are reported in the order that they
are added to the sections
The output can be written to anything that supports a write method,
e.g. a file or the console (default)
The output is intended to in markdown format
Parameters
----------
output : stream
file or console stream - anything that support write - where
the output will be written
"""
if not hasattr(output, 'write'):
raise TypeError("output must have a write method")
if self.header:
output.write("{}\n{}\n".format(self.header,
("="*len(self.header))))
for name, sect in self.sections.items():
sect.report(output)
if self.total_time is not None:
output.write("\nSummary\n-------\n")
output.write("{}\t solver total time\n".format(
_format_time(self.total_time)))
[docs] def clear(self):
"""
Clear counts, timings and messages from all sections
"""
for sect in self.sections.values():
sect.clear()
self.total_time = None
class _StatsSection(object):
"""
Not intended to be directly instantiated
This is the type for the SolverStats.sections values
The method parameter descriptions are the same as for those the parent
with the same method name
Parameters
----------
name : string
key for the parent sections dictionary
will also be used as the header
parent : `class` : SolverStats
The container for all the sections
Attributes
----------
name : string
key for the parent sections dictionary
will also be used as the header
parent : `class` : SolverStats
The container for all the sections
header : string
Used as heading for section in report
counts : OrderedDict
The integer type statistics for the stats section
timings : OrderedDict
The timing type statistics for the stats section
Expected to contain float values representing values in seconds
messages : OrderedDict
Text type output to be reported
total_time : float
Total time for processing in the section
Can be None, meaning that section timing percentages will be reported
"""
def __init__(self, name, parent):
self.parent = parent
self.header = str(name)
self.name = name
self.counts = OrderedDict()
self.timings = OrderedDict()
self.messages = OrderedDict()
self.total_time = None
def add_count(self, key, value):
"""
Add value to count. If key does not already exist in section then
it is created with this value.
If key already exists it is increased by the given value
value is expected to be an integer
"""
if not isinstance(value, int):
try:
value = int(value)
except:
raise TypeError("value is expected to be an integer")
if key in self.counts:
self.counts[key] += value
else:
self.counts[key] = value
def add_timing(self, key, value):
"""
Add value to timing. If key does not already exist in section then
it is created with this value.
If key already exists it is increased by the give value
value is expected to be a float, and given in seconds.
"""
if not isinstance(value, float):
try:
value = float(value)
except:
raise TypeError("value is expected to be a float")
if key in self.timings:
self.timings[key] += value
else:
self.timings[key] = value
def add_message(self, key, value, sep=";"):
"""
Add value to message. If key does not already exist in section then
it is created with this value.
If key already exists the value is added to the message
The value will be converted to a string
"""
value = str(value)
if key in self.messages:
if sep is not None:
try:
value = sep + value
except:
TypeError("It is not possible to concatenate the value with "
"the given seperator")
self.messages[key] += value
else:
self.messages[key] = value
def report(self, output=sys.stdout):
"""
Report the counts, timings and messages for this section.
Note the percentage of the section and solver total times will be
given if the parent and or section total_time is set
"""
if self.header:
output.write("\n{}\n{}\n".format(self.header,
("-"*len(self.header))))
# TODO: Make the timings and counts ouput in a table format
# Generally make more pretty
# Report timings
try:
ttt = self.parent.total_time
except:
ttt = None
tt = self.total_time
output.write("### Timings:\n")
for key, value in self.timings.items():
l = " - {}\t{}\n".format(_format_time(value, tt, ttt), key)
output.write(l)
if tt is not None:
output.write(" - {}\t{} total time\n".format(_format_time(tt),
self.name))
# Report counts
output.write("### Counts:\n")
for key, value in self.counts.items():
l = " - {}\t{}\n".format(value, key)
output.write(l)
# Report messages
output.write("### Messages:\n")
for key, value in self.messages.items():
l = " - {}:\t{}\n".format(key, value)
output.write(l)
def clear(self):
"""
Clear counts, timings and messages from this section
"""
self.counts.clear()
self.timings.clear()
self.messages.clear()
self.total_time = None
def _solver_safety_check(H, state, c_ops=[], e_ops=[], args={}):
# Input is std Qobj (Hamiltonian or Liouvillian)
if isinstance(H, Qobj):
Hdims = H.dims
Htype = H.type
_structure_check(Hdims, Htype, state)
# Input H is function
elif isinstance(H, (FunctionType, BuiltinFunctionType)):
Hdims = H(0,args).dims
Htype = H(0,args).type
_structure_check(Hdims, Htype, state)
# Input is td-list
elif isinstance(H, list):
if isinstance(H[0], Qobj):
Hdims = H[0].dims
Htype = H[0].type
elif isinstance(H[0], list):
Hdims = H[0][0].dims
Htype = H[0][0].type
else:
raise Exception('Invalid td-list element.')
# Check all operators in list
for ii in range(len(H)):
if isinstance(H[ii], Qobj):
_temp_dims = H[ii].dims
_temp_type = H[ii].type
elif isinstance(H[ii], list):
_temp_dims = H[ii][0].dims
_temp_type = H[ii][0].type
else:
raise Exception('Invalid td-list element.')
_structure_check(_temp_dims,_temp_type,state)
else:
raise Exception('Invalid time-dependent format.')
for ii in range(len(c_ops)):
if isinstance(c_ops[ii], Qobj):
_temp_state = c_ops[ii]
elif isinstance(c_ops[ii], list):
_temp_state = c_ops[ii][0]
else:
raise Exception('Invalid td-list element.')
_structure_check(Hdims, Htype, _temp_state)
for ii in range(len(e_ops)):
if isinstance(e_ops[ii], Qobj):
_temp_state = e_ops[ii]
elif isinstance(e_ops[ii], list):
_temp_state = e_ops[ii][0]
else:
raise Exception('Invalid td-list element.')
_structure_check(Hdims,Htype,_temp_state)
def _structure_check(Hdims, Htype, state):
# Input state is a ket vector
if state.type == 'ket':
# Input is Hamiltonian
if Htype == 'oper':
if Hdims[1] != state.dims[0]:
raise Exception('Input operator and ket do not share same structure.')
# Input is super and state is ket
elif Htype == 'super':
if Hdims[1][1] != state.dims[0]:
raise Exception('Input operator and ket do not share same structure.')
else:
raise Exception('Invalid input operator.')
# Input state is a density matrix
elif state.type == 'oper':
# Input is Hamiltonian and state is density matrix
if Htype == 'oper':
if Hdims[1] != state.dims[0]:
raise Exception('Input operators do not share same structure.')
# Input is super op. and state is density matrix
elif Htype == 'super':
if Hdims[1] != state.dims:
raise Exception('Input operators do not share same structure.')
#
# create a global instance of the SolverConfiguration class
#
config = SolverConfiguration()
# for backwards compatibility
Odeoptions = Options
Odedata = Result