# Source code for qutip.control.propcomp

# -*- coding: utf-8 -*-
# This file is part of QuTiP: Quantum Toolbox in Python.
#
#    Copyright (c) 2014 and later, Alexander J G Pitchford
#
#    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: Alexander Pitchford
# @email1: agp1@aber.ac.uk
# @email2: alex.pitchford@gmail.com
# @organization: Aberystwyth University
# @supervisor: Daniel Burgarth

"""
Propagator Computer
Classes used to calculate the propagators,

Note the methods in the _Diag class was inspired by:
DYNAMO - Dynamic Framework for Quantum Optimal Control
See Machnes et.al., arXiv.1011.4874
"""

# import os
import warnings
import numpy as np
import scipy.linalg as la
import scipy.sparse as sp
# QuTiP
from qutip.qobj import Qobj
# QuTiP logging
import qutip.logging_utils as logging
logger = logging.get_logger()
# QuTiP control modules
from qutip.control import errors

def _func_deprecation(message, stacklevel=3):
"""
Issue deprecation warning
Using stacklevel=3 will ensure message refers the function
calling with the deprecated parameter,
"""
warnings.warn(message, DeprecationWarning, stacklevel=stacklevel)

[docs]class PropagatorComputer(object):
"""
Base for all  Propagator Computer classes
that are used to calculate the propagators,
Note: they must be instantiated with a Dynamics object, that is the
container for the data that the functions operate on
This base class cannot be used directly. See subclass descriptions
and choose the appropriate one for the application

Attributes
----------
log_level : integer
level of messaging output from the logger.
Options are attributes of qutip_utils.logging,
in decreasing levels of messaging, are:
DEBUG_INTENSE, DEBUG_VERBOSE, DEBUG, INFO, WARN, ERROR, CRITICAL
Anything WARN or above is effectively 'quiet' execution,
assuming everything runs as expected.
The default NOTSET implies that the level will be taken from
the QuTiP settings file, which by default is WARN

indicates whether the computer class instance is capable
of computing propagator gradients. It is used to determine
whether to create the Dynamics prop_grad array
"""
def __init__(self, dynamics, params=None):
self.parent = dynamics
self.params = params
self.reset()

[docs]    def reset(self):
"""
reset any configuration data
"""
self.id_text = 'PROP_COMP_BASE'
self.log_level = self.parent.log_level

[docs]    def apply_params(self, params=None):
"""
Set object attributes based on the dictionary (if any) passed in the
instantiation, or passed as a parameter
This is called during the instantiation automatically.
The key value pairs are the attribute name and value
Note: attributes are created if they do not exist already,
and are overwritten if they do.
"""
if not params:
params = self.params

if isinstance(params, dict):
self.params = params
for key in params:
setattr(self, key, params[key])

@property
def log_level(self):
return logger.level

@log_level.setter
def log_level(self, lvl):
"""
Set the log_level attribute and set the level of the logger
that is call logger.setLevel(lvl)
"""
logger.setLevel(lvl)

def compute_propagator(self, k):
_func_deprecation("'compute_propagator' has been replaced "
"by '_compute_propagator'")
return self._compute_propagator(k)

def _compute_propagator(self, k):
"""
calculate the progator between X(k) and X(k+1)
Uses matrix expm of the dyn_gen at that point (in time)
Assumes that the dyn_gen have been been calculated,
i.e. drift and ctrls combined
Return the propagator
"""
dyn = self.parent
dgt = dyn._get_phased_dyn_gen(k)*dyn.tau[k]
if dyn.oper_dtype == Qobj:
prop = dgt.expm()
else:
prop = la.expm(dgt)
return prop

def compute_diff_prop(self, k, j, epsilon):
_func_deprecation("'compute_diff_prop' has been replaced "
"by '_compute_diff_prop'")
return self._compute_diff_prop( k, j, epsilon)

def _compute_diff_prop(self, k, j, epsilon):
"""
Calculate the propagator from the current point to a trial point
a distance 'epsilon' (change in amplitude)
in the direction the given control j in timeslot k
Returns the propagator
"""
raise errors.UsageError("Not implemented in the baseclass."
" Choose a subclass")

"""
Calculate the gradient of propagator wrt the control amplitude
in the timeslot.
"""
raise errors.UsageError("Not implemented in the baseclass."
" Choose a subclass")

"""
This subclass can be used when the propagator is calculated simply
by expm of the dynamics generator, i.e. when gradients will be calculated
using approximate methods.
"""

[docs]    def reset(self):
"""
reset any configuration data
"""
PropagatorComputer.reset(self)
self.id_text = 'APPROX'
self.apply_params()

def _compute_diff_prop(self, k, j, epsilon):
"""
Calculate the propagator from the current point to a trial point
a distance 'epsilon' (change in amplitude)
in the direction the given control j in timeslot k
Returns the propagator
"""
dyn = self.parent
dgt_eps = (dyn._get_phased_dyn_gen(k) +
epsilon*dyn._get_phased_ctrl_dyn_gen(k, j))*dyn.tau[k]

if dyn.oper_dtype == Qobj:
prop_eps = dgt_eps.expm()
else:
prop_eps = la.expm(dgt_eps)

return prop_eps

[docs]class PropCompDiag(PropagatorComputer):
"""
Coumputes the propagator exponentiation using diagonalisation of
of the dynamics generator
"""
[docs]    def reset(self):
"""
reset any configuration data
"""
PropagatorComputer.reset(self)
self.id_text = 'DIAG'
self.apply_params()

def _compute_propagator(self, k):
"""
Calculates the exponentiation of the dynamics generator (H)
As part of the calc the the eigen decomposition is required, which
is reused in the propagator gradient calculation
"""
dyn = self.parent
dyn._ensure_decomp_curr(k)

