Retrieval

The retrieve() entry point, the shared solver base, the ten algorithms (croak.retrieve.ALGORITHMS), the result object, the live progress reporter and the trace-error metrics. See Retrieval algorithms for how each solver works and Choosing a solver for when to use which.

The retrieve wrapper

The retrieve() convenience wrapper and the algorithm registry.

croak.retrieve.retrieve(trace, omega, delays, interaction, *, algorithm='warm-lbfgs', guess=None, weights=None, rng=None, omega0=None, callback=None, progress='auto', **kwargs)

Retrieve a pulse from a measured trace with the chosen algorithm.

A thin convenience wrapper around the solver classes: it constructs the solver named by algorithm from kwargs and calls its run(). See Choosing a solver for the supported combinations.

Parameters:

• trace (array_like) – Measured trace intensity (Nomega, Ndelay).

• omega (array_like) – Centred angular-frequency axis (rad/s).

• delays (array_like) – Delay axis (s).

• interaction (str) – "shg", "sd" or "pg".

• algorithm (str, optional) – "warm-lbfgs" (the default: lbfgs-ad warm-started from COPRA’s local sweep), "copra", "copra-jax", "lbfgs", "lbfgs-hand", "lbfgs-ad", "lbfgs-optx", "lm", "lm-optx" or "cma-es" — the keys of ALGORITHMS.

• guess (None, array_like or Pulse, optional) – Initial guess.

• weights (array_like, optional) – Per-frequency residual weights (length Nomega), squared in the least-squares loss. Use inverse standard deviations for independent errors with known variance.

• rng (numpy.random.Generator, optional) – Random generator for a random initial guess.

• omega0 (float, optional) – Carrier angular frequency (rad/s). For a dispersive retrieval it sets the carrier at which the slab dispersion is built (forwarded to the solver constructor) as well as labelling the result’s wavelength axes; for a thin local model it only labels the axes. The chromatic focal and collection models also need the correct carrier.

• callback (callable, optional) – callback(iteration, R, best_R, snapshot=...) progress hook, where snapshot is an optional zero-argument in-progress-result builder (or None) for a live full-plot preview. When supplied (e.g. by the GUI) the caller is assumed to render its own progress, so the built-in progress bar/report is suppressed.

• progress ({"auto", True, False}, optional) – Built-in progress reporting when no callback is given. "auto" (default) shows a live bar and a final summary report (error, retrieved and transform-limited FWHM, iterations, time) in a terminal or notebook, and stays silent when output is non-interactive (piped/CI/tests). True forces the report on; False disables it.

• **kwargs – Solver-specific options forwarded to the algorithm constructor (e.g. maxiters, alpha for COPRA; material, thickness, npoints for dispersive propagation — paired with omega0 above).

Returns:

RetrievalResult

Return type:

RetrievalResult

Notes

The penalty-capable solvers (every family except COPRA) default to the default regularisation weights — reg_spectrum = 0.01 and reg_amp = 0.03 for the gradient objective, mapped to 1e-5/3e-5 for the Levenberg–Marquardt residual form (croak.solver.REG_DEFAULTS). The spectral term binds only when a spectrum_target is supplied; pass reg_amp=0, reg_spectrum=0 for an unregularised fit.

croak.retrieve.ALGORITHMS: dict[str, type[Retriever]] = {'cma-es': <class 'croak.cmaes.CMAES'>, 'copra': <class 'croak.copra.COPRA'>, 'copra-jax': <class 'croak.copra_jax.COPRAJax'>, 'lbfgs': <class 'croak.lbfgs.LBFGS'>, 'lbfgs-ad': <class 'croak.lbfgs_ad.LBFGSAD'>, 'lbfgs-hand': <class 'croak.lbfgs_hand.LBFGSHand'>, 'lbfgs-optx': <class 'croak.optimistix_lbfgs.OptxLBFGS'>, 'lm': <class 'croak.lm.LM'>, 'lm-optx': <class 'croak.optimistix_lm.OptxLM'>, 'warm-lbfgs': <class 'croak.warm_lbfgs.WarmLBFGS'>}

Registry mapping algorithm names to solver classes. Dispersive propagation is a parameter of the solvers (material/thickness), not a separate algorithm. "lbfgs" uses hand-derived analytic gradients (pure NumPy); "lbfgs-hand" uses the same hand-derived gradients but jitted in JAX (vmap over delays, no AD); "lbfgs-ad" and "lm" use JAX automatic differentiation (gradient / Jacobian). The "-optx" variants are JAX-native twins backed by Optimistix (the whole optimisation loop runs on-device). "copra-jax" is the JAX twin of "copra" (vmap/jit global step, lax.scan local sweep) — the same algorithm, vectorised over delays.

croak.retrieve.algorithm_params(name)

Return the keyword names accepted by a registered solver constructor.

The solver’s __init__ signature is the single source of truth for what each algorithm supports ("copra"/"copra-jax" lack reg_amp etc., so those are absent from the returned set); the set includes inherited options for wrappers such as warm-lbfgs. The pipeline dispatch and the GUI’s option greying use it to forward only accepted options and to disable controls an algorithm cannot use. Unknown names raise ValueError.

Parameters:

name (str)

Return type:

set[str]

Solver base

Shared retrieval-solver base class and helpers.

Defines Retriever, the common base for every solver registered in croak.retrieve.ALGORITHMS. It standardises the run entry point — trace normalisation, grid construction, initial-guess resolution and weight handling — so every algorithm shares an identical calling convention and return type (RetrievalResult).

class croak.solver.Retriever(*, maxiters=100, verbose=False)

Bases: object

Base class for retrieval algorithms.

Subclasses implement _solve. Algorithm and forward-model hyperparameters are set on the constructor; run() takes the measurement and returns a RetrievalResult.

Parameters:

• maxiters (int)

• verbose (bool)

name: str = 'retriever'

Algorithm name recorded in the result.

run(trace, omega, delays, interaction, *, guess=None, weights=None, rng=None, omega0=None, callback=None)

Retrieve a pulse from a measured FROG trace.

Parameters:

• trace (array_like) – Measured trace intensity, shape (Nomega, Ndelay). Copied and normalised internally.

• omega (array_like) – Centred angular-frequency axis (rad/s).

• delays (array_like) – Delay axis (s).

• interaction (str) – "shg", "sd" or "pg".

• guess (None, array_like or Pulse, optional) – Initial guess (see resolve_guess()).

• weights (array_like, optional) – Per-frequency weights (length Nomega); uniform if omitted.

• rng (numpy.random.Generator, optional) – Random generator for a random initial guess.

• omega0 (float, optional) – Carrier angular frequency (rad/s) recorded on the result for wavelength axes; defaults to the solver’s omega0 (dispersive solvers) or 0.

• callback (callable, optional) – Called as callback(iteration, R, best_R, snapshot=...) after each iteration (or objective evaluation), e.g. for live GUI progress. snapshot is an optional zero-argument builder of the current in-progress RetrievalResult (for a live full-plot preview); solvers that cannot cheaply build one pass None. A callback may ignore it.

Returns:

RetrievalResult

Return type:

RetrievalResult

croak.solver.resolve_guess(grid, guess, *, rng=None)

Turn a flexible guess argument into a complex spectrum.

Parameters:

• grid (Grid) – Time/frequency grid.

• guess (None, array_like or Pulse) – None draws a random Gaussian-spectrum initial guess; an array is used directly as the spectrum; a Pulse contributes its spectrum.

• rng (numpy.random.Generator, optional) – Random generator for the None case.

Returns:

numpy.ndarray – A complex spectrum of length grid.n (centred order).

Return type:

NDArray[complex128]

croak.solver.build_weights(n, weights)

Return a length-n per-frequency weight vector (uniform if None).

Parameters:

• n (int)

• weights (ArrayLike | None)

Return type:

NDArray[float64]

croak.solver.spectral_target_amplitude(intensity)

Peak-normalised spectral amplitude from a measured spectral intensity.

Shared by the gradient/LM solvers’ spectral-match regularisation: converts a measured spectral intensity \(I(\omega)\) (centred order, on the retrieval grid) to the peak-normalised amplitude \(\sqrt{I}\) compared against the retrieved \(|\tilde E(\omega)|\). Returns None when no spectrum is supplied.

Parameters:

intensity (ArrayLike | None)

Return type:

NDArray[float64] | None

croak.solver.assemble_result(*, spectrum, t_sim, grid, t_meas, delays, weights, interaction, algorithm, errors, R_omega=False, thickness=None, tau0=0.0, smear_scale=None, smear_scale_delta=None, omega0=0.0)

Assemble a RetrievalResult from a candidate solution.

Every retriever, at convergence, takes its retrieved spectrum and the corresponding unscaled simulated trace t_sim and computes the optimal intensity scale \(\mu\), the scaled trace and the FROG error \(R\) identically. This helper is that shared assembly, so the final result and any in-progress snapshot (for the GUI’s live preview) are built the same way — a single source of truth for the mu/error maths.

Parameters:

• spectrum (numpy.ndarray) – Retrieved complex spectrum \(\tilde E(\omega)\) (centred order).

• t_sim (numpy.ndarray) – Unscaled simulated trace \(T_\mathrm{sim}\), shape (Nomega, Ndelay) (i.e. before the optimal \(\mu\) scaling is applied).

• grid (Grid) – Retrieval time/frequency grid.

• t_meas (numpy.ndarray) – Normalised measured trace, same shape as t_sim.

• delays (numpy.ndarray) – Delay axis (s).

• weights (numpy.ndarray) – Per-frequency weights (length Nomega).

• interaction (str) – "shg"/"sd"/"pg" (normalised via get_interaction()).

• algorithm (str) – Algorithm name recorded on the result.

• errors (sequence of float) – FROG error per iteration so far (copied onto the result).

• R_omega (bool, optional) – Per-frequency adaptive scaling (a \(\mu\) vector) instead of a single scalar.

• thickness (float or None, optional) – Fitted slab thickness (m), or None when not fitted.

• tau0 (float, optional) – Fitted delay-zero offset (s); 0.0 when not fitted.

• smear_scale (float or None, optional) – Fitted geometric-smearing multiplier, or None when not fitted. For a split fit this is the gate-shape (p) channel.

• smear_scale_delta (float or None, optional) – Fitted delay-channel multiplier for a split smearing fit (fit_smearing_split), or None for a joint fit.

• omega0 (float, optional) – Carrier angular frequency (rad/s) recorded for the wavelength axes.

Returns:

RetrievalResult

Return type:

RetrievalResult

croak.solver.adapt_callback(callback)

Return callback wrapped so it always accepts a snapshot= keyword.

The documented progress signature is callback(iteration, R, best_R, snapshot=None), but the shorter callback(iteration, R, best_R) form is also supported. Every solver can now supply an in-progress snapshot, so calling a three-argument callback with the keyword would fail with a TypeError deep inside the iteration. The adapter drops the snapshot for callbacks that cannot receive it, so the same callback works with every solver.

