# 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__ = ['eseries', 'esval', 'esspec', 'estidy']
import numpy as np
import scipy.sparse as sp
from qutip.qobj import Qobj
[docs]class eseries():
"""
Class representation of an exponential-series expansion of
time-dependent quantum objects.
Attributes
----------
ampl : ndarray
Array of amplitudes for exponential series.
rates : ndarray
Array of rates for exponential series.
dims : list
Dimensions of exponential series components
shape : list
Shape corresponding to exponential series components
Methods
-------
value(tlist)
Evaluate an exponential series at the times listed in tlist
spec(wlist)
Evaluate the spectrum of an exponential series at frequencies in wlist.
tidyup()
Returns a tidier version of the exponential series
"""
__array_priority__ = 101
def __init__(self, q=None, s=np.array([])):
if isinstance(s, (int, float, complex)):
s = np.array([s])
if q is None:
self.ampl = np.array([])
self.rates = np.array([])
self.dims = [[1, 1]]
self.shape = [1, 1]
elif (len(s) == 0):
if isinstance(q, eseries):
self.ampl = q.ampl
self.rates = q.rates
self.dims = q.dims
self.shape = q.shape
elif isinstance(q, (np.ndarray, list)):
q = np.asarray(q, dtype=object)
ind = np.shape(q)
num = ind[0] # number of elements in q
if any([Qobj(x).shape != Qobj(q[0]).shape for x in q]):
raise TypeError('All amplitudes must have same dimension.')
self.ampl = np.array([x for x in q], dtype=object)
self.rates = np.zeros(ind)
self.dims = self.ampl[0].dims
self.shape = self.ampl[0].shape
elif isinstance(q, Qobj):
qo = Qobj(q)
self.ampl = np.array([qo], dtype=object)
self.rates = np.array([0])
self.dims = qo.dims
self.shape = qo.shape
else:
self.ampl = np.array([q])
self.rates = np.array([0])
self.dims = [[1, 1]]
self.shape = [1, 1]
elif len(s) != 0:
if isinstance(q, (np.ndarray, list)):
q = np.asarray(q, dtype=object)
ind = np.shape(q)
num = ind[0]
if any([Qobj(x).shape != Qobj(q[0]).shape for x in q]):
raise TypeError('All amplitudes must have same dimension.')
self.ampl = np.array([Qobj(q[x]) for x in range(0, num)],
dtype=object)
self.dims = self.ampl[0].dims
self.shape = self.ampl[0].shape
else:
num = 1
self.ampl = np.array([Qobj(q)], dtype=object)
self.dims = self.ampl[0].dims
self.shape = self.ampl[0].shape
if isinstance(s, (int, complex, float)):
if num != 1:
raise TypeError('Number of rates must match number ' +
'of members in object array.')
self.rates = np.array([s])
elif isinstance(s, (np.ndarray, list)):
if len(s) != num:
raise TypeError('Number of rates must match number ' +
' of members in object array.')
self.rates = np.array(s)
if len(self.ampl) != 0:
# combine arrays so that they can be sorted together
zipped = list(zip(self.rates, self.ampl))
zipped.sort() # sort rates from lowest to highest
rates, ampl = list(zip(*zipped)) # get back rates and ampl
self.ampl = np.array(ampl, dtype=object)
self.rates = np.array(rates)
def __str__(self): # string of ESERIES information
self.tidyup()
s = "ESERIES object: " + str(len(self.ampl)) + " terms\n"
s += "Hilbert space dimensions: " + str(self.dims) + "\n"
for k in range(0, len(self.ampl)):
s += "Exponent #" + str(k) + " = " + str(self.rates[k]) + "\n"
if isinstance(self.ampl[k], sp.spmatrix):
s += str(self.ampl[k]) + "\n"
else:
s += str(self.ampl[k]) + "\n"
return s
def __repr__(self):
return self.__str__()
# Addition with ESERIES on left (ex. ESERIES+5)
def __add__(self, other):
right = eseries(other)
if self.dims != right.dims:
raise TypeError("Incompatible operands for ESERIES addition")
out = eseries()
out.dims = self.dims
out.shape = self.shape
out.ampl = np.append(self.ampl, right.ampl)
out.rates = np.append(self.rates, right.rates)
return out
# Addition with ESERIES on right(ex. 5+ESERIES)
def __radd__(self, other):
return self + other
# define negation of ESERIES
def __neg__(self):
out = eseries()
out.dims = self.dims
out.shape = self.shape
out.ampl = -self.ampl
out.rates = self.rates
return out
# Subtraction with ESERIES on left (ex. ESERIES-5)
def __sub__(self, other):
return self + (-other)
# Subtraction with ESERIES on right (ex. 5-ESERIES)
def __rsub__(self, other):
return other + (-self)
# Multiplication with ESERIES on left (ex. ESERIES*other)
def __mul__(self, other):
if isinstance(other, eseries):
out = eseries()
out.dims = self.dims
out.shape = self.shape
for i in range(len(self.rates)):
for j in range(len(other.rates)):
out += eseries(self.ampl[i] * other.ampl[j],
self.rates[i] + other.rates[j])
return out
else:
out = eseries()
out.dims = self.dims
out.shape = self.shape
out.ampl = self.ampl * other
out.rates = self.rates
return out
# Multiplication with ESERIES on right (ex. other*ESERIES)
def __rmul__(self, other):
out = eseries()
out.dims = self.dims
out.shape = self.shape
out.ampl = other * self.ampl
out.rates = self.rates
return out
#
# todo:
# select_ampl, select_rate: functions to select some terms given the ampl
# or rate. This is done with {ampl} or (rate) in qotoolbox. we should use
# functions with descriptive names for this.
#
#
# evaluate the eseries for a list of times
#
[docs] def value(self, tlist):
"""
Evaluates an exponential series at the times listed in ``tlist``.
Parameters
----------
tlist : ndarray
Times at which to evaluate exponential series.
Returns
-------
val_list : ndarray
Values of exponential at times in ``tlist``.
"""
if self.ampl is None or len(self.ampl) == 0:
# no terms, evalue to zero
return np.zeros(np.shape(tlist))
if isinstance(tlist, float) or isinstance(tlist, int):
tlist = [tlist]
if isinstance(self.ampl[0], Qobj):
# amplitude vector contains quantum objects
val_list = []
for j in range(len(tlist)):
exp_factors = np.exp(np.array(self.rates) * tlist[j])
val = 0
for i in range(len(self.ampl)):
val += self.ampl[i] * exp_factors[i]
val_list.append(val)
val_list = np.array(val_list, dtype=object)
else:
# the amplitude vector contains c numbers
val_list = np.zeros(np.size(tlist), dtype=complex)
for j in range(len(tlist)):
exp_factors = np.exp(np.array(self.rates) * tlist[j])
val_list[j] = np.sum(np.dot(self.ampl, exp_factors))
if all(np.imag(val_list) == 0):
val_list = np.real(val_list)
if len(tlist) == 1:
return val_list[0]
else:
return val_list
[docs] def spec(self, wlist):
"""
Evaluate the spectrum of an exponential series at frequencies
in ``wlist``.
Parameters
----------
wlist : array_like
Array/list of frequenies.
Returns
-------
val_list : ndarray
Values of exponential series at frequencies in ``wlist``.
"""
val_list = np.zeros(np.size(wlist))
for i in range(len(wlist)):
val_list[i] = 2 * np.real(
np.dot(self.ampl, 1. / (1.0j * wlist[i] - self.rates)))
return val_list
[docs] def tidyup(self, *args):
""" Returns a tidier version of exponential series.
"""
#
# combine duplicate entries (same rate)
#
rate_tol = 1e-10
ampl_tol = 1e-10
ampl_dict = {}
unique_rates = {}
ur_len = 0
for r_idx in range(len(self.rates)):
# look for a matching rate in the list of unique rates
idx = -1
for ur_key in unique_rates.keys():
if abs(self.rates[r_idx] - unique_rates[ur_key]) < rate_tol:
idx = ur_key
break
if idx == -1:
# no matching rate, add it
unique_rates[ur_len] = self.rates[r_idx]
ampl_dict[ur_len] = [self.ampl[r_idx]]
ur_len = len(unique_rates)
else:
# found matching rate, append amplitude to its list
ampl_dict[idx].append(self.ampl[r_idx])
# create new amplitude and rate list with only unique rates, and
# nonzero amplitudes
self.rates = np.array([])
self.ampl = np.array([])
for ur_key in unique_rates.keys():
total_ampl = np.sum(np.asarray(ampl_dict[ur_key], dtype=object))
if (isinstance(total_ampl, float) or
isinstance(total_ampl, complex)):
if abs(total_ampl) > ampl_tol:
self.rates = np.append(self.rates, unique_rates[ur_key])
self.ampl = np.append(self.ampl, total_ampl)
else:
if abs(total_ampl.full()).max() > ampl_tol:
self.rates = np.append(self.rates, unique_rates[ur_key])
self.ampl = np.append(self.ampl,
np.asarray([total_ampl],
dtype=object))
return self
# -----------------------------------------------------------------------------
#
# wrapper functions for accessing the class methods (for compatibility with
# quantum optics toolbox)
#
def esval(es, tlist):
"""
Evaluates an exponential series at the times listed in ``tlist``.
Parameters
----------
tlist : ndarray
Times at which to evaluate exponential series.
Returns
-------
val_list : ndarray
Values of exponential at times in ``tlist``.
"""
return es.value(tlist)
def esspec(es, wlist):
"""Evaluate the spectrum of an exponential series at frequencies
in ``wlist``.
Parameters
----------
wlist : array_like
Array/list of frequenies.
Returns
-------
val_list : ndarray
Values of exponential series at frequencies in ``wlist``.
"""
return es.spec(wlist)
def estidy(es, *args):
"""
Returns a tidier version of exponential series.
"""
return es.tidyup()