5. Modules

5.1. hamitonians

This module holds the class that holds the Hamiltonian.

class hamiltonians.hamiltonian(H0, V, A=None, g=None, gradg=None)[source]

A class to hold and manipulate the Hamiltonian and dissipation terms.

The Hamiltonian should contain a static part, a list of possible perturbations, and also a list of dissipative terms for the Lindblad equation (which is not, of course, a proper term of the Hamiltonian)

There are two options for the user to pass the Hamiltonian to qocttools: as a list of Qobj objects, or as a Qobj-valued user-defined function.

In the first case, the Hamiltonian has the form:

\[H(t) = H_0 + \sum_{i=1}^N f_i(t) V_i\]

The pulses \(f_i(t)\) are not described in the hamiltonian object. The user must then pass H0 and a list of perturbations V (it may also be just one). In addition, for open systems, the hamiltonian object should also contain the Lindblad operators, which should be supplied by the user as a list of Qobj objects.
In the second case, the user should supply Qobj-valued functions. In this case, the H0 function contains the full Hamiltonian, and depends on time and on the set of “pulses” \(f_i(t)\):

\[H_0 = H_0(t, f_1(t), \dots, f_N(t))\]

The user must then also supply, in the argument V, a list with the derivatives of \(H0\) with respect to each \(f_i\) argument:

\[V_i(t, f_1(t), \dots, f_N(t)) = \frac{\partial H_o}{\partial f_i}(t, f_1(t), \dots, f_N(t))\]

The function H0 (and the functions in the V list) should have the following inteface:
```
def H0(t, args):
    f = args["f"]
    # f[0], f[1], f[2], ..., f[N-1] are the time-dependent functions
```
The arguments are a float t that means time, and a dictionary args with just one element with name “f”, and whose value is a list of functions of time.

Parameters:

H0 (Qobj or Qobj-valued function) – The qutip Qobj instance holding the static part of the Hamiltonian.
V (Qobj or list of Qobj or Qobj-valued function or list of Qobj-valued functions.) – One or various Qobj instances holding external perturbations.
A (Qobj or list of Qobj, default = None) – One or various Qobj instances holding the jump operators defining the Lindbladian.
g (still undocumnted, default = None) –
gradg (still undocumented, default = None) –

H0

The qutip Qobj instance holding the static part of the Hamiltonian.

Type:: Qobj

V

One or various Qobj instances holding external perturbations.

Type:: list of Qobj

A

One or various Qobj instances holding the jump operators defining the Lindblad operator.

Type:: list of Qobj

has_dissipative_terms()[source]

Returns True if the dissipative terms are not None

Returns:: True if the attribute holding the dissipative terms is not None, False otherwise.
Return type:: bool

hamiltonians.toliouville(h)[source]

This function transforms a Hamiltonian, defined in Hilbert space, together with the set of Lindblad operators, to the corresponding Liouvillian in Liouville space. It is also returned as a hamiltonian object

Parameters:: h (hamiltonian) – The input Hamiltonian, to be transformed.
Returns:: A hamiltonian object with the transformed Hamiltonian.
Return type:: hamiltonian

5.2. qoct

class qoct.Qoct(H, T, ntsteps, tg, f, y0, interaction_picture=False, solve_method='cfmagnus4', equality_constraints=None, sparse=False, output_file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>, new_parametrization=False, mpi_comm=None, floquet_mode='qoct')[source]

The class that is used to setup the optimization.

It stores the Hamiltonian, the definition of the target, the pulses, the initial state, and some parameters that determine how the optimization is done.

Parameters:

H (hamiltonian or list of hamiltonian) – A hamiltonian class instance, or list of instances for multitarget optimizations.
T (float) – total propagation time
ntsteps (int) – number of time steps
tg (Target) – The object that defines the target functional.
f (pulse or list of pulse) – The pulse object, or list of pulse objects (in the case the Hamiltonian(s) have more than one perturbations)
y0 (Qobj or list of Qobj) – The initial state(s) for the optimization. It should be a list for multitarget optimizations. It can be a Hilbert space state, a density matrix, or a propagator.
interaction_picture (bool, default = False) – The calculations can be done in the interaction picture, or in Schrödinger’s picture.
solve_method (str, default = 'cfmagnus4') – The method to use for the propagations.
equality_constraints (list, default = None) – The optimization can be done with constraints, that should be specified as a list of [func, tol] objects: func is a function that defines the constraint, and tol is a float that sets the tolerance accepted for the fulfillment of that constraint.
sparse (bool, default = False) – It forces the use of sparse algebra only (otherwise, some parts of the calculation are done with dense matrices).
of (file, default = sys.stdout) – One may ask the internal procedures of the qoct object to output messages to a file different from standard output