Parameters:

callback (callable or None) – A user progress callback, or None.

Returns:

callable or None – A callable accepting snapshot=, or None if none was given.

croak.solver.REG_DEFAULTS: dict[str, tuple[float, float]] = {'gradient': (0.01, 0.03), 'lm': (1e-05, 3e-05)}

Default regularisation weights (reg_spectrum, reg_amp) per solver family.

The two families do not minimise the same objective, so the same number does not mean the same thing to them. The gradient solvers minimise \(R + \lambda P\); the Levenberg-Marquardt solvers append the penalty as residual rows, hence minimise \(R^2 + \lambda' P\). Matching the stationary conditions at the same field pairs them as \(\lambda' = 2 R \lambda\), i.e. at the \(R \approx 2\times10^{-4}\) of a good dispersive TG-FROG retrieval the LM weight is some three orders of magnitude smaller.

The difference matters in practice. On a simulated 2 fs pulse through a 9.5 um dispersive TG-FROG medium, sweeping the LM weight over four decades gave:

LM (spectrum, amp)	R (%)	eps	dFWHM (%)	rough.
1e-2, 3e-2 (as LBFGS)	0.0237	0.305	+1.11	1.00
1e-5, 3e-5 (mapped)	0.0193	0.112	-0.11	1.01
0, 0	0.0174	0.466	+3.12	2.57

(eps is the complex-field error against the known pulse and rough. the spectral roughness relative to the truth.) Both failure modes are visible: the gradient weights over-regularise LM (a worse pulse and a worse trace error), and no penalty at all overfits (the best trace error of any arm but the worst pulse, with 2.6x the truth’s spectral roughness and +3.1 % on duration). The measured optimum, 1e-5, sits within 30 % of the \(2 R \lambda\) prediction, which supports the mapping as the rule rather than a coincidence of one trace.

The LM entry therefore scales with the trace error actually reached: multiply it by \(R/2\times10^{-4}\) on a noisier measurement. See Regularisation.

COPRA

COPRA step-based retrieval.

Implements the Common Pulse Retrieval Algorithm (COPRA), reformulated so that the frequency-domain signal \(\psi(\omega)\) produced by ForwardModel is the working quantity. The single reverse-mode adjoint kernel of the forward model provides every gradient, so the same engine serves both the thin-medium and the dispersive case (the only difference is the forward model the engine is handed).

COPRA alternates two modes:

• local — sweep the delays in random order, taking a per-delay gradient step on \(\tilde E\) towards the amplitude-projected signal;

• global — once the local mode stalls, take one joint gradient-descent step in the signal and then in \(\tilde E\) over all delays at once.

class croak.copra.COPRA(*, maxiters=100, alpha=0.25, reltol=0.0, abstol=0.0, material=None, thickness=0.0, npoints=20, omega0=0.0, quadrature='gausslegendre', phase_only=False, R_omega=False, verbose=False)

Bases: Retriever

COPRA step-based retrieval solver (optionally with dispersive propagation).

The same engine handles the thin medium and the dispersive case — set material/thickness to propagate the field through a slab (depth integral over npoints quadrature nodes). With material=None (the default) it is the plain thin-medium algorithm.

Parameters:

• maxiters (int, optional) – Maximum iterations (default 100).

• alpha (float, optional) – Global-mode step-size scaling (default 0.25).

• reltol (float, optional) – Relative-improvement convergence tolerance (default 0 = run all maxiters); see run_copra().

• abstol (float, optional) – Absolute target FROG error (default 0 = off); see run_copra().

• material (str or None, optional) – Slab material for dispersive propagation (e.g. "SiO2"); None for a thin medium.

• thickness (float, optional) – Slab thickness (m).

• npoints (int, optional) – Number of depth quadrature points (default 20).

• omega0 (float, optional) – Carrier angular frequency (rad/s) used to build the dispersion.

• quadrature (str, optional) – Quadrature method ("gausslegendre" or "uniform").

• phase_only (bool, optional) – Retrieve only the spectral phase, holding the guess amplitude fixed.

• R_omega (bool, optional) – Use per-frequency adaptive scaling factors instead of a single global mu (the Rω mode).

• verbose (bool, optional) – Print per-iteration progress.

name: str = 'copra'

Algorithm name recorded in the result.

croak.copra.run_copra(model, grid, t_meas, ew0, weights, *, maxiters=100, alpha=0.25, reltol=0.0, abstol=0.0, verbose=False, rng=None, callback=None, fixed_amplitude=None, R_omega=False, snapshot_fn=None)

Run the COPRA iteration. Returns (spectrum, errors, trace, mu, best_R).

snapshot_fn, when given, is (ew, t_sim, errors) -> RetrievalResult: it turns the current iterate into an in-progress result, handed to callback as its snapshot= argument so the GUI can draw a live full plot. Both arguments are already computed each iteration for the error bookkeeping, so a snapshot costs no extra forward pass.

If fixed_amplitude is given, the spectral amplitude is held fixed to it after every update (phase-only retrieval): E ← fixed_amplitude · e^{iφ}.

If R_omega is set, the optimal scale mu is computed per frequency (one factor per spectral row) instead of as a single global value, matching the Rω adaptive-scaling mode of the gradient solvers. The returned mu is then a length-Nomega vector.

Parameters:

• model (ForwardModel) – Configured forward model (thin or dispersive).

• grid (Grid) – Time/frequency grid.

• t_meas (numpy.ndarray) – Measured trace (Nomega, Ndelay), normalised to unit maximum.

• ew0 (numpy.ndarray) – Initial complex spectrum (centred order).

• weights (numpy.ndarray) – Per-frequency weights (length Nomega).

• maxiters (int, optional) – Maximum number of iterations.

• alpha (float, optional) – Global-mode step-size scaling.

• reltol (float, optional) – Relative-improvement convergence tolerance (0 disables). Iteration stops once the best FROG error has improved by less than reltol (relative) over the last _COPRA_PATIENCE iterations — i.e. the run has plateaued (including after the global-mode switch).

• abstol (float, optional) – Absolute target FROG error (0 disables): stop once the best error drops to or below abstol.

• verbose (bool, optional) – Print per-iteration FROG error.

• rng (Generator | None)

• fixed_amplitude (NDArray[float64] | None)

• R_omega (bool)

Return type:

tuple[NDArray[complex128], list[float], NDArray[float64], NDArray[float64] | float, float]

COPRA (JAX)

The copra-jax twin: identical algorithm, per-iteration numerics vectorised with vmap/jit and lax.scan.

JAX-accelerated COPRA (copra-jax).

A drop-in twin of croak.copra.COPRA that runs the per-iteration work on the JAX forward/adjoint kernel (croak.forward_jax.make_forward_adjoint_fns()) instead of the per-delay NumPy loops. The update rules follow croak.copra.run_copra() — same two modes, same step sizes, same local→global switch and convergence tests — only the two stage bodies are vectorised:

• the global step does one vmap forward (keeping the (test, gate) tape) and one vmap adjoint over all delays inside a single jax.jit;

• the local sweep, which is inherently sequential (each delay updates \(\tilde E\) for the next), is expressed as a jax.lax.scan over the delay-index permutation. The permutation is drawn host-side with the same NumPy rng as the reference, so the two implementations track each other (up to floating-point differences in the vectorised arithmetic).

The outer iteration loop (mode switching, best-error tracking, callbacks, abstol/reltol stops) stays in Python, exactly as croak.lbfgs_hand.LBFGSHand keeps its loop in Python and delegates the numerics to jitted JAX.

class croak.copra_jax.COPRAJax(*, maxiters=100, alpha=0.25, reltol=0.0, abstol=0.0, material=None, thickness=0.0, npoints=20, omega0=0.0, quadrature='gausslegendre', phase_only=False, R_omega=False, verbose=False, smearing=None, stall_patience=5)

Bases: Retriever

JAX-accelerated COPRA retriever (copra-jax).

Same parameters and behaviour as croak.copra.COPRA (maxiters, alpha, reltol/abstol, the dispersive-slab material/thickness/npoints/omega0/quadrature options, phase_only and R_omega); see that class for the full description. The per-iteration numerics run on the JAX forward/adjoint kernel (vmap/jit global step, lax.scan local sweep). See run_copra_jax().

Notes

Like croak.copra.COPRA, it takes none of the gradient-solver regularisation penalties (reg_amp/reg_phase/reg_spectrum/ reg_time) or the B-spline phase basis.

Parameters:

• maxiters (int)

• alpha (float)

• reltol (float)

• abstol (float)

• material (str | None)

• thickness (float)

• npoints (int)

• omega0 (float)

• quadrature (str)

• phase_only (bool)

• R_omega (bool)

• verbose (bool)

• stall_patience (int)

name: str = 'copra-jax'

Algorithm name recorded in the result.

stall_patience

Iterations without improvement before the local projection sweep hands over to the global descent. Sets how much of a run is projection and how much is gradient descent, which is the whole difference between the two stages – and, with a smeared model, the difference between the stage that can represent the physics and the one that cannot. Geib’s value is 5; a very large value never leaves the local sweep.

croak.copra_jax.run_copra_jax(signal_tape, adjoint, delays, t_meas, ew0, weights, *, maxiters=100, alpha=0.25, reltol=0.0, abstol=0.0, verbose=False, rng=None, callback=None, fixed_amplitude=None, R_omega=False, smeared_trace_fn=None, snapshot_fn=None, stall_patience=5)

JAX clone of croak.copra.run_copra().

signal_tape / adjoint are the pair from make_forward_adjoint_fns(). Returns (spectrum, errors, trace, mu, best_R) as numpy, matching run_copra.

snapshot_fn, when given, is (ew, t_sim, errors) -> RetrievalResult: it turns the current iterate into an in-progress result, which is handed to callback as its snapshot= argument so the GUI can draw a live full plot. The sweep already computes both arguments each iteration for its own error bookkeeping, so a snapshot costs no extra forward pass.

Parameters:

• delays (NDArray[float64])

• t_meas (NDArray[float64])

• ew0 (NDArray[complex128])

• weights (NDArray[float64])

• maxiters (int)

• alpha (float)

• reltol (float)

• abstol (float)

• verbose (bool)

• rng (Generator | None)

• fixed_amplitude (NDArray[float64] | None)

• R_omega (bool)

• stall_patience (int)

Return type:

tuple[NDArray[complex128], list[float], NDArray[float64], NDArray[float64] | float, float]

L-BFGS (analytic gradients)

L-BFGS retrieval via NLopt with hand-derived gradients.

Minimises the FROG error \(R(u)\) over a pulse parameterisation u (ArrayPulse) using NLopt’s LD_LBFGS. The gradient is analytic: the FROG-error cotangent on the trace is propagated through adjoint_single() and then through the pulse’s ew_vjp() — no automatic differentiation.

