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
- 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:
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.
- 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.
- 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.
- gradf(t)[source]
Grdient of the function at time t
“Gradient of the function” means the gradient with respect to all the control parameters.
- 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:
- 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:
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
5.5. target
- class target.Target(targettype, Fyu=None, dFdy=None, dFdu=None, Pu=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: …
- 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
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 optimice <ψ(T)|O|ψ(T)>
S – Penalty function
tolerance – Precision for the algorithm
alpha – Penalty factor
- Return type:
- 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