y0

Type:: list of Qobj

H

Type:: list of hamiltonian

almost_converged(fraction=0.99)[source]

Returns the iteration number at which the convergence was almost achieved

It analyzes the convergence history of an optimization, and looks for the iteration at which the value of the functional was already a fraction (0.99 by default) of the maximum.

Parameters:

fraction (float, default = 0.99) – The fraction of the maximum value at which the process is considered to be almost converged. If, at iteration 76, the value of the functional is already fraction * maximum, then this function will return 76 (along with the number of propagations necessary to reach iteration 76.

Returns:

float – The iteration at which the process was almost converged
float – The number of propagations done until that iteration.

check_grad(u, m=None, delta=0.01)[source]

A check on the accuracy of the gradient calculation

It computes one component of the gradient of the target functional using both the QOCT expression, and a finite-difference algorithm (the Ridders algorithm).

Parameters:

u (ndarray) – A numpy array holding all the control parameters
m (int, default = None) – Which gradient component to compute. If None, it will use the one with the largest absolute value
delta (float, default = 0.01) – A parameter determining the starting finite difference in the Ridders algorithm.

Returns:

float – The gradient component as computed with the QOCT formula.
float – The gradient component as computed with the finite-difference formula.
float – An estimation of the gradient accuracy computed with the finite difference formula.

dgfunc(u)[source]

Returns the value and gradient of the target functional for a set of parameters u.

It just calls the dgfunction function, using all the information contained in the class object.

Parameters:

u (ndarray) – A numpy array holding all the control parameters

Returns:

float – The value of the target parameter.
ndarray – The gradient of the target functional.

gfunc(u)[source]

Returns the value of the target functional for a set of parameters u.

It just calls the gfunction function, using all the information contained in the class object.

Parameters:: u (ndarray) – A numpy array holding all the control parameters
Returns:: The value of the target parameter.
Return type:: float

maximize(maxeval=100, stopval=None, verbose=False, algorithm=11, local_algorithm=None, upper_bounds=None, lower_bounds=None, tolerance=1e-06)[source]

Performs the QOCT maximization

This is the main procedure of the class, as it is the one that launches the optimization process. There are a number of parameters controlling the way in which this optimization is done, that can be controlled here.

The starting guess parameters for the optimizatio are read from the values of the pulses, that were associated to the class object when it was created.

5.3. pulses

This module includes the ‘pulse’ class, that is used to contain the time-dependent real function that is used to define the perturbation part of the Hamiltonian.

Also, it contains definitions of ‘typical’ pulse forms, such as the pi-pulse.

pulses.pi_pulse(t, u)[source]

A pi-pulse value at time t, parameters specifiey by the ndarray u

Returns the pi-pulse value at a specific time

Input: A: pi-pulse amplitude w: pi-pulse frecuency t0: time in which start the propagation of the pulse t: time in which the pulse is calculated t_duration: pi-pulse duration
Output:: if the t value in which the pulse calculation is requested is between t0 and t0 + t_duration, return the value of a pi-pulse with frecuency w and amplitude A, otherwise return 0.

pulses.pi_pulse_chain(t, u)[source]: Returns the value of a chain of pi pulse at a specific time

class pulses.pulse(type, T, u, bound=None)[source]

Definition of a class to hold control functions, a.k.a. pulses.

The “pulses” are real time-dependent functions defined in the interval \([0, T]\), defined as functions \(f(u_1, u_2, ..., u_P, T)\), where \(u = u_1, ..., u_P\) are parameters (the control parameters). These pulses are the control functions used by qocttools.

Depending on the type of the pulse, the parametrized form of these pulses is:

fourier

In this case, the pulse is a simple Fourier expansion:

\[f(u, t) = \frac{u_0}{\sqrt{T}} + \frac{2}{\sqrt{T}} \sum_i^M \left[ u_{2i}\cos(\omega_i t) + u_{2i-1}\sin(\omega_i t)\right]\]

The pulse therfore requires an odd number of parameters, \(P = 2M + 1\). The frequencies are :math`iomega_0 = ifrac{2pi}{T}` for \(i=1, \dots, M\). The cutoff M will be decided by the number of parameters that are pussed when initializing the object.
bound_fourier

This is a modification of the normal Fourier expansion, that ensures that the absolute value of the pulse is never larger than a certain bound. The definition is:

\[f(u, t) = \Phi(g(u, t))\]

where \(g(u, t)\) is a normal Fourier expansion as the one given above, and:

\[\Phi(x) = \kappa\frac{x}{1+\vert x\vert}\]

The value of \(\kappa\) must be supplied by the bound argument.
realtime

The (…)
enveloped_fourier

This is a normal Fourier expansion, but multiplied by a function \(S(t)\).
user_defined

This is the most general possibility: the user defines the pulse as a a python function. The function must be communicated to the pulse object using the assign_user_defined_function().

Parameters:

type (string) – The type of function: ‘fourier’, ‘bound_fourier’, ‘realtime’, ‘user_defined’, or ‘enveloped’ fourier
T (float) – The duration of the pulse
u (ndarray) – The parameters
bound (float, default = None) – For the ‘bound_fourier’ type, the bound

type

Type:: string

T

Type:: float

u

Type:: ndarray

nu

Type:: int

bound

Type:: float

assign_envelope_function(efunc)[source]: Sets the envelope function that may be multiplied by the pulse itself.

assign_user_defined_function(fu, dfu)[source]: For user-defined functions, set the function definition, and the definition of the gradient of the function.

dfu(t, m, u=None)[source]: Derivative of the pulse with respect to the m-th parameter, at time t

envelope(t)[source]: Returns the value of the envelope function at time t

f(t, args)[source]: The value of the pulse at time t

fitparams(fref, times, u0)[source]: Fit the parameters of a function to the best mpossible much wrt a reference function

fu(t, u=None)[source]

Returns the value of pulse at time t

Optionally, one can pass the control parameters u and in that way those are updated before the computation.

fw(t, u, w=None, explicit_formula=False)[source]: Fourier transform of the pulse

gradf(t)[source]

Grdient of the function at time t

“Gradient of the function” means the gradient with respect to all the control parameters.

print(filename)[source]

Prints a pulse to a file called ‘filename’.

It does not print the possible constraint functions.

Parameters:: filename (str) – Name of the file where the info about the pulse is written.

set_constraint(which, val=0.0)[source]: Sets a constraint on the shape of the pulse

set_parameters(u)[source]: Sets the parameters \(u = u_1, \dots, u_P\) on the pulse.

pulses.pulse_collection_get_parameters(f)[source]

Sets the parameters of a collection of pulses

It receives a numpy array of parameters u, whose dimension should be equal to the summ of parameters of the pulses in the collection f.

pulses.pulse_collection_parameter_range(f, j)[source]

Returns the ‘parameter range’ of a given pulse in a collection

Given a list of pulses ‘f’, the full control parameter list will be an array joining all the control parameters of each of of them. This function returns the indexes that would correspond, in that array, to one of the pulses.

Thus, for example, if we have two pulses with nu1 and nu2 parameters each, the range of the first one would be (0, nu1), whereas the range of the second would be (nu1, nu1+nu2).

Parameters:

f (list of pulse) – A list with pulse objects
j (int) – The pulse for which we want to get the parameter range.

Returns:

a list of two integer numbers, with the starting index of the parameters corresponding to the first pulse, and the final index.

Return type:

list of int

pulses.pulse_collection_set_parameters(f, u)[source]

Sets the parameters of a collection of pulses

It receives a numpy array of parameters u, whose dimension should be equal to the summ of parameters of the pulses in the collection f.

pulses.pulse_constraint_functions(f, param_range)[source]

Given a pulse, returns a list with all its constraint functions.

The constraint functions contained in the pulse object are functions of the nu control parameters on which the pulse depends. In contrast, the functions returned by this function are functions of all the control parameters in a list of pulses. Therefore, we need as an input a ‘parameter range’, i.e. the starting and finishin indexes of the parameters of this particular pulse in the list.

pulses.read_pulse(filename)[source]

Reads a pulse from the info contained in file.

Parameters:: filename (str) – Name of the file where the info about the pulse is written.
Returns:: The pulse whose information was previously written to filename.
Return type:: pulse

pulses.rotation(theta, n, dim=2, i=None, j=None)[source]

Returns a rotation operator between two levels, given an angle an axis.

The rotation affects level i and j (or the only two levels, if dim = 2, which is the default). The angle of rotation is given by theta, and the axis is given by the unitary vector n.

Parameters:

theta (float) – The rotation angle.
n (ndarray) – A three-dimensional float array containing a unit vector.
dim (int, default = 2) – The dimension of the rotation operator that will be created.
i (int, default = None) – One of the states affected by the rotation.
j (int, default = None) – One of the states affected by the rotation.

Returns:

A Qobj operator with the rotation operator.

Return type:

Qobj

pulses.rotationpulse(axis, A, omega, mu0, theta)[source]

Returns a pulse that implements a rotation betwen two states.

Given an amplitude for a pulse A, and two levels characterized by an energy difference omega and a coupling mu0, creates a pulse object that implements a rotation between those two levels in direction given by axis (0 => x, 1 => y, 2 => z), and angle theta.

Parameters:

axis (int) – The rotation index (0 => x, 1 => y, 2 => z).
A (float) – The amplitude of the pulse.
omega (float) – The frequency difference between the levels.
mu0 (float) – The coupling matrix element
theta (float) – The rotation angle

Returns:

A pulse object with the pulse that produces the rotation.

Return type:

pulse

5.4. solvers

Solvers for the Schrödinger/Lindblad equation

The code may use the qutip solvers, or use some internal solvers that may be faster in some circumstances.

solvers.cfmagnus2solver(H, f, y0, time, returnQoutput=True, interaction_picture=False, cops=None)[source]

Implementation of the exponential midpoint rule.

It is also the second-order commutator-free Magnus method.

solvers.cfmagnus4solver(H, f, psi0, time, returnQoutput=True, interaction_picture=False, cops=None)[source]: Implementation of the fourth order commutator-free Magnus method.

solvers.solve(solve_method, H, f, u, y0, time)[source]

This function propagate the object y0 using either the Schrödinger equation with H = H0 + f(t)*V as Hamiltonian, or Lindblad’s equation. If an inverse time array is given, the solver take y0 as the state at final time, making a backward propagation.

The problem may be solved in the interaction picture if interaction_picture = True. Note that in that case H0 should be diagonal on entry (the system should be represented in the eigenbasis of H0). In this case, the output will also be in the interaction picture.

Parameters:

solve_method (string indicating the propagation method) –
H0 (hamiltonian) – Hamiltonian’s time independet part
V (list of functions) – Hamiltonian’s perturbation component
f (list of pulse) – pulse class defined in typical_pulse.py based on Fourier expansion parametrization
y0 – initial state
time (ndarray) – array that contain each time step.
u (ndarray) – array that contain the control parameters of the pulse, needed to build the args diccionary defined in pulse class
Returns –
....... –
qutip.result – qutiip.result type object with the propagation data

solvers.solve_h_is_function(solve_method, H, f, y0, time, returnQoutput=True, options=None, interaction_picture=False)[source]: Solver for the case in which the hamiltonian is given as a function

5.5. target

This module contains the Target class, and associated procedures.

class target.Target(targettype, Fyu=None, dFdy=None, dFdu=None, operator=None, Utarget=None, alpha=None, S=None, targeteps=None, T=None, fepsilon=None, dfdepsilon=None)[source]

The class that holds the definition of the target

Parameters:

targettype (string) –
The code admits the following types of characters:
1. generic
  
  The user may input any function as target, defined by the driver program, using the Fyu parameter (see below).
2. expectationvalue
  
  The target is the expectation value of some operator:
  
  \[F(y, u) = \langle \Psi(T)\vert O \vert \Psi(T)\rangle\]
  
  or, if the multitarget option is used,
  
  \[F(y, u) = \frac{1}{Q} \sum_i \langle \Psi_i(T)\vert O_i \vert \Psi_i(T)\rangle\]
  
  Here, \(Q\) is the number of targets. The operator(s) should be provided using the operator argument (below).
3. evolutionoperator
  
  The goal now is to find a set of control parameters that leads to an evolution of the system characterized by a fixed target evolution operator. This target evolution operator (or operators, if a multitarget setup is used), is provided by the argument Utarget. The target definition is:
  
  \[F(U, u) = \frac{1}{Q} \sum_i \vert U(T) \cdot U^{(i)}_{\rm target} \vert^2\]
4. floquet
Fyu (function) –
If targettype == “generic”, the user should pass the function used to define the target. The interface should be:
```
def Fyu(y: qutip.Qobj or list of qutip.Qobj, u: ndarray) -> float:
```
The function Fyu should receive as argument a system state (be it a Hilbert state, evolution operator, or density matrix), that is supposed to be the state at the end of the propagation, and a control parameters set. It should output the target functional value. It may also take as argument a list of states, if one is running the optimizaion in multitarget mode.
dFdy (function) –
If targettype == “generic”, the user should pass the function used to define the target (Fyu), and the derivative of that function with respect to the state, which would be given by this dFdy parameter. The interface should be:
```
def dFdy(y: qutip.Qobj or list of qutip.Qobj, u: ndarray) -> list of float:
```
dFdu (function) –
If targettype == “generic”, the user should pass the function used to define the target (Fyu), and, if it is not zero, the derivative of that function with respect to the control parameters, which would be given by this dFdu parameter. The interface should be:
```
def dFdu(u: ndarray, m: int) -> float:
```

Fyu(y, u)[source]

The functional F of the trayectory y that is to be maximized

It may also be an explicit function of the control parameters u.

Parameters:

y (qutip.result) – The trajectory of the quantum system
u (ndarray) – The control parameters that were used to generate the trajectory.

Returns:

The value of the target functional.

Return type:

float

dFdy(y, u)[source]

Derivative of the target functional wrt the quantum trajectory.

Right now, it in fact assumes that the functional only depends on the final state of the quantum trajectory. This functional derivative is used to determine the boundary condition used in the definition of the costate.

Parameters:

y (qutip.Qobj or list of qutip.Qobj) – The state of the system at the final time of the propagation
u (ndarray) – The control parameters

Returns:

The derivative of F with respect to the quantum state.

Return type:

qutip.Qobj or list of qutip.Qobj

5.6. floquet

This module contains all the procedures needed to do Floquet optimization.

(…)

class floquet.Nessopt(target_operator, T, nts, L, glist)[source]: Optimization of NESSs

floquet.dmtovector(rho)[source]

Transforms a Qobj density matrix object into a Qobj vector

Parameters:: rho (Qobj) – Density matrix
Returns:: A Qobj vector containing the vectorized version of the input density matrix
Return type:: Qobj

floquet.epsi(UT, T)[source]: (…)

floquet.epsilon(H, T, args=None)[source]: (…)

floquet.epsilon2(H0, H1, H2, Ax, Ay, u, T)[source]: (…)

floquet.epsilon3(H, f, u, T)[source]: (…)

floquet.fourier_steady_state2(T, ntsteps, L, f, nu=None, dgt=None)[source]: (…)

floquet.gradepsilon(H, f, u, T)[source]: (…)

floquet.vectortodm(rhovec)[source]: (…)

5.7. krotov

krotov.Coestate_eq(solve_method, H, f, psiF, times, O, obj, interaction_picture=False)[source]

This function calculate the coestate equation with the final condition described by

\[\vert\chi(T)\rangle = O\vert \Psi(t=T)\rangle\]

In the case of having a state propagating, and

\[\vert B(T)\rangle = Tr [U(T)+*O]*O/d^2\]

in the case of having an operator propagating.

The main idea is to solve it as an linear schrödinger equation but backwards, so we have to take care of the time and the pulse in order to make this correctly (with realtime parametrized pulses only is needed to invert the time). In principle this should work with any solve_method: rk4, cfm2 or cmf4.

krotov.Krotov(solve_method, H, f_initial, psi0, times, O, S, alpha=2, tolerance=0.001, verbose=False, maxeval=100, interaction_picture=False)[source]

This is the main function on the module.

We can divide it in 3 main steps (after some initializing):

Resolve a non-linear Schrödingers equation
Checking if the pulse have converged enough
Resolve a linear co-estate equation using the final condition of Schrödinger state (basically going backwards)

Parameters:

solve_method – string indicating the propagation method
H0 – Hamiltonian’s time independet part.
V – Hamiltonian’s perturbation component.
f_initial – pulse class defined in typical_pulse.py, specifically ‘Realtime’ parametrization
psi0 – initial state.
time – array that contain each time step.
O – Target operator that we want to optimize \(\langle\Psi(T)\vert O\vert\Psi(T)\rangle\)
S – Penalty function
tolerance – Precision for the algorithm
alpha – Penalty factor

Return type:

pulse

krotov.cfmagnus2_nonlinear(H0, V, f, chi, psi, times, alpha=2)[source]: This function solves the nonlinear schordinger equation using cfmagnus2 as the integrator algorithm.

krotov.cfmagnus4_nonlinear(H0, V, f, chi, psi0, time, S, alpha=2, interaction_picture=False, cops=None)[source]: This function solves the nonlinear schordinger equation using cfmagnus4 as the integrator algorithm.

krotov.check_pulse(f, times)[source]

This function checks if the pulse parametrization is realtime.

If it is not, it convert the pulse into realtime parametrization.

krotov.coestate_interpolator(chi, times, coestate_flag, solve_method)[source]

This function interpolates the value of the coestate.

In this way, we can solve properly the non-linear Schrodinger equation with the rk4 and cfmagnus4 integrators. This function returns the value of the coestate in the time that´s required

krotov.have_converged(old, new, times, epsilon=1e-06)[source]

This function checks if the algorithm have converged enough.

The idea is to check every time of the pulse.

Parameters:

old – The pulse before solving the non-linear schrodinger
new – The pulse after solving the non-linear schrodinger
epsilon – The tolerance we have so we can say the pulse have converged enough

Returns:

A boolean which is ‘True’ if every amplitude have converged, and is ‘False’ if any
of the amplitudes have not converged enough

krotov.non_linear_schr(solve_method, H0, V, f, chi, psi0, times, S, alpha=2, interaction_picture=False)[source]

This function solves the non-linear schrödinger equation that make the step 1 of the krotov´s algorithm.

The main idea is to solve the ecuation time by time, calculating the pulse every time-step. (Must check if is good)

For now, im just sticking to solve by rk4, but later i think is a good idea that this should be generaliced for cfmagnus4 and cfmagnus2.

krotov.rk4_nonlinear(H0, V, f, chi, psi0, times, S, alpha=2, interaction_picture=False)[source]

This function solves the nonlinear schordinger equation using rk4 as the integrator algorithm.

It is a modification of the function rk4solver in the solvers module of qocttools, so it can solve a slightly different equation and it can return the new state and the new pulse.

5.8. math_extra

math_extra.diff_ridders(function, origin, stepsize, con=1.4, safe=2.0, maxiter=15)[source]

Estimate first-order derivative with Ridders’ finite difference method.

This implementation is based on the one from the book Numerical Recipes. The code is pythonized and no longer using fixed-size arrays. Also, the output of the function can be an array.

Parameters:

function (function) – The function to be differentiated.
origin (float) – The point at which must be differentiated.
stepsize (float) – The initial step size.
con (float) – The rate at which the step size is decreased (contracted). Must be larger than one.
safe (float) – The safety check used to terminate the algorithm. If Errors between successive orders become larger than safe times the error on the best estimate, the algorithm stop. This happens due to round-off errors.
maxiter (int) – The maximum number of iterations, equals the maximum number of function calls and also the highest polynomial order in the Neville method.

Returns:

estimate (float) – The best estimate of the first-order derivative.
error (float) – The (optimistic) estimate of the error on the derivative.

math_extra.frobenius_product(A, B)[source]: return the Frobenius product between the operators A and B

math_extra.infidelity(U, dmi, dmo)[source]

Returns a measure of the infidelity of a quantum process.

(…)

math_extra.timegrid(H0, T, delta)[source]: Generates a time grid in the interval [0, T]