Optional features: per-frequency adaptive scaling (R_omega), phase-only retrieval, and second-difference smoothness regularisation of the spectral amplitude and phase.

class croak.lbfgs.LBFGS(*, maxiters=300, material=None, thickness=0.0, npoints=20, omega0=0.0, quadrature='gausslegendre', phase_only=False, R_omega=False, reg_amp=None, reg_phase=0.0, reg_spectrum=None, spectrum_target=None, reltol=0.0001, abstol=1e-08, targeterr=0.0, verbose=False)

Bases: Retriever

NLopt L-BFGS solver with analytic gradients.

Parameters:

• maxiters (int, optional) – Maximum number of objective evaluations (default 300).

• material (str or None, optional) – Dispersive slab material; None (default) is a thin medium.

• thickness (float, optional) – Slab thickness (m).

• npoints (int, optional) – Depth quadrature points (default 20).

• omega0 (float, optional) – Carrier angular frequency (rad/s) for the dispersion.

• quadrature (str, optional) – Quadrature method.

• phase_only (bool, optional) – Retrieve only the spectral phase, holding the guess amplitude fixed.

• R_omega (bool, optional) – Use per-frequency adaptive scaling factors.

• reg_amp (float or None, optional) – Amplitude second-difference smoothness weight. None (default) resolves to the package default weight 0.03 (croak.solver.REG_DEFAULTS); pass 0 to disable.

• reg_phase (float, optional) – Phase second-difference smoothness weight (default 0).

• reg_spectrum (float or None, optional) – Weight of the spectral-match penalty pulling the retrieved spectral amplitude towards spectrum_target (full mode only). None (default) resolves to the default 0.01 and is inert without a spectrum_target; pass 0 to disable.

• spectrum_target (array_like or None, optional) – Measured spectral intensity on the retrieval grid (centred order, e.g. TraceData.Iomega) used by reg_spectrum; None disables the penalty.

• reltol (float, optional) – NLopt relative/absolute function tolerances (defaults 1e-4 / 1e-8).

• abstol (float, optional) – NLopt relative/absolute function tolerances (defaults 1e-4 / 1e-8).

• targeterr (float, optional) – Stop early when the FROG error drops below this value. Two robustness guards apply during optimisation: a non-finite objective or gradient is replaced by a large finite penalty so the line search backtracks (guard_nonfinite), and an internal NLopt failure (line-search breakdown, roundoff limit) returns the best iterate found instead of raising; both warn once.

• verbose (bool, optional) – Print per-evaluation progress.

name: str = 'lbfgs'

Algorithm name recorded in the result.

L-BFGS (hand-gradient JAX)

The lbfgs-hand twin: the same analytic Wirtinger gradient as LBFGS, but jitted and vmap-ed in JAX (no autodiff).

Retrieve with analytic Wirtinger gradients evaluated in JAX (lbfgs-hand).

A twin of croak.lbfgs.LBFGS: the same pulse parameterisation, NLopt LD_LBFGS optimiser and options (phase-only retrieval, R_omega adaptive scaling, amplitude/phase/spectral penalties, dispersive or thin propagation). It does not support the additional spatial and fitted-parameter options of croak.lbfgs_ad.LBFGSAD.

The gradient is the same analytic Wirtinger adjoint as croak.lbfgs.LBFGS (adjoint_single → ew_vjp), but expressed entirely in JAX and evaluated by a single jax.jit()-compiled function that vmaps the forward and adjoint over all delays; no automatic differentiation is used. This gives the hand-derived gradients the same vectorisation/JIT speedup that makes lbfgs-ad faster than the pure-NumPy lbfgs (the gap is the trace evaluation, not the gradient maths). It is thus a fast, AD-free reference path: an independent derivative route against which the AD gradient is checked in tests/test_ad_vs_analytic.py.

class croak.lbfgs_hand.LBFGSHand(*, maxiters=300, material=None, thickness=0.0, npoints=20, omega0=0.0, quadrature='gausslegendre', phase_only=False, R_omega=False, reg_amp=None, reg_phase=0.0, reg_spectrum=None, spectrum_target=None, reltol=0.0001, abstol=1e-08, targeterr=0.0, verbose=False)

Bases: Retriever

NLopt L-BFGS solver with jitted JAX hand-derived (non-AD) gradients.

Parameters mirror croak.lbfgs.LBFGS exactly (maxiters, material/thickness/npoints/omega0/quadrature, phase_only, R_omega, reg_amp/reg_phase/reg_spectrum + spectrum_target, reltol/abstol/targeterr, verbose); see that class for the full description. The forward model, FROG error and its gradient are all evaluated with the JAX twins, the gradient using the explicit Wirtinger adjoint (not jax.grad()).

Notes

Like croak.lbfgs.LBFGS — and unlike croak.lbfgs_ad.LBFGSAD — it does not expose the B-spline phase basis (phase_basis/n_nodes) or the temporal penalty (reg_time/time_window).

Parameters:

• maxiters (int)

• material (str | None)

• thickness (float)

• npoints (int)

• omega0 (float)

• quadrature (str)

• phase_only (bool)

• R_omega (bool)

• reg_amp (float | None)

• reg_phase (float)

• reg_spectrum (float | None)

• reltol (float)

• abstol (float)

• targeterr (float)

• verbose (bool)

name: str = 'lbfgs-hand'

Algorithm name recorded in the result.

L-BFGS (autodiff)

L-BFGS retrieval with automatic-differentiation gradients (lbfgs-ad).

The most feature-complete solver in the package. It shares the NLopt LD_LBFGS optimiser and the amplitude/phase parameterisation with croak.lbfgs.LBFGS, but takes its gradient from jax.value_and_grad() of the JAX forward model and error metric (croak.forward_jax, croak.metrics_jax) rather than from the hand-derived adjoint (croak.forward.ForwardModel.adjoint_single() + croak.pulses.Pulse.ew_vjp()).

Differentiating the model automatically is what makes the extras possible, and they are only available here and on the other autodiff solvers:

• smearing — the geometric time-smearing kernel (croak.smearing),

• fit_thickness / fit_tau0 / fit_smearing — extra fitted scalars alongside the field,

• phase_basis="bspline" — the low-dimensional spline phase basis,

• reg_time / time_window — the temporal energy penalty.

Everything the plain lbfgs supports (phase-only, R_omega adaptive scaling, amplitude/phase/spectral regularisation, thin or dispersive propagation) works here too, so the two can also be run head-to-head on the same trace. The gradients themselves are compared directly in tests/test_ad_vs_analytic.py.

class croak.lbfgs_ad.LBFGSAD(*, maxiters=300, material=None, thickness=0.0, npoints=20, omega0=0.0, quadrature='gausslegendre', smearing=None, focal=None, depth_weight=None, depth_transverse=False, phase_only=False, phase_basis='pointwise', n_nodes=20, R_omega=False, reg_amp=None, reg_phase=0.0, reg_spectrum=None, spectrum_target=None, reg_time=0.0, time_window=None, tau0=0.0, smear_scale=1.0, smear_scale_delta=1.0, fit_thickness=False, fit_tau0=False, fit_smearing=False, fit_smearing_split=False, polish=False, delay_origin='coincidence', tau0_bound=0.0, reltol=0.0001, abstol=1e-08, targeterr=0.0, verbose=False)

Bases: Retriever

NLopt L-BFGS solver with JAX automatic-differentiation gradients.

The forward model and FROG error are evaluated with the JAX twins (croak.forward_jax, croak.metrics_jax) and the gradient comes from jax.value_and_grad(). That makes this the most capable solver in the package: geometric smearing, the fitted extras (fit_thickness, fit_tau0, fit_smearing), the B-spline phase basis and the temporal penalty are all reachable here and not from the analytic-gradient solvers.

The optimiser, parameterisation and objective otherwise match croak.lbfgs.LBFGS, so the two solvers can be compared directly on the same trace.

Parameters:

• maxiters (int, optional) – Maximum number of objective evaluations (default 300).

• material (str | None) – Dispersive-slab parameters, identical to croak.lbfgs.LBFGS; material=None (default) is a thin medium. See croak.forward.ForwardModel.

• thickness (float) – Dispersive-slab parameters, identical to croak.lbfgs.LBFGS; material=None (default) is a thin medium. See croak.forward.ForwardModel.

• npoints (int) – Dispersive-slab parameters, identical to croak.lbfgs.LBFGS; material=None (default) is a thin medium. See croak.forward.ForwardModel.

• omega0 (float) – Dispersive-slab parameters, identical to croak.lbfgs.LBFGS; material=None (default) is a thin medium. See croak.forward.ForwardModel.

• quadrature (str) – Dispersive-slab parameters, identical to croak.lbfgs.LBFGS; material=None (default) is a thin medium. See croak.forward.ForwardModel.

• smearing (SmearingKernel or None, optional) – Geometric time-smearing kernel of a non-collinear BOXCARS geometry (default None = off). PG (TG) and SD only. See croak.smearing and Modelling geometric time smearing.

• focal (FocalMixture or None, optional) –

  Chromatic focal mixture (croak.focal) in place of smearing: the focal average is carried explicitly, each node applying its own wavelength-dependent amplitude filter, where the reduced kernel evaluates a single \(|A|^6\) weight at the carrier. Mutually exclusive with smearing, and incompatible with fit_smearing (an explicit mixture has no kernel width to scale).

  Worth its cost for a broadband pulse. The reduced kernel’s widths scale as \(\sigma\propto\lambda\), so freezing them at the carrier is a -16 %/+24 % spread of width across a 1 fs deep-UV band; in our earlier instrument scans that cost a mean 12 % of retrieved duration there, against 3.5 % at 2 fs. The mixture needs one field build per focal node rather than per Gauss-Hermite node, so it is one to two orders of magnitude slower; check quadrature convergence for the chosen pulse and geometry.

• depth_weight (array_like or None, optional) – Optional longitudinal weighting for crossing beams – focal-spot evolution, beamlet walk-through and Gouy phase – as a smooth per-depth factor on the generated signal; needed only when the slab is not much thinner than the beams’ Rayleigh range. See croak.forward_jax.make_param_trace_fn().

• depth_transverse (bool, optional) – Model the transverse decoherence of the depth integral (default False): each depth node’s signal carries its own transverse phase \([k_z(\omega,k_\perp)-k_z(\omega,0)](L-z_q)\) at the collection aperture’s wavevectors before the coherent depth sum. Parameter-free. Requires focal with a collection and a dispersive slab. See croak.forward_jax.make_param_trace_fn() and Modelling the collection aperture.

• phase_only (bool, optional) – Retrieve only the spectral phase, holding the guess amplitude fixed.

• phase_basis ({"pointwise", "bspline"}, optional) – Phase parameterisation. "pointwise" (default) optimises one phase value per frequency bin. "bspline" optimises the n_nodes control points of a cubic B-spline over the spectral support, a strong low-dimensional regulariser that implies phase-only retrieval (the amplitude is held at the guess). See croak._jax_pulse.make_spline_phase_basis().

• n_nodes (int, optional) – Number of B-spline control points when phase_basis="bspline" (default 20); ignored otherwise.

• R_omega (bool, optional) – Use per-frequency adaptive scaling factors instead of a single scalar.

• reg_amp (float or None, optional) – Amplitude second-difference smoothness weight. None (default) resolves to the package default weight 0.03 (croak.solver.REG_DEFAULTS); pass 0 to disable.

• reg_phase (float, optional) – Phase second-difference smoothness weight (default 0). See Regularisation.

• reg_spectrum (float or None, optional) – Weight of the spectral-match penalty pulling the retrieved spectral amplitude towards spectrum_target (full mode only). None (default) resolves to the default 0.01 and is inert without a spectrum_target; pass 0 to disable.

• spectrum_target (array_like or None, optional) – Measured spectral intensity on the retrieval grid (centred order, e.g. TraceData.Iomega) used by reg_spectrum; None disables it.

• reg_time (float, optional) – Weight of a temporal penalty on the fraction of pulse energy falling outside time_window (0 = off). Suppresses spurious satellite pulses and pedestals far from the main feature. Not available in croak.lbfgs.LBFGS.

• time_window (tuple of float or None, optional) – (t_lo, t_hi) in seconds bounding the region kept free of penalty for reg_time. None (default) uses (delays.min(), delays.max()). The window edges are Planck-tapered. See croak._jax_pulse.time_window_mask().

• reltol (float, optional) – NLopt relative / absolute function tolerances (defaults 1e-4 / 1e-8).

• abstol (float, optional) – NLopt relative / absolute function tolerances (defaults 1e-4 / 1e-8).

• tau0 (float, optional) – Starting delay-zero offset \(\tau_0\) in seconds (default 0) — the centre fit_tau0 optimises around, e.g. the previous fit’s value when warm-starting. Held at this value when fit_tau0 is off.

• smear_scale (float, optional) – Starting multiplier on the smearing-kernel widths (default 1.0 = the kernel as supplied) — the centre fit_smearing optimises around. Held at this value when fit_smearing is off.

• fit_thickness (bool, optional) – Also fit the dispersive-slab thickness as an extra parameter (default False). Requires a dispersive slab (material set, thickness > 0) — only meaningful for PG/SD, where the depth-resolved interaction makes the trace depend on thickness beyond a pure input spectral phase. The prior thickness is the starting value. See Fitting the medium thickness and delay-zero.

• fit_tau0 (bool, optional) – Also fit a delay-zero offset \(\tau_0\) (default False), correcting a mis-set experimental time zero. Works for every interaction; enters the forward model as exp(iω(τ − τ₀)).

• fit_smearing (bool, optional) – Also fit the geometric-smearing strength as an extra parameter (default False): a dimensionless multiplier on both kernel widths, equivalently the mask ratio d/D. Requires smearing. Do not enable it together with fit_thickness — the smearing width is a nuisance parameter that can absorb dispersion-model error. See Modelling geometric time smearing.

• fit_smearing_split (bool, optional) – Fit the gate-shape (p) and delay (delta) kernel widths as two independent multipliers instead of one joint one (default False) — a diagnostic for which channel deviates from the geometric prediction. Requires fit_smearing; the fitted values land in smear_scale (p) and smear_scale_delta (delta) on the result.

• smear_scale_delta (float, optional) – Starting delay-channel multiplier for a split fit (default 1.0); used only when fit_smearing_split is set.

• polish (bool, optional) – Two-phase polish (default False): first retrieve the pulse with the extra parameters held fixed, then free them for a joint final phase. Only meaningful when fit_thickness or fit_tau0 is set.

• delay_origin ({"coincidence", "marginal_peak"}, optional) – Where the model’s delay zero sits. "coincidence" (default): gate–probe coincidence, the forward model’s native convention. "marginal_peak": every model trace is re-centred on its own delay-marginal peak by the same sub-sample rule the preprocessing applies to the data, so the two sides share one convention with no free parameter (croak.delay_origin). Use it when the data were centred by their marginal peak, as for coherently collected single-cycle traces, where the marginal peak sits tens of attoseconds from coincidence and a fixed origin cost ~10 % of duration in our earlier instrument scans; unlike fit_tau0 it stays put on chirped pulses. Requires uniform delays; the row-gain weighting of the centring is differentiated through.

• tau0_bound (float, optional) – Box bound on the fitted delay offset, |tau0| <= tau0_bound (s); 0 (default) leaves it unbounded. Only meaningful with fit_tau0.

• targeterr (float, optional) – Stop early when the FROG error drops below this value. Two robustness guards apply during optimisation: a non-finite objective or gradient is replaced by a large finite penalty so the line search backtracks (guard_nonfinite), and an internal NLopt failure (line-search breakdown, roundoff limit) returns the best iterate found instead of raising; both warn once.

• verbose (bool, optional) – Print per-evaluation progress.

name: str = 'lbfgs-ad'

Algorithm name recorded in the result.

L-BFGS (Optimistix)

The lbfgs-optx twin: the same JAX objective as LBFGSAD, driven by Optimistix’s on-device L-BFGS instead of NLopt.

Retrieve a pulse with Optimistix L-BFGS and JAX gradients (lbfgs-optx).

A JAX-native twin of croak.lbfgs_ad.LBFGSAD: it minimises the same scalar objective (FROG error plus the same amplitude/phase/spectral penalties) over the same real parameter vector (croak._jax_pulse), with the same JAX forward model and autodiff gradient. The difference is the optimiser: instead of NLopt’s LD_LBFGS driving a host-side loop (with the gradient transferred device→host every evaluation), Optimistix’s L-BFGS runs the limited-memory BFGS loop on-device, each step jitted. It does not expose the temporal penalty or the jointly fitted scalar parameters of the NLopt-backed solver.

It is the counterpart to OptxLM for head-to-head comparison against the NLopt-backed lbfgs-ad.

class croak.optimistix_lbfgs.OptxLBFGS(*, maxiters=300, material=None, thickness=0.0, npoints=20, omega0=0.0, quadrature='gausslegendre', smearing=None, focal=None, phase_only=False, phase_basis='pointwise', n_nodes=20, R_omega=False, reg_amp=None, reg_phase=0.0, reg_spectrum=None, spectrum_target=None, reltol=0.0001, abstol=1e-08, targeterr=0.0, verbose=False)

Bases: Retriever

Optimistix limited-memory BFGS solver with JAX autodiff gradients.

The objective and parameterisation match croak.lbfgs_ad.LBFGSAD; the optimiser is Optimistix’s on-device L-BFGS in place of NLopt’s host-side LD_LBFGS. The temporal penalty and the jointly fitted scalar parameters of LBFGSAD are not available.

Parameters:

• maxiters (int, optional) – Maximum number of solver steps (default 300).

• material (str | None) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• thickness (float) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• npoints (int) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• omega0 (float) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• quadrature (str) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• smearing (SmearingKernel or None, optional) – Geometric time-smearing kernel of a non-collinear BOXCARS geometry (default None = off). PG (TG) and SD only. See croak.smearing and Modelling geometric time smearing.

• focal (FocalMixture or None, optional) –

  Chromatic focal mixture (croak.focal) in place of smearing: the focal average is carried explicitly, each node applying its own wavelength-dependent amplitude filter, where the reduced kernel evaluates a single \(|A|^6\) weight at the carrier. Mutually exclusive with smearing, and incompatible with fit_smearing (an explicit mixture has no kernel width to scale).

  Worth its cost for a broadband pulse. The reduced kernel’s widths scale as \(\sigma\propto\lambda\), so freezing them at the carrier is a -16 %/+24 % spread of width across a 1 fs deep-UV band; in our earlier instrument scans that cost a mean 12 % of retrieved duration there, against 3.5 % at 2 fs. The mixture needs one field build per focal node rather than per Gauss-Hermite node, so it is one to two orders of magnitude slower; check quadrature convergence for the chosen pulse and geometry.

• phase_only (bool, optional) – Hold the initial spectral amplitude fixed and retrieve phase only.

• phase_basis ({"pointwise", "bspline"}, optional) – Phase parameterisation; "bspline" uses n_nodes cubic-spline control points (implies phase-only). See croak._jax_pulse.make_spline_phase_basis().

• n_nodes (int, optional) – Number of B-spline control points when phase_basis="bspline".

• R_omega (bool, optional) – Use per-frequency adaptive scaling factors.

• reg_amp (float | None) – Smoothness / spectral-match penalties, identical to croak.lbfgs_ad.LBFGSAD. See Regularisation.

• reg_phase (float) – Smoothness / spectral-match penalties, identical to croak.lbfgs_ad.LBFGSAD. See Regularisation.

• reg_spectrum (float | None) – Smoothness / spectral-match penalties, identical to croak.lbfgs_ad.LBFGSAD. See Regularisation.

• spectrum_target – Smoothness / spectral-match penalties, identical to croak.lbfgs_ad.LBFGSAD. See Regularisation.

• reltol (float, optional) – Convergence tolerances (defaults 1e-4 / 1e-8). abstol is the Optimistix atol; reltol is passed as the Optimistix rtol and drives a relative-improvement plateau stop in run_optx (Optimistix’s own Cauchy termination is atol-dominated and effectively ignores rtol, so the plateau stop is what makes reltol take effect).

• abstol (float, optional) – Convergence tolerances (defaults 1e-4 / 1e-8). abstol is the Optimistix atol; reltol is passed as the Optimistix rtol and drives a relative-improvement plateau stop in run_optx (Optimistix’s own Cauchy termination is atol-dominated and effectively ignores rtol, so the plateau stop is what makes reltol take effect).

• targeterr (float, optional) – Stop early when the FROG error drops below this value.

• verbose (bool, optional) – Print per-step progress.

Notes

Unlike croak.lbfgs_ad.LBFGSAD, the temporal penalty (reg_time / time_window) is not exposed here.

name: str = 'lbfgs-optx'

Algorithm name recorded in the result.

Levenberg–Marquardt

Retrieve by nonlinear least squares with SciPy (lm).

A least-squares retriever that drives the flattened FROG residual \((T_\mathrm{meas}-\mu T_\mathrm{sim})\,w/\sqrt{\mathrm{denom}}\) to zero. It shares the JAX forward model (croak.forward_jax), the residual builder (croak.metrics_jax) and the pulse parameterisation (croak._jax_pulse) with LBFGSAD.

Choice of method

The default method="trf" is SciPy’s trust-region reflective least squares, driven by the analytic JAX Jacobian (forward mode, jax.jacfwd(), since there are far fewer parameters than residuals). Explicit method="lm" selects MINPACK’s classic Levenberg–Marquardt with a finite-difference Jacobian; "dogbox" is also available. MINPACK is not the default because its Fortran driver calls the user Jacobian back from compiled code, and evaluating the JAX/XLA Jacobian program from inside that callback crashes the interpreter (an XLA re-entrancy fault; the JAX residual is fine). "trf" is pure-Python SciPy, accepts the analytic Jacobian, and for these unconstrained problems reaches the same minima; the numeric Jacobian keeps the MINPACK route from crashing.

The normalised trace residual has Euclidean norm \(R\). Regularisation is appended as residual rows, so the squared objective is \(R^2 + \sum_k\lambda_k P_k\), rather than the gradient solvers’ \(R + \sum_k\lambda_k P_k\). Equal penalty weights therefore do not define equal regularised problems. See Retrieval algorithms.

The model supports thin local SHG/PG/SD and dispersive PG/SD, with optional spatial models, phase constraints and fitted scalar parameters as documented below. Bounds depend on the selected least-squares method.

class croak.lm.LM(*, maxiters=300, method='trf', material=None, thickness=0.0, npoints=20, omega0=0.0, quadrature='gausslegendre', smearing=None, focal=None, phase_only=False, phase_basis='pointwise', n_nodes=20, R_omega=False, reg_amp=None, reg_phase=0.0, reg_spectrum=None, spectrum_target=None, tau0=0.0, smear_scale=1.0, smear_scale_delta=1.0, fit_thickness=False, fit_tau0=False, fit_smearing=False, fit_smearing_split=False, polish=False, reltol=0.0001, abstol=1e-08, verbose=False)

Bases: Retriever

Retrieve with SciPy nonlinear least squares and a JAX Jacobian by default.

Parameters:

• maxiters (int, optional) – Maximum number of residual evaluations (max_nfev).

• method ({"trf", "lm", "dogbox"}, optional) – SciPy least_squares method. Default "trf" uses the analytic JAX Jacobian; "lm" falls back to a numeric Jacobian (see the module docstring for why).

• material (str | None) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• thickness (float) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• npoints (int) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• omega0 (float) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• quadrature (str) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• smearing (SmearingKernel or None, optional) – Geometric time-smearing kernel of a non-collinear BOXCARS geometry (default None = off). PG (TG) and SD only. See croak.smearing and Modelling geometric time smearing.

• focal (FocalMixture or None, optional) –

  Chromatic focal mixture (croak.focal) in place of smearing: the focal average is carried explicitly, each node applying its own wavelength-dependent amplitude filter, where the reduced kernel evaluates a single \(|A|^6\) weight at the carrier. Mutually exclusive with smearing, and incompatible with fit_smearing (an explicit mixture has no kernel width to scale).

  Worth its cost for a broadband pulse. The reduced kernel’s widths scale as \(\sigma\propto\lambda\), so freezing them at the carrier is a -16 %/+24 % spread of width across a 1 fs deep-UV band; in our earlier instrument scans that cost a mean 12 % of retrieved duration there, against 3.5 % at 2 fs. The mixture needs one field build per focal node rather than per Gauss-Hermite node, so it is one to two orders of magnitude slower; check quadrature convergence for the chosen pulse and geometry.

• phase_only (bool, optional) – Hold the initial spectral amplitude fixed and retrieve phase only.

• phase_basis ({"pointwise", "bspline"}, optional) – Phase parameterisation. "pointwise" (default) optimises one phase value per frequency bin. "bspline" optimises the n_nodes control points of a cubic B-spline over the spectral support — a strong low-dimensional regulariser that implies phase-only retrieval. See croak._jax_pulse.make_spline_phase_basis().

• n_nodes (int, optional) – Number of B-spline control points when phase_basis="bspline" (default 20); ignored otherwise.

• R_omega (bool, optional) – Use per-frequency adaptive scaling factors.

• reg_amp (float or None, optional) – Amplitude second-difference smoothness weight, added as extra residual rows. None (default) resolves to the default LM-family weight 3e-5 (croak.solver.REG_DEFAULTS; the residual-rows objective needs far smaller weights than the gradient solvers); pass 0 to disable.

• reg_phase (float, optional) – Phase second-difference smoothness weight (default 0).

• reg_spectrum (float or None, optional) – Spectral-match weight pulling the retrieved spectral amplitude towards spectrum_target (full mode only; added as extra residual rows). None (default) resolves to the default LM-family 1e-5 and is inert without a spectrum_target; pass 0 to disable.

• spectrum_target (array_like or None, optional) – Measured spectral intensity on the retrieval grid (e.g. TraceData.Iomega) for reg_spectrum; None disables it.

• reltol (float, optional) – Convergence tolerances (defaults 1e-4 / 1e-8, matching the reference LM solver). reltol is the relative cost-reduction tolerance (SciPy ftol); abstol bounds termination by the step and gradient (SciPy xtol and gtol).

• abstol (float, optional) – Convergence tolerances (defaults 1e-4 / 1e-8, matching the reference LM solver). reltol is the relative cost-reduction tolerance (SciPy ftol); abstol bounds termination by the step and gradient (SciPy xtol and gtol).

• tau0 (float, optional) – Starting delay-zero offset \(\tau_0\) in seconds (default 0) — the centre fit_tau0 optimises around, e.g. the previous fit’s value when warm-starting. Held at this value when fit_tau0 is off.

• smear_scale (float, optional) – Starting multiplier on the smearing-kernel widths (default 1.0 = the kernel as supplied) — the centre fit_smearing optimises around. Held at this value when fit_smearing is off.

• fit_thickness (bool, optional) – Also fit the dispersive-slab thickness as an extra parameter (default False). Requires a dispersive slab (material set, thickness > 0) — PG/SD only. The prior thickness is the starting value. See Fitting the medium thickness and delay-zero.

• fit_tau0 (bool, optional) – Also fit a delay-zero offset \(\tau_0\) (default False), correcting a mis-set experimental time zero (all interactions).

• fit_smearing (bool, optional) – Also fit the geometric-smearing strength as an extra parameter (default False): a dimensionless multiplier on both kernel widths, equivalently the mask ratio d/D. Requires smearing. Do not enable it together with fit_thickness — the smearing width is a nuisance parameter that can absorb dispersion-model error. See Modelling geometric time smearing.

• fit_smearing_split (bool, optional) – Fit the gate-shape (p) and delay (delta) kernel widths as two independent multipliers instead of one joint one (default False) — a diagnostic for which channel deviates from the geometric prediction. Requires fit_smearing; the fitted values land in smear_scale (p) and smear_scale_delta (delta) on the result.

• smear_scale_delta (float, optional) – Starting delay-channel multiplier for a split fit (default 1.0); used only when fit_smearing_split is set.

• polish (bool, optional) – Two-phase polish (default False): first retrieve the pulse with the extra parameters held fixed, then free them for a joint final phase.

• verbose (bool, optional) – Print per-evaluation progress.

name: str = 'lm'

Algorithm name recorded in the result.

Levenberg–Marquardt (Optimistix)

The lm-optx solver uses Optimistix Levenberg–Marquardt on device. The lm interface uses SciPy trust-region reflective least squares; see the algorithm guide for constraints and objectives.

Retrieve with on-device Optimistix Levenberg–Marquardt (lm-optx).

A JAX-native twin of croak.lm.LM: it drives the same flattened FROG residual (croak.metrics_jax.residual_vector(), with the same amplitude/phase/spectral regularisation rows appended) to zero, over the same real parameter vector, with the same JAX forward model. Instead of SciPy least_squares driving a host-side trust-region loop with a JAX Jacobian materialised and transferred every iteration, Optimistix’s Levenberg–Marquardt runs the whole damped Gauss–Newton loop on-device (each step jitted) with the linear algebra handled by lineax. Being pure JAX, it is a genuine Levenberg–Marquardt with the analytic Jacobian, which the MINPACK route of croak.lm.LM cannot offer (see croak.lm).

The implementation is unconstrained: it does not enforce a non-negative fitted thickness, so check fitted values.

See Retrieval algorithms for the objective, linear-system choices and comparison with SciPy’s trust-region solver.

class croak.optimistix_lm.OptxLM(*, maxiters=300, material=None, thickness=0.0, npoints=20, omega0=0.0, quadrature='gausslegendre', smearing=None, focal=None, phase_only=False, phase_basis='pointwise', n_nodes=20, R_omega=False, reg_amp=None, reg_phase=0.0, reg_spectrum=None, spectrum_target=None, tau0=0.0, smear_scale=1.0, smear_scale_delta=1.0, fit_thickness=False, fit_tau0=False, fit_smearing=False, fit_smearing_split=False, polish=False, linear_solver='normal', reltol=0.0001, abstol=1e-08, verbose=False)

Bases: Retriever

Optimistix Levenberg–Marquardt solver with a JAX Jacobian.

The residual conventions match croak.lm.LM; the optimiser runs on-device and is unconstrained, so bounds are not enforced.

Parameters:

• maxiters (int, optional) – Maximum number of solver steps (default 300).

• material (str | None) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• thickness (float) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• npoints (int) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• omega0 (float) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• quadrature (str) – Dispersive-slab parameters; see croak.forward.ForwardModel.

• smearing (SmearingKernel or None, optional) – Geometric time-smearing kernel of a non-collinear BOXCARS geometry (default None = off). PG (TG) and SD only. See croak.smearing and Modelling geometric time smearing.

• focal (FocalMixture or None, optional) –

  Chromatic focal mixture (croak.focal) in place of smearing: the focal average is carried explicitly, each node applying its own wavelength-dependent amplitude filter, where the reduced kernel evaluates a single \(|A|^6\) weight at the carrier. Mutually exclusive with smearing, and incompatible with fit_smearing (an explicit mixture has no kernel width to scale).

  Worth its cost for a broadband pulse. The reduced kernel’s widths scale as \(\sigma\propto\lambda\), so freezing them at the carrier is a -16 %/+24 % spread of width across a 1 fs deep-UV band; in our earlier instrument scans that cost a mean 12 % of retrieved duration there, against 3.5 % at 2 fs. The mixture needs one field build per focal node rather than per Gauss-Hermite node, so it is one to two orders of magnitude slower; check quadrature convergence for the chosen pulse and geometry.

• phase_only (bool, optional) – Hold the initial spectral amplitude fixed and retrieve phase only.

• phase_basis ({"pointwise", "bspline"}, optional) – Phase parameterisation; "bspline" uses n_nodes cubic-spline control points (implies phase-only). See croak._jax_pulse.make_spline_phase_basis().

• n_nodes (int, optional) – Number of B-spline control points when phase_basis="bspline".

• R_omega (bool, optional) – Use per-frequency adaptive scaling factors.

• reg_amp (float | None) – Smoothness / spectral-match penalties, appended as extra residual rows exactly as in croak.lm.LM. See Regularisation.

• reg_phase (float) – Smoothness / spectral-match penalties, appended as extra residual rows exactly as in croak.lm.LM. See Regularisation.

• reg_spectrum (float | None) – Smoothness / spectral-match penalties, appended as extra residual rows exactly as in croak.lm.LM. See Regularisation.

• spectrum_target – Smoothness / spectral-match penalties, appended as extra residual rows exactly as in croak.lm.LM. See Regularisation.

• linear_solver ({"normal", "qr"}, optional) – How to solve the damped Gauss–Newton system each iteration (default "normal"). "normal" forms J^T J + lambda I and takes its Cholesky factor; "qr" factorises the stacked [J; sqrt(lambda) I] operator directly, which is Optimistix’s own default. "normal" was roughly 2x faster on a full-size pointwise retrieval in our benchmarks, because the Gram matrix is a single BLAS gemm while the tall-skinny QR is panel-bound. It squares the condition number, so switch to "qr" if a badly scaled problem stalls or returns a non-finite result. See Retrieval algorithms.

• reltol (float, optional) – Convergence tolerances (defaults 1e-4 / 1e-8). abstol is the Optimistix atol; reltol is passed as the Optimistix rtol and drives a relative-improvement plateau stop in run_optx (stop once the best FROG error improves by less than reltol over the last few steps). The plateau stop is needed because Optimistix’s own Cauchy termination is atol-dominated near convergence and effectively ignores rtol — so without it reltol would have no practical effect.

• abstol (float, optional) – Convergence tolerances (defaults 1e-4 / 1e-8). abstol is the Optimistix atol; reltol is passed as the Optimistix rtol and drives a relative-improvement plateau stop in run_optx (stop once the best FROG error improves by less than reltol over the last few steps). The plateau stop is needed because Optimistix’s own Cauchy termination is atol-dominated near convergence and effectively ignores rtol — so without it reltol would have no practical effect.

• tau0 (float, optional) – Starting (and held) centres for the delay-zero offset (s) and the smearing-width multiplier, exactly as in croak.lm.LM.

• smear_scale (float, optional) – Starting (and held) centres for the delay-zero offset (s) and the smearing-width multiplier, exactly as in croak.lm.LM.

• fit_thickness (bool, optional) – Also fit the dispersive-slab thickness as an extra parameter (default False). Requires a dispersive slab (material set, thickness > 0) — PG/SD only. The prior thickness is the starting value. See Fitting the medium thickness and delay-zero.

• fit_tau0 (bool, optional) – Also fit a delay-zero offset \(\tau_0\) (default False), correcting a mis-set experimental time zero (all interactions).

• fit_smearing (bool, optional) – Also fit the geometric-smearing strength as an extra parameter (default False): a dimensionless multiplier on both kernel widths, equivalently the mask ratio d/D. Requires smearing. Do not enable it together with fit_thickness — the smearing width is a nuisance parameter that can absorb dispersion-model error. See Modelling geometric time smearing.

• fit_smearing_split (bool, optional) – Fit the gate-shape (p) and delay (delta) kernel widths as two independent multipliers instead of one joint one (default False) — a diagnostic for which channel deviates from the geometric prediction. Requires fit_smearing; the fitted values land in smear_scale (p) and smear_scale_delta (delta) on the result.

• smear_scale_delta (float, optional) – Starting delay-channel multiplier for a split fit (default 1.0); used only when fit_smearing_split is set.

• polish (bool, optional) – Two-phase polish (default False): first retrieve the pulse with the extra parameters held fixed, then free them for a joint final phase.

• verbose (bool, optional) – Print per-step progress.

Notes

There is no method argument: Optimistix Levenberg–Marquardt has no trust-region variants to select between.

Optimistix Levenberg–Marquardt is unconstrained, so unlike SciPy’s bounded trf it cannot hold the fitted thickness non-negative — it ignores the lower bound croak._jax_pulse.augment() supplies and relies on starting from a good prior (\(\delta_L = 0\)), exactly as the SciPy method="lm" path does. Fit the thickness as a polish from a good prior.

name: str = 'lm-optx'

Algorithm name recorded in the result.

CMA-ES (global)

The cma-es evolution-strategy retriever (CMA-ES / sep-CMA / DE via evosax) for global search from poor initial guesses, typically paired with the B-spline phase basis. Requires the optional croak[evo] extra.

Global retrieval by evolution strategy (cma-es), powered by evosax.

Every other retriever (except COPRA) is a local optimiser, so a poor initial guess can leave it in the wrong basin. This retriever runs a global population search — CMA-ES (default), Separable CMA-ES or Differential Evolution — over the same differentiable parameterisation the gradient solvers use (full amplitude+phase, phase-only, or the reduced cubic B-spline phase basis). Because the forward model is JAX, the whole population is evaluated in a single jitted vmap device call, so a generation costs roughly one batched forward pass (scaling with the population size and the forward-model quadrature).

The intended workflow is a global search here followed by a local refinement: run cma-es (typically with the low-dimensional phase_basis="bspline"), then feed its result as the guess to lbfgs-ad (the GUI’s “Reuse previous result” initial guess). CMA-ES suits the ~20–30-parameter spline basis; Separable CMA-ES / DE scale better if used on the full-dimensional parameterisation.

The same smoothness / spectral-match regularisation as the gradient solvers is available (reg_amp/reg_phase/reg_spectrum): the penalties are added to the FROG error to form the total the evolution strategy minimises, while selection uses the regularised objective and reporting exposes the bare FROG error separately. The penalties self-gate by mode/basis, so on the common B-spline path only the phase penalty is active (amplitude/spectral are no-ops).

evosax is an optional dependency — pip install "croak[evo]" (or uv add "croak[evo]").

class croak.cmaes.CMAES(*, maxiters=200, strategy='cma', popsize=0, std_init=0.5, material=None, thickness=0.0, npoints=20, omega0=0.0, quadrature='gausslegendre', smearing=None, focal=None, phase_only=False, phase_basis='pointwise', n_nodes=20, R_omega=False, reg_amp=None, reg_phase=0.0, reg_spectrum=None, spectrum_target=None, reg_time=0.0, time_window=None, tau0=0.0, smear_scale=1.0, smear_scale_delta=1.0, fit_thickness=False, fit_tau0=False, fit_smearing=False, fit_smearing_split=False, seed=None, verbose=False)

Bases: Retriever

Global evolution-strategy retriever (cma-es); see the module docstring.

Parameters:

• maxiters (int) – Number of generations.

• strategy ({"cma", "sep-cma", "de"}) – Evolution strategy (CMA-ES, Separable CMA-ES, Differential Evolution).

• popsize (int) – Population size (0 → evosax default 4 + ⌊3 ln N⌋).

• std_init (float) – Initial search std (CMA mean spread / DE population scatter).

• phase_only (bool) – Mode options shared with the gradient solvers (phase_basis="bspline" gives the reduced-dimensional phase search).

• phase_basis (str) – Mode options shared with the gradient solvers (phase_basis="bspline" gives the reduced-dimensional phase search).

• n_nodes (int) – Mode options shared with the gradient solvers (phase_basis="bspline" gives the reduced-dimensional phase search).

• R_omega (bool) – Mode options shared with the gradient solvers (phase_basis="bspline" gives the reduced-dimensional phase search).

• material (str | None) – Forward-model options shared with the gradient solvers.

• thickness (float) – Forward-model options shared with the gradient solvers.

• npoints (int) – Forward-model options shared with the gradient solvers.

• omega0 (float) – Forward-model options shared with the gradient solvers.

• quadrature (str) – Forward-model options shared with the gradient solvers.

• smearing (SmearingKernel or None, optional) – Geometric time-smearing kernel of a non-collinear BOXCARS geometry (default None = off). PG (TG) and SD only. See croak.smearing and Modelling geometric time smearing.

• reg_amp (float | None) – Smoothness / spectral-match regularisation weights, identical to the gradient solvers (see croak.lbfgs_ad.LBFGSAD). The penalties are added to the FROG error to form the total the search minimises and selects on; the reported trace error excludes them. The penalties self-gate by mode/basis, so amplitude/spectral terms are no-ops in phase-only or B-spline mode; reg_spectrum additionally needs an independent spectrum_target (full mode).

• reg_phase (float) – Smoothness / spectral-match regularisation weights, identical to the gradient solvers (see croak.lbfgs_ad.LBFGSAD). The penalties are added to the FROG error to form the total the search minimises and selects on; the reported trace error excludes them. The penalties self-gate by mode/basis, so amplitude/spectral terms are no-ops in phase-only or B-spline mode; reg_spectrum additionally needs an independent spectrum_target (full mode).

• reg_spectrum (float | None) – Smoothness / spectral-match regularisation weights, identical to the gradient solvers (see croak.lbfgs_ad.LBFGSAD). The penalties are added to the FROG error to form the total the search minimises and selects on; the reported trace error excludes them. The penalties self-gate by mode/basis, so amplitude/spectral terms are no-ops in phase-only or B-spline mode; reg_spectrum additionally needs an independent spectrum_target (full mode).

• spectrum_target – Smoothness / spectral-match regularisation weights, identical to the gradient solvers (see croak.lbfgs_ad.LBFGSAD). The penalties are added to the FROG error to form the total the search minimises and selects on; the reported trace error excludes them. The penalties self-gate by mode/basis, so amplitude/spectral terms are no-ops in phase-only or B-spline mode; reg_spectrum additionally needs an independent spectrum_target (full mode).

• reg_time (float) – Temporal (pedestal) penalty weight and (lo, hi) time window in seconds. Penalises the fraction of E(t) = ifft(spectrum) energy falling outside the window — a self-consistent, mode-agnostic alternative to the post-retrieval temporal filter (it keeps spectrum and field an exact Fourier pair, so it introduces no spectral artefacts). Active in all modes including B-spline. time_window=None defaults to the measurement delay range.

• time_window (tuple[float, float] | None) – Temporal (pedestal) penalty weight and (lo, hi) time window in seconds. Penalises the fraction of E(t) = ifft(spectrum) energy falling outside the window — a self-consistent, mode-agnostic alternative to the post-retrieval temporal filter (it keeps spectrum and field an exact Fourier pair, so it introduces no spectral artefacts). Active in all modes including B-spline. time_window=None defaults to the measurement delay range.

• tau0 (float, optional) – Search centres for the delay-zero offset (s, default 0) and the smearing-width multiplier (default 1.0 = the kernel as supplied), e.g. the previous fit’s values when warm-starting. Each is held at its centre when the matching fit_* flag is off.

• smear_scale (float, optional) – Search centres for the delay-zero offset (s, default 0) and the smearing-width multiplier (default 1.0 = the kernel as supplied), e.g. the previous fit’s values when warm-starting. Each is held at its centre when the matching fit_* flag is off.

• focal (FocalMixture or None, optional) –

  Chromatic focal mixture (croak.focal) in place of smearing: the focal average is carried explicitly, each node applying its own wavelength-dependent amplitude filter, where the reduced kernel evaluates a single \(|A|^6\) weight at the carrier. Mutually exclusive with smearing, and incompatible with fit_smearing (an explicit mixture has no kernel width to scale).

  Worth its cost for a broadband pulse. The reduced kernel’s widths scale as \(\sigma\propto\lambda\), so freezing them at the carrier is a -16 %/+24 % spread of width across a 1 fs deep-UV band; in our earlier instrument scans that cost a mean 12 % of retrieved duration there, against 3.5 % at 2 fs. The mixture needs one field build per focal node rather than per Gauss-Hermite node, so it is one to two orders of magnitude slower; check quadrature convergence for the chosen pulse and geometry.

• seed (int or None) – Optional explicit PRNG seed (otherwise derived from the framework rng).

• smear_scale_delta (float)

• fit_thickness (bool)

• fit_tau0 (bool)

• fit_smearing (bool)

• fit_smearing_split (bool)

• verbose (bool)

Notes

The extra parameters are appended to the search vector as dimensionless, O(1)-scaled offsets (croak._jax_pulse.augment()) so the isotropic std_init spread treats them sensibly. The search is unconstrained, so — as for lm-optx — the thickness lower bound is not enforced; it relies on the prior as the search centre. There is no polish option (the global search frees the extras from the start); fit the thickness here only with a decent prior, or refine a pulse-only cma-es result with lbfgs-ad (polish=True) as in the usual global-then-local workflow.

name: str = 'cma-es'

Algorithm name recorded in the result.

croak.cmaes.run_cmaes(objective, u0, *, strategy='cma', popsize=0, std_init=0.5, generations=200, key, callback=None, make_snapshot=None, err_log=None, verbose=False)

Minimise objective(u) with an evosax strategy; return (best_u, R_log).

objective must be a batched map (P, N) -> ((P,), (P,)) (i.e. already jit(vmap(...))) returning, per population member, the (total, frog_error) pair: the total is the regularised objective the evolution strategy is driven by (ask/tell), while the bare frog_error is reported separately. The returned best_u is the member with the lowest total ever evaluated (consistent with the gradient solvers, which return the minimiser of the regularised objective), and R_log is that member’s bare FROG error after each generation.

callback is called as callback(generation, gen_R, best_R, snapshot=...) where snapshot is a zero-argument builder of the current best RetrievalResult (for the GUI live preview) when make_snapshot(best_u) -> RetrievalResult is supplied, else None. Pass err_log to accumulate into an existing list a snapshot builder can read.

Parameters:

• u0 (NDArray[float64])

• strategy (str)

• popsize (int)

• std_init (float)

• generations (int)

• make_snapshot (Callable[[NDArray[float64]], object] | None)

• err_log (list[float] | None)

• verbose (bool)

Return type:

tuple[NDArray[float64], list[float]]

Warm-started L-BFGS

The default warm-lbfgs retriever: LBFGSAD on the full extended model, started from a COPRA local-projection sweep rather than from the raw initial guess. The final stage uses the lbfgs-ad objective; the COPRA start makes it much less sensitive to the initial guess than lbfgs-ad alone.

Retrieve with L-BFGS-AD after a COPRA local-sweep warm start.

COPRA has two stages that behave very differently. Its local sweep visits delays in random order and replaces each one’s signal modulus with the measured value, which injects far more measured information per iteration than a gradient step; its global stage is a gradient descent, and on the cases we measured engaging it made the result worse. The local sweep alone therefore produces a good, cheap estimate, but one that cannot represent geometrical smearing: a smeared trace is an incoherent mixture with no single signal whose modulus to replace.

This solver uses that estimate as the starting point for the regularised extended model under L-BFGS. The final objective and options are those of croak.lbfgs_ad.LBFGSAD; only the initial guess changes. The warm-up uses the core dispersive model (material, thickness, quadrature, and phase_only when set) and, when supplied, the reduced smearing kernel to score its iterates. It does not include the chromatic focal or collection model, fitted extras, regularisation or marginal-peak centring of the final stage, and its local projection is the unsmeared single-field update.

Measured effect

On simulated TG-FROG traces (the Gaussian-beam control and the apertured mask, each at transform limit and ±2 fs², plus a mask ±0.625 fs² pair; eight cases at a 9.5 µm substrate), compared with a cold lbfgs-ad start:

• Accuracy: mean complex error 0.333 → 0.208 (38 % lower), median 0.408 → 0.150 (63 % lower); mean |dFWHM| 17.7 % → 8.5 %, worst 69.3 % → 20.8 %. The duration improved in seven of the eight cases and the complex error degraded by more than 20 % in none. It rescues bad basins rather than refining good solutions.

• Variance over starting points is the largest effect: the spread in retrieved duration collapsed five-fold to twelve-fold on every chirped arm (e.g. ±36.9 → ±3.1 points), while the transform-limited arms, which were never start-sensitive, were unchanged.

• Known regression: on the Gaussian transform-limited arm at depths away from 9.5 µm the warm complex error reached 0.63 against a cold maximum of 0.18, i.e. where the cold start was already reliable.

Noise overfitting

The COPRA warm-up has no spectral penalty. On a noisy trace whose grid spans far more bandwidth than the signal, any solver without a spectral target can fit noise into signal-free spectral bins near the grid edges; the warm-up does so before the penalised stage starts. On the noisy trace of examples/example_thickness_uncertainty.py, both warm-lbfgs and lbfgs-ad without a target leave ~3 % of the spectral energy at the edges and return durations from about 1.9 to 3.2 fs for a 3.0 fs pulse. With the measured spectrum as target, lbfgs-ad recovers 2.97 fs while warm-lbfgs still ends at 3.16 fs. Check croak.edge_energy_fraction() on noisy data and supply the measured spectrum.

The result combines both convergence histories. COPRA sweeps and L-BFGS function evaluations are different units of work; stage_boundaries marks the handover.

class croak.warm_lbfgs.WarmLBFGS(*args, warmup_iters=300, **kwargs)

Bases: LBFGSAD

Retrieve with a COPRA local-sweep initialisation followed by L-BFGS-AD.

Accepts every LBFGSAD argument for the final stage. The warm-up inherits the dispersive model and the smearing kernel, so it judges its iterates against the same physics: its projection step cannot use the kernel, but the error by which it selects the returned iterate is the smeared one. Per-frequency factors are off during the warm-up (COPRA’s projection has no use for them) and free to move afterwards. See croak.warm_lbfgs for the remaining differences.

Parameters:

warmup_iters (int, optional) – Maximum COPRA warm-up sweeps (default 300); the sweep normally converges well inside this. Separate from the final L-BFGS maxiters budget.

name: str = 'warm-lbfgs'

Name recorded in results and used by uncertainty replicates when no alternative solver is requested. It must not inherit LBFGSAD’s label: a warm-started result would then claim to be a plain lbfgs-ad (contradicting its stage_boundaries), and a bootstrap would re-run its replicates with the cold solver.

Result

The RetrievalResult returned by every retrieval algorithm.

class croak.result.RetrievalResult(spectrum, grid, delays, interaction, algorithm, error, errors=<factory>, trace=None, mu=1.0, omega0=0.0, thickness=None, tau0=0.0, smear_scale=None, smear_scale_delta=None, stage_boundaries=<factory>)

Bases: object

Outcome of a pulse retrieval, uniform across all algorithms.

Parameters:

• spectrum (NDArray[complex128])

• grid (Grid)

• delays (NDArray[float64])

• interaction (str)

• algorithm (str)

• error (float)

• errors (list[float])

• trace (NDArray[float64] | None)

• mu (float | NDArray[float64])

• omega0 (float)

• thickness (float | None)

• tau0 (float)

• smear_scale (float | None)

• smear_scale_delta (float | None)

• stage_boundaries (list[int])

spectrum

Retrieved complex spectrum \(\tilde E(\omega)\) (centred order).

Type:

numpy.ndarray

grid

Time/frequency grid the retrieval ran on.

Type:

Grid

delays

Delay axis (s).

Type:

numpy.ndarray

interaction

Interaction name ("shg"/"sd"/"pg").

Type:

str

algorithm

Algorithm name — one of the keys of croak.retrieve.ALGORITHMS ("copra", "copra-jax", "lbfgs", "lbfgs-hand", "lbfgs-ad", "warm-lbfgs", "lbfgs-optx", "lm", "lm-optx", "cma-es").

Type:

str

error

Final FROG error R.

Type:

float

errors

FROG error at each iteration.

Type:

list of float

trace

Retrieved (simulated) trace, (Nomega, Ndelay), optimally scaled.

Type:

numpy.ndarray

mu

Final intensity scale factor (scalar, or per-frequency in Rω mode).

Type:

float or numpy.ndarray

thickness

Fitted dispersive-slab thickness (m) when the retriever fitted it as an extra parameter (lbfgs-ad/lm with fit_thickness); None when thickness was held fixed.

Type:

float or None

tau0

Fitted delay-zero offset (s) when the retriever fitted it (fit_tau0); 0.0 otherwise. Sign convention: the true delay-zero of the measured trace lies at +tau0, i.e. the simulated gate uses the delay phase exp(iω(τ − τ0)).

Type:

float

smear_scale

Fitted multiplier on the geometric-smearing kernel widths when the retriever fitted it (fit_smearing); None when the kernel was held fixed or absent. 1.0 means the fit kept the nominal widths. For a split fit (fit_smearing_split) this is the gate-shape (p) channel multiplier.

Type:

float or None

smear_scale_delta

Fitted multiplier on the delay (delta) kernel width when the two smearing channels were fitted separately (fit_smearing_split); None for a joint fit (where smear_scale covers both channels) or when smearing was not fitted.

Type:

float or None

stage_boundaries

Indices into errors where one stage of a multi-stage solver handed over to the next, so a convergence plot can mark the join. Empty for the single-stage solvers. A spliced curve needs this: warm-lbfgs counts one COPRA iteration (a full local sweep over every delay) before the handover and one L-BFGS function evaluation after it, and those are not the same unit of work — an unmarked kink there reads as convergence behaviour rather than as a change of algorithm.

Type:

list of int

property omega: NDArray[float64]

Angular-frequency axis (rad/s, centred).

property wavelength: NDArray[float64]

\(2\pi c/(\omega+\omega_0)\).

Type:

Absolute wavelength axis (m)

property t: NDArray[float64]

Time axis (s).

property field: NDArray[complex128]

Retrieved temporal field \(E(t)\).

property intensity_omega: NDArray[float64]

Spectral intensity \(|\tilde E(\omega)|^2\).

property phase_omega: NDArray[float64]

Spectral phase \(\arg \tilde E(\omega)\), unwrapped.

property intensity_t: NDArray[float64]

Temporal intensity \(|E(t)|^2\).

property phase_t: NDArray[float64]

Temporal phase \(\arg E(t)\), unwrapped.

Progress reporting

The live progress bar and end-of-retrieval summary, installed automatically by retrieve() (controlled by its progress argument) unless you pass your own callback.

Terminal / notebook progress reporting for croak.retrieve().

When the caller of retrieve() does not pass its own callback, a ProgressReporter is installed automatically. It renders a live progress bar and prints a one-line summary report when retrieval finishes (final FROG error, retrieved and transform-limited FWHM, iteration count and wall-clock time).

The live bar is drawn with tqdm.auto, which picks the right widget for the environment: an ipywidgets bar in a Jupyter notebook (or a plain text bar if ipywidgets is not installed) and a text bar on the terminal. If tqdm is somehow unavailable, a small built-in carriage-return bar is used instead so the module stays self-contained.

The GUI passes its own callback and so is left untouched: it draws its own convergence curve and status line (see The wizard GUI).

Environments

detect_environment distinguishes three cases the wrapper cares about:

"notebook"

A Jupyter / IPython kernel (ZMQInteractiveShell); tqdm.auto renders an inline bar and the report is printed below it.

"terminal"

An interactive TTY. The bar is written to stderr so it does not pollute anything piped from stdout.

"plain"

Non-interactive (output redirected to a file, a CI log, a pytest capture, …). Under the default progress="auto" no live bar or report is shown, so library calls in scripts and test suites stay silent; pass progress=True to force the report anyway.

class croak.progress.ProgressReporter(*, total=None, env=None, live=True, label='retrieving', width=30, min_interval=0.1, use_tqdm=True)

Bases: object

Live progress bar + final summary report for a single retrieval.

Its callback() matches the solver hook signature callback(iteration, R, best_R, snapshot=None) and is what croak.retrieve.retrieve() forwards to the solver (the trailing snapshot is ignored here). finish() is called once the solver returns.

Parameters:

• total (int | None)

• env (str | None)

• live (bool)

• label (str)

• width (int)

• min_interval (float)

• use_tqdm (bool)

callback(iteration, R, best_R, snapshot=None)

Solver progress hook: record the latest iteration and redraw the bar.

snapshot (an optional zero-argument in-progress-result builder some solvers pass for the GUI live preview) is ignored here — the terminal reporter only needs the scalar progress.

Parameters:

• iteration (int)

• R (float)

• best_R (float)

Return type:

None

finish(result, elapsed=None, *, processed=None)

Close the live bar (if any) and print the summary report.

Parameters:

• result (RetrievalResult)

• elapsed (float | None)

• processed (ProcessedResult | None)

Return type:

None

croak.progress.detect_environment()

Return "notebook", "terminal" or "plain" (see module docs).

Return type:

str

croak.progress.make_reporter(progress, *, total)

Build a ProgressReporter for the progress setting, or None.

progress is the croak.retrieve.retrieve() argument:

• "auto" — report only in a terminal or notebook ("plain" is silent),

• True — always report, forcing a live bar even when non-interactive,

• False — never report (returns None).

Parameters:

total (int | None)

Return type:

ProgressReporter | None

croak.progress.format_report(result, elapsed, *, processed=None)

Format the end-of-retrieval summary (error, FWHMs, iterations, time).

Parameters:

• result (RetrievalResult)

• elapsed (float)

• processed (ProcessedResult | None)

Return type:

str

Metrics

Trace-error metrics for nonlinear-process pulse retrieval (FROG, TDP, …).

All retrieval algorithms compare a measured trace T_meas against a simulated trace T_sim (both intensity) through three quantities:

• the optimal scale factor \(\mu\) minimising the weighted residual,

• the residual sum of squares \(r\),

• the normalised, dimensionless trace error \(R\) (the standard “FROG error”, but the same metric for any PNPS trace).

Traces may be in either orientation for the global metrics (which sum over all elements); compute_mu_per_freq() is the exception and takes an explicit frequency axis.

croak.metrics.compute_mu(t_meas, t_sim, weights)

Optimal global scale factor \(\mu\).

Minimises \(\sum (T_\mathrm{meas}-\mu T_\mathrm{sim})^2 w^2\), giving \(\mu = \sum(T_\mathrm{meas}T_\mathrm{sim}w^2)/\sum(T_\mathrm{sim}^2 w^2)\).

Parameters:

• t_meas (ArrayLike)

• t_sim (ArrayLike)

• weights (ArrayLike)

Return type:

float

croak.metrics.compute_mu_per_freq(t_meas, t_sim, weights, *, freq_axis=0)

Per-frequency scale factors (for the Rω adaptive-scaling mode).

Each factor is the per-row least-squares optimum sum(t_sim * t_meas) / sum(t_sim**2). The weights multiply numerator and denominator identically and therefore cancel: the factors are weight-independent. Their only effect is that a zero-weight frequency slice has a vanishing denominator and falls back to 1.0, so rows excluded from the fit are left unscaled rather than scaled by a meaningless ratio.

Parameters:

• t_meas (array_like) – Measured and simulated traces (2-D).

• t_sim (array_like) – Measured and simulated traces (2-D).

• weights (array_like) – Per-frequency weight vector. Only its zero pattern matters (see above).

• freq_axis (int, optional) – Axis corresponding to frequency (default 0).

Returns:

numpy.ndarray – One scale factor per frequency slice (1.0 where the denominator vanishes).

Return type:

NDArray[float64]

croak.metrics.compute_r(t_meas, t_sim, weights, mu)

Weighted residual sum of squares.

\(r = \sum[(T_\mathrm{meas}-\mu T_\mathrm{sim})w]^2\).

mu is the optimal scale: a scalar (global) or a per-frequency column vector broadcast against the trace (the Rω mode).

Parameters:

• t_meas (ArrayLike)

• t_sim (ArrayLike)

• weights (ArrayLike)

• mu (float | NDArray[float64])

Return type:

float

croak.metrics.compute_R(r, t_meas, weights)

Normalised FROG error \(R = \sqrt{r/(MN\,\max(T_\mathrm{meas}w)^2)}\).

Dimensionless and independent of the global intensity scale; typically 0.001–0.05 for a good retrieval. Values still depend on grid coverage, weights, preprocessing and noise, and a small trace error does not by itself certify the retrieved field.

Parameters:

• r (float)

• t_meas (ArrayLike)

• weights (ArrayLike)

Return type:

float

croak.metrics.trace_error(t_meas, t_sim, weights=None)

General name for the normalised trace error (the “FROG error” generalised to any nonlinear-process trace). Alias of frog_error().

Parameters:

• t_meas (ArrayLike)

• t_sim (ArrayLike)

• weights (ArrayLike | None)

Return type:

float

croak.metrics.frog_error(t_meas, t_sim, weights=None)

Return the FROG error R via compute_mu() / r / R.

Parameters:

• t_meas (array_like) – Measured and simulated traces.

• t_sim (array_like) – Measured and simulated traces.

• weights (array_like, optional) – Weight matrix (same shape as the traces) or per-frequency vector (broadcast along axis 0). Defaults to uniform weights.

Return type:

float

Delay-origin convention

Match the model’s delay origin to preprocessing. See weighted gains and centring for the distinction between marginal-peak centring and fitting a free delay offset.

Match a model trace to the selected delay-origin convention.

A measured trace does not know where its delay zero is. The preprocessing (croak.preprocess.regrid()) therefore centres the data on their delay-marginal peak, to sub-sample precision. The forward model instead places zero delay at gate–probe coincidence. For a delay-symmetric trace the two conventions agree. For a tightly apertured single-cycle instrument they do not: coherent collection makes the trace asymmetric and puts its marginal peak tens of attoseconds from coincidence, growing with the chirp accumulated in the medium. A retrieval with the delay zero held at coincidence can pay for that misalignment only by broadening the pulse; in our earlier scans of the simulated TG-FROG instrument this accounted for a ~10 % duration over-report at a 0.5 mm collection hole.

Fitting a free delay offset (fit_tau0) works for transform-limited pulses but not for chirped ones: the chirp itself displaces the marginal peak, so tau0 wanders by hundreds of attoseconds and trades against the retrieved chirp. This module implements the alternative: with origin="marginal_peak", centre the model trace on its own marginal peak by the same parabolic rule the preprocessing applies to the data, at every evaluation, including the row-gain weighting when measured data are supplied. Both sides then share one convention, for any pulse, with no free parameter. This is a coordinate convention, distinct from fitting an unknown delay offset.

The shift is an exact Fourier translation along the uniform delay axis, so it costs two FFTs per evaluation. The translation is periodic, so the trace must have negligible support at both delay ends. The parabolic peak position is differentiated while its discrete peak index (argmax) is held fixed; the map is only piecewise differentiable where competing peaks exchange order. See Retrieval as nonlinear least squares for gains and centring assumptions.

croak.delay_origin.DELAY_ORIGINS = ('coincidence', 'marginal_peak')

delay zero at gate–probe coincidence (the model’s native convention; the default). "marginal_peak": the model trace is shifted so its delay-marginal peak sits at delay zero, matching the data’s preprocessing.

Type:

"coincidence"

croak.delay_origin.marginal_peak_delay_jax(delays, marginal)

Estimate the marginal-peak delay with a local parabola.

Parameters:

• delays (jax.Array) – Uniform delay samples in seconds, with at least three points.

• marginal (jax.Array) – One marginal value per delay sample.

Returns:

jax.Array – Scalar peak delay in seconds. The fitted vertex is clamped to the neighbouring samples. For an edge maximum, the current implementation uses the adjacent interior delay sample; avoid clipped traces. A non-concave interior triplet uses its central sample. The peak index is not differentiated; the vertex location within that branch is.

Return type:

Array

croak.delay_origin.recentre_trace(trace, delays, row_weights=None)

Shift a trace so its estimated delay-marginal peak sits at zero.

Parameters:

• trace (jax.Array) – Trace with shape (Nomega, Ndelay).

• delays (jax.Array) – Uniform delay samples in seconds, with at least three points.

• row_weights (jax.Array or None, optional) – Length-Nomega intensity gains applied only when forming the marginal. Pass the row gains mapping the prediction onto the measured trace when the data have a different spectral response; see recentring().

Returns:

jax.Array – Translated trace, with its original shape. The Fourier translation is periodic, so check support at both delay boundaries. The supplied row gains set the centring position but are not applied to the returned intensity. This low-level function assumes a uniform delay axis.

Return type:

Array

croak.delay_origin.recentring(trace_fn, delays, origin='coincidence', t_meas=None)

Wrap a trace function so its output honours origin.

Parameters:

• trace_fn (callable) – trace_fn(*args) -> (Nomega, Ndelay) array (e.g. the trace_param of croak.forward_jax.make_param_trace_fn()).

• delays (array_like) – The (uniform) delay axis the trace is evaluated on (s).

• origin ({"coincidence", "marginal_peak"}) – "coincidence" returns trace_fn unchanged.

• t_meas (array_like, optional) – The measured trace on the same grid. When given, the model’s marginal is formed with the per-row least-squares factors \(\mu_\omega = \langle T_{meas}, T\rangle_\omega / \langle T, T\rangle_\omega\) that map the model rows onto the data rows, so the two marginals carry the same row weighting whatever generation-response exponent the data were corrected with. Recommended; the factors are recomputed at every evaluation from the current model, and differentiated through.

Returns:

callable