if dyn.oper_dtype == Qobj:

prop = (dyn._dyn_gen_eigenvectors[k]*dyn._prop_eigen[k]*
else:
prop = dyn._dyn_gen_eigenvectors[k].dot(
dyn._prop_eigen[k]).dot(

return prop

"""
Calculate the gradient of propagator wrt the control amplitude
in the timeslot.

Returns:
"""
dyn = self.parent
dyn._ensure_decomp_curr(k)

if compute_prop:
prop = self._compute_propagator(k)

if dyn.oper_dtype == Qobj:
# put control dyn_gen in combined dg diagonal basis
dyn._get_phased_ctrl_dyn_gen(k, j)*
dyn._dyn_gen_eigenvectors[k])
# multiply (elementwise) by timeslice and factor matrix
cdg = Qobj(np.multiply(cdg.full()*dyn.tau[k],
dyn._dyn_gen_factormatrix[k]), dims=dyn.dyn_dims)
else:
# put control dyn_gen in combined dg diagonal basis
dyn._get_phased_ctrl_dyn_gen(k, j)).dot(
dyn._dyn_gen_eigenvectors[k])
# multiply (elementwise) by timeslice and factor matrix
cdg = np.multiply(cdg*dyn.tau[k], dyn._dyn_gen_factormatrix[k])

if compute_prop:
else:

class PropCompAugMat(PropagatorComputer):
"""
Augmented Matrix (deprecated - see _Frechet)

It should work for all systems, e.g. open, symplectic
There will be other PropagatorComputer subclasses that are more efficient
The _Frechet class should provide exactly the same functionality
more efficiently.

Note the propagator gradient calculation using the augmented matrix
is taken from:
'Robust quantum gates for open systems via optimal control:
Markovian versus non-Markovian dynamics'
Frederik F Floether, Pierre de Fouquieres, and Sophie G Schirmer
"""
def reset(self):
PropagatorComputer.reset(self)
self.id_text = 'AUG_MAT'
self.apply_params()

def _get_aug_mat(self, k, j):
"""
Generate the matrix [[A, E], [0, A]] where
A is the overall dynamics generator
E is the control dynamics generator
for a given timeslot and control
returns this augmented matrix
"""
dyn = self.parent
dg = dyn._get_phased_dyn_gen(k)

if dyn.oper_dtype == Qobj:
A = dg.data*dyn.tau[k]
E = dyn._get_phased_ctrl_dyn_gen(k, j).data*dyn.tau[k]
Z = sp.csr_matrix(dg.data.shape)
aug = Qobj(sp.vstack([sp.hstack([A, E]), sp.hstack([Z, A])]))
elif dyn.oper_dtype == np.ndarray:
A = dg*dyn.tau[k]
E = dyn._get_phased_ctrl_dyn_gen(k, j)*dyn.tau[k]
Z = np.zeros(dg.shape)
aug = np.vstack([np.hstack([A, E]), np.hstack([Z, A])])
else:
A = dg*dyn.tau[k]
E = dyn._get_phased_ctrl_dyn_gen(k, j)*dyn.tau[k]
Z = dg*0.0
aug = sp.vstack([sp.hstack([A, E]), sp.hstack([Z, A])])
return aug

"""
Calculate the gradient of propagator wrt the control amplitude
in the timeslot using the exponentiation of the the augmented
matrix.
The propagtor is calculated for 'free' in this method
and hence it is returned if compute_prop==True
Returns:
"""
dyn = self.parent
dg = dyn._get_phased_dyn_gen(k)
aug = self._get_aug_mat(k, j)

if dyn.oper_dtype == Qobj:
aug_exp = aug.expm()
dims=dyn.dyn_dims)
if compute_prop:
prop = Qobj(aug_exp.data[:dg.shape[0], :dg.shape[1]],
dims=dyn.dyn_dims)
else:
aug_exp = la.expm(aug)
if compute_prop:
prop = aug_exp[:dg.shape[0], :dg.shape[1]]

if compute_prop:
else:

[docs]class PropCompFrechet(PropagatorComputer):
"""
Frechet method for calculating the propagator:
exponentiating the combined dynamics generator
It should work for all systems, e.g. unitary, open, symplectic
There are other PropagatorComputer subclasses that may be more efficient
"""
[docs]    def reset(self):
PropagatorComputer.reset(self)
self.id_text = 'FRECHET'
self.apply_params()

"""
Calculate the gradient of propagator wrt the control amplitude
in the timeslot using the expm_frechet method
The propagtor is calculated (almost) for 'free' in this method
and hence it is returned if compute_prop==True
Returns:
"""
dyn = self.parent

if dyn.oper_dtype == Qobj:
A = dyn._get_phased_dyn_gen(k).full()*dyn.tau[k]
E = dyn._get_phased_ctrl_dyn_gen(k, j).full()*dyn.tau[k]
if compute_prop:
prop = Qobj(prop_dense, dims=dyn.dyn_dims)
dims=dyn.dyn_dims)
else:
dims=dyn.dyn_dims)
elif dyn.oper_dtype == np.ndarray:
A = dyn._get_phased_dyn_gen(k)*dyn.tau[k]
E = dyn._get_phased_ctrl_dyn_gen(k, j)*dyn.tau[k]
if compute_prop:
else:
compute_expm=False)
else:
# Assuming some sparse matrix
spcls = dyn._dyn_gen[k].__class__
A = (dyn._get_phased_dyn_gen(k)*dyn.tau[k]).toarray()
E = (dyn._get_phased_ctrl_dyn_gen(k, j)*dyn.tau[k]).toarray()
if compute_prop: