magnus.magnus

magnus.py

Compute the time-evolution operator using the Magnus expansion.

This module contains the numerical core of Magnus: routines to compute the matrix exponential of the Magnus expansion of a (possibly time-dependent) matrix function \(A(t)\), i.e.,

\[U(t_1, t_0) = \exp\!\left[\Omega_1 + \Omega_2 + \cdots + \Omega_k\right] ,\]

where the terms \(\Omega_k\) are built from time-ordered integrals of nested commutators of \(A(t)\). For neutrino oscillations, \(A(t) = -i H(t)\), with \(H(t)\) the Hamiltonian, but the routines below work for arbitrary matrix-valued \(A(t)\).

The terms are generated with the standard recursion based on Bernoulli numbers [1] (in the \(B_1 = -1/2\) convention):

\[\begin{split}\Omega_1(t) &= \int_0^t A(s)\, ds \\ \Omega_n(t) &= \sum_{j=1}^{n-1} \frac{B_j}{j!} \int_0^t S_n^{(j)}(s)\, ds ,\end{split}\]

with \(S_n^{(j)}\) the sums of nested commutators of the lower-order terms with \(A\). Orders 1–6 are implemented (\(B_3 = B_5 = 0\), so those groups vanish identically).

Two families of methods are available, selected via integration_method:

'gl' (the default): Gauss-Legendre commutator-free collocation [1] [2]. For a slab of width \(h\) it needs only 1, 2, or 3 evaluations of \(A\) to reach order 2, 4, or 6, respectively, with quadrature error matched to the truncation order. n_tpts is ignored. Both faster and more accurate than the alternatives whenever \(A(t)\) is smooth within each slab, which is the common case – and, for the Earth, is what aligning slab edges with the PREM layer boundaries is for.
'trapezoid' / 'simpson': sample \(A(t)\) on a uniform grid of n_tpts points and evaluate the nested integrals with cumulative quadrature. Fully general, and so the safer choice if \(A(t)\) has a kink or a discontinuity inside a slab, where Gauss-Legendre loses its order advantage. The quadrature error (\(\mathcal{O}(h^2)\) or \(\mathcal{O}(h^4)\)) can dominate the Magnus truncation error at high orders unless n_tpts grows accordingly.

References

Routine listings

commutator - Returns [X, Y] = XY - YX

probe_eval_mode - Determines how a matrix function can be evaluated

suggest_n_slabs - Suggests a starting number of time slabs

magnus_expansion - Computes \(\exp(\Omega)\) for a single time slab

evolution_operators_from_samples - Evolution operators of a chain
of slabs from precomputed samples of A

gl_nodes - Returns the Gauss-Legendre nodes used by the ‘gl’ method

magnus_expansion_multislab - Computes the evolution operators of
all time slabs at once, from A directly

MagnusConvergenceWarning - Warning class for slabs too wide for
guaranteed Magnus convergence

Attributes

`B`
`F1`
`F2`
`f1`
`f2`
`MAGNUS_EXP_ORDER_MAX`
`valid_integration_methods`
`valid_expm_backends`	The accepted values of `EXPM_BACKEND` and of every
`EXPM_BACKEND`	Module-level switch selecting how \(\exp(\Omega)\) is computed.
`USE_PALINDROME`	Module-level switch.

Exceptions

`MagnusHighOrderCostWarning`	Warns that a Magnus order above 6 costs substantially more per slab.
`ScalarHamiltonianWarning`	Warns that `H_func` accepts only one position at a time.
`MagnusConvergenceWarning`	Warns that a time slab may be too wide for the Magnus series.

Functions

`commutator`(→ numpy.ndarray)	Returns the commutator [X, Y] = X Y - Y X.
`cached_eval_mode`(→ str)	`probe_eval_mode` for a callable that will be probed more than once.
`probe_eval_mode`(→ str)	Determine how the matrix function A can be evaluated.
`suggest_n_slabs`(→ int)	Suggest a starting number of time slabs for [t0, t1].
`ordered_product`(→ numpy.ndarray)	Time-ordered product of a stack of slab operators, earliest slab first.
`magnus_expansion`(→ numpy.ndarray)	Compute \(\exp(\Omega_1 + \cdots + \Omega_\text{order})\) of \(A(t)\) from
`evolution_operators_from_samples`(→ numpy.ndarray)	Evolution operators of a chain of slabs from precomputed samples.
`gl_nodes`(→ numpy.ndarray)	Returns the Gauss-Legendre nodes on [0, 1] used by the 'gl' method.
`palindromic`(→ bool)	Returns whether every array given reads the same both ways.
`magnus_expansion_multislab`(→ numpy.ndarray)	Compute the evolution operators of all time slabs at once.

Module Contents

exception magnus.magnus.MagnusHighOrderCostWarning[source]

Bases: UserWarning

Warns that a Magnus order above 6 costs substantially more per slab.

Orders 1-6 are written out inline with their shared subexpressions named and reused. Above that the terms are generated from the recursion, and their number roughly doubles per order (9 terms at order 6; 17, 33, 65, 129 at orders 7-10), so the work per slab grows with it – measured at roughly 2.7x order 6 at order 7, rising to about 17x at order 10, for the same grid.

Higher order buys a genuinely faster convergence rate in the slab width, so this is a trade rather than a mistake. But it is often the worse side of the trade: narrowing the slabs at order 4 or 6 usually reaches a given accuracy for less total work, and beyond the Magnus series’ convergence radius no order helps at all (see MagnusConvergenceWarning).

Added in version 1.0.0.

exception magnus.magnus.ScalarHamiltonianWarning[source]

Bases: UserWarning

Warns that H_func accepts only one position at a time.

The engine evaluates the Hamiltonian at every quadrature node of every slab – often a few hundred positions for a single probability, and the adaptive refinement repeats that at each level. _evaluate_A therefore tries a single vectorized call, A(times), and uses the result if it has the right shape and agrees with a scalar spot-check. If that fails it falls back to a Python loop, one call per position.

That fallback is correct but typically several times slower, and it is silent: nothing about a scalar-only H_func looks wrong, so the slow path is easy to sit on indefinitely. Measured on a three-flavor exponential-density profile, making the same H_func array-capable cut the time per magnus.oscprob.osc_prob() call from 7.8 ms to 1.7 ms, a factor of 4.6, with bit-identical output.

Making a Hamiltonian array-capable usually means no more than writing its position dependence with NumPy and letting the matrix part broadcast:

# slow: one position at a time
def H_func(l):
    VCC = matter.VCC_func(l, num_density_e_func)
    return (1.0/energy)*h_vac + hamiltonians.hamiltonian_3nu_matter(VCC)

# fast: the same physics, evaluated for all positions at once
e00 = np.diag([1.0, 0.0, 0.0])
def H_func(l):
    l = np.asarray(l, dtype=float)
    VCC = vcc_of(l)                      # returns an array
    return (1.0/energy)*h_vac + VCC[..., None, None]*e00

The trailing [..., None, None] is what lets one potential per position multiply a stack of matrices. A Hamiltonian that ignores its argument entirely is detected separately and costs nothing, so constant-density cases never trigger this.

Pass A_eval_mode='scalar' to magnus_expansion() (or accept the warning) when a scalar-only Hamiltonian is genuinely unavoidable.

Added in version 1.0.0.

exception magnus.magnus.MagnusConvergenceWarning[source]

Bases: UserWarning

Warns that a time slab may be too wide for the Magnus series.

What was detected. The Magnus series is guaranteed to converge when \(\int_{t_0}^{t_1} \lVert A(t)\rVert_2\, dt < \pi\). \(\lVert\Omega\rVert_2 \geq \pi\) is used as a cheap proxy for that integral – it comes free from the eigenvalues already computed for the matrix exponential – so this fires when a sufficient condition for convergence was not met on at least one slab. The message says how far past \(\pi\), in three buckets, which is the one quantity this check actually knows.

What it means for the answer: unknown, and that is the honest answer. This is a statement about the slab width, not about the error. The condition is sufficient, not necessary, so exceeding it does not imply a wrong answer – and it fires on results accurate to 1.6e-06 (docs/dev/DECISION_DISPATCH_ORDER.md §5) as well as on results seven times outside a requested 1e-3. Anything that claims to tell you which of those you have is claiming more than this check can support; magnus.oscprob.ToleranceNotAchievedWarning is the one that reports a failed convergence test.

