# Source code for qutip.nonmarkov.transfertensor

# -*- coding: utf-8 -*-
# This file is part of QuTiP: Quantum Toolbox in Python.
#
#    Copyright (c) 2015 and later, Arne L. Grimsmo
#
#    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.
###############################################################################

# @author: Arne L. Grimsmo
# @email1: arne.grimsmo@gmail.com
# @organization: University of Sherbrooke

"""
This module contains an implementation of the non-Markovian transfer tensor
method (TTM), introduced in [1].

[1] Javier Cerrillo and Jianshu Cao, Phys. Rev. Lett 112, 110401 (2014)
"""

import numpy as np

from qutip import (Options, spre, vector_to_operator, operator_to_vector,
ket2dm, isket)
from qutip.solver import Result
from qutip.expect import expect_rho_vec

[docs]class TTMSolverOptions:
"""Class of options for the Transfer Tensor Method solver.

Attributes
----------
dynmaps : list of :class:qutip.Qobj
List of precomputed dynamical maps (superoperators),
or a callback function that returns the
superoperator at a given time.

times : array_like
List of times :math:t_n at which to calculate :math:\\rho(t_n)

learningtimes : array_like
List of times :math:t_k to use as learning times if argument
dynmaps is a callback function.

thres : float
Threshold for halting. Halts if  :math:||T_{n}-T_{n-1}|| is below
treshold.

options : :class:qutip.solver.Options
Generic solver options.
"""

def __init__(self, dynmaps=None, times=[], learningtimes=[],
thres=0.0, options=None):

if options is None:
options = Options()

self.dynmaps = dynmaps
self.times = times
self.learningtimes = learningtimes
self.thres = thres
self.store_states = options.store_states

[docs]def ttmsolve(dynmaps, rho0, times, e_ops=[], learningtimes=None, tensors=None,
**kwargs):
"""
Solve time-evolution using the Transfer Tensor Method, based on a set of
precomputed dynamical maps.

Parameters
----------
dynmaps : list of :class:qutip.Qobj
List of precomputed dynamical maps (superoperators),
or a callback function that returns the
superoperator at a given time.

rho0 : :class:qutip.Qobj
Initial density matrix or state vector (ket).

times : array_like
list of times :math:t_n at which to compute :math:\\rho(t_n).
Must be uniformily spaced.

e_ops : list of :class:qutip.Qobj / callback function
single operator or list of operators for which to evaluate
expectation values.

learningtimes : array_like
list of times :math:t_k for which we have knowledge of the dynamical
maps :math:E(t_k).

tensors : array_like
optional list of precomputed tensors :math:T_k

kwargs : dictionary
Optional keyword arguments. See
:class:qutip.nonmarkov.ttm.TTMSolverOptions.

Returns
-------
output: :class:qutip.solver.Result
An instance of the class :class:qutip.solver.Result.
"""

opt = TTMSolverOptions(dynmaps=dynmaps, times=times,
learningtimes=learningtimes, **kwargs)

diff = None

if isket(rho0):
rho0 = ket2dm(rho0)

output = Result()
e_sops_data = []

if callable(e_ops):
n_expt_op = 0
expt_callback = True

else:
try:
tmp = e_ops[:]
del tmp

n_expt_op = len(e_ops)
expt_callback = False

if n_expt_op == 0:
# fall back on storing states
opt.store_states = True

for op in e_ops:
e_sops_data.append(spre(op).data)
if op.isherm and rho0.isherm:
output.expect.append(np.zeros(len(times)))
else:
output.expect.append(np.zeros(len(times), dtype=complex))
except TypeError:
raise TypeError("Argument 'e_ops' should be a callable or" +
"list-like.")

if tensors is None:
tensors, diff = _generatetensors(dynmaps, learningtimes, opt=opt)

if rho0.isoper:
# vectorize density matrix
rho0vec = operator_to_vector(rho0)
else:
# rho0 might be a super in which case we should not vectorize
rho0vec = rho0

K = len(tensors)
states = [rho0vec]
for n in range(1, len(times)):
states.append(None)
for k in range(n):
if n-k < K:
states[-1] += tensors[n-k]*states[k]
for i, r in enumerate(states):
if opt.store_states or expt_callback:
if r.type == 'operator-ket':
states[i] = vector_to_operator(r)
else:
states[i] = r
if expt_callback:
# use callback method
e_ops(times[i], states[i])
for m in range(n_expt_op):
if output.expect[m].dtype == complex:
output.expect[m][i] = expect_rho_vec(e_sops_data[m], r, 0)
else:
output.expect[m][i] = expect_rho_vec(e_sops_data[m], r, 1)

output.solver = "ttmsolve"
output.times = times

output.ttmconvergence = diff

if opt.store_states:
output.states = states

return output

def _generatetensors(dynmaps, learningtimes=None, **kwargs):
"""
Generate the tensors :math:T_1,\dots,T_K from the dynamical maps
:math:E(t_k).

A stationary process is assumed, i.e., :math:T_{n,k} = T_{n-k}.

Parameters
----------
dynmaps : list of :class:qutip.Qobj
List of precomputed dynamical maps (superoperators) at the times
specified in learningtimes, or a callback function that returns the
superoperator at a given time.

learningtimes : array_like
list of times :math:t_k to use if argument dynmaps is a callback
function.

kwargs : dictionary
Optional keyword arguments. See
:class:qutip.nonmarkov.ttm.TTMSolverOptions.

Returns
-------
Tlist: list of :class:qutip.Qobj.
A list of transfer tensors :math:T_1,\dots,T_K
"""

# Determine if dynmaps is callable or list-like
if callable(dynmaps):
if learningtimes is None:
raise TypeError("Argument 'learnintimes' required when 'dynmaps'" +
"is a callback function.")

def dynmapfunc(n): return dynmaps(learningtimes[n])
Kmax = len(learningtimes)
else:
try:
tmp = dynmaps[:]
del tmp

def dynmapfunc(n): return dynmaps[n]
Kmax = len(dynmaps)
except TypeError:
raise TypeError("Argument 'dynmaps' should be a callable or" +
"list-like.")

if "opt" not in kwargs:
opt = TTMSolverOptions(dynmaps=dynmaps, learningtimes=learningtimes,
**kwargs)
else:
opt = kwargs['opt']

Tlist = []
diff = [0.0]
for n in range(Kmax):
T = dynmapfunc(n)
for m in range(1, n):
T -= Tlist[n-m]*dynmapfunc(m)
Tlist.append(T)
if n > 1:
diff.append((Tlist[-1]-Tlist[-2]).norm())
if diff[-1] < opt.thres:
# Below threshold for truncation
print('breaking', (Tlist[-1]-Tlist[-2]).norm(), n)
break
return Tlist, diff