Code documentation

system_builder

system_builder.diffusion_op(dim, C_vector, triple_prods, basis_norms_sq, basis, **kwargs)

Return a matrix \(D\) such that when \(\rho\) is vectorized the expression

\[\frac{d\rho}{dt}=\mathcal{D}[c]\rho=c\rho c^\dagger- \frac{1}{2}(c^\dagger c\rho+\rho c^\dagger c)\]

can be calculated by:

\[\frac{d\vec{\rho}}{dt}=D\vec{\rho}\]

Vectorization is done according to the order prescribed in basis, with the component proportional to identity in the last place.

Parameters:

• coupling_op (numpy.array) – The operator \(c\) in matrix form
• basis (list(numpy.array)) – An almost complete (minus identity), Hermitian, traceless, orthogonal basis for the operators (does not need to be normalized).

Returns:

The matrix \(D\) operating on a vectorized density operator

Return type:

numpy.array

system_builder.double_comm_op(dim, C_vector, triple_prods, M_sq, basis_norms_sq, basis, **kwargs)

Return a matrix \(E\) such that when \(\rho\) is vectorized the expression

\[\frac{d\rho}{dt}=\left(\frac{M^*}{2}[c,[c,\rho]]+ \frac{M}{2}[c^\dagger,[c^\dagger,\rho]]\right)\]

can be calculated by:

\[\frac{d\vec{\rho}}{dt}=E\vec{\rho}\]

Vectorization is done according to the order prescribed in basis, with the component proportional to identity in the last place.

Parameters:

• coupling_op (numpy.array) – The operator \(c\) in matrix form
• M_sq (complex) – Complex squeezing parameter \(M\) defined by \(\langle dB(t)dB(t)\rangle=Mdt\).
• basis (list(numpy.array)) – An almost complete (minus identity), Hermitian, traceless, orthogonal basis for the operators (does not need to be normalized).

Returns:

The matrix \(E\) operating on a vectorized density operator

Return type:

numpy.array

system_builder.dualize(operator, basis)

Take an operator to its dual vectorized form in a particular operator basis.

Designed to work in conjunction with vectorize so that, given an orthogonal basis \(\{\Lambda^m\}\) where \(\operatorname{Tr}[{\Lambda^m}^\dagger\Lambda^n]\propto\delta_{mn}\), the dual action of an operator \(A\) on another operator \(B\) interpreted as \(\operatorname{Tr}[A^\dagger B]\) can be easily calculated by sum([a*b for a, b in zip(dualize(A), vectorize(B))] (in other words it becomes an ordinaty dot product in this particular representation).

Parameters:

• operator (numpy.array) – The operator to vectorize
• basis (list(numpy.array)) – The basis to vectorize the operator in

Returns:

The vector components

Return type:

numpy.array

system_builder.hamiltonian_op(dim, H_vector, double_prods, basis_norms_sq, basis, **kwargs)

Return a matrix \(F\) such that when \(\rho\) is vectorized the expression

\[\frac{d\rho}{dt}=-i[H,\rho]\]

can be calculated by:

\[\frac{d\vec{\rho}}{dt}=F\vec{\rho}\]

Vectorization is done according to the order prescribed in basis, with the component proportional to identity in the last place.

Parameters:

• hamiltonian (numpy.array) – The Hamiltonian \(H\) in matrix form
• basis (list(numpy.array)) – An almost complete (minus identity), Hermitian, traceless, orthogonal basis for the operators (does not need to be normalized).

Returns:

The matrix \(F\) operating on a vectorized density operator

Return type:

numpy.array

system_builder.norm_squared(operator)

Returns the square of the Frobenius norm of the operator.

Parameters:operator (numpy.array) – The operator for which to calculate the squared norm Returns:The square of the norm of the operator Return type:Positive real

system_builder.op_calc_setup(coupling_op, M_sq, N, H, partial_basis)

Handle the repeated tasks performed every time a superoperator matrix is computed.

system_builder.recur_dot(mats)

Perform numpy.dot on a list in a right-associative manner.

system_builder.vectorize(operator, basis)

Vectorize an operator in a particular operator basis.

Parameters:

• operator (numpy.array) – The operator to vectorize
• basis (list(numpy.array)) – The basis to vectorize the operator in

Returns:

The vector components

Return type:

numpy.array

system_builder.weiner_op(dim, C_vector, double_prods, basis_norms_sq, basis, **kwargs)

Return a the matrix-vector pair \((G,\vec{k})\) such that when \(\rho\) is vectorized the expression

\[d\rho=dW\,\left(c\rho+\rho c^\dagger- \rho\operatorname{Tr}[(c+c^\dagger)\rho]\right)\]

can be calculated by:

\[d\vec{\rho}=dW\,(G+\vec{k}\cdot\vec{\rho})\vec{\rho}\]

Vectorization is done according to the order prescribed in basis, with the component proportional to identity in the last place.

Parameters:

• coupling_op (numpy.array) – The operator \(c\) in matrix form
• basis (list(numpy.array)) – An almost complete (minus identity), Hermitian, traceless, orthogonal basis for the operators (does not need to be normalized).

Returns:

The matrix-vector pair \((G,\vec{k})\) operating on a vectorized density operator (k is returned as a row-vector)

Return type:

tuple(numpy.array)

gellmann

gellmann.gellmann(j, k, d)

Returns a generalized Gell-Mann matrix of dimension d. According to the convention in Bloch Vectors for Qubits by Bertlmann and Krammer (2008), returns \(\Lambda^j\) for \(1\leq j=k\leq d-1\), \(\Lambda^{kj}_s\) for \(1\leq k<j\leq d\), \(\Lambda^{jk}_a\) for \(1\leq j<k\leq d\), and \(I\) for \(j=k=d\).

Parameters:

• j (positive integer) – First index for generalized Gell-Mann matrix
• k (positive integer) – Second index for generalized Gell-Mann matrix
• d (positive integer) – Dimension of the generalized Gell-Mann matrix

Returns:

A genereralized Gell-Mann matrix.

Return type:

numpy.array

gellmann.get_basis(d)

Return a basis of orthogonal Hermitian operators on a Hilbert space of dimension d, with the identity element in the last place.

gramschmidt

gramschmidt.orthonormalize(A)

Return an orthonormal basis where A has support only on the first three elements. The first element is guaranteed to be proportional to the identity.

integrate

sde

sde.euler(drift_fn, diffusion_fn, X0, ts, Us)

Integrate a system of ordinary stochastic differential equations subject to scalar noise:

\[d\vec{X}=\vec{a}(\vec{X},t)\,dt+\vec{b}(\vec{X},t)\,dW_t\]

Uses the Euler method:

\[\vec{X}_{i+1}=\vec{X}_i+\vec{a}(\vec{X}_i,t_i)\Delta t_i+ \vec{b}(\vec{X}_i,t_i)\Delta W_i\]

where \(\Delta W_i=U_i\sqrt{\Delta t}\), \(U\) being a normally distributed random variable with mean 0 and variance 1.

Parameters:

• drift_fn (callable(X, t)) – Computes the drift coefficient \(\vec{a}(\vec{X},t)\)
• diffusion_fn (callable(X, t)) – Computes the diffusion coefficient \(\vec{b}(\vec{X},t)\)
• X0 (array) – Initial condition on X
• ts (array) – A sequence of time points for which to solve for X. The initial value point should be the first element of this sequence.
• Us (array, shape=(len(t) - 1)) – Normalized Weiner increments for each time step (i.e. samples from a Gaussian distribution with mean 0 and variance 1).

Returns:

Array containing the value of X for each desired time in t, with the initial value X0 in the first row.

Return type:

numpy.array, shape=(len(ts), len(X0))

sde.faulty_milstein(drift, diffusion, b_dx_b, X0, ts, Us)

Integrate a system of ordinary stochastic differential equations subject to scalar noise:

\[d\vec{X}=\vec{a}(\vec{X},t)\,dt+\vec{b}(\vec{X},t)\,dW_t\]

Uses a faulty Milstein method (i.e. missing the factor of 1/2 in the term added to Euler integration):

\[\vec{X}_{i+1}=\vec{X}_i+\vec{a}(\vec{X}_i,t_i)\Delta t_i+ \vec{b}(\vec{X}_i,t_i)\Delta W_i+ \left(\vec{b}(\vec{X}_i,t_i)\cdot\vec{\nabla}_{\vec{X}}\right) \vec{b}(\vec{X}_i,t_i)\left((\Delta W_i)^2-\Delta t_i\right)\]

where \(\Delta W_i=U_i\sqrt{\Delta t}\), \(U\) being a normally distributed random variable with mean 0 and variance 1.

Parameters:

• drift (callable(X, t)) – Computes the drift coefficient \(\vec{a}(\vec{X},t)\)
• diffusion (callable(X, t)) – Computes the diffusion coefficient \(\vec{b}(\vec{X},t)\)
• b_dx_b (callable(X, t)) – Computes the correction coefficient \(\left(\vec{b}(\vec{X},t)\cdot \vec{\nabla}_{\vec{X}}\right)\vec{b}(\vec{X},t)\)
• X0 (array) – Initial condition on X
• ts (array) – A sequence of time points for which to solve for X. The initial value point should be the first element of this sequence.
• Us (array, shape=(len(t) - 1)) – Normalized Weiner increments for each time step (i.e. samples from a Gaussian distribution with mean 0 and variance 1).

Returns:

Array containing the value of X for each desired time in t, with the initial value X0 in the first row.

Return type:

numpy.array, shape=(len(t), len(X0))

sde.meas_euler(drift_fn, diffusion_fn, dW_fn, X0, ts, dMs)

Integrate a system of ordinary stochastic differential equations conditioned on an incremental measurement record:

\[d\vec{X}=\vec{a}(\vec{X},t)\,dt+\vec{b}(\vec{X},t)\,dW_t\]

Uses the Euler method:

\[\vec{X}_{i+1}=\vec{X}_i+\vec{a}(\vec{X}_i,t_i)\Delta t_i+ \vec{b}(\vec{X}_i,t_i)\Delta W_i\]

where \(\Delta W_i=f(\Delta M_i,\vec{X}, t)\), \(\Delta M_i\) being the incremental measurement record being used to drive the SDE.

Parameters:

• drift_fn (callable(X, t)) – Computes the drift coefficient \(\vec{a}(\vec{X},t)\)
• diffusion_fn (callable(X, t)) – Computes the diffusion coefficient \(\vec{b}(\vec{X},t)\)
• dW_fn (callable(dM, dt, X, t)) – The function that converts the incremental measurement and current state to the Wiener increment.
• X0 (array) – Initial condition on X
• ts (array) – A sequence of time points for which to solve for X. The initial value point should be the first element of this sequence.
• dMs (array, shape=(len(t) - 1)) – Incremental measurement outcomes used to drive the SDE.

Returns:

Array containing the value of X for each desired time in t, with the initial value X0 in the first row.

Return type:

numpy.array, shape=(len(ts), len(X0))

sde.meas_milstein(drift_fn, diffusion_fn, b_dx_b_fn, dW_fn, X0, ts, dMs)

Integrate a system of ordinary stochastic differential equations conditioned on an incremental measurement record:

\[d\vec{X}=\vec{a}(\vec{X},t)\,dt+\vec{b}(\vec{X},t)\,dW_t\]

Uses the Milstein method:

\[\vec{X}_{i+1}=\vec{X}_i+\vec{a}(\vec{X}_i,t_i)\Delta t_i+ \vec{b}(\vec{X}_i,t_i)\Delta W_i+ \frac{1}{2}\left(\vec{b}(\vec{X}_i,t_i)\cdot\vec{\nabla}_{\vec{X}}\right) \vec{b}(\vec{X}_i,t_i)\left((\Delta W_i)^2-\Delta t_i\right)\]

where \(\Delta W_i=f(\Delta M_i,\vec{X}, t)\), \(\Delta M_i\) being the incremental measurement record being used to drive the SDE.

Parameters:

• drift_fn (callable(X, t)) – Computes the drift coefficient \(\vec{a}(\vec{X},t)\)
• diffusion_fn (callable(X, t)) – Computes the diffusion coefficient \(\vec{b}(\vec{X},t)\)
• b_dx_b_fn (callable(X, t)) – Computes the correction coefficient \(\left(\vec{b}(\vec{X},t)\cdot \vec{\nabla}_{\vec{X}}\right)\vec{b}(\vec{X},t)\)
• dW_fn (callable(dM, dt, X, t)) – The function that converts the incremental measurement and current state to the Wiener increment.
• X0 (array) – Initial condition on X
• ts (array) – A sequence of time points for which to solve for X. The initial value point should be the first element of this sequence.
• dMs (array, shape=(len(t) - 1)) – Incremental measurement outcomes used to drive the SDE.

Returns:

Array containing the value of X for each desired time in t, with the initial value X0 in the first row.

Return type:

numpy.array, shape=(len(ts), len(X0))

sde.milstein(drift, diffusion, b_dx_b, X0, ts, Us)

Integrate a system of ordinary stochastic differential equations subject to scalar noise:

\[d\vec{X}=\vec{a}(\vec{X},t)\,dt+\vec{b}(\vec{X},t)\,dW_t\]

Uses the Milstein method:

\[\vec{X}_{i+1}=\vec{X}_i+\vec{a}(\vec{X}_i,t_i)\Delta t_i+ \vec{b}(\vec{X}_i,t_i)\Delta W_i+ \frac{1}{2}\left(\vec{b}(\vec{X}_i,t_i)\cdot\vec{\nabla}_{\vec{X}}\right) \vec{b}(\vec{X}_i,t_i)\left((\Delta W_i)^2-\Delta t_i\right)\]

where \(\Delta W_i=U_i\sqrt{\Delta t}\), \(U\) being a normally distributed random variable with mean 0 and variance 1.

Parameters:

• drift (callable(X, t)) – Computes the drift coefficient \(\vec{a}(\vec{X},t)\)
• diffusion (callable(X, t)) – Computes the diffusion coefficient \(\vec{b}(\vec{X},t)\)
• b_dx_b (callable(X, t)) – Computes the correction coefficient \(\left(\vec{b}(\vec{X},t)\cdot \vec{\nabla}_{\vec{X}}\right)\vec{b}(\vec{X},t)\)
• X0 (array) – Initial condition on X
• ts (array) – A sequence of time points for which to solve for X. The initial value point should be the first element of this sequence.
• Us (array, shape=(len(t) - 1)) – Normalized Weiner increments for each time step (i.e. samples from a Gaussian distribution with mean 0 and variance 1).

Returns:

Array containing the value of X for each desired time in t, with the initial value X0 in the first row.

Return type:

numpy.array, shape=(len(ts), len(X0))

sde.time_ind_taylor_1_5(drift, diffusion, b_dx_b, b_dx_a, a_dx_b, a_dx_a, b_dx_b_dx_b, b_b_dx_dx_b, b_b_dx_dx_a, X0, ts, U1s, U2s)

Integrate a system of ordinary stochastic differential equations with time-independent coefficients subject to scalar noise:

\[d\vec{X}=\vec{a}(\vec{X})\,dt+\vec{b}(\vec{X})\,dW_t\]

Uses an order 1.5 Taylor method:

\[\begin{split}\begin{align} \rho^\mu_{i+1}&=\rho^\mu_i+a^\mu_i\Delta t_i+ b^\mu_i\Delta W_i+\frac{1}{2}b^\nu_i\partial_\nu b^\mu_i\left( (\Delta W_i)^2-\Delta t_i\right)+ \\ &\quad b^\nu_i\partial_\nu a^\mu_i\Delta Z_i +\left(a^\nu_i\partial_\nu +\frac{1}{2}b^\nu_ib^\sigma_i\partial_\nu\partial_\sigma\right) b^\mu_i\left(\Delta W_i\Delta t_i-\Delta Z_i\right)+ \\ &\quad\frac{1}{2}\left(a^\nu_i\partial_\nu +\frac{1}{2}b^\nu_ib^\sigma_i\partial_\nu\partial_\sigma\right) a^\mu_i\Delta t_i^2 +\frac{1}{2}b^\nu_i\partial_\nu b^\sigma_i\partial_\sigma b^\mu_i\left( \frac{1}{3}(\Delta W_i)^2-\Delta t_i\right)\Delta W_i \end{align}\end{split}\]

Parameters:

• drift (callable(X)) – Computes the drift coefficient \(\vec{a}(\vec{X})\)
• diffusion (callable(X)) – Computes the diffusion coefficient \(\vec{b}(\vec{X})\)
• b_dx_b (callable(X)) – Computes the coefficient \(\left(\vec{b}(\vec{X})\cdot \vec{\nabla}_{\vec{X}}\right)\vec{b}(\vec{X})\)
• b_dx_a (callable(X)) – Computes the coefficient \(\left(\vec{b}(\vec{X})\cdot \vec{\nabla}_{\vec{X}}\right)\vec{a}(\vec{X})\)
• a_dx_b (callable(X)) – Computes the coefficient \(\left(\vec{a}(\vec{X})\cdot \vec{\nabla}_{\vec{X}}\right)\vec{b}(\vec{X})\)
• a_dx_a (callable(X)) – Computes the coefficient \(\left(\vec{a}(\vec{X})\cdot \vec{\nabla}_{\vec{X}}\right)\vec{a}(\vec{X})\)
• b_dx_b_dx_b (callable(X)) – Computes the coefficient \(\left(\vec{b}(\vec{X})\cdot \vec{\nabla}_{\vec{X}}\right)^2\vec{b}(\vec{X})\)
• b_b_dx_dx_b – Computes \(b^\nu b^\sigma\partial_\nu\partial_\sigma b^\mu\hat{e}_\mu\).
• b_b_dx_dx_a – Computes \(b^\nu b^\sigma\partial_\nu\partial_\sigma a^\mu\hat{e}_\mu\).
• X0 (array) – Initial condition on X
• ts (array) – A sequence of time points for which to solve for X. The initial value point should be the first element of this sequence.
• U1s (array, shape=(len(t) - 1)) – Normalized Weiner increments for each time step (i.e. samples from a Gaussian distribution with mean 0 and variance 1).
• U2s (array, shape=(len(t) - 1)) – Normalized Weiner increments for each time step (i.e. samples from a Gaussian distribution with mean 0 and variance 1).

Returns:

Array containing the value of X for each desired time in t, with the initial value X0 in the first row.

Return type:

numpy.array, shape=(len(t), len(X0))

grid_conv

Functions for testing convergence rates using grid convergence

grid_conv.calc_rate(integrator, rho_0, times, U1s=None, U2s=None)

Calculate the convergence rate for some integrator.

Parameters:

• integrator – An Integrator object.
• rho_0 (numpy.array) – The initial state of the system
• times – Sequence of times (assumed to be evenly spaced, defining a number of increments divisible by 4).
• U1s (numpy.array(len(times) - 1)) – Samples from a standard-normal distribution used to construct Wiener increments \(\Delta W\) for each time interval. If not provided will be generated by the function.
• U2s (numpy.array(len(times) - 1)) – Samples from a standard-normal distribution used to construct multiple-Ito increments \(\Delta Z\) for each time interval. If not provided will be generated by the function.

Returns:

The convergence rate as a power of \(\Delta t\).

Return type:

float

grid_conv.double_increments(times, U1s, U2s=None)

Take a list of times (assumed to be evenly spaced) and standard-normal random variables used to define the Ito integrals on the intervals and return the equivalent lists for doubled time intervals. The new standard-normal random variables are defined in terms of the old ones by

Parameters:

• times (numpy.array) – List of evenly spaced times defining an even number of time intervals.
• U1s (numpy.array(N, len(times) - 1)) – Samples from a standard-normal distribution used to construct Wiener increments \(\Delta W\) for each time interval. Multiple rows may be included for independent trajectories.
• U2s (numpy.array(N, len(times) - 1)) – Samples from a standard-normal distribution used to construct multiple-Ito increments \(\Delta Z\) for each time interval. Multiple rows may be included for independent trajectories.

Returns:

Times sampled at half the frequency and the modified standard-normal-random-variable samples for the new intervals. If U2s=None, only new U1s are returned.

Return type:

(numpy.array(len(times)//2 + 1), numpy.array(len(times)//2)[, numpy.array(len(times)//2]))