What to change. More, narrower slabs: request a smaller rtol/atol, or raise n_slabs. Raising magnus_exp_order does not help in this regime – beyond the series’ radius no order converges. If the profile has a density jump or a kink, pass t_breakpoints there as well: a slab straddling one is never fixed by more slabs, only narrowed.

When it is safe to ignore. When the answer has been checked another way – a tighter tolerance giving the same result, or magnus.oscprob.cross_check_strategies() showing a different engine agreeing. Not merely because a tolerance was requested. That advice used to be in this message and it is false in exactly the cases where the warning matters: measured on a sawtooth density with rtol=atol=1e-3 explicitly requested, under both strategy='auto' and strategy='magnus', the adaptive refinement ran and the answer was still 7.484e-03, seven times outside the tolerance asked for, with this warning showing.

Measured rates (docs/dev/adversarial_batteries/warn_fp.py, 168 configurations across the profile families this package serves, d = 2-5, scored against solve_ivp or, for piecewise profiles, against expm): fired 70 times, of which 17 true positives and 53 false positives – a 76 % false-positive rate, the highest of any warning here. That is the price of reporting a sufficient condition, and it is why the text above refuses to translate the condition into a claim about the error.

Where that noise comes from, and what would fix it. Of 66 single-point calls, some refinement level exceeded \(\pi\) in 46 – but the level whose answer was actually returned did so in only 7. So 39 of 46 firings, 85 %, describe an intermediate grid that nobody receives: the ladder started coarse, said so, then refined and never retracted it. Keying the warning to the returned level alone would cut false alarms from 31 to 5 at a similar rate (67 % against 71 %). That change is mechanical – capture the norm per level and emit once the loop has decided – and is deliberately not made here, because it touches the refinement loop and the warning plumbing several tests depend on. It is written down with its numbers so it can be made deliberately rather than rediscovered.

Added in version 1.0.0.

magnus.magnus.B[source]

magnus.magnus.F1 = 0.08333333333333333[source]

magnus.magnus.F2 = -0.001388888888888889[source]

magnus.magnus.f1 = 0.08333333333333333[source]

magnus.magnus.f2 = -0.001388888888888889[source]

magnus.magnus.MAGNUS_EXP_ORDER_MAX = 10[source]

magnus.magnus.valid_integration_methods = ['gl', 'trapezoid', 'simpson'][source]

magnus.magnus.commutator(X: numpy.ndarray, Y: numpy.ndarray) → numpy.ndarray[source]

Returns the commutator [X, Y] = X Y - Y X.

Works on single matrices and on stacks of matrices (the matrix product broadcasts over all leading axes).

Added in version 1.0.0.

Parameters:

X (np.ndarray) – Left matrix (or stack of matrices).
Y (np.ndarray) – Right matrix (or stack of matrices), broadcastable against X.

Returns:

The commutator X @ Y - Y @ X.

Return type:

np.ndarray

Examples

import numpy as np

from magnus import magnus

X = np.array([[0.0, 1.0], [0.0, 0.0]])
Y = np.array([[0.0, 0.0], [1.0, 0.0]])

print(magnus.commutator(X, Y))
print('antisymmetric:',
      np.array_equal(magnus.commutator(X, Y), -magnus.commutator(Y, X)))

[[ 1.  0.]
 [ 0. -1.]]
antisymmetric: True

magnus.magnus.cached_eval_mode(A: Callable, t0: float, t1: float, key=None) → str[source]

probe_eval_mode for a callable that will be probed more than once.

Returns the same value probe_eval_mode() would for this interval, and remembers it against (key or A, t0, t1) so a repeated call on the same Hamiltonian over the same span does not evaluate it three more times. Falls straight through for anything that cannot be weakly referenced or hashed.

key exists because callers often have to wrap the object they want cached: probe_eval_mode needs \(A = -iH\), and a fresh lambda t: -1j*H(t) per call would miss every time. Passing key=H_func caches against the thing whose signature is actually being described. Multiplying by a constant cannot change whether a function accepts an array, so the two share a verdict.

Added in version 1.0.0.

Parameters:

A (Callable) – The matrix function to probe.
t0 (float) – The interval to probe over. Part of the cache key: see the comment above _EVAL_MODE_CACHE for the wrong answer that omitting it produced.
t1 (float) – The interval to probe over. Part of the cache key: see the comment above _EVAL_MODE_CACHE for the wrong answer that omitting it produced.
key (optional) – Object to cache against instead of A. Defaults to A.

Returns:

‘vector’, ‘scalar’ or ‘constant’; see probe_eval_mode().

Return type:

str

magnus.magnus.probe_eval_mode(A: Callable, t0: float, t1: float, n_probe: int | None = 5) → str[source]

Determine how the matrix function A can be evaluated.

Returns ‘vector’ if A accepts an array of times (fast path), ‘constant’ if A ignores its argument, and ‘scalar’ otherwise. Use the result as the A_eval_mode argument of magnus_expansion() and magnus_expansion_multislab() to avoid re-probing A on every call.

Added in version 1.0.0.

Parameters:

A (Callable) – Matrix function of time; see magnus_expansion().
t0 (float) – Interval over which A is probed (t1 >= t0).
t1 (float) – Interval over which A is probed (t1 >= t0).
n_probe (int, optional) – Number of sample times used for the probe. Default: 5.

Returns:

‘vector’, ‘constant’, or ‘scalar’.

Return type:

str

Examples

import numpy as np

from magnus import magnus

def vectorised(t):
    return -1j*np.eye(2)*np.asarray(t)[..., None, None]

print('array-capable :', magnus.probe_eval_mode(vectorised, 0.0, 1.0))
print('constant      :', magnus.probe_eval_mode(lambda t: -1j*np.eye(2),
                                                0.0, 1.0))

array-capable : vector
constant      : constant

magnus.magnus.suggest_n_slabs(A: Callable, t0: float, t1: float, A_eval_mode: str | None = None, n_probe: int | None = 17, phase_per_slab: float | None = 2.0 * np.pi) → int[source]

Suggest a starting number of time slabs for [t0, t1].

Estimates the accumulated phase \(\lVert\Omega_1\rVert_2\) over the whole interval from a coarse sample of A (with the trace removed, since a global phase does not affect the probabilities) and suggests enough slabs to keep roughly phase_per_slab (radians) of phase per slab. Starting an adaptive refinement from this estimate skips most of the geometric ladder that would otherwise climb from a single slab.

The default of \(2\pi\) radians per slab is deliberately looser than the Magnus convergence guarantee (\(\pi\)): empirically, for smooth profiles, order-4 methods reach ~1e-3 accuracy already at this slab width, and the adaptive refinement loop – which remains the sole arbiter of accuracy – grows the slab count from here when the requested tolerance demands it.

Added in version 1.0.0.

Parameters:

A (Callable) – Matrix function of time; see magnus_expansion().
t0 (float) – Interval over which the phase is estimated (t1 >= t0).
t1 (float) – Interval over which the phase is estimated (t1 >= t0).
A_eval_mode (str, optional) – Skip probing how A can be evaluated; see probe_eval_mode().
n_probe (int, optional) – Number of sample points used to estimate the accumulated phase. Default: 17.
phase_per_slab (float, optional) – Target accumulated phase per slab, in radians. Default: \(2\pi\).

Returns:

Suggested starting number of slabs (at least 1).

Return type:

int

Examples

import numpy as np

from magnus import magnus

H = np.array([[0.0, 1.0], [1.0, 0.0]])
print('slabs suggested:',
      magnus.suggest_n_slabs(lambda t: -1j*20.0*H, 0.0, 1.0))

slabs suggested: 4

magnus.magnus.ordered_product(U: numpy.ndarray) → numpy.ndarray[source]

Time-ordered product of a stack of slab operators, earliest slab first.

Returns \(U_{n-1} \cdots U_1 U_0\) for U of shape (n, d, d) – the same quantity as functools.reduce(np.matmul, U[::-1]) and, because matrix multiplication is associative, the same value.

The difference is how it gets there. reduce walks the stack one matrix at a time, which is \(n-1\) separate Python-level calls into NumPy for matrices of size 3; the array was already materialised as a single (n, d, d) block, so nearly all of that time is call overhead rather than arithmetic. Multiplying adjacent pairs instead collapses the stack in \(\lceil\log_2 n\rceil\) batched matmuls, each of which does its whole level in one call.

Measured on unitary 3x3 stacks: 92 -> 29 us at n = 108, 1764 -> 370 us at n = 2048, agreeing to 8e-16 with no systematic loss of unitarity. The gain grows with the slab count, which is the direction the adaptive refinement moves in.

Associativity is what makes this legitimate; commutativity is not required and is not assumed. Adjacent pairs are combined in order, so the operator ordering is preserved exactly – an odd element is carried forward untouched rather than being folded in out of turn.

Added in version 1.0.0.

Parameters:: U (np.ndarray) – Stack of operators, shape (n, d, d), ordered earliest slab first.
Returns:: The ordered product, shape (d, d).
Return type:: np.ndarray

magnus.magnus.valid_expm_backends = ['auto', 'numba', 'eigh'][source]

The accepted values of EXPM_BACKEND and of every expm_backend parameter.

Type:: list of str

magnus.magnus.EXPM_BACKEND = 'auto'[source]

Module-level switch selecting how \(\exp(\Omega)\) is computed.

Which routine exponentiates each slab. This is not a correctness switch: the two backends agree to about 1e-15 wherever the kernel is used, which is the accuracy either one has – and where it would not, it is not used: the kernel reports the conditioning of its own characteristic cubic and eigh answers instead. See magnus.expmkernels.SEV_TOL.

'auto' (the default): the compiled Cayley-Hamilton kernel of magnus.expmkernels for 2x2 and 3x3 matrices when numba is installed, and numpy.linalg.eigh for everything else. Never fails: without numba, or at dimension 4 and above, it is silently 'eigh'.
'numba': the same, except that a missing numba is an error rather than a fallback – for a caller who means to be sure the fast path is the one running. Dimensions 4 and 5 still use eigh even here, because there is no practical closed form for a 4x4 or 5x5 Hermitian eigenproblem; 4nu and 5nu stay correct and are simply not accelerated.
'eigh': numpy.linalg.eigh always, ignoring numba. The reference route, and what to set when comparing the two.

eigh costs about 1.25 us per 3x3 whatever the stack size, because it loops over LAPACK internally instead of vectorising, which makes it roughly a quarter of a 108-slab Magnus pass. The kernel removes that.

Setting this is the way to reach the whole package, including every magnus.oscprob wrapper; the expm_backend parameter on magnus_expansion(), evolution_operators_from_samples() and magnus_expansion_multislab() overrides it for one call.

That includes n_jobs != 1, but only because it is carried across deliberately: a module global does not survive a process boundary, and loky re-imports magnus in each worker with this back at its default. oscprob.osc_prob_energy_baseline reads the value in the parent and re-applies it inside the worker. Anything that adds a second parallel entry point has to do the same, or that path silently runs 'auto' whatever this says.

Added in version 1.0.0.

Type:: str

Compute \(\exp(\Omega_1 + \cdots + \Omega_\text{order})\) of \(A(t)\) from t0 to t1.

Added in version 1.0.0.

Parameters:

A (Callable) – Matrix function of time; must return a (d, d) NumPy array for a scalar time. If it also accepts an array of times (returning a (n, d, d) stack), the vectorized form is used automatically for speed; this is detected silently and verified against a scalar evaluation.
t0 (float) – Integration limits (t1 >= t0).
t1 (float) – Integration limits (t1 >= t0).
n_tpts (int, optional) – Number of uniformly spaced time points used to evaluate the integrals (‘trapezoid’/’simpson’ methods only; >= 2).
order (int, optional) – Highest Magnus order (1 to 6).
integration_method (str, optional) – ‘gl’ (Gauss-Legendre collocation; ignores n_tpts and uses 1, 2, or 3 nodes for orders <= 2, <= 4, <= 6, respectively), ‘trapezoid’, or ‘simpson’. Default: ‘gl’.
return_magnus_terms (bool, optional) – If True, also return the individual Magnus terms. For the ‘gl’ method the terms are not separable, and a single-element list containing the total \(\Omega\) is returned instead.
validate_input (bool, optional) – If True, validate order and integration_method (raises ValueError on invalid input).
A_eval_mode (str, optional) – Skip probing how A can be evaluated by declaring it up front (‘vector’, ‘scalar’, or ‘constant’); see probe_eval_mode(). If None (default), it is probed once and detected automatically.
expm_backend (str, optional) – Which routine exponentiates the slab: 'auto', 'numba' or 'eigh'. If None (default), the module-level EXPM_BACKEND decides.

Returns:

The evolution operator \(U = \exp(\sum_k \Omega_k)\); if return_magnus_terms is True, also the stacked terms.

Return type:

np.ndarray, or (np.ndarray, np.ndarray)

Examples

import numpy as np

from magnus import magnus

H = np.array([[0.0, 1.0], [1.0, 0.0]])
U = magnus.magnus_expansion(lambda t: -1j*H, 0.0, np.pi/4, order=4)

print(np.round(U, 6))
print('unitary to %.1e' % np.max(np.abs(U.conj().T @ U - np.eye(2))))

[[0.707107+0.j       0.      -0.707107j]
 [0.      -0.707107j 0.707107+0.j      ]]
unitary to 3.3e-16

Unitary to rounding, and that is structural rather than lucky: the truncated series is anti-Hermitian at any order.

Evolution operators of a chain of slabs from precomputed samples.

Mid-level entry point for callers that build the samples of A themselves – e.g., to batch extra axes (such as the neutrino energy) in front of the slab axis, which this routine broadcasts through all operations.

Added in version 1.0.0.

Parameters:

At (np.ndarray) – Samples of A, shape (…, n_slabs, m, d, d). For the quadrature methods (‘trapezoid’/’simpson’), the m samples of each slab lie on the uniform grid spanning the slab (endpoints included). For ‘gl’, they lie on the Gauss-Legendre nodes (m = 1, 2, or 3 for orders <= 2, <= 4, <= 6; see gl_nodes()).
widths (list or np.ndarray) – Slab widths, shape (n_slabs,) (or broadcastable to the leading axes of At without the last three).
order (int, optional) – Highest Magnus order (1 to 6).
integration_method (str, optional) – ‘gl’, ‘trapezoid’, or ‘simpson’. Default: ‘gl’.
A_is_const (bool, optional) – Set to True if A is constant in time to skip the (inapplicable) slab-width convergence warning.
validate_input (bool, optional) – If True, validate order and integration_method.
expm_backend (str, optional) – Which routine exponentiates each slab: 'auto', 'numba' or 'eigh'. If None (default), the module-level EXPM_BACKEND decides.

Returns:

Evolution operators, shape (…, n_slabs, d, d).

Return type:

np.ndarray

magnus.magnus.gl_nodes(order: int) → numpy.ndarray[source]

Returns the Gauss-Legendre nodes on [0, 1] used by the ‘gl’ method.

Added in version 1.0.0.

Parameters:: order (int) – Requested Magnus order; mapped to the smallest GL scheme with at least that order (1-2 -> 1 node, 3-4 -> 2 nodes, 5-6 -> 3 nodes).
Returns:: GL nodes on [0, 1] (1, 2, or 3 of them).
Return type:: np.ndarray

Examples

import numpy as np

from magnus import magnus

for order in (2, 4, 6):
    print('order %d -> %s' % (order, np.round(magnus.gl_nodes(order), 6)))

order 2 -> [0.5]
order 4 -> [0.211325 0.788675]
order 6 -> [0.112702 0.5      0.887298]

One, two or three nodes: the scheme uses the fewest that reach the order.

magnus.magnus.USE_PALINDROME = True[source]

Module-level switch.

Whether a slab chain whose profile reads the same from either end may be built by evaluating \(A\) on its first half only, the mirrored half following by reversal. True by default: every Earth chord qualifies, because a chord through a spherically symmetric Earth meets every radius twice.

Set it to False to evaluate every slab in full. This is not a correctness switch, but neither is it a no-op: the two routes agree to a few times 1e-15 rather than bitwise, because the mirrored slab’s nodes are reached as (L - b) + h*s on one route and a + h*s on the other, which are different floating-point expressions for the same real number. Set it to False to ask for the plain per-slab evaluation when a comparison needs one.

The saving is halved evaluations of the caller’s Hamiltonian, so it is worth most where that Hamiltonian is expensive: with f the share of slab time spent inside it, the speed-up is about \(1/(1 - f/2)\).

Added in version 1.0.0.

Type:: bool

magnus.magnus.palindromic(*arrays: numpy.ndarray) → bool[source]

Returns whether every array given reads the same both ways.

Added in version 1.0.0.

The comparison is exact, deliberately. The saving relies on the mirrored slab’s inputs being identical to the reversal of its partner’s, which follows from identical inputs and from nothing weaker; a tolerance here would silently return a different answer for a nearly-symmetric profile, which is the one thing an optimisation must never do. It is the producer’s business to make a profile exactly symmetric rather than nearly so.

This mirrors fastkernels.palindromic in NuOscProbExact, deliberately, down to treating an empty call and any array shorter than two entries as trivially palindromic.

Parameters:: arrays (np.ndarray) – Arrays to test, given as separate arguments and each reversed along its first axis.
Returns:: Whether every array equals its own reverse exactly.
Return type:: bool

Examples

import numpy as np

from magnus import magnus

print(magnus.palindromic(np.array([1.0, 2.0, 1.0])))
print(magnus.palindromic(np.array([1.0, 2.0, 3.0])))

True
False

Compute the evolution operators of all time slabs at once.

Vectorized (batched) version of magnus_expansion() for a chain of time slabs: A is evaluated for all slabs in a single call (when it supports array input), and the quadrature, commutator algebra, and matrix exponentials are evaluated as batched NumPy operations with the slab axis leading. This is much faster than calling magnus_expansion() slab by slab.

Added in version 1.0.0.

Parameters:

A (Callable) – Matrix function of time (see magnus_expansion()).
t_slab_edges (list or np.ndarray) – Slab edges, shape (n_slabs, 2): [[t0, t1], [t1, t2], …]. Slabs of zero width yield identity operators.
n_tpts_per_slab (int, optional) – Number of time points per slab (‘trapezoid’/’simpson’ only).
order (int, optional) – Highest Magnus order (1 to 6).
integration_method (str, optional) – ‘gl’, ‘trapezoid’, or ‘simpson’. Default: ‘gl’.
validate_input (bool, optional) – If True, validate input (raises ValueError on invalid input).
A_eval_mode (str, optional) – Skip probing how A can be evaluated by declaring it up front (‘vector’, ‘scalar’, or ‘constant’); see probe_eval_mode(). If None (default), it is probed once and detected automatically.
symmetric_over (tuple, optional) –
Caller’s declaration that A(t) == A(lo + hi - t) on (lo, hi). When given, and when the slab chain is found to span exactly that interval with exactly palindromic widths, A is evaluated on the first half of the slabs only and the rest follows by reversal – halving the calls to the caller’s Hamiltonian. Ignored when USE_PALINDROME is False.

This is a declaration, not a test: it is not checked, and cannot be cheaply, since testing it would require the evaluations it exists to avoid. Declaring it of a profile that is not symmetric returns a silently wrong answer – measured at 3.3e-01 on a monotonic profile. It is therefore not a user-facing knob: it is set by the Earth entry points, where the symmetry is a fact of chord geometry rather than a claim. See docs/dev/PLAN_PALINDROMIC_PROFILES.md section 3d(ii).
expm_backend (str, optional) – Which routine exponentiates each slab: 'auto', 'numba' or 'eigh'. If None (default), the module-level EXPM_BACKEND decides.

Returns:

Stack of evolution operators, shape (n_slabs, d, d), ordered like t_slab_edges (i.e., earliest slab first).

Return type:

np.ndarray

Notes

The time-ordered product over the chain is U_total = U[n_slabs-1] @ ... @ U[1] @ U[0], i.e., the last slab is the leftmost factor.