Workflow

Loading files, cleaning and regridding measured traces, the high-level pipeline, post-processing, dispersion tuning (including chirped mirrors, back-propagated and beam-path coating mirrors) and saving. The dispersion sources here are described in Materials and mirrors.

I/O

Loading experimental FROG/spectrum files and unit handling.

A thin, GUI-agnostic layer over HDF5 (h5py), NumPy .npz and delimited text files. It enumerates the numeric datasets in a file, auto-matches the FROG trace / wavelength / scan-axis arrays by name, loads selected datasets, and converts the loaded axes to SI units.

The actual file choosing is the GUI’s job (QFileDialog); this module only reads paths it is given, so the same code serves the library and the GUI.

class croak.io.DatasetInfo(name, ndim, shape)

Bases: object

Metadata for one numeric dataset within a file.

Parameters:

• name (str)

• ndim (int)

• shape (tuple[int, ...])

class croak.io.FileDatasets(path, filetype, datasets)

Bases: object

The numeric datasets discovered in a file, with lazy loading.

Parameters:

• path (str)

• filetype (str)

• datasets (list[DatasetInfo])

path

File path.

Type:

str

filetype

"hdf5" or "npz".

Type:

str

datasets

All numeric datasets, in discovery order.

Type:

list of DatasetInfo

property names: list[str]

Dataset names, in discovery order.

property matrix_names: list[str]

Names of the 2-D datasets (candidate FROG traces).

property vector_names: list[str]

Names of the 1-D datasets (candidate wavelength/scan axes).

property option_labels: list[str]

name  (shape) labels for every dataset (for spectrum pickers).

load(name)

Load one dataset by name as a float array.

Parameters:

name (str)

Return type:

NDArray[float64]

auto_index(kind)

Index of the best dataset for kind (ifrog/lambda/scan).

Parameters:

kind (str)

Return type:

int | None

croak.io.list_datasets(path)

Enumerate the numeric datasets in an HDF5 or NPZ file.

Parameters:

path (str) – Path to a .h5/.hdf5 or .npz file.

Returns:

FileDatasets

Return type:

FileDatasets

croak.io.load_dataset(path, name, filetype=None)

Load a single dataset from an HDF5 or NPZ file as float64.

Parameters:

• path (str) – File path.

• name (str) – Dataset name/path within the file.

• filetype (str, optional) – "hdf5" or "npz"; inferred from the extension if omitted.

Return type:

NDArray[float64]

croak.io.load_csv(path, *, skip=None, delimiter='auto')

Load a delimited numeric text file (CSV or whitespace).

Comment (#) and blank lines are skipped; the header is auto-detected when skip is omitted (first all-numeric row).

Parameters:

• path (str) – File path.

• skip (int, optional) – Number of leading data rows to skip. Auto-detected if omitted.

• delimiter (str or None, optional) – ",", None (whitespace), or "auto" (default) to detect.

Returns:

numpy.ndarray – A 2-D (rows, columns) float array.

Return type:

NDArray[float64]

croak.io.match_dataset(names, pattern)

Index of the shortest dataset name matching pattern (ties: first).

Parameters:

• names (list[str])

• pattern (Pattern)

Return type:

int | None

croak.io.guess_scanaxis_type(name)

Guess whether a scan-axis name denotes a "delay" or "position".

Parameters:

name (str)

Return type:

str | None

croak.io.guess_wavelength_unit(values)

Guess the unit ("m"/"µm"/"nm") of a wavelength axis.

Parameters:

values (NDArray[float64])

Return type:

str

croak.io.guess_scanaxis_unit(values, scan_type)

Guess the unit of a scan axis given its type ("delay" or "position").

Parameters:

• values (NDArray[float64])

• scan_type (str)

Return type:

str

croak.io.csv_wavelength_unit(path)

Guess the wavelength unit of a 2-column CSV’s first column (nm on error).

Parameters:

path (str)

Return type:

str

croak.io.to_1d(arr, name='')

Flatten a near-1-D array (e.g. (N,1)) to 1-D; error if genuinely 2-D.

Used for wavelength axes loaded from HDF5/NPZ where the dataset may carry a singleton dimension.

Parameters:

• arr (NDArray)

• name (str)

Return type:

NDArray[float64]

croak.io.reduce_to_1d(arr, name, match_len)

Reduce an intensity array to length match_len (averaging multi-shot axes).

A 2-D (N_lambda, N_shots) array is averaged over the non-matching axis. Square arrays default to averaging axis 1.

Parameters:

• arr (NDArray)

• name (str)

• match_len (int)

Return type:

NDArray[float64]

croak.io.unit_to_si(unit)

Multiplicative factor converting unit to SI (metres or seconds).

Parameters:

unit (str)

Return type:

float

croak.io.position_to_delay(z)

Convert stage position (m) to double-pass delay (s), \(\tau = 2z/c\).

Parameters:

z (NDArray[float64])

Return type:

NDArray[float64]

croak.io.UNIT_FACTORS: dict[str, float] = {'fs': 1e-15, 'm': 1.0, 'mm': 0.001, 'nm': 1e-09, 'ps': 1e-12, 's': 1.0, 'um': 1e-06, 'µm': 1e-06}

Multiplicative factors converting a unit string to SI (metres or seconds).

class croak.io.MaskWindowSpec(window_type, hole_x, hole_y, hole_diameter, z_mask, apod, apod_param, delta_k, reference_wavelength, omega_reference, weighting='amplitude')

Bases: object

The collection window a simulated scan was recorded through.

scansave files store the signal-extraction window as flattened /grid/window_def_* scalars rather than the (gigabyte-scale) window array, so the aperture can be rebuilt losslessly. This is the raw record; turning it into a quadrature is croak.collection.aperture_from_scan().

Parameters:

• window_type (str)

• hole_x (float)

• hole_y (float)

• hole_diameter (float)

• z_mask (float)

• apod (str)

• apod_param (float | None)

• delta_k (float)

• reference_wavelength (float)

• omega_reference (float)

• weighting (str)

window_type

The simulator’s window class, e.g. "PhysicalMaskWindow".

Type:

str

hole_x, hole_y

Hole centre in the mask plane (m).

Type:

float

hole_diameter

Hole diameter (m).

Type:

float

z_mask

Mask-plane distance from the focus (m).

Type:

float

apod

Apodisation form: "hard", "supergauss", "tanh" or "rcos" (the compact raised-cosine edge pnps windows use).

Type:

str

apod_param

Apodisation parameter, or None when the file recorded the simulator’s "default" — which depends on the transverse k-grid and so must be re-derived (croak.collection.resolve_tanh_width()) from delta_k, omega and reference_wavelength.

Type:

float or None

weighting

How the edge profile enters the collected energy: "amplitude" (a field filter, squared in the energy — ModelPNPS files, and the default when the record has no weighting scalar) or "quadrature" (a weight over a hard hole, entering once — pnps files). See croak.collection.CollectionAperture.

Type:

str

delta_k

Transverse k-grid spacing of the simulation (rad/m), from /grid/kx.

Type:

float

reference_wavelength

The carrier wavelength (m) the simulator’s defaults were evaluated at.

Type:

float

omega_reference

The simulation’s own frequency grid point nearest 2 pi c / reference_wavelength (rad/s) — which is what the simulator evaluates its default apodisation width at, and on a coarse spectral grid is not the same number.

Type:

float

class croak.io.SimulatedScan(omega, omega0, tau, trace, Iomega, Iomega_beamlet, t, It, tau_fwhm, window_key, z_index, z_positions, delay_convention='legacy', Eomega_beamlet=None, Eomega=None, Iomega_beamlet_reimaged=None, Eomega_beamlet_reimaged=None, mask_diameter=None, mask_spacing=None, f_foc=None, geometry=None, mask_window=None)

Bases: object

A numerically simulated FROG scan read from a scansave HDF5 file.

The arrays are returned in croak’s conventions: the angular-frequency axis is absolute and sorted ascending, and the trace is oriented (Nomega, Ndelay). Spectral quantities (trace, Iomega, Iomega_beamlet) are angular-frequency densities \(|\tilde E|^2\), not the wavelength densities a spectrometer records — converting between the two is the caller’s job (see croak.preprocess.omega_to_lambda_density()).

Parameters:

• omega (NDArray[float64])

• omega0 (float)

• tau (NDArray[float64])

• trace (NDArray[float64])

• Iomega (NDArray[float64])

• Iomega_beamlet (NDArray[float64] | None)

• t (NDArray[float64])

• It (NDArray[float64])

• tau_fwhm (float)

• window_key (str)

• z_index (int)

• z_positions (NDArray[float64] | None)

• delay_convention (str)

• Eomega_beamlet (NDArray[complex128] | None)

• Eomega (NDArray[complex128] | None)

• Iomega_beamlet_reimaged (NDArray[float64] | None)

• Eomega_beamlet_reimaged (NDArray[complex128] | None)

• mask_diameter (float | None)

• mask_spacing (float | None)

• f_foc (float | None)

• geometry (str | None)

• mask_window (MaskWindowSpec | None)

omega

Absolute angular-frequency axis (rad/s), ascending. Shape (Nomega,).

Type:

numpy.ndarray

omega0

Carrier angular frequency (rad/s).

Type:

float

tau

Delay axis (s). Shape (Ndelay,).

Type:

numpy.ndarray

trace

Signal spectral intensity (Nomega, Ndelay) (ω-density), as simulated.

Type:

numpy.ndarray

Iomega

Reference (input) pulse spectrum (Nomega,) (ω-density).

Type:

numpy.ndarray

Iomega_beamlet

Vignetted gate-beam spectrum (Nomega,) for the mask correction, or None when the file does not provide it.

Type:

numpy.ndarray or None

t

Truth temporal axis (s); the oversampled To if present.

Type:

numpy.ndarray

It

Truth temporal intensity for the chosen truth_source (the post-mask It_beamlet by default, else the source It); the oversampled Ito* if present.

Type:

numpy.ndarray

tau_fwhm

Intensity FWHM of the input pulse (s).

Type:

float

window_key

The trace-window dataset that was read.

Type:

str

z_index

The resolved (non-negative) propagation-slice index that was loaded. The selected propagation distance is z_positions[z_index] when z_positions is available.

Type:

int

z_positions

The saved propagation distances (m) of the slices stacked along the window’s nz axis (the file’s /grid/zsave), strictly increasing from the entrance (0) to the substrate exit (zmax). None for the older single-thickness format that does not store zsave.

Type:

numpy.ndarray or None

delay_convention

The file’s delay-axis convention: "gate" when the simulation stored the trace directly in the gate-delay (paper) frame (the file carries /grid/delay_convention), else "legacy" — the probe-delayed frame of older files, whose delay axis must be negated on loading (SimulatedLoadParams.reverse_trace).

Type:

str

Eomega_beamlet

The complex vignetted gate-beam spectrum (Eω_beamlet_re/_im in newer files): a stored reference including its spectral phase. This is the transverse-integrated beamlet convention; convert to a common field reference before comparing with a retrieved pulse. Enables complex-field retrieval-error metrics and direct truth-GDD measurement, the temporal/spectral truth-phase overlays, and truth-seeded retrievals. None for files that store only intensities.

Type:

numpy.ndarray or None

Eomega

The complex pre-mask input (source) spectrum, same encoding; None when absent.

Type:

numpy.ndarray or None

Iomega_beamlet_reimaged: NDArray[float64] | None = None

On-axis beamlet spectrum and field, when stored. In the declared diffraction convention this differs from the transverse-integrated beamlet by one power of absolute frequency in field amplitude. It is the entrance field represented by the chromatic focal model; the detector’s collection mode is a separate choice.

mask_diameter: float | None = None

Mask geometry recorded in the scan (m). Use these values to keep the model geometry consistent with the generating simulation. The reduced kernel uses d/D = (spacing + D)/(2D). None if not recorded.

f_foc: float | None = None

Focusing focal length (m), when recorded. The reduced smearing kernel never needs it (f cancels between spot size and tilt), but the explicit chromatic focal mixture (croak.focal) does.

geometry: str | None = None

Interaction geometry the file was generated with ("tg" / "sd"), when recorded. It selects the smearing layout: a two-beam SD has no gate-shape channel at all, so the boxcars kernel is wrong for it.

mask_window: MaskWindowSpec | None = None

The collection window the trace was recorded through, when the file stores it. Needed to model the finite collection aperture (croak.collection); None for files written before the simulator recorded the window.

croak.io.SIMULATED_WINDOW_KEYS: tuple[str, ...] = ('Iω_win', 'Iω_win_reimaged', 'Iω_full')

the aperture-windowed signal integrated over all k (Iω_win) or at the on-axis re-imaged pixel (Iω_win_reimaged), and the full signal-beam collection without the aperture crop (Iω_full, newer files only). The unicode ω is part of the on-disk key string. read_simulated_scan reads any of these as the trace; not every file stores every one (use read_simulated_window_keys() to list those actually present).

Type:

HDF5 dataset names of the simulated-scan trace windows

croak.io.SIMULATED_TRUTH_SOURCES: tuple[str, ...] = ('beamlet', 'beamlet_reimaged', 'source')

the post-mask gate beamlet (transverse-integrated, It_beamlet), its ON-AXIS counterpart (It_beamlet_reimaged — one power of ω bluer in amplitude, the field the chromatic focal mixture’s ew represents; written by pnps), and the pre-mask source (It).

Type:

Time-domain truth sources a scan may store, in preference order

croak.io.read_simulated_scan(path, *, window_key='Iω_win', z_index=-1, z_thickness=None, truth_source='beamlet')

Read a numerically simulated FROG scan from a scansave HDF5 file.

The file layout (a Luna scansave output) is a /grid group with the frequency/time axes and reference pulse, a /scanvariables/τ delay axis, and one or more top-level trace windows (Iω_win, Iω_win_reimaged). This reader only reads and orients; the wavelength conversion, Jacobian, third-order and mask corrections live in croak.session.pipeline.assemble_simulated_load_data().

Parameters:

• path (str) – Path to the .h5/.hdf5 scan file.

• window_key (str, optional) – Trace-window dataset to read: "Iω_win" (first collection window, default) or "Iω_win_reimaged" (on-axis re-imaged).

• z_index (int, optional) – Propagation-slice index along the window’s nz axis. Default -1 (the last slice, i.e. the substrate exit), matching the Julia :end. Ignored when z_thickness is given.

• z_thickness (float, optional) – Material thickness (m) to load from a multi-thickness scan. When given, the slice whose saved position /grid/zsave is nearest this value is selected (argmin|zsave − z_thickness|), taking precedence over z_index. Requires the file to store /grid/zsave (the newer multi-thickness format); a value is invalid on the older format.

• truth_source ({"beamlet", "beamlet_reimaged", "source"}, optional) – Which stored reference pulse to return as the time-domain truth (t, It): the post-mask gate beamlet (It_beamlet/Ito_beamlet, default — the integrated beamlet reference) or the pre-mask ideal source input (It/Ito). The oversampled grid (To/Ito*) is used when present. Falls back to the other source when the requested one is absent (e.g. a Gaussian-beam run with no beamlet, or an older file with only It). See read_simulated_truth_keys().

Returns:

SimulatedScan – Axes and trace in croak conventions (absolute ω ascending, (Nomega, Ndelay)), with the resolved z_index and (when present) the z_positions thickness axis.

Raises:

• KeyError – If the /grid group, the /scanvariables/τ axis, or window_key is missing (the error lists the available top-level datasets), or if z_thickness is requested but the file has no /grid/zsave.

• IndexError – If z_index is out of range for the window’s nz axis.

• ValueError – If path is not an HDF5 file or truth_source is not one of SIMULATED_TRUTH_SOURCES.

Return type:

SimulatedScan

croak.io.read_simulated_z_positions(path)

Read just the propagation-thickness axis of a simulated scan.

A lightweight metadata reader for populating a thickness selector without loading the (large) trace windows: it opens the scansave HDF5 file and returns only /grid/zsave.

Parameters:

path (str) – Path to the .h5/.hdf5 scan file.

Returns:

numpy.ndarray or None – The saved propagation distances (m), strictly increasing from the entrance (0) to the substrate exit (zmax), or None for the older single-thickness format that does not store zsave.

Raises:

KeyError – If the /grid group is missing (not a scansave file).

Return type:

NDArray[float64] | None

croak.io.read_simulated_window_keys(path)

List the trace windows actually stored in a simulated scan.

A lightweight metadata reader for populating a trace-window selector. The canonical trio (SIMULATED_WINDOW_KEYS — Iω_win, Iω_win_reimaged, Iω_full) comes first in that order; any further trace datasets follow — in particular the numbered windows of a multi-aperture scan (Iω_win_2 … with their _reimaged partners), where one propagation was reduced through a whole series of collection holes. Every returned key is loadable via read_simulated_scan()’s window_key.

Parameters:

path (str) – Path to the .h5/.hdf5 scan file.

Returns:

list of str – The available trace-window dataset names (possibly empty).

Return type:

list[str]

croak.io.read_simulated_mask_window(path, window=None)

Read just the collection-window record of a simulated scan.

A lightweight metadata reader, in the same spirit as read_simulated_z_positions(): it opens the scansave file and returns the /grid/window_def_* record without touching the (large) trace windows.

Parameters:

• path (str) – Path to the .h5/.hdf5 scan file.

• window (str or int, optional) – Which collection hole of a multi-aperture scan: the trace-window dataset name ("Iω_win_5", _reimaged accepted) or the bare index; None is the first window (see window_def_prefix()).

Returns:

MaskWindowSpec or None – None when the file does not record a window (or records it only partly).

Raises:

KeyError – If the /grid group is missing (not a scansave file).

Return type:

MaskWindowSpec | None

croak.io.window_def_prefix(window)

Map a trace-window selector to its window_def key prefix.

Multi-aperture scans store one record per collection hole: window_def_* for the first window (Iω_win/Iω_win_reimaged) and window_def_N_* for Iω_win_N (and its _reimaged partner). Accepts the dataset name ("Iω_win_5", "Iω_win_5_reimaged"), the bare index (5), or None/"Iω_win" for the first window.

Parameters:

window (str | int | None)

Return type:

str

croak.io.read_simulated_truth_keys(path)

List the time-domain truth sources a simulated scan provides.

A lightweight metadata reader for populating the truth selector (and choosing its default): it returns the subset of SIMULATED_TRUTH_SOURCES the file stores, in preference order ("beamlet" first). "beamlet" is listed when /grid/It_beamlet or /grid/Ito_beamlet is present; "beamlet_reimaged" when /grid/It_beamlet_reimaged or /grid/Ito_beamlet_reimaged is (pnps files); "source" when /grid/It or /grid/Ito is present.

Parameters:

path (str) – Path to the .h5/.hdf5 scan file.

Returns:

tuple of str – The available truth sources (possibly empty for a malformed file).

Raises:

KeyError – If the /grid group is missing (not a scansave file).

Return type:

tuple[str, …]

croak.io.read_simulated_wavelength_range(path)

Return the full positive-frequency wavelength span of a simulated scan.

A lightweight metadata reader (reads only /grid/ω) for defaulting the load band to the file’s full data range. The simulation grid is an FFT axis spanning positive and negative absolute frequencies; only the positive bins carry a wavelength, and λ = 2πc/ω maps the largest ω to the shortest λ and the smallest positive ω to the longest λ.

Parameters:

path (str) – Path to the .h5/.hdf5 scan file.

Returns:

tuple of float – (lam_min_nm, lam_max_nm) — the shortest and longest wavelengths (nm) of the positive-frequency bins.

Raises:

• KeyError – If the /grid group or its ω axis is missing.

• ValueError – If the grid has no positive-frequency bins.

Return type:

tuple[float, float]

Preprocessing

Cleaning, filtering and regridding of measured FROG data.

The orchestrator load_and_clean() turns a measured FROG trace (on a wavelength axis and a delay/position axis, with an optional independent fundamental spectrum) into a normalised trace on a uniform retrieval grid:

filter (fringe/DC/low-pass) -> spectral window -> regrid (lambda->omega resample + anti-alias) -> marginal-correct -> threshold

Background subtraction, calibration and third-order scaling are not performed here: apply them to the trace before calling load_and_clean(), or fold a per-wavelength response into its scalecurve argument. The standalone helpers background_subtract(), calibration_curve_scale() and third_order_scale() are provided for that.

The delay filter uses a non-uniform FFT (FINUFFT) so the fringe band-stop and dtau_min low-pass remain correctly calibrated even when the delay axis is not uniformly sampled; its stop bands roll off with a steep Planck taper rather than a brick-wall cut. Thresholding is applied once, at the very end, and is off by default.

Arrays use the croak convention (Nlambda, Ndelay) (wavelength on axis 0, delay on axis 1), matching the retrieval trace orientation (Nomega, Ndelay). load_and_clean() returns a TraceData.

class croak.preprocess.TraceData(interaction, grid, delays, trace, omega0_pulse, omega0_trace, Iomega, lam_min, lam_max, lamm_lims, tau_lims, lam_frog, tau_meas, Ifrog_meas, tau_filt, Ifrog_filt, lam_spec, Ilam_spec_meas, Ilam_spec_filt, scalecurve, taper_loss=0.0)

Bases: object

A cleaned, regridded FROG trace ready for retrieval, plus provenance.

Parameters:

• interaction (str)

• grid (Grid)

• delays (NDArray[float64])

• trace (NDArray[float64])

• omega0_pulse (float)

• omega0_trace (float)

• Iomega (NDArray[float64] | None)

• lam_min (float)

• lam_max (float)

• lamm_lims (tuple[float, float] | None)

• tau_lims (tuple[float, float] | None)

• lam_frog (NDArray[float64])

• tau_meas (NDArray[float64])

• Ifrog_meas (NDArray[float64])

• tau_filt (NDArray[float64])

• Ifrog_filt (NDArray[float64])

• lam_spec (NDArray[float64] | None)

• Ilam_spec_meas (NDArray[float64] | None)

• Ilam_spec_filt (NDArray[float64] | None)

• scalecurve (NDArray[float64])

• taper_loss (float)

interaction

"shg"/"sd"/"pg".

Type:

str

grid

Uniform retrieval grid.

Type:

Grid

delays

Regridded delay axis (s).

Type:

numpy.ndarray

trace

Regridded, normalised trace (Nomega, Ndelay).

Type:

numpy.ndarray

omega0_pulse, omega0_trace

Pulse and trace carrier frequencies (rad/s).

Type:

float

Iomega

Independent fundamental spectrum resampled on the grid, or None.

Type:

numpy.ndarray or None

lam_min, lam_max

Retrieval wavelength range (m) — the output grid extent.

Type:

float

lamm_lims

Measurement wavelength window λm (m): where the trace carries signal. Use this (not lam_min/lam_max) for any post-retrieval spectral windowing — see post_filter().

Type:

tuple or None

tau_lims

Measurement delay window τm (s), or None. The companion to lamm_lims for post-retrieval temporal windowing.

Type:

tuple or None

lam_frog, tau_meas, Ifrog_meas, tau_filt, Ifrog_filt

Measured and filtered trace in the wavelength domain (for plotting).

lam_spec, Ilam_spec_meas, Ilam_spec_filt

Independent spectrum (measured/filtered) or None.

scalecurve

Per-wavelength scaling applied during filtering.

Type:

numpy.ndarray

taper_loss

Fraction of the measured trace removed by regrid()’s edge taper (see TAPER_COLLAR_BINS). A band chosen with room to spare leaves this well under 0.01; a large value means lam_min/lam_max are clipping live signal, which both discards data and frees the retrieved field at the affected frequencies. 0.0 when the trace never went through the regrid.

Type:

float

property omega: NDArray[float64]

Centred angular-frequency axis (rad/s).

croak.preprocess.filter_spectrum(lam, Ilam, lam_min, lam_max, *, rolloff=10)

Planck-taper window applied to a spectrum in the frequency domain.

The window passes [lam_min, lam_max] and rolls off symmetrically in angular frequency (the same frequency width on the blue and red edges, rolloff source samples wide), which is the natural choice for a spectrum that will be resampled onto a uniform frequency grid. Normalised to the in-band peak; the input is not modified.

Parameters:

• lam (ArrayLike)

• Ilam (ArrayLike)

• lam_min (float)

• lam_max (float)

• rolloff (int)

Return type:

NDArray[float64]

croak.preprocess.filter_frog_trace(Ifrog, lam, tau, *, tau_lims=None, filter_fringes=False, filter_dc=False, dtau_min=None, defringe=False, defringe_fraction=0.5, on_alias='raise', baseline=False, baseline_smoothness=100.0, baseline_ratio=0.001, baseline_clip=False, baseline_tau_exclude=None, baseline_downsample=1, baseline_skip_below=0.0, lamm_lims=None, subsample=None, scalecurve=None, roll_bins=3)

Filter a measured FROG trace along the delay axis.

Removes interferometric fringes, the DC delay component, and delay structure faster than dtau_min (low-pass) via a non-uniform FFT (see _delay_filter); applies a spectral window (lamm_lims), centres tau=0 on the delay-marginal peak, applies an optional delay window (tau_lims) and per-wavelength scalecurve, then normalises and optionally subsamples.

Fringes can be removed two ways — pick one. filter_fringes is the band- stop notch at the carrier; defringe is the Takeda-FTSI low-pass that keeps the slow DC band and discards the ±f_c sidebands (see defringe_carrier()). They are two solutions to the same problem and are mutually exclusive — enabling both raises ValueError.

Parameters:

• Ifrog (array_like) – Trace (Nlambda, Ndelay).

• lam (array_like) – Wavelength axis (m), length Nlambda.

• tau (array_like) – Delay axis (s), length Ndelay (need not be uniformly sampled).

• filter_fringes (bool or tuple, optional) – Band-stop notch at the fringe carrier. Default False — the filters here all remove real signal along with their target, so none is applied unless asked for. A tuple gives an explicit stop band.

• filter_dc (bool, optional) – Remove the delay-independent component of each wavelength row. Default False; enable only for data with a genuine delay-independent background, since on a clean trace it subtracts real signal.

• dtau_min (float, optional) – Low-pass the delay axis at 1/dtau_min (Hz). None (the default) or 0 applies no low-pass.

• defringe (bool, optional) – Remove fringes by the Takeda-FTSI carrier low-pass instead of the filter_fringes notch. Default False.

• defringe_fraction (float, optional) – De-fringe low-pass cutoff as a fraction of the carrier f_c = c/λ; must lie in (0, 1). Default 0.5. Ignored unless defringe.

• on_alias ({"raise", "warn", "clamp"}, optional) – De-fringe behaviour when the carrier aliases (delay step too coarse); see defringe_carrier(). Default "raise".

• baseline (bool, optional) – Remove the broad, slowly varying leakage background along delay per wavelength row by asymmetric least squares (see arpls_baseline()), applied right after de-fringing. Default False.

• baseline_smoothness (optional) – Stiffness, convergence tolerance and an optional |tau| < baseline_tau_exclude hold-out, forwarded to arpls_baseline(). Ignored unless baseline.

• baseline_ratio (optional) – Stiffness, convergence tolerance and an optional |tau| < baseline_tau_exclude hold-out, forwarded to arpls_baseline(). Ignored unless baseline.

• baseline_tau_exclude (optional) – Stiffness, convergence tolerance and an optional |tau| < baseline_tau_exclude hold-out, forwarded to arpls_baseline(). Ignored unless baseline.

• baseline_downsample (optional) – Speed controls forwarded to arpls_baseline() (coarse-grid solve and near-zero-row skip). Ignored unless baseline.

• baseline_skip_below (optional) – Speed controls forwarded to arpls_baseline() (coarse-grid solve and near-zero-row skip). Ignored unless baseline.

• baseline_clip (bool, optional) – Clip the trace to be non-negative after baseline subtraction. Default False (small negative values are retained).

• roll_bins (int, optional) – Width (in delay-frequency modes) of the steep Planck-taper roll-off on the fringe notch and dtau_min low-pass. Small ⇒ steep edges.

• tau_lims (tuple[float, float] | None)

• lamm_lims (tuple[float, float] | None)

• subsample (int | None)

• scalecurve (ArrayLike | None)

Returns:

tuple – (tau_filtered, Ifrog_filtered).

Raises:

ValueError – If both filter_fringes and defringe are enabled, or if defringe_fraction is outside (0, 1).

Return type:

tuple[NDArray[float64], NDArray[float64]]

Notes

No thresholding is applied here; small negative values produced by the filters are preserved, avoiding additional censoring of the noise. load_and_clean() applies an optional threshold once, at the very end.

croak.preprocess.defringe_carrier(Ifrog, lam, tau, *, carrier_fraction=0.5, dtau_min=None, roll_bins=3, on_alias='raise', eps=1e-09)

Remove interferometric fringes by keeping the DC band (Takeda FTSI).

A measured TG/PG-FROG trace decomposes per wavelength row as

\[I(\lambda,\tau) = a(\lambda,\tau) + c(\lambda,\tau)\,e^{+i\omega(\lambda)\tau} + c^{*}(\lambda,\tau)\,e^{-i\omega(\lambda)\tau},\]

where the slow term \(a = |E_s|^2 + |E_b|^2\) is the fringe-free magnitude trace and the cross term (the fringes) oscillates at the optical carrier \(\omega(\lambda)=2\pi c/\lambda\), i.e. ordinary fringe frequency \(f_c(\lambda)=c/\lambda\) — known a priori from the wavelength calibration, not fitted. The fringes are diagonal because \(f_c\) tilts with \(\lambda\).

This keeps only the DC band with a smooth Planck-taper low-pass whose cutoff is carrier_fraction * f_c(lambda) per row — below the carrier and above the trace’s own slow delay-bandwidth — and discards the \(\pm f_c\) sidebands, recovering a smooth, fringe-free \(a(\lambda,\tau)\). Smooth windows avoid the Gibbs ringing of a brick-wall notch. Separation requires the signal and fringe sidebands to occupy distinct delay-frequency bands. The transform is a non-uniform FFT, so the cutoff stays correctly calibrated on a jittered or unevenly stepped delay axis.

Unlike DC removal (filter_dc), this preserves the slow signal content including its zero-frequency mean and low-frequency wings; the broad leakage background \(|E_b|^2\) survives here and is removed by a separate baseline step.

Parameters:

• Ifrog (array_like) – Measured trace (Nlambda, Ndelay) (delay on axis 1).

• lam (array_like) – Wavelength axis (m), length Nlambda, used to compute the carrier.

• tau (array_like) – Delay axis (s), length Ndelay (need not be uniformly sampled).

• carrier_fraction (float, optional) – Cutoff as a fraction of the per-row carrier f_c = c/lambda. Must satisfy 0 < carrier_fraction < 1; default 0.5 (mid-way between the trace’s slow DC band and the fringe carrier).

• dtau_min (float, optional) – Optional absolute low-pass at 1/dtau_min (Hz) applied in addition; the more restrictive of the two cutoffs wins per row. None (default) uses only the carrier-relative cutoff.

• roll_bins (int, optional) – Width (delay-frequency modes) of the Planck-taper roll-off. Small ⇒ steep edges, less passband smoothing, while still avoiding brick-wall ringing. Default 3.

• on_alias ({"raise", "warn", "clamp"}, optional) – Behaviour when the carrier f_c(lambda) exceeds the delay Nyquist 1/(2*dtau) (the fringe aliases) at a signal-bearing wavelength. Default "raise".

• eps (float, optional) – NUFFT accuracy. Default 1e-9.

Returns:

numpy.ndarray – The smooth fringe-free trace a(lambda, tau), shape (Nlambda, Ndelay).

Raises:

ValueError – If carrier_fraction is outside (0, 1), or if on_alias="raise" and the carrier aliases at a signal-bearing wavelength.

Return type:

NDArray[float64]

Examples

>>> import numpy as np
>>> from croak.constants import C
>>> tau = np.linspace(-60e-15, 60e-15, 256)
>>> lam = np.array([800e-9])
>>> a = np.exp(-0.5 * (tau / 8e-15) ** 2)            # slow envelope
>>> fringe = a * np.cos(2 * np.pi * (C / lam[0]) * tau)
>>> clean = defringe_carrier((a + fringe)[None, :], lam, tau)
>>> bool(np.allclose(clean[0, 120:136], a[120:136], atol=2e-2))
True

croak.preprocess.arpls_baseline(Ifrog, *, smoothness=100.0, ratio=0.001, max_iter=50, tau=None, tau_exclude=None, downsample=1, skip_below=0.0)

Estimate a smooth, asymmetric per-row baseline (arPLS, Baek et al. 2015).

For each wavelength row y of a trace (Nlambda, Ndelay) (intensity vs delay), fit a free-form smooth baseline z that rides under the broad, slowly varying leakage background \(|E_b|^2\) while downweighting the compact signal peak near tau = 0. Where background and signal overlap the baseline can still remove signal; suggested_tau_exclude() gives a hold-out that prevents this. The baseline minimises

\[F(z) = \sum_i w_i (y_i - z_i)^2 + \lambda \sum_i (\Delta^2 z_i)^2,\]

where \(\Delta^2\) is the discrete second difference (a curvature penalty) and \(\lambda\) (smoothness) sets the stiffness. The normal equations \((W + \lambda D^{\mathsf{T}}D)\,z = W y\) are solved repeatedly; after each solve the weights are updated from the residual \(d = y - z\) using only its negative part (assumed background/noise): with \(m=\mathrm{mean}(d<0)\), \(s=\mathrm{std}(d<0)\),

\[w_i = \frac{1}{1 + \exp\!\big(2\,(d_i - (2s - m))/s\big)},\]

so points well above the baseline (signal) are driven toward zero weight and the curve settles onto the lower envelope.

Because z is a free-form smooth curve fit to that lower envelope — no symmetry or shape prior — it follows a background that is asymmetric about tau = 0 (e.g. a broad leakage shelf strong at negative delay and weak at positive delay, as in DUV TG-FROG) just as readily as a flat one (the zero-curvature limit). Subtract the result to clean the trace: Ifrog - arpls_baseline(Ifrog).

Parameters:

• Ifrog (array_like) – Trace (Nlambda, Ndelay) (delay on axis 1), or a single row (Ndelay,). Rows are treated independently.

• smoothness (float, optional) – Curvature penalty \(\lambda\). Larger ⇒ stiffer (flatter) baseline; smaller ⇒ the baseline bends more and risks eating the signal peak. The penalty is internally normalised by (256/Ndelay)**4 and each row by its peak, so this value is roughly portable across delay-sample counts and trace amplitudes. Default 1e2; useful range ~1–1e5. Must be positive.

• ratio (float, optional) – Convergence tolerance on the relative change of the baseline between iterations, ||z_new - z|| / ||z_new||. Default 1e-3. Must be positive. (Converging on the baseline rather than the weights reaches the same fixed point in far fewer iterations.)

• max_iter (int, optional) – Maximum reweighting iterations per row. Default 50. Must be >= 1.

• tau (array_like, optional) – Delay axis (s), length Ndelay. Required only when tau_exclude is given (to locate |tau| < tau_exclude).

• tau_exclude (float, optional) – Hold out the delay window |tau| < tau_exclude (s) from the fit — those points keep w = 0 every iteration and are excluded from the residual statistics, so the compact signal peak never pulls the baseline up. None (default) lets the reweighting find the peak on its own. Requires tau.

• downsample (int, optional) – Solve on every downsample-th delay (block-averaging the row onto the coarse grid) and linearly interpolate the baseline back to full resolution. The baseline is smooth by construction, so a coarse grid captures it at a fraction of the cost; smoothness is grid-normalised so the same value transfers. Default 1 (exact, full-resolution). Must be >= 1. A value of 4 is a good speed/accuracy trade for interactive use (≈2.5 % baseline difference; the test test_arpls_baseline_downsample_matches_full bounds it below 5 % of the baseline peak); use 1 for a final calculation.

• skip_below (float, optional) – Skip rows whose peak is below skip_below times the global peak, returning a zero baseline for them (no leakage background is fitted where there is no signal). Default 0.0 (process every row). Range [0, 1).

Returns:

numpy.ndarray – The estimated baseline, same shape as Ifrog.

Raises:

ValueError – If smoothness <= 0, ratio <= 0, max_iter < 1, downsample < 1, skip_below is outside [0, 1), the input is not 1-D or 2-D, the delay axis has fewer than 3 samples (after downsampling), tau_exclude is given without a matching-length tau, or tau_exclude covers the whole delay axis.

Return type:

NDArray[float64]

Notes

The broad leakage background survives defringe_carrier() (which removes only the interferometric cross term); this is the companion step that removes it. The intended order is de-fringe → arpls_baseline subtract → optional clip. No clipping is done here — small negatives in the subtracted trace are retained to avoid additional noise censoring. See Cleaning and regridding measured data for the statistical assumptions and the effects of clipping.

References

Baek, Park, Ahn & Choo, “Baseline correction using asymmetrically reweighted penalized least squares smoothing”, Analyst 140, 250 (2015).

Examples

>>> import numpy as np
>>> tau = np.linspace(-300e-15, 300e-15, 600)
>>> shelf = 0.4 * np.exp(np.minimum(-tau / 150e-15, 0.0))  # high at negative delay
>>> peak = np.exp(-0.5 * (tau / 8e-15) ** 2)               # compact signal
>>> z = arpls_baseline(shelf + peak, smoothness=1e2)
>>> wings = np.abs(tau) > 150e-15                          # away from the peak
>>> bool(np.allclose(z[wings], shelf[wings], atol=5e-2))
True

croak.preprocess.suggested_tau_exclude(Ifrog, tau, *, factor=2.0, max_fraction=0.4)

Suggest an arpls_baseline() tau_exclude from the delay marginal.

Without a hold-out the arPLS baseline bends up into the signal at τ≈0 and removes part of the peak (“peak suck-up”). How much depends on how well each wavelength row constrains the fit, so it varies sharply from row to row: on one of our measured deep-UV TG-FROG traces the worst rows kept 23–36 % of their peak while their immediate neighbours kept ~70 %, which reads as a thin stripe at one wavelength that appears only when baseline removal is on.

Holding out the delay range where the signal lives fixes it. The suggested half-width is factor times the delay-marginal FWHM, limited by max_fraction of the delay half-span. Twice the FWHM worked best on that trace: it took the worst row from 0.23 to 0.87 of its peak and left no row below 0.87 (thirteen rows were below 0.8 unguarded), and it was faster, because the fit settles sooner. Once (too narrow) and three times (too little wing left to fit) were both worse.

Parameters:

• Ifrog (array_like) – Trace, (Nlambda, Ndelay), or a single delay marginal (Ndelay,). Any stage of processing will do: the width is insensitive to it.

• tau (array_like) – Delay axis (s), matching the last dimension of Ifrog.

• factor (float, optional) – Multiple of the marginal FWHM to hold out (default 2.0).

• max_fraction (float, optional) – Cap, as a fraction of the delay half-span (default 0.4), so the hold-out can never approach the whole axis — which arpls_baseline() rejects.

Returns:

float – Hold-out half-width in seconds, for arpls_baseline(tau_exclude=...).

Raises:

ValueError – If the shapes disagree, or the delay marginal has no discernible peak (its half maximum is never crossed), so there is nothing to hold out.

Return type:

float

Examples

>>> tau = np.linspace(-50e-15, 50e-15, 401)
>>> marginal = np.exp(-((tau / 4e-15) ** 2))  # 6.66 fs FWHM
>>> round(suggested_tau_exclude(marginal, tau) / 1e-15, 2)
13.32

croak.preprocess.resample_spectrum(lam, Ilam, omega, *, prefilter=True)

Resample a spectrum from a wavelength grid onto an (absolute) frequency grid.

Applies the Jacobian Ilam * lam**2 for the domain change, interpolates with a cubic spline, anti-aliases by low-passing in the frequency domain when the target grid is coarser than the source, and tapers the edges.

Parameters:

• lam (array_like) – Source wavelength axis (m).

• Ilam (array_like) – Source spectral intensity.

• omega (array_like) – Target absolute angular-frequency grid (rad/s), e.g. g.omega + omega0.

• prefilter (bool, optional) – Apply the anti-alias low-pass before downsampling.

Return type:

NDArray[float64]

croak.preprocess.regrid(tau, lam, trace, lam_spec, Ilam_spec, interaction, *, lam_min, lam_max, trange=None, resample_t=False, marginal_correct=False, prefilter=True, taper_warn_level=0.01)

Regrid a filtered FROG trace onto a uniform frequency/delay grid.

Parameters:

• tau (array_like) – Filtered delay axis (s).

• lam (array_like) – Wavelength axis (m).

• trace (array_like) – Filtered trace (Nlambda, Ndelay).

• lam_spec (array_like or None) – Independent fundamental spectrum, or None.

• Ilam_spec (array_like or None) – Independent fundamental spectrum, or None.

• interaction (str) – "shg"/"sd"/"pg" (sets the trace carrier scaling).

• lam_min (float) – Fundamental-pulse retrieval wavelength range (m), also for SHG.

• lam_max (float) – Fundamental-pulse retrieval wavelength range (m), also for SHG.

• trange (float, optional) – Output temporal window (s); default twice the delay span.

• resample_t (bool, optional) – Resample the delay axis onto the grid time axis.

• marginal_correct (bool, optional) – Apply SHG autoconvolution marginal correction (needs a spectrum).

• prefilter (bool, optional) – Anti-alias before downsampling in resample_spectrum().

• taper_warn_level (float, optional) – Warn when the edge taper removes more than this fraction of the measured trace, i.e. when lam_min/lam_max are too tight. Default 0.01 (1 %); pass 0.0 to silence the check.

Returns:

tuple – (tau, grid, omega0_pulse, omega0_trace, trace_regridded, Iomega, taper_loss).

Warns:

UserWarning – If the spectral taper removes more than taper_warn_level of the trace.

croak.preprocess.background_subtract(intensity, background, *, clamp=False)

Subtract a background; optionally clamp negatives to zero.

Parameters:

• intensity (ArrayLike)

• background (ArrayLike)

• clamp (bool)

Return type:

NDArray[float64]

croak.preprocess.calibration_curve_scale(lam, calib_lam_nm, calib_curve)

Per-wavelength calibration factors from a (λ_nm, scale) curve.

The curve is spline-interpolated (constant extrapolation) onto lam (m).

Parameters:

• lam (ArrayLike)

• calib_lam_nm (ArrayLike)

• calib_curve (ArrayLike)

Return type:

NDArray[float64]

croak.preprocess.third_order_scale(lam, exponent=4.0)

TG-FROG third-order efficiency scaling \((\lambda/\mu m)^{exp}\).

Parameters:

• lam (ArrayLike)

• exponent (float)

Return type:

NDArray[float64]

croak.preprocess.omega_to_lambda_density(intensity_omega, lam)

Convert an angular-frequency spectral density to a wavelength density.

A measured (spectrometer) FROG trace is a wavelength density \(I_\lambda\), and regrid() multiplies it by \(\lambda^2\) to recover the angular-frequency density \(I_\omega \propto I_\lambda\,\lambda^2\) the retrieval forward model produces. Numerically simulated traces are already \(I_\omega\) (uniform in \(\omega\)), so they must be divided by \(\lambda^2\) before entering the wavelength pipeline; the subsequent \(\times\lambda^2\) then returns the original \(I_\omega\) undistorted.

The Jacobian is \(|\mathrm{d}\omega/\mathrm{d}\lambda| = 2\pi c/\lambda^2 \propto 1/\lambda^2\); the constant \(2\pi c\) is irrelevant because the trace is renormalised downstream, so only the \(1/\lambda^2\) factor is applied.

Parameters:

• intensity_omega (array_like) – Spectral intensity as an angular-frequency density. A 1-D spectrum (Nlambda,) or a 2-D trace (Nlambda, Ndelay).

• lam (array_like) – Wavelength axis (m), aligned with axis 0 of intensity_omega.

Returns:

numpy.ndarray – The corresponding wavelength density (same shape as intensity_omega).

Return type:

NDArray[float64]

croak.preprocess.marginal_peak_delay(delays, marginal)

Locate the peak of a delay marginal to sub-sample precision.

Time zero is not known a priori on a measured scan; it is taken to be where the delay marginal peaks. Reading that off as the largest sample quantises the answer to the delay step, leaving up to half a step of offset, which a retrieval then shows as a red/blue dipole straddling the trace in the residual, because the measurement and the model sit at slightly different delays.

This refines the estimate by fitting a parabola through the peak sample and its two neighbours and returning the vertex, which is exact for a locally quadratic peak and costs nothing: the delay axis is only translated, never resampled. The marginal peak need not coincide with gate–probe coincidence for an asymmetric trace; use the same origin convention for data and model (see croak.delay_origin).

The fit is done in local coordinates and the vertex clamped to the bracketing samples, so noise cannot throw it outside the interval its three points span. A peak sitting on either end of the axis, or a fit that does not curve downward, falls back to the sample position.

Parameters:

• delays (array_like) – Delay axis (s). Need not be uniformly sampled.

• marginal (array_like) – Delay marginal (the trace summed over frequency), same length.

Returns:

float – Estimated delay (s) of the marginal’s peak.

Return type:

float

Examples

>>> tau = np.linspace(-10e-15, 10e-15, 21)          # 1 fs steps
>>> marg = np.exp(-(((tau - 0.4e-15) / 3e-15) ** 2))  # true peak at 0.4 fs
>>> round(float(marginal_peak_delay(tau, marg)) / 1e-15, 3)  # parabolic estimate
0.392

croak.preprocess.uniform_core_indices(x, rtol=0.001)

Find the longest contiguous uniformly spaced run of x.

Several simulated campaigns extend a uniform delay core with coarser wing points (e.g. a 200-point ±25 fs core at 0.2513 fs plus 1 fs-step wings to ±40 fs for the Raman wake). The wings are valuable data, but the smearing kernel’s delay-axis convolution is evaluated on a uniform delay grid, so a retrieval that models smearing must restrict itself to the core. This finds it: the maximal contiguous stretch whose consecutive spacings agree with the stretch’s first spacing to rtol (relative).

Parameters:

• x (array_like) – Sample positions, sorted ascending. Not modified.

• rtol (float, optional) – Relative tolerance on the spacing within a run.

Returns:

numpy.ndarray of intp – Indices (into x) of the selected run, contiguous and ascending. The whole axis when it is already uniform; ties go to the earliest run.

Return type:

NDArray[int64]

Examples

>>> import numpy as np
>>> core = np.arange(0.0, 10.0, 0.5)
>>> wings = np.array([-3.0, -2.0, -1.0])
>>> x = np.concatenate([wings, core, [11.0, 13.0, 15.0]])
>>> idx = uniform_core_indices(x)
>>> np.allclose(x[idx], core)
True

croak.preprocess.load_and_clean(Ifrog, lam_frog, scanaxis, interaction, *, lam_min, lam_max, lamm_lims=None, input_unit='position', trange=None, tau_crop=None, tau_lims=None, filter_fringes=False, filter_dc=False, dtau_min=None, defringe=False, defringe_fraction=0.5, on_alias='raise', baseline=False, baseline_smoothness=100.0, baseline_ratio=0.001, baseline_clip=False, baseline_tau_exclude=None, baseline_downsample=1, baseline_skip_below=0.0, resample_t=False, subsample=None, scalecurve=None, threshold=None, roll_bins=3, marginal_correct=False, prefilter=True, taper_warn_level=0.01, lam_spec=None, Ilam_spec=None)

Clean and regrid a measured FROG trace from pre-loaded arrays.

Parameters:

• Ifrog (array_like) – Measured trace (Nlambda, Ndelay), a wavelength density \(I_\lambda\) (see Notes).

• lam_frog (array_like) – Wavelength axis (m).

• scanaxis (array_like) – Stage position (m) if input_unit="position" (converted via τ = 2z/c) or delay (s) if input_unit="delay".

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

• lam_min (float) – Fundamental-pulse retrieval wavelength range (m), also for SHG.

• lam_max (float) – Fundamental-pulse retrieval wavelength range (m), also for SHG.

• lamm_lims (tuple, optional) – Useful measurement wavelength window λm (m); defaults to the full extent. Applied as a spectral window before regridding, and stored on the returned TraceData for optional post-retrieval windowing.

• input_unit ({"position", "delay"}, optional) – Whether scanaxis is a stage position (m) or a delay (s). Default "position", converted with τ = 2z/c.

• trange (float, optional) – Time-window length (s) of the retrieval grid. None (default) uses twice the delay span. Together with lam_min/lam_max this fixes the grid: the number of points is trange × the spectral bandwidth, while the time step is set by the wavelength window alone (dt ≈ 1/Δν). Widen the window to resolve the pulse more finely; lengthen trange to get more points at the same dt.

• tau_crop (tuple, optional) – (τ_min, τ_max) in s: discard delays outside this range after filtering, regridding and optional thresholding. None (default) keeps the whole scan.

• tau_lims (tuple, optional) – Useful measurement delay window τm (s), applied as a delay taper during filtering and stored on the returned TraceData for optional post-retrieval windowing.

• filter_fringes (bool or tuple, optional) – Notch out interferometric fringes in the delay-frequency domain. Off by default — every filter here removes real signal along with whatever it targets, so none is applied unless asked for. Turn this on for an interferometric measurement that actually carries fringes; a tuple gives an explicit stop band. Mutually exclusive with defringe, which is the better method when the carrier is known a priori (it assumes separated delay-frequency sidebands).

• filter_dc (bool, optional) – Remove the delay-independent (zero delay-frequency) component of each wavelength row. Default False. Only switch it on when the measurement genuinely carries a delay-independent background: on a clean trace it subtracts real signal, which shows up as large negative excursions in the filtered trace, a percent-level floor on the achievable trace error and bias in the retrieved field.

• dtau_min (float, optional) – Low-pass the delay axis at 1/dtau_min (Hz), suppressing structure faster than the delay stage can resolve. None (the default) or 0 applies no low-pass.

• defringe (bool, optional) – Remove fringes by the Takeda-FTSI carrier low-pass (see defringe_carrier()) instead of the filter_fringes notch — the two are mutually exclusive. Default False.

• defringe_fraction (float, optional) – De-fringe cutoff as a fraction of the per-row carrier f_c = c/λ, in (0, 1). Default 0.5, midway between the trace’s own slow delay bandwidth and the fringe carrier.

• on_alias ({"raise", "warn", "clamp"}, optional) – What to do when the carrier exceeds the delay Nyquist frequency, i.e. the scan is too coarsely stepped to resolve the fringes. Default "raise".

• baseline (bool, optional) – Remove the broad, slowly varying leakage background |E_b|² along delay (see arpls_baseline()), applied right after de-fringing. Default False.

• baseline_smoothness (float, optional) – arPLS curvature penalty. Larger is stiffer (flatter); too small and the baseline bends up into the signal peak. Default 1e2.

• baseline_ratio (float, optional) – arPLS convergence tolerance on the relative change of the baseline between iterations. Default 1e-3.

• baseline_clip (bool, optional) – Clip the baseline-subtracted trace at zero. Default False, which keeps negatives visible rather than hiding a bad subtraction.

• baseline_tau_exclude (float, optional) – Hold out |τ| < this value (s) when fitting the baseline, so a strong τ≈0 peak cannot drag it up. None (default) holds out nothing; see suggested_tau_exclude() for a sensible width.

• baseline_downsample (int, optional) – Solve the baseline on every n-th delay and interpolate back — a speed/accuracy trade for interactive use. Default 1 (exact).

• baseline_skip_below (float, optional) – Skip baseline fitting on rows whose peak is below this fraction of the trace maximum. Default 0.0 (fit every row).

• resample_t (bool, optional) – Resample the delay axis onto the retrieval grid’s uniform time step. Default False, which keeps the measured delays — the forward model evaluates at arbitrary delays, so resampling is rarely needed.

• subsample (int, optional) – Keep every subsample-th delay, for thinning an oversampled scan. None (default) keeps them all.

• scalecurve (array_like, optional) – Per-wavelength multiplicative response correction (a calibration curve), applied before filtering. None (default) applies none.

• threshold (float, optional) – Noise floor relative to the peak. None (the default) keeps all values; otherwise applied once, at the very end (after regridding). Thresholding biases a maximum-likelihood fit — see the preprocessing guide — so prefer leaving it off and letting the weighting handle the noise.

• roll_bins (int, optional) – Steepness of the delay-filter stop-band roll-offs; see filter_frog_trace().

• marginal_correct (bool, optional) – Rescale the trace so its frequency marginal matches the one predicted from the independent spectrum (SHG only; see croak.marginal_checks). Default False.

• prefilter (bool, optional) – Apply the anti-alias low-pass before downsampling in regrid(). Default True, and it should normally stay on: without it, spectral structure finer than the retrieval grid folds back as aliasing. It is a sharp cut in the Fourier domain, so it rings — the filtered trace carries small negative excursions, which croak keeps rather than clamping so they stay visible.

• taper_warn_level (float, optional) – Warn when regrid’s edge taper removes more than this fraction of the measured trace, i.e. when lam_min/lam_max clip live signal. Default 0.01; pass 0.0 to silence it. See TraceData.taper_loss.

• lam_spec (array_like, optional) – Independent fundamental spectrum, on its own wavelength axis. Used as a spectral reference and carried on the result for plotting; it can initialise the pulse and provides the target for spectral regularisation, which the high-level retrieval pipeline enables by default when the spectrum is present.

• Ilam_spec (array_like, optional) – Independent fundamental spectrum, on its own wavelength axis. Used as a spectral reference and carried on the result for plotting; it can initialise the pulse and provides the target for spectral regularisation, which the high-level retrieval pipeline enables by default when the spectrum is present.

Returns:

TraceData

Return type:

TraceData

Notes

The trace must be a wavelength density \(I_\lambda\), which is what a spectrometer records. regrid() multiplies by \(\lambda^2\) to recover the \(I_\omega\) the forward model produces. A numerically synthesised trace is already \(I_\omega\), so pass it through omega_to_lambda_density() first, or the conversion is applied twice.

croak.preprocess.TAPER_COLLAR_BINS = 10

Width, in grid bins, of the Planck-taper roll-off that regrid() applies at each end of the retrieval grid’s frequency axis. The taper is not optional: the forward model FFTs the trace on a periodic grid, so the resampled trace has to reach zero at the grid boundary or the wrap-around discontinuity rings back into the fit. Its cost is that the outermost TAPER_COLLAR_BINS rows at each end carry little or no measurement, which leaves the field at those frequencies almost unconstrained. Check TraceData.taper_loss (enlarge the grid if appreciable signal is removed) and, on the retrieved field, croak.processing.edge_energy_fraction().

Pipeline

High-level retrieval from cleaned TraceData.

One call that takes a cleaned trace all the way to a retrieved field: it builds a physically sensible initial guess from the trace (and any independent spectrum), selects the solver and pulse mode, runs the retrieval, and returns a RetrievalResult carrying the carrier frequency.

croak.pipeline.retrieve_from_tracedata(td, *, algorithm='warm-lbfgs', full=True, R_omega=False, guess='auto', perfect_fwhm=1e-14, truth=None, material=None, thickness=0.0, npoints=20, smearing=None, focal=None, depth_weight=None, depth_transverse=False, reg_amp=None, reg_phase=0.0, reg_spectrum=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, phase_basis='pointwise', n_nodes=20, strategy='cma', popsize=0, std_init=0.5, seed=None, maxiters=100, alpha=0.25, stall_patience=5, reltol=0.0001, abstol=1e-08, rng=None, callback=None, verbose=False)

Retrieve a pulse from a cleaned TraceData.

Parameters:

• td (TraceData) – Cleaned, regridded trace.

• algorithm ({"warm-lbfgs", "copra", "copra-jax", "lbfgs", "lbfgs-hand", "lbfgs-ad", "lbfgs-optx", "lm", "lm-optx", "cma-es"}, optional) – Retrieval algorithm — any key in croak.retrieve.ALGORITHMS. Default "warm-lbfgs". "warm-lbfgs" (the default) is "lbfgs-ad" warm-started from COPRA’s local sweep; "copra" and "lbfgs" use hand-derived analytic gradients (pure NumPy); "copra-jax" and "lbfgs-hand" run the same hand-derived gradients jitted in JAX (no AD); "lbfgs-ad" and "lm" use JAX automatic differentiation. The "-optx" variants are JAX-native (Optimistix) twins of "lbfgs-ad" and "lm" that run the whole optimisation loop on-device, and "cma-es" is the derivative-free evolutionary strategy (evosax).

• full (bool, optional) – Retrieve amplitude and phase (True) or phase only (False, L-BFGS with the guess amplitude held fixed).

• R_omega (bool, optional) – Per-frequency adaptive scaling (supported by every algorithm).

• guess (str, array_like or Pulse, optional) – "auto"/"tl"/"perfect"/"random"/"truth" (see initial_guess()) or an explicit spectrum/pulse passed straight through.

• perfect_fwhm (float, optional) – FWHM (s) for the "perfect" guess.

• truth (TruthPulse, optional) – Known complex pulse used when guess="truth".

• material (str | None) – Dispersive-slab parameters (dispersive propagation through a material).

• thickness (float) – Dispersive-slab parameters (dispersive propagation through a material).

• npoints (int) – Dispersive-slab parameters (dispersive propagation through a material).

• smearing (SmearingKernel or None, optional) – Geometric time-smearing kernel of a non-collinear BOXCARS geometry (None = off). Only the autodiff solvers can model it; asking for it with any other algorithm raises rather than silently dropping it. See croak.smearing and Modelling geometric time smearing.

• focal (FocalMixture or None, optional) – Explicit chromatic focal-plane mixture, optionally carrying a collection aperture (None = off). Replaces the reduced smearing kernel with the node-resolved mixture of croak.focal; like smearing it is a physical model term, so requesting it with a solver that cannot represent it raises rather than silently dropping it. See Modelling the collection aperture.

• reg_amp (float or None, optional) – Spectral-amplitude smoothness weight. None (default) resolves to the package default weight for the solver’s objective family — 0.03 for the gradient solvers, 3e-5 for the LM residual form (see croak.solver.REG_DEFAULTS). Pass 0 to disable.

• reg_phase (float, optional) – Spectral-phase smoothness weight (default 0).

• reg_spectrum (float or None, optional) – Spectral-match regularisation weight (full mode): pulls the retrieved spectral amplitude towards the measured independent spectrum td.Iomega. None (default) resolves to the package default weight (0.01 gradient / 1e-5 LM) and silently does nothing when no independent spectrum is available; an explicit reg_spectrum > 0 with no spectrum is ignored with a warning. Pass 0 to disable.

• reg_time (float, optional) – Temporal (pedestal) regularisation weight (LBFGS-AD / cma-es): penalises the fraction of E(t) energy outside time_window. A self-consistent alternative to the temporal post-filter — it keeps spectrum and field a Fourier pair, so introduces no spectral artefacts, though the window is a prior on the pulse’s support. Active in all modes.

• time_window (tuple of float, optional) – (lo, hi) window in seconds for reg_time; defaults to the measurement delay range (td.delays.min(), td.delays.max()).

• tau0 (float, optional) – Starting centres for the two fitted extras that have no nominal argument of their own: the delay-zero offset (s, default 0) and the multiplier on the smearing-kernel widths (default 1.0 = the kernel as supplied). Pass the previous fit’s values to warm start a re-run, as the GUI’s “reuse previous result” does; the slab thickness is warm-started through thickness itself. Each is held at its centre when the matching fit_* flag is off.

• smear_scale (float, optional) – Starting centres for the two fitted extras that have no nominal argument of their own: the delay-zero offset (s, default 0) and the multiplier on the smearing-kernel widths (default 1.0 = the kernel as supplied). Pass the previous fit’s values to warm start a re-run, as the GUI’s “reuse previous result” does; the slab thickness is warm-started through thickness itself. Each is held at its centre when the matching fit_* flag is off.

• fit_thickness (bool, optional) – Fit the dispersive-slab thickness and/or a delay-zero offset as extra parameters alongside the pulse (lbfgs-ad / warm-lbfgs / lm / lm-optx / cma-es; fit_thickness needs a dispersive slab, PG/SD). See Fitting the medium thickness and delay-zero.

• fit_tau0 (bool, optional) – Fit the dispersive-slab thickness and/or a delay-zero offset as extra parameters alongside the pulse (lbfgs-ad / warm-lbfgs / lm / lm-optx / cma-es; fit_thickness needs a dispersive slab, PG/SD). See Fitting the medium thickness and delay-zero.

• fit_smearing (bool, optional) – Fit the geometric-smearing strength (one dimensionless multiplier on the kernel widths) alongside the pulse. Requires smearing. 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. Requires fit_smearing; smear_scale_delta (default 1.0) is the delay channel’s starting centre. See Modelling geometric time smearing.

• depth_weight (array_like or None, optional) – Complex generation-envelope weights, one per depth-quadrature node (lbfgs-ad and warm-lbfgs): the transverse-geometry factor a 1D depth integral cannot see. See croak.forward_jax.make_param_trace_fn().

• depth_transverse (bool, optional) – Model the transverse decoherence of the depth integral (lbfgs-ad and warm-lbfgs): a parameter-free per-depth transverse phase inside the collection transform. Requires focal with a collection and a dispersive slab. See croak.forward_jax.make_param_trace_fn().

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

• delay_origin ({"coincidence", "marginal_peak"}, optional) – The model’s delay-zero convention (AD solvers): gate–probe coincidence, or the model trace re-centred on its own marginal peak like the data. See croak.delay_origin.

• tau0_bound (float, optional) – |tau0| <= tau0_bound (s) for a fitted delay offset; 0 = unbounded.

• maxiters (int, optional) – Maximum iterations.

• reltol (float, optional) – Convergence tolerances forwarded to every solver (defaults 1e-4 / 1e-8). For the gradient/LM solvers these are the optimiser function tolerances; for COPRA they enable relative-plateau / absolute-target early stopping (off by default in the bare COPRA class).

• abstol (float, optional) – Convergence tolerances forwarded to every solver (defaults 1e-4 / 1e-8). For the gradient/LM solvers these are the optimiser function tolerances; for COPRA they enable relative-plateau / absolute-target early stopping (off by default in the bare COPRA class).

• rng (numpy.random.Generator, optional) – Generator (random guess / COPRA permutation).

• callback (callable, optional) – callback(iteration, R, best_R, snapshot=...) progress hook; the optional snapshot is a zero-argument in-progress-result builder (or None) for a live full-plot preview.

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

• smear_scale_delta (float)

• phase_basis (str)

• n_nodes (int)

• strategy (str)

• popsize (int)

• std_init (float)

• seed (int | None)

• alpha (float)

• stall_patience (int)

Returns:

RetrievalResult

Return type:

RetrievalResult

croak.pipeline.initial_guess(td, *, mode='auto', perfect_fwhm=1e-14, truth=None, rng=None)

Build an initial complex spectrum from a cleaned trace.

Parameters:

• td (TraceData) – Cleaned trace.

• mode ({"auto", "tl", "perfect", "random", "truth"}) –

  The first four start from a spectral amplitude and differ in the phase, which is what selects the basin the optimiser falls into:

  "tl"

  the measured amplitude with flat phase — the transform limit of the spectrum actually measured. The natural experimental guess: it uses the real amplitude and assumes nothing about the duration.

  "perfect"

  a synthetic transform-limited Gaussian of perfect_fwhm at the spectral peak. Note this discards the measured amplitude and encodes an assumed duration; "tl" is the better-justified version of the same flat-phase start and performs equivalently.

  "auto"

  the measured amplitude with a centroid-delay phase, which already encodes the group-delay slope and so starts a chirped pulse close to its answer.

  "random"

  the measured amplitude with a uniform random phase.

  "truth"

  the known complex spectrum carried by truth, aligned to the retrieval grid. ModelPNPS native-grid data are copied exactly; regridded data interpolate amplitude and unwrapped phase.

  The flat-phase starts ("tl", "perfect") suit near-transform- limited pulses; a strongly chirped pulse started flat must build its quadratic phase from nothing and tends to retain temporal structure, where "auto" does not.

• perfect_fwhm (float, optional) – Temporal FWHM for the "perfect" seed (s).

• truth (TruthPulse, optional) – Known pulse carrying a complex spectrum; required for mode="truth".

• rng (numpy.random.Generator, optional) – Generator for "random".

Returns:

numpy.ndarray – Complex spectrum (centred order).

Return type:

NDArray[complex128]

Processing

Post-retrieval processing for analysis and display.

Turns a RetrievalResult into a ProcessedResult carrying everything the plots need: oversampled temporal intensity/phase (linear ramp removed), the transform-limited reference, wavelength-domain spectral intensity/phase with fitted GDD/TOD, trace marginals, residuals and a Gabor spectrogram helper.

croak.processing.marginals(trace)

Normalised frequency and delay marginals of a trace (Nomega, Ndelay).

Returns (omega_marginal, delay_marginal); each is floored at zero, has its minimum removed and is scaled to unit maximum.

Parameters:

trace (ArrayLike)

Return type:

tuple[NDArray[float64], NDArray[float64]]

croak.processing.residuals(measured, retrieved, *, weights=None)

Weighted residual (measured - mu*retrieved) with optimal scale mu.

Parameters:

• measured (ArrayLike)

• retrieved (ArrayLike)

• weights (ArrayLike | None)

Return type:

NDArray[float64]

croak.processing.oversample(t, field, factor=8)

Smoothly oversample a complex field by zero-padding in the frequency domain.

Parameters:

• t (array_like) – Time axis (uniform).

• field (array_like) – Complex field samples.

• factor (int, optional) – Oversampling factor (default 8).

Returns:

tuple – (t_oversampled, field_oversampled).

Return type:

tuple[NDArray[float64], NDArray[complex128]]

croak.processing.shiftnorm(grid, field)

Shift a time field so its peak sits at t=0 and normalise its peak to 1.

Parameters:

• grid (Grid)

• field (ArrayLike)

Return type:

NDArray[complex128]

croak.processing.normalize_spectral_phase(omega, ew, mask_w, Iw)

Spectral phase in the display gauge: unwrapped, de-tilted, peak-zeroed.

The gauge matters because spectral phase is only defined up to a constant and a linear term (an overall phase and a delay), so two phases can only be compared after both are pinned the same way. Shared by the retrieved spectrum and by a known truth phase precisely so the two cannot drift apart: plotting a truth normalised any other way would show a spurious tilt or offset that is pure gauge.

croak.processing.normalize_temporal_phase(t, field, mask_t)

Return temporal phase in the display gauge used for retrieved pulses.

The constant phase is not observable in the envelope retrieval; the linear temporal ramp is a carrier-frequency offset. Removing both from a known field and a retrieved field makes their temporal phases directly comparable. The ramp removal is a display choice: use the complex fields when a calibrated carrier offset matters.

Parameters:

• t (array_like) – Time axis (s).

• field (array_like) – Complex temporal envelope sampled on t.

• mask_t (array_like) – Boolean support mask used for the linear fit and final centring.

Returns:

numpy.ndarray – Unwrapped, de-ramped temporal phase (rad).

croak.processing.spectrogram(result, *, n_t=None, fw=None, lam_min=None, lam_max=None)

Gabor spectrogram of the retrieved pulse (for the dispersion view).

Returns (t_centers, wavelength, S) with S of shape (Nwavelength, Nt), converted to wavelength density and cropped to [lam_min, lam_max].

Parameters:

• result (RetrievalResult)

• n_t (int | None)

• fw (float | None)

• lam_min (float | None)

• lam_max (float | None)

croak.processing.post_filter(result, *, lam_lims=None, tau_lims=None)

Window a retrieved field in time and/or wavelength.

The temporal window is applied first (a Planck taper over tau_lims, transformed back to the spectrum and renormalised), then the spectral window (a Planck taper over the frequencies corresponding to lam_lims). Either may be None. Order matters: windowing in time redistributes spectral content, so the spectral taper must see the already-windowed field.

lam_lims/tau_lims should be the measurement windows λm/τm — the ranges over which the trace actually carries signal (available as lamm_lims / tau_lims). Passing the full retrieval-grid extent instead would taper only the outermost few bins and do essentially nothing.

Returns a new RetrievalResult with the filtered spectrum (the original is unchanged), so the effect can be toggled without re-running the retrieval. The stored trace, error and scale factor are not recomputed; they still describe the unfiltered fit.

Parameters:

• result (RetrievalResult)

• lam_lims (tuple[float, float] | None)

• tau_lims (tuple[float, float] | None)

Return type:

RetrievalResult

croak.processing.peak_wavelength(wavelength, intensity)

Return the wavelength (m) of an intensity spectrum’s maximum.

The result depends on the spectral density convention: pass the wavelength density (e.g. ProcessedResult.Ilam, not ProcessedResult.Iw) to match what the spectrum panel plots.

Parameters:

• wavelength (array_like) – Wavelength axis (m), ascending or descending. Non-positive and non-finite entries are dropped (see the grid’s zero crossing).

• intensity (array_like) – Spectral intensity on wavelength; scaling is irrelevant.

Returns:

float – Wavelength of the peak (m), or nan when no usable bin remains.

Return type:

float

Examples

>>> lam = np.linspace(700e-9, 900e-9, 201)
>>> round(float(peak_wavelength(lam, np.exp(-(((lam - 800e-9) / 20e-9) ** 2)))), 12)
8e-07

croak.processing.spectral_fwhm(wavelength, intensity, *, level=0.5)

Return the full width at half maximum (m) of an intensity spectrum.

Thin wrapper over croak.maths.fwhm() that first restricts the axis to its physical branch (see peak_wavelength() for the density convention).

Parameters:

• wavelength (array_like) – Wavelength axis (m), ascending or descending.

• intensity (array_like) – Spectral intensity on wavelength.

• level (float, optional) – Fraction of the peak at which to measure the width (default 0.5).

Returns:

float – Width in metres, or nan when the level is not crossed on both sides (a spectrum clipped by the grid) or no usable bin remains.

Return type:

float

Examples

>>> lam = np.linspace(700e-9, 900e-9, 2001)
>>> width = spectral_fwhm(lam, np.exp(-(((lam - 800e-9) / 20e-9) ** 2)))
>>> round(width / 1e-9, 1)
33.3

croak.processing.edge_energy_fraction(spectrum, *, collar_bins=10)

Fraction of the retrieved pulse energy sitting in the outermost grid bins.

A retrieval-quality flag. croak.preprocess.regrid() tapers the outermost collar_bins frequency rows at each end of the grid gradually to zero, so the trace carries little or no measurement there. The nonlinear mixing still couples those frequencies to the live rows, but only weakly (a probe spike at the grid edge moves the trace error by ~100× less than the same spike mid-band), so the field is nearly free there. Because \(|E(\omega)|^2 \ge 0\) and the true value sits on that bound, an unconstrained edge bin can only drift upward, which is why the failure always looks like spectral energy appearing at the band edges.

A well-conditioned retrieval on a generously chosen band returns well under 0.01 here. Values of a few percent mean part of the reported spectrum, and therefore part of the reported pulse duration, is an artefact. The usual causes, in order of how much they matter:

1. R_omega=True without reg_spectrum or phase_only. Per-frequency scaling discards exactly the information that pins the spectral amplitude, and it lowers the reported trace error while doing so.

2. A lam_min/lam_max band tight enough that the taper cuts live signal (croak.preprocess.regrid() warns about this at load time and records the discarded fraction as croak.preprocess.TraceData.taper_loss).

3. Noise fitted into signal-free bins on a noisy trace when no measured spectrum is supplied as a target; the COPRA warm-up of warm-lbfgs has no spectral penalty even when one is (see croak.warm_lbfgs).

A low value does not by itself establish field accuracy. See Post-processing and plotting results.

Parameters:

• spectrum (array_like) – Complex spectrum on the retrieval grid, centred order.

• collar_bins (int, optional) – Number of bins at each end to count; defaults to croak.preprocess.TAPER_COLLAR_BINS, matching the taper.

Returns:

float – Energy fraction in [0, 1]; 0.0 for a zero-energy spectrum, and for a grid too short for two non-overlapping collars.

Return type:

float

Examples

>>> import numpy as np
>>> ew = np.zeros(64, dtype=complex)
>>> ew[28:36] = 1.0                       # all energy mid-grid
>>> float(edge_energy_fraction(ew))
0.0
>>> ew[0] = 2.0                           # ... plus a spike on the edge
>>> round(float(edge_energy_fraction(ew)), 3)
0.333

class croak.processing.TruthPulse(t, It, fwhm, lam=None, phi_w=None, Iw=None, phi_t=None, omega=None, Eomega=None)

Bases: object

A known ground-truth pulse, in physical units, for overlaying and saving.

Carries the temporal intensity envelope and, when available, the complex spectrum of a known pulse so the retrieval screen can show it beneath the retrieved curves, the retrieval pipeline can use it as an initial guess, and the saver can record it. Two workflows produce one: the synthetic generator (from the generating complex field, via from_spectrum()) and the numerically simulated-scan loader (from the scan’s stored reference intensity, via from_intensity()).

The arrays are peak-centred (t=0 at the intensity peak) and peak-normalised so the overlay sits naturally against the normalised retrieved pulse.

Parameters:

• t (NDArray[float64])

• It (NDArray[float64])

• fwhm (float)

• lam (NDArray[float64] | None)

• phi_w (NDArray[float64] | None)

• Iw (NDArray[float64] | None)

• phi_t (NDArray[float64] | None)

• omega (NDArray[float64] | None)

• Eomega (NDArray[complex128] | None)

t

Temporal axis (s), peak-centred on the intensity maximum.

Type:

numpy.ndarray

It

Temporal intensity (power) envelope, normalised to unit peak.

Type:

numpy.ndarray

fwhm

Intensity FWHM of the temporal envelope (s).

Type:

float

lam

Wavelength axis (m), ascending, for the spectral intensity (or None).

Type:

numpy.ndarray or None

Iw

Spectral intensity vs lam (wavelength density), unit peak (or None).

Type:

numpy.ndarray or None

phi_t

Temporal phase (rad) on t, in the same display gauge as a retrieved phase, or None when no complex field is available.

Type:

numpy.ndarray or None

omega

Absolute angular-frequency axis (rad/s), ascending, for Eomega.

Type:

numpy.ndarray or None

Eomega

Known complex spectrum on omega in croak’s Fourier convention.

Type:

numpy.ndarray or None

phi_w: NDArray[float64] | None = None

Known spectral phase (rad) on lam, when the source provides a complex field. Normalised as the retrieved phase is — unwrapped, linear term removed, zeroed at the spectral peak — so the two are directly comparable on the same axes. None when only intensities are known.

classmethod from_spectrum(grid, ew, omega0, *, oversampling=8)

Build a TruthPulse from a complex spectrum (synthetic pulse).

The temporal intensity is oversampled (matching process_result()) for a smooth overlay and an accurate FWHM; the spectral intensity is the wavelength density ∝ |E(ω)|²·ω² (the coordinate Jacobian) mapped to ascending wavelength.

Parameters:

• grid (Grid) – Grid on which ew is sampled.

• ew (array_like) – Complex spectrum in centred order.

• omega0 (float) – Carrier angular frequency (rad/s).

• oversampling (int, optional) – Temporal zero-padding factor.

Returns:

TruthPulse – Known pulse retaining ew for truth-seeded retrievals.

Return type:

TruthPulse

classmethod from_intensity(t, It, *, lam=None, Iw=None, phi_w=None, phi_t=None, omega=None, Eomega=None, fwhm=None)

Build a TruthPulse from a sampled temporal intensity.

Used when the known pulse is available only as an intensity envelope (no phase), e.g. a Luna scansave reference pulse. Optional phase and complex-spectrum arrays preserve the full known field when the source file supplies it. fwhm defaults to the measured intensity FWHM of (t, It).

Parameters:

• t (array_like) – Temporal axis (s) and matching intensity envelope.

• It (array_like) – Temporal axis (s) and matching intensity envelope.

• lam (array_like, optional) – Matching wavelength axis (m) and spectral intensity.

• Iw (array_like, optional) – Matching wavelength axis (m) and spectral intensity.

• phi_w (array_like, optional) – Spectral phase (rad) sampled on lam.

• phi_t (array_like, optional) – Temporal phase (rad) sampled on t.

• omega (array_like, optional) – Matching absolute angular-frequency axis (rad/s) and complex spectrum in croak’s Fourier convention. Supply both or neither.

• Eomega (array_like, optional) – Matching absolute angular-frequency axis (rad/s) and complex spectrum in croak’s Fourier convention. Supply both or neither.

• fwhm (float, optional) – Known temporal intensity FWHM (s); measured from t/It when omitted.

Returns:

TruthPulse – Peak-centred, peak-normalised known pulse.

Raises:

ValueError – If paired axes/fields have inconsistent shapes.

Return type:

TruthPulse

spectrum_on_grid(grid, omega0)

Return the known complex spectrum on a retrieval grid.

Parameters:

• grid (Grid) – Retrieval grid whose frequencies are relative to omega0.

• omega0 (float) – Pulse carrier angular frequency (rad/s).

Returns:

numpy.ndarray – Complex spectrum in centred order. An identical ModelPNPS native grid is copied exactly; a regridded workflow interpolates amplitude and unwrapped phase and sets frequencies outside the stored range to zero.

Raises:

ValueError – If this truth carries no complex spectrum or its spectrum is empty.

Return type:

NDArray[complex128]

peak_power(energy)

Peak power (W) of this pulse when it carries energy joules.

It is peak-normalised, so its integral is the pulse’s effective duration and the peak power is energy / ∫(I/I_peak) dt — the same definition process_result() uses for the retrieved and transform-limited pulses, evaluated over the truth’s own time axis.

Parameters:

energy (float or None) – Pulse energy (J). None or non-positive means the absolute scale is unknown.

Returns:

float or None – Peak power in watts, or None when energy is unusable or the envelope integrates to zero.

Return type:

float | None

class croak.processing.ProcessedResult(t_over, It_retr, t_tl, It_tl, phi_t, mask_t, fwhm_retr, fwhm_tl, omega, omega0, wavelength, Iw, Ilam, Iw_meas, phi_w, phi_w_fit, mask_w, gdd_fs2, tod_fs3, fod_fs4, edge_energy, delays, trace_retr, trace_meas, residual, omega_marg_retr, omega_marg_meas, tau_marg_retr, tau_marg_meas, errors, error, result, energy=None, peak_power=None, peak_power_tl=None)

Bases: object

Display quantities derived from a RetrievalResult.

Parameters:

• t_over (NDArray[float64])

• It_retr (NDArray[float64])

• t_tl (NDArray[float64])

• It_tl (NDArray[float64])

• phi_t (NDArray[float64])

• mask_t (NDArray[bool])

• fwhm_retr (float)

• fwhm_tl (float)

• omega (NDArray[float64])

• omega0 (float)

• wavelength (NDArray[float64])

• Iw (NDArray[float64])

• Ilam (NDArray[float64])

• Iw_meas (NDArray[float64] | None)

• phi_w (NDArray[float64])

• phi_w_fit (NDArray[float64])

• mask_w (NDArray[bool])

• gdd_fs2 (float)

• tod_fs3 (float)

• fod_fs4 (float)

• edge_energy (float)

• delays (NDArray[float64])

• trace_retr (NDArray[float64])

• trace_meas (NDArray[float64] | None)

• residual (NDArray[float64] | None)

• omega_marg_retr (NDArray[float64])

• omega_marg_meas (NDArray[float64] | None)

• tau_marg_retr (NDArray[float64])

• tau_marg_meas (NDArray[float64] | None)

• errors (list[float])

• error (float)

• result (RetrievalResult)

• energy (float | None)

• peak_power (float | None)

• peak_power_tl (float | None)

edge_energy: float

Fraction of the retrieved pulse energy sitting in the outermost grid bins, where the trace carries little or no measurement. A quality flag, not a physical quantity; see edge_energy_fraction().

croak.processing.process_result(result, *, measured=None, Iomega_meas=None, oversampling=8, energy=None)

Compute display-ready temporal/spectral/trace quantities.

Parameters:

• result (RetrievalResult) – A completed retrieval.

• measured (array_like, optional) – The measured trace (Nomega, Ndelay) for marginals/residuals.

• Iomega_meas (array_like, optional) – Independent measured spectrum on the grid (for the spectral panel).

• oversampling (int, optional) – Temporal oversampling factor.

• energy (float, optional) – Measured pulse energy in joules. When given (and positive), the temporal intensity is turned into absolute instantaneous power and the peak powers are reported. With the intensity normalised to unit peak, instantaneous power is \(P(t)=E\,I(t)/\!\int I(t)\,dt\), so the peak power is \(E/\!\int I\,dt\) (the energy divided by the effective duration). The transform-limited reference carries the same energy, so it gives the maximum achievable peak power for the same spectral amplitude (its FWHM need not be the shortest for a structured pulse). None or a non-positive value leaves the plots in normalised units.

Returns:

ProcessedResult

Return type:

ProcessedResult

croak.processing.resolve_time_direction(result, truth, *, tolerance=0.001)

Pick the time direction that matches a known pulse, for ambiguous geometries.

An SHG-FROG trace is unchanged by \(E(t) \to E^*(-t)\): the two pulses fit the measurement equally well and a solver returns whichever branch it landed in. When the true pulse is known — a synthetic test, a simulated scan, or a cross-check against another diagnostic — the ambiguity can be settled by picking the branch that matches, which is what this does.

The choice is made on the temporal intensity profiles, compared by their best normalised cross-correlation over all relative lags. Using intensity rather than the field keeps the comparison independent of the absolute phase (also ambiguous), and maximising over lag keeps it independent of the time origin (likewise). Neither profile has to be a good match for the comparison to work — only for the right one to be the better of the two — so this stays reliable on an imperfect retrieval. A time-symmetric intensity profile cannot distinguish the two branches.

Reversal is applied by conjugating the spectrum, which is exactly equivalent and leaves the stored trace, mu and error valid, since the trace is invariant under it. Note that the retrieved GDD, TOD and the whole spectral phase change sign, so resolving the direction is what makes those signs comparable with the truth.

Parameters:

• result (RetrievalResult) – A completed retrieval.

• truth (TruthPulse) – The known pulse to match against.

• tolerance (float, optional) – Minimum improvement in correlation required to flip. A pulse whose intensity is close to time-symmetric scores almost equally either way; requiring a margin keeps the answer stable against noise rather than letting it toggle on a rounding difference. Default 1e-3.

Returns:

(RetrievalResult, bool) – The result oriented to match truth — the input object itself when no flip was applied — and whether it was flipped.

Return type:

tuple[RetrievalResult, bool]

Notes

Returns the result unchanged for interactions that fix the direction of time (SD, PG), where the retrieved orientation is already meaningful and flipping it would be wrong. See time_reversal_ambiguous.

Known-truth errors

Errors of a retrieval against a known pulse — available for the synthetic and simulated workflows, where the answer is known in advance. Distinct from croak.metrics, which scores the retrieval against the measured trace: a small trace error does not imply a correct pulse. See Loading simulated traces.

Retrieval errors against a known pulse.

croak.metrics answers “how well does the retrieval reproduce the measured trace” — the FROG error \(R\), which is all one can compute from a real measurement. When the true pulse is known (a synthetic generation, or a ModelPNPS scan carrying its complex field), a second and much sharper question becomes available: how close is the retrieved pulse to the true one. The two are not the same. A trace error can be small while the pulse is wrong, because a FROG trace does not determine the pulse uniquely and because a good optimiser can fit noise.

Three errors, all dimensionless and all zero for a perfect retrieval:

eps_It

temporal intensity \(|E(t)|^2\), compared after optimal intensity scaling, on the truth’s own time axis, minimised over the retrieval’s time-translation gauge.

eps_Iw

spectral intensity \(|\tilde E(\omega)|^2\) over the truth’s support band, compared after optimal intensity scaling so the arbitrary scale cancels.

eps_Ew

the complex field \(\tilde E(\omega)\), minimised over the three quantities a FROG measurement genuinely cannot determine — amplitude scale, absolute phase and delay. This is Geib’s \(\varepsilon\), the metric the PNPS literature quotes, and the only one of the three that cannot be satisfied by getting amplitude or phase right alone.

Each is a least-squares distance after removing exactly the gauge freedoms that are unobservable, and no others: a metric that quotients out more would flatter the retrieval, one that quotients out less would score a pure relabelling as error. Compare the same field reference and spectral density in both arguments. For SHG, resolve the conjugate time-reversal branch separately with croak.processing.resolve_time_direction().

class croak.truth_metrics.TruthErrors(eps_It, eps_Iw=None, eps_Ew=None)

Bases: object

The three known-truth errors of one retrieval.

Parameters:

• eps_It (float)

• eps_Iw (float | None)

• eps_Ew (float | None)

eps_It

Temporal-intensity error (eps_temporal_intensity()).

Type:

float

eps_Iw

Spectral-intensity error (eps_spectral_intensity()); None when the truth carries no complex spectrum.

Type:

float or None

eps_Ew

Complex-field error (eps_complex_field()); None when the truth carries no complex spectrum.

Type:

float or None

croak.truth_metrics.eps_temporal_intensity(t_retrieved, I_retrieved, t_truth, I_truth)

Least-squares relative error between two temporal intensities.

Both intensities are brought onto the truth’s time axis; the compared quantity is the two-norm relative error minimised over the free scale, \(\min_c \lVert c I - I_{\text{true}}\rVert_2 / \lVert I_{\text{true}}\rVert_2 = \sin\theta\), so the comparison is of pulse shape, not of scale. The retrieval’s time origin is a gauge freedom — a FROG trace is unchanged by translating the pulse — so it is removed in two steps: the peaks are aligned first (a coarse shift, without which a pulse parked elsewhere in the wider retrieval window can be interpolated away entirely and the error saturates near 1), then the cross-correlation peak refines it to sub-sample precision via croak.preprocess.marginal_peak_delay().

The refinement matters. Aligning only to the nearest sample leaves up to half a step of residual offset, which a normalised rms scores as roughly dt / 2 sigma. On an 8x-oversampled 20 fs Gaussian (0.125 fs step), a half-step delay scores 5.2e-3 under integer-bin alignment and 2.3e-5 with the sub-sample fit: a floor 220x lower.

Parameters:

• t_retrieved (array_like) – Retrieved time axis (s) and its temporal intensity (any scale).

• I_retrieved (array_like) – Retrieved time axis (s) and its temporal intensity (any scale).

• t_truth (array_like) – Known time axis (s), which must be uniformly sampled, and its intensity. The comparison is made on this axis.

• I_truth (array_like) – Known time axis (s), which must be uniformly sampled, and its intensity. The comparison is made on this axis.

Returns:

float – \(\varepsilon_{I_t} \ge 0\); 0 for shapes identical up to energy scale and a time translation.

Return type:

float

croak.truth_metrics.eps_spectral_intensity(spectrum, truth, *, threshold=0.02)

Least-squares relative error between two spectral intensities.

Compares \(|\tilde E(\omega)|^2\) on a shared frequency grid. Both are restricted to where the truth actually carries signal; the compared quantity is the two-norm relative error minimised over the free scale (\(\sin\theta\)), which removes the retrieval’s arbitrary overall scale; spectral intensity is already blind to the constant-phase and delay gauges. Like eps_complex_field(), the residual is formed directly, so identical shapes score zero to machine precision.

Parameters:

• spectrum (array_like) – Retrieved and known complex spectra on the same grid (centred order).

• truth (array_like) – Retrieved and known complex spectra on the same grid (centred order).

• threshold (float, optional) – Band cut as a fraction of the truth’s peak spectral intensity. Bins below it are excluded: they carry no information but do carry the retrieval’s noise floor, which would otherwise dominate a normalised distance.

Returns:

float – \(\varepsilon_{I_\omega} \ge 0\); 0 for identical shapes.

Raises:

ValueError – If the two spectra differ in shape, or the truth is empty in band.

Return type:

float

croak.truth_metrics.eps_complex_field(spectrum, truth)

Complex-field error, minimised over the unobservable gauges.

Geib’s \(\varepsilon\): the normalised distance between two complex spectra after minimising over amplitude scale, absolute phase and delay — exactly the three continuous transformations a PNPS measurement cannot resolve. Resolve any SHG conjugate time-reversal ambiguity separately. Written as an overlap,

\[\varepsilon = \sqrt{1 - \max_\tau \frac{|\langle \tilde E_\mathrm{truth}, \tilde E\,e^{i\omega\tau} \rangle|^2}{\|\tilde E_\mathrm{truth}\|^2\,\|\tilde E\|^2}},\]

where maximising the overlap magnitude over \(\tau\) handles the delay and taking the modulus handles the absolute phase, while the normalisation handles the scale. The delay scan is the DFT of \(\tilde E^*_\mathrm{truth}\tilde E\), zero-padded 32x so the optimum is found between grid steps rather than snapped to one.

Parameters:

• spectrum (array_like) – Retrieved and known complex spectra on the same grid (centred order).

• truth (array_like) – Retrieved and known complex spectra on the same grid (centred order).

Returns:

float – \(\varepsilon_{E_\omega} \in [0, 1]\); 0 for fields identical up to scale, absolute phase and delay, 1 for uncorrelated ones.

Raises:

ValueError – If the two spectra differ in shape, or either is identically zero.

Return type:

float

Notes

The displayed overlap formula is evaluated through its algebraically equivalent least-squares projection residual. This avoids the catastrophic cancellation in sqrt(1 - x**2) as the normalised overlap x approaches one, so fields identical up to the removed gauges score zero to machine precision rather than acquiring a platform-dependent ~sqrt(eps) floor.

croak.truth_metrics.truth_errors(result, truth, *, frame_p=0.0)

Compute every known-truth error available for one retrieval.

The spectral errors need the truth’s complex field, so they are reported only when it has one (a synthetic generation always does; a ModelPNPS scan does when the file stores the complex Eω datasets) and it overlaps the retrieval grid. The temporal error needs only the truth’s intensity envelope and is always available.

Parameters:

• result (RetrievalResult) – The retrieval to score.

• truth (TruthPulse) – The known pulse, from the synthetic or simulated load.

• frame_p (float, optional) – The spectrum_frame_p the retrieval ran with. A frame-corrected retrieval carries its amplitude in the on-axis frame, while the stored integrated beamlet truth lives in the measurement frame. This function divides the retrieved complex spectrum by ((omega + omega0) / omega0)**frame_p and recomputes its temporal field, so all three errors compare like with like: a phase-only retrieval then scores eps_Iw = 0 instead of the fixed frame footprint (~0.107 for a 1 fs spectrum at p = 0.5). frame_p=1 is the declared on-axis convention; 0.5 an empirical effective one. Use zero when no conversion was applied. See Conventions.

Returns:

TruthErrors – With eps_Iw/eps_Ew set to None for an intensity-only truth.

Return type:

TruthErrors

Examples

>>> errs = truth_errors(result, truth)
>>> errs.eps_It < 0.01
True

Marginal consistency checks

Pre-retrieval consistency checks from the frequency marginal: the marginal predicted from the independent spectrum (exact for SHG, transform-limited for PG/TG/SD), its centroid/width diagnostics, and the centroid-matching efficiency-exponent picker. The predictions assume the thin local model; propagation and collection change the marginal. The background is in FROG trace marginals as consistency checks.

Pre-retrieval consistency checks from FROG trace frequency marginals.

The frequency marginal of a FROG-type trace is its integral over delay,

\[M(\omega) = \int T(\omega, \tau)\,\mathrm{d}\tau ,\]

a one-dimensional function of optical frequency. A thin local FROG signal factorises into an undelayed part \(P(t)\) and a delayed gate \(g(t-\tau)\), and integrating over delay collapses the delay phase to a delta function, leaving a convolution of spectral intensities

\[M(\omega) \;\propto\; \bigl(|\tilde P|^2 * |\tilde g|^2\bigr)(\omega).\]

Two consequences hold for any pulse (intensity-weighted moments add under convolution): the marginal centroid is the sum of the factor centroids and its variance is the sum of the factor variances. From an independently measured fundamental spectrum \(S(\omega) = |\tilde E(\omega)|^2\) this gives, per geometry:

• SHG — marginal \(S * S\) (exact); centroid \(2\bar\omega_S\); RMS width \(\sqrt2\,\sigma_S\). Both anchors hard.

• PG / TG — marginal \(S * |\widehat{|E|^2}|^2\) (transform-limited); centroid \(\bar\omega_S\) (exact); width in \([\sigma_S, \sigma_{TL}]\).

• SD — marginal \(|\widehat{E^2}|^2 * S(-\omega)\) (transform-limited); centroid \(\approx\bar\omega_S\) (TL/symmetric only); width \(>\sigma_S\).

SHG is the strongest check: for a thin local medium with full delay support its marginal is fully phase-independent, so a mismatch points to the trace (calibration, spectral response, detector nonlinearity, scatter, delay support) rather than to the pulse. PG/TG keep one rigorous anchor (the centroid) and one rigorous lower bound (the width); SD has only soft, transform-limited (TL) anchors because its phase-dependent factor is the SHG spectrum.

This module computes those predictions (predicted_marginal(), predicted_centroid()), the moments of a measured marginal (centroid(), rms_width()), and an exponent picker that tilts a trace marginal onto the predicted centroid (auto_third_order_exponent()). The full physics is in the manual’s FROG trace marginals explanation page.

The predictions assume a thin local interaction and matching field references; spatial collection and propagation in a thick medium modify them. See FROG trace marginals as consistency checks.

class croak.marginal_checks.MarginalPrediction(interaction, omega, marginal, centroid, rms_width, width_bracket, centroid_is_hard, width_is_hard, sigma_spectrum)

Bases: object

A predicted FROG marginal and its consistency anchors.

Parameters:

• interaction (str)

• omega (NDArray[float64])

• marginal (NDArray[float64])

• centroid (float)

• rms_width (float)

• width_bracket (tuple[float, float | None] | None)

• centroid_is_hard (bool)

• width_is_hard (bool)

• sigma_spectrum (float)

interaction

Resolved geometry: "shg", "pg" or "sd" ("tg" maps to "pg").

Type:

str

omega

Absolute angular-frequency axis of the predicted marginal (rad/s).

Type:

numpy.ndarray

marginal

Predicted marginal, peak-normalised to unit maximum.

Type:

numpy.ndarray

centroid

Predicted centroid (rad/s): \(2\bar\omega_S\) for SHG, else \(\bar\omega_S\).

Type:

float

rms_width

Predicted RMS width (rad/s): the exact \(\sqrt2\,\sigma_S\) for SHG, otherwise the transform-limited width \(\sigma_\mathrm{TL}\).

Type:

float

width_bracket

(lower, upper) RMS-width bracket. For PG/TG the measured width should lie in [\sigma_S, \sigma_\mathrm{TL}] (lower bound rigorous, upper soft). For SD only the lower bound \(\sigma_S\) is rigorous (upper is None). None for SHG (the width is exact).

Type:

tuple of float or None

centroid_is_hard

Whether the centroid is a rigorous, phase-independent check (SHG, PG/TG) or a soft TL/symmetric anchor (SD).

Type:

bool

width_is_hard

Whether the RMS width is an exact check (SHG only).

Type:

bool

sigma_spectrum

RMS width of the fundamental spectrum \(\sigma_S\) (rad/s).

Type:

float

croak.marginal_checks.centroid(omega, marginal)

Intensity-weighted centroid \(\bar\omega = \sum\omega M/\sum M\).

Parameters:

• omega (array_like) – Coordinate axis (typically angular frequency, rad/s).

• marginal (array_like) – Non-negative weights (the marginal intensity).

Returns:

float – The first moment of marginal over omega.

Return type:

float

Notes

This is a discrete weighted sum. For non-uniform axes, include the integration weights in marginal before calling.

croak.marginal_checks.rms_width(omega, marginal)

Root-mean-square width of a distribution about its centroid.

Computes \(\sigma = \sqrt{\langle\omega^2\rangle - \bar\omega^2}\).

Parameters:

• omega (array_like) – Coordinate axis (typically angular frequency, rad/s).

• marginal (array_like) – Non-negative weights (the marginal intensity).

Returns:

float – The standard deviation of marginal over omega (floored at zero before the square root to absorb discretisation noise in the tails).

Return type:

float

Notes

This is a discrete weighted sum. For non-uniform axes, include the integration weights in marginal before calling.

croak.marginal_checks.spectral_autoconvolution(spectrum, *, circular=False)

Autoconvolution \(S * S\) of a spectral intensity.

This is the exact, phase-independent SHG-FROG marginal shape for a thin local medium with full delay support. Two boundary conventions are offered:

• circular=False (default): the linear convolution, which doubles the support to 2N - 1 samples so there is no wrap-around. Use this for the true predicted marginal shape (it lands on a doubled, twice-as-wide axis).

• circular=True: the circular convolution via the convolution theorem (\(S*S = \mathcal F\{(\mathcal F^{-1}S)^2\}\) on a centred array), returning N samples on the same grid as spectrum. Use this to scale a trace’s frequency rows onto the predicted SHG marginal in place (see croak.preprocess.load_and_clean()’s marginal_correct); it aliases if the autoconvolution support exceeds the grid.

Parameters:

• spectrum (array_like) – Spectral intensity \(S(\omega)\) on a uniform grid (>= 0). For circular=True it must be a centred array (zero frequency in the middle), matching croak’s grid convention.

• circular (bool, optional) – Select the circular, grid-matched convolution instead of the linear one.

Returns:

numpy.ndarray – The autoconvolution: length 2N - 1 (linear) or N (circular).

Return type:

NDArray[float64]

croak.marginal_checks.apply_marginal_correction(trace, Iomega)

Rescale a SHG trace’s frequency rows onto the spectrum-predicted marginal.

For a thin local medium the SHG frequency marginal is the spectral autoconvolution \(S * S\) of the independent fundamental spectrum Iomega (both on the same centred frequency grid as the trace). Each frequency row of trace is rescaled so its delay-integrated marginal matches that prediction (rows with no marginal power are zeroed), which corrects a measured/efficiency tilt along the spectral axis without touching the delay structure.

Shared by croak.preprocess.load_and_clean() (the marginal_correct option) and the marginal-check stage’s live preview, so the correction is defined in exactly one place.

Parameters:

• trace (array_like) – SHG trace (Nomega, Ndelay) on the centred frequency grid.

• Iomega (array_like) – Independent fundamental spectral intensity on the same grid (centred).

Returns:

numpy.ndarray – A new corrected trace (Nomega, Ndelay); the input is not mutated.

Return type:

NDArray[float64]

croak.marginal_checks.predicted_centroid(omega, spectrum, interaction)

Predicted marginal centroid (rad/s) from the fundamental spectrum.

Returns \(2\bar\omega_S\) for SHG (the marginal sits on the doubled axis) and \(\bar\omega_S\) for PG/TG/SD (the marginal sits at the fundamental). For SHG and PG/TG this is a hard, phase-independent check; for SD it holds only for a transform-limited or temporally symmetric pulse.

Parameters:

• omega (array_like) – Absolute angular-frequency axis of the spectrum (rad/s).

• spectrum (array_like) – Fundamental spectral intensity \(S(\omega)\) (>= 0).

• interaction (str) – "shg", "sd", "pg" or "tg" (case-insensitive).

Return type:

float

croak.marginal_checks.predicted_marginal(omega, spectrum, interaction)

Predict the FROG marginal and its checks from the fundamental spectrum.

Computes the predicted marginal shape by spectral convolution (exact for SHG, transform-limited for PG/TG/SD), then places it on the analytically predicted centroid and reports the centroid/width consistency anchors. The shape comes from the numerics; the placement from the exact centroid identity \(\bar\omega[M] = \bar\omega[|\tilde P|^2] + \bar\omega[|\tilde g|^2]\), which avoids fragile absolute-carrier bookkeeping through the FFTs.

Parameters:

• omega (array_like) – Absolute angular-frequency axis of the spectrum (rad/s), uniform and ascending.

• spectrum (array_like) – Fundamental spectral intensity \(S(\omega)\) (>= 0).

• interaction (str) – "shg", "sd", "pg" or "tg" (case-insensitive); "tg" is identical to "pg" (the transient-grating gate is \(|E|^2\)).

Returns:

MarginalPrediction – The predicted marginal on its own (doubled, for SHG) frequency axis and the centroid/width anchors.

Raises:

ValueError – If interaction is unknown or omega/spectrum are degenerate.

Return type:

MarginalPrediction

croak.marginal_checks.auto_third_order_exponent(lam, marginal, lam_spec, Ilam_spec, interaction, *, window=(2.1e-07, 1e-06), bounds=(0.0, 6.0), target_centroid=None)

Pick the (\lambda/\mu m)^n exponent that matches the marginal centroid.

The TG/PG third-order efficiency correction multiplies the trace marginal by a power law \((\lambda/\mu\mathrm{m})^n\) (croak.preprocess.third_order_scale()). This chooses n so the centroid of the corrected marginal (computed as an angular-frequency density, within a wavelength window) equals the centroid predicted from the independent spectrum (predicted_centroid(): \(\bar\omega_S\) at the fundamental for PG/TG/SD, \(2\bar\omega_S\) for SHG).

target_centroid (rad/s) overrides the spectrum-derived target. Pass the centroid of an already-computed prediction so the matched centroid is exactly the one being displayed — otherwise small differences in how the spectrum is windowed/resampled between the prediction and this fit leave the corrected marginal a few nm off the drawn prediction line.

Increasing n shifts weight to long wavelengths (low frequency), so the centroid decreases monotonically with n for a positive marginal and the match has at most one root, found by scipy.optimize.brentq() and clamped to bounds. If the target lies outside the achievable centroid range the nearer bound is returned.

Parameters:

• lam (array_like) – Wavelength grid of the trace marginal (m).

• marginal (array_like) – The trace’s wavelength marginal before the power-law correction.

• lam_spec (array_like) – Independent fundamental spectrum and its wavelength grid (m).

• Ilam_spec (array_like) – Independent fundamental spectrum and its wavelength grid (m).

• interaction (str) – "shg", "sd", "pg" or "tg" (case-insensitive).

• window (tuple of float, optional) – (lo, hi) trust window in metres (default (210e-9, 1000e-9)); only wavelengths inside it enter the centroid.

• bounds (tuple of float, optional) – (lo, hi) clamp on the returned exponent (default (0.0, 6.0)).

• target_centroid (float, optional) – Predicted centroid to match (rad/s). When None (default) it is computed from lam_spec/Ilam_spec over the window.

Returns:

float – The centroid-matching exponent, clamped to bounds.

Raises:

ValueError – If fewer than three wavelengths carry signal in the window, or the spectrum centroid cannot be formed.

Return type:

float

Notes

This fits an empirical tilt over the chosen window; it does not identify which physical factor (source or detector response) causes it.

Dispersion tuning

Interactive residual-dispersion tuning of a retrieved pulse.

Applies Taylor (GDD/TOD/FOD) and material dispersion to a retrieved complex spectrum and provides auto-tuners that maximise the temporal peak power (i.e. compress the pulse).

Sign convention matches croak.pulses.gaussian_pulse(): a positive GDD adds the phase \(+\tfrac12\,\mathrm{GDD}\,\omega^2\), so applying a material’s GDD and then its negative cancels exactly, and material_gdd is positive for normal dispersion.

Chirped-mirror dispersion is data-dependent; MIRRORS is an empty, extensible registry (use register_mirror()).

croak.dispersion.apply_dispersion(omega, omega0, ew, *, gdd=0.0, tod=0.0, fod=0.0, material_thicknesses=None, mirror_bounces=None)

Apply Taylor + material + mirror dispersion to a spectrum (returns a copy).

Parameters:

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

• omega0 (float) – Carrier angular frequency (rad/s).

• ew (array_like) – Complex spectrum (centred order).

• gdd (float, optional) – Taylor coefficients in SI (s², s³, s⁴). Phase added: GDD·ω²/2 + TOD·ω³/6 + FOD·ω⁴/24.

• tod (float, optional) – Taylor coefficients in SI (s², s³, s⁴). Phase added: GDD·ω²/2 + TOD·ω³/6 + FOD·ω⁴/24.

• fod (float, optional) – Taylor coefficients in SI (s², s³, s⁴). Phase added: GDD·ω²/2 + TOD·ω³/6 + FOD·ω⁴/24.

• material_thicknesses (dict, optional) – {material_name: thickness_m} of material to propagate through.

• mirror_bounces (dict, optional) – {mirror_name: n_bounces} for registered mirrors. A negative count back-propagates (removes) the mirror. Mirrors that also have a registered amplitude (MIRROR_AMPLITUDES) have their reflectivity applied too (removal divides by it, clamped by _MIRROR_AMPLITUDE_CAP).

Returns:

numpy.ndarray – The dispersed complex spectrum.

Return type:

NDArray[complex128]

Notes

Negative material lengths and bounce counts invert the supplied transfer function (back-propagation to an upstream plane). Phase-only propagation conserves spectral energy; mirror amplitude functions change it, and removing a low-reflectivity mirror amplifies noise in those bands.

croak.dispersion.dispersion_phase(omega, omega0, *, gdd=0.0, tod=0.0, fod=0.0, material_thicknesses=None, mirror_bounces=None)

Spectral phase that apply_dispersion() adds (Taylor/material/mirror).

This is exactly the argument of the exp(i·φ) factor applied to the spectrum, with the same sign convention; useful for plotting the applied dispersion separately from the pulse’s own retrieved phase.

Parameters:

• omega (ArrayLike)

• omega0 (float)

• gdd (float)

• tod (float)

• fod (float)

• material_thicknesses (dict[str, float] | None)

• mirror_bounces (dict[str, int] | None)

Return type:

NDArray[float64]

croak.dispersion.material_gdd(material, omega0, thickness)

GDD (fs²) of thickness of material at the carrier omega0.

Parameters:

• material (str)

• omega0 (float)

• thickness (float)

Return type:

float

croak.dispersion.material_tod(material, omega0, thickness)

TOD (fs³) of thickness of material at the carrier omega0.

Estimated from the third derivative of β via two second-derivative samples.

Parameters:

• material (str)

• omega0 (float)

• thickness (float)

Return type:

float

croak.dispersion.auto_taylor(grid, ew, omega0, *, bounds=((-10000.0, 10000.0), (-100000.0, 100000.0), (-1000000.0, 1000000.0)), method='differential_evolution')

Find (GDD, TOD, FOD) in SI that maximise the absolute peak power.

Parameters:

• grid – Grid, complex spectrum and carrier frequency.

• ew (ArrayLike) – Grid, complex spectrum and carrier frequency.

• omega0 (float) – Grid, complex spectrum and carrier frequency.

• bounds (sequence of (lo, hi)) – Search bounds for (GDD, TOD, FOD) in fs² / fs³ / fs⁴ (use the slider bounds).

• method (str) – One of OPTIMIZERS.

Returns:

tuple – (GDD, TOD, FOD) in SI (s², s³, s⁴).

Return type:

tuple[float, float, float]

croak.dispersion.auto_material(grid, ew, omega0, material, *, bounds=(-0.02, 0.02), method='differential_evolution')

Find the material thickness (m) maximising the absolute peak power.

Parameters:

• ew (ArrayLike)

• omega0 (float)

• material (str)

• bounds (tuple[float, float])

• method (str)

Return type:

float

croak.dispersion.register_mirror(name, phase_fn, amplitude_fn=None)

Register a mirror’s per-bounce spectral-phase (and optional amplitude).

Parameters:

• name (str) – Registry key used in mirror_bounces.

• phase_fn (callable) – phase(omega_abs) -> ndarray, the per-bounce spectral phase (rad).

• amplitude_fn (callable, optional) – r(omega_abs) -> ndarray, the per-bounce amplitude reflectivity (0–1), returning 1 outside its band. When given, apply_dispersion() applies it (removal — a negative bounce count — divides by it). When omitted the mirror is lossless (phase-only) and any stale amplitude under name is dropped.

Return type:

None

croak.dispersion.unregister_mirror(name)

Remove a mirror’s phase and amplitude functions from the registries.

Parameters:

name (str)

Return type:

None

croak.dispersion.OPTIMIZERS = ('differential_evolution', 'dual_annealing', 'direct', 'nelder-mead')

Available global/local optimisers for the dispersion auto-tuners.

croak.dispersion.MIRRORS: dict[str, Callable[[NDArray[float64]], NDArray[float64]]] = {}

Registry of chirped-mirror dispersion data, keyed by name. Each entry is a callable phase(omega_abs) -> ndarray giving the per-bounce spectral phase.

croak.dispersion.MIRROR_AMPLITUDES: dict[str, Callable[[NDArray[float64]], NDArray[float64]]] = {}

Optional per-bounce amplitude (reflectivity) registry, keyed by the same name as MIRRORS. Each entry is a callable r(omega_abs) -> ndarray returning the amplitude reflectivity (0–1), with r = 1 outside its valid band so it is a no-op there. A name present here has its reflectivity applied by apply_dispersion() (removal divides by it); a name absent here is treated as lossless (phase-only) — the convention for chirped compressors.

Mirrors

Built-in and custom mirror reflectivity / dispersion — chirped compressors (BUILTIN_MIRRORS) and beam-path coatings (BUILTIN_COATINGS, which carry reflectivity and are usually back-propagated) — registered into the croak.dispersion.MIRRORS table via register_mirror().

Chirped-mirror dispersion (reflectivity + group-delay data).

A port of the lookup_mirror / process_mirror_data routines of Luna.jl (PhysData.jl), an open-source (MIT) nonlinear pulse-propagation package, along with its bundled mirror data. Each mirror is described by measured reflectivity R(λ) and group-delay dispersion GDD(λ) (or group delay GD(λ)); the dispersion is integrated to a spectral phase, its overall group delay removed, and a Planck-taper window rolls it smoothly to zero outside the design band.

The public products are per-bounce spectral-phase functions phase(omega_abs) -> ndarray in croak’s sign convention (exp(+iφ), matching croak.dispersion), so they drop straight into croak.dispersion.MIRRORS via croak.dispersion.register_mirror().

Both the chirped compressor mirrors (BUILTIN_MIRRORS) and the beam-path coatings (BUILTIN_COATINGS) return an amplitude (reflectivity) function alongside the phase, so the measured loss per bounce is modelled: applying a mirror attenuates by r(λ) and removing it divides by r(λ) (clamped). For good chirped compressors R is ≳99 % in band, so this barely shifts the peak-power auto-compressor; for the coatings (bare Si, MgF₂-coated Al, exported by jaxmirror with per-polarisation reflectivity, phase and GDD) it matters, as they are usually removed (back-propagated) to recover the pulse upstream, where a small reflectivity amplifies noise. See croak.dispersion.register_mirror().

croak.mirrors.BUILTIN_COATINGS: dict[str, str] = {'MgF2/Al @ 0°': 'MgF2Al_00deg.csv', 'MgF2/Al @ 45°': 'MgF2Al_45deg.csv', 'Si @ 0°': 'Si_00deg.csv', 'Si @ 45°': 'Si_45deg.csv', 'Si @ 74°': 'Si_74deg.csv'}

Built-in beam-path coating mirrors (jaxmirror exports), keyed by a friendly "<coating> @ <angle>" name → CSV filename under data/mirrors/. These files share the column layout wavelength_nm, R_s, R_p, phase_s_rad, phase_p_rad, GDD_s_fs2, GDD_p_fs2 (phase convention sign +1, matching croak). Loaded by load_coating(), which builds the per-bounce phase from the GDD column (so the wrapped phase columns and the s/p ambiguity at normal incidence never bite) and the amplitude from the reflectivity column for the chosen polarisation.

croak.mirrors.BUILTIN_MIRRORS: dict[str, MirrorSpec] = {'HD120': MirrorSpec(r_file='HD120.csv', r_lam_col=0, r_cols=(1,), r_combine='single', d_file='HD120.csv', d_lam_col=2, d_cols=(3,), mode='gdd', lam0_nm=1030.0, taper_nm=(880.0, 900.0, 1200.0, 1220.0), fitorder=5, fit_window_fs=None, delim=','), 'HD59': MirrorSpec(r_file='', r_lam_col=0, r_cols=(), r_combine='single', d_file='HD59_GDD.dat', d_lam_col=0, d_cols=(1,), mode='gdd', lam0_nm=1030.0, taper_nm=(993.0, 1000.0, 1060.0, 1075.0), fitorder=5, fit_window_fs=None, delim=None), 'PC147': MirrorSpec(r_file='PC147.txt', r_lam_col=0, r_cols=(1,), r_combine='single', d_file='PC147.txt', d_lam_col=2, d_cols=(3,), mode='gdd', lam0_nm=1030.0, taper_nm=(640.0, 650.0, 1350.0, 1360.0), fitorder=5, fit_window_fs=None, delim=None), 'PC1611': MirrorSpec(r_file='PC1611.txt', r_lam_col=0, r_cols=(1,), r_combine='single', d_file='PC1611.txt', d_lam_col=2, d_cols=(3,), mode='gdd', lam0_nm=1030.0, taper_nm=(850.0, 855.0, 1195.0, 1200.0), fitorder=5, fit_window_fs=None, delim=None), 'PC1821': MirrorSpec(r_file='PC1821.txt', r_lam_col=0, r_cols=(1,), r_combine='single', d_file='PC1821.txt', d_lam_col=2, d_cols=(3,), mode='gdd', lam0_nm=1030.0, taper_nm=(800.0, 805.0, 1345.0, 1350.0), fitorder=5, fit_window_fs=None, delim=None), 'PC70': MirrorSpec(r_file='PC70_R.csv', r_lam_col=0, r_cols=(1, 2), r_combine='geom2', d_file='PC70_GDD.csv', d_lam_col=0, d_cols=(1, 2), mode='gdd', lam0_nm=800.0, taper_nm=(400.0, 450.0, 1200.0, 1300.0), fitorder=5, fit_window_fs=(2.0, 4.0), delim=','), 'ThorlabsUMC': MirrorSpec(r_file='UCxx-15FS_R.csv', r_lam_col=0, r_cols=(2,), r_combine='single', d_file='UCxx-15FS_GD.csv', d_lam_col=0, d_cols=(2,), mode='gd', lam0_nm=800.0, taper_nm=(640.0, 650.0, 1050.0, 1100.0), fitorder=3, fit_window_fs=(2.0, 4.0), delim=',')}

Built-in mirrors ported from Luna.jl, keyed by name.

class croak.mirrors.MirrorSpec(r_file, r_lam_col, r_cols, r_combine, d_file, d_lam_col, d_cols, mode, lam0_nm, taper_nm, fitorder, fit_window_fs, delim)

Bases: object

Description of a built-in mirror’s data files and processing constants.

Reflectivity (%) and dispersion come from one or two CSV/whitespace files; mode is "gdd" (fs², double-integrated) or "gd" (fs, single-integrated). taper_nm are the four Planck-taper corners (left0, left1, right1, right0) in nm; the inner pair is the valid range.

Parameters:

• r_file (str)

• r_lam_col (int)

• r_cols (tuple[int, ...])

• r_combine (str)

• d_file (str)

• d_lam_col (int)

• d_cols (tuple[int, ...])

• mode (str)

• lam0_nm (float)

• taper_nm (tuple[float, float, float, float])

• fitorder (int)

• fit_window_fs (tuple[float, float] | None)

• delim (str | None)

property valid_range_m: tuple[float, float]

the inner Planck-taper corners.

Type:

Valid wavelength range (m)

croak.mirrors.load_builtin(name)

Load a built-in mirror by name.

Returns (phase_fn, amplitude_fn, valid_range_m). amplitude_fn is the per-bounce amplitude reflectivity ready for application (floored, and 1 outside the design band where the phase is also tapered to zero); it is 1 everywhere for the lossless mirrors that ship without R data.

Parameters:

name (str)

Return type:

tuple[Callable[[ArrayLike], NDArray[float64]], Callable[[ArrayLike], NDArray[float64]], tuple[float, float]]

croak.mirrors.load_coating(name, pol='s')

Load a built-in beam-path coating mirror for one polarisation.

Unlike the chirped compressors of load_builtin(), these are simple optics whose dispersion is usually removed: the returned amplitude function therefore carries the measured reflectivity (applied — and divided out on removal — by croak.dispersion.apply_dispersion()).

The per-bounce phase is built from the GDD column (double-integrated and group-delay-removed by mirror_from_arrays()), which is wrap-free and, at normal incidence where GDD_s == GDD_p, sidesteps the ~π s/p ambiguity of the tabulated phase columns.

Parameters:

• name (str) – A key of BUILTIN_COATINGS (e.g. "MgF2/Al @ 45°").

• pol ({"s", "p"}) – Polarisation; selects the reflectivity and GDD columns. At 0° incidence the two are identical.

Returns:

(phase_fn, amplitude_fn, valid_range_m) – phase_fn(omega_abs) is the per-bounce spectral phase (croak’s exp(+iφ) convention); amplitude_fn(omega_abs) is the amplitude reflectivity r = √R (1 outside the data band, floored in-band); valid_range_m is the usable wavelength range (m).

Raises:

ValueError – If name is not a built-in coating or pol is not "s"/"p".

Return type:

tuple[Callable[[ArrayLike], NDArray[float64]], Callable[[ArrayLike], NDArray[float64]], tuple[float, float]]

croak.mirrors.mirror_from_arrays(lam_m, values, *, mode, lam0_m=None, taper_m=None, fitorder=5)

Build a per-bounce phase function from a custom mirror file’s arrays.

Parameters:

• lam_m (array_like) – Wavelength samples (m).

• values (array_like) – Either spectral phase (rad, mode="phase") or GDD (s², mode="gdd").

• mode ({"phase", "gdd"}) – Interpretation of values. "phase" is used as-is (after removing the overall group delay); "gdd" is double-integrated first.

• lam0_m (float, optional) – Reference wavelength for group-delay removal (defaults to the intensity-unweighted band centre).

• taper_m (tuple, optional) – Planck-taper corners (m). Defaults to a soft 2 %-margin taper just inside the data extents.

• fitorder (int) – Degree of the group-delay-removal polynomial.

Returns:

(phase_fn, valid_range_m)

Return type:

tuple[Callable[[ArrayLike], NDArray[float64]], tuple[float, float]]

croak.mirrors.process_mirror_data(lam_R, R, lam_disp, disp, *, mode, lam0_m, taper_m, fitorder=5, fit_window_fs=None)

Build per-bounce phase(omega_abs) and r(lambda_m) from mirror data.

Parameters:

• lam_R (array_like or None) – Reflectivity wavelength (m) and amplitude reflectivity (0–1), or None for a lossless mirror (r = 1).

• R (array_like or None) – Reflectivity wavelength (m) and amplitude reflectivity (0–1), or None for a lossless mirror (r = 1).

• lam_disp (array_like) – Dispersion wavelength (m) and the dispersion value: GDD (s²) if mode == "gdd" or group delay (s) if mode == "gd".

• disp (array_like) – Dispersion wavelength (m) and the dispersion value: GDD (s²) if mode == "gdd" or group delay (s) if mode == "gd".

• mode ({"gdd", "gd"}) – Whether disp is GDD (double-integrated to phase) or group delay (single-integrated).

• lam0_m (float) – Central wavelength used to reference the removed group delay.

• taper_m (tuple) – Planck-taper corners (left0, left1, right1, right0) in m.

• fitorder (int) – Degree of the group-delay-removal polynomial.

• fit_window_fs (tuple or None) – (lo, hi) band in ωfs over which to fit the group delay (edge kinks otherwise corrupt the fit). None fits all points.

Returns:

(phase_fn, r_fn) – phase_fn(omega_abs) is the per-bounce spectral phase in croak’s exp(+iφ) convention; r_fn(lambda_m) is the amplitude reflectivity (metadata only — not applied by the dispersion engine).

Return type:

tuple[Callable[[ArrayLike], NDArray[float64]], Callable[[ArrayLike], NDArray[float64]]]

Saving and loading

Saving retrieval results (HDF5) and GUI session options (TOML).

save_result() writes a RetrievalResult (optionally enriched with a ProcessedResult) to an HDF5 file using a flat key schema — one dataset per quantity, so the file is readable from any HDF5 binding without knowing this package. The same schema can be written into a named group instead of the root, which is how a series of retrievals (a parameter sweep, say) share one file. save_options() / load_options() persist the GUI’s per-stage parameters as options.toml.

croak.save.save_result(result, path, *, processed=None, measured=None, truth=None, group=None, force=False)

Write a retrieval result to an HDF5 file.

Parameters:

• result (RetrievalResult) – The retrieval to save.

• path (str) – Output .h5 path.

• processed (ProcessedResult, optional) – Pre-computed processing to include; if omitted and measured is given, it is computed.

• measured (numpy.ndarray, optional) – Measured trace, used to compute processed and stored.

• truth (TruthPulse, optional) – Known ground-truth pulse (synthetic / simulated workflows); its temporal intensity It_truth and, when present, spectral intensity, temporal / spectral phase, and the complex truth spectrum are stored.

• group (str, optional) – Write the (unchanged) flat schema into this group instead of the file root, appending to an existing file rather than truncating it. Use it to collect a series of retrievals — a parameter sweep, one group per run — in a single file; a slash-separated name ("full/z00") nests. The default (None) writes at the root, one result per file.

• force (bool, optional) – Overwrite an existing file — or, with group, an existing group of that name (other groups are left alone). Otherwise raise.

Returns:

str – The path written.

Raises:

FileExistsError – If the file (or, with group, the group) already exists and force is false.

Return type:

str

croak.save.load_result(path)

Read a saved result HDF5 file back into a (possibly nested) dict.

Top-level attrs and datasets become keys; any subgroup (e.g. the uncertainty group written by save_uncertainty()) becomes a nested dict.

Parameters:

path (str)

Return type:

dict[str, Any]

croak.save.result_from_saved(saved, *, grid=None)

Rebuild a RetrievalResult from a saved-result dict.

Reconstructs the retrieval (the retrieved spectrum, trace, error history and fitted extras) from a load_result() dict, so a previously saved retrieval can be reloaded without re-running the solver.

Parameters:

• saved (mapping) – A load_result() dict (the flat top-level entries of result.h5).

• grid (Grid, optional) – Grid to attach. Defaults to one matched to the saved omega axis (croak.grid.Grid.from_omega()); pass the live preprocess grid when available so the rehydrated result shares it exactly.

Returns:

RetrievalResult

Raises:

ValueError – If grid is supplied but its size does not match the saved spectrum.

Return type:

RetrievalResult

croak.save.save_uncertainty(results, path, *, force=False)

Write FWHM-uncertainty estimates into an uncertainty group of an HDF5 file.

Appends to an existing file (e.g. a result.h5 from save_result()) when present, otherwise creates it. Each estimate becomes a subgroup uncertainty/<method> holding the point estimate, 68 %/95 % intervals and ± (femtoseconds), the raw FWHM samples, and — when the estimate carries them — the drawn thicknesses (µm) and the 68 %/95 % temporal-intensity confidence bands (median and lower/upper envelopes) shown shaded in the plot.

Parameters:

• results (mapping of str to UncertaintyResult, or UncertaintyResult) – The estimates to save, keyed by a label (e.g. "parametric", "thickness", "combined"). A single result is wrapped under its method.

• path (str) – Output .h5 path (appended to if it already exists).

• force (bool, optional) – Overwrite an existing uncertainty group (otherwise raise).

Returns:

str – The path written.

Raises:

• ValueError – If results is empty.

• FileExistsError – If an uncertainty group already exists and force is false.

Return type:

str

croak.save.save_options(options, path)

Write GUI session options to a TOML file.

A croak_version / options_version header is written first so older files can be reloaded (and migrated, if ever needed) gracefully. Nested dicts (one per stage) become TOML tables; numpy values are coerced.

Parameters:

• options (dict[str, Any])

• path (str)

Return type:

str

croak.save.load_options(path)

Read GUI session options from a TOML file.

Parameters:

path (str)

Return type:

dict[str, Any]

Uncertainty analysis

Statistical and systematic error bars on the retrieved FWHM (the resampling bootstraps and the substrate-thickness systematic), and the fast analytic linearised Gauss–Newton covariance. Each estimator takes an optional propagation transfer function (from transfer_function()) to report the bar at a different beamline point. The background and caveats are in Uncertainty of the retrieved pulse duration.

Statistical uncertainty of a retrieved pulse from a single trace.

This module estimates a statistical (precision) error bar on the retrieved FWHM intensity duration by resampling, following the analysis in docs/explanation/uncertainty_estimation.md. Two statistical estimators are provided, alongside propagation of an external thickness uncertainty:

• Method A — parametric Monte-Carlo bootstrap (parametric_bootstrap()): simulate synthetic traces from the best-fit clean trace plus a NoiseModel, re-retrieve each (warm-started), and take the spread of the per-replicate FWHM.

• Method B — data-resampling bootstrap (resampling_bootstrap()): resample the measured trace over delay columns (a block bootstrap, recommended for delay-scanned data) or over frequency rows, and re-retrieve. Needs no noise model.

The FWHM is gauge-invariant under all the trivial FROG ambiguities (translation, CEP, and — for the TG/PG kernel — there is no direction-of-time ambiguity), so it is a well-defined scalar to bootstrap with no gauge alignment required.

The interval is a statistical bound conditional on the forward and noise models; it excludes systematic error (calibration, geometry, model mismatch), which often dominates. coverage_calibration() validates, on synthetic ground truth, that the nominal interval has the coverage it claims.

class croak.uncertainty.NoiseModel(sigma, gain=0.0)

Bases: object

Additive (and optionally shot-dependent) Gaussian trace noise.

The per-pixel variance is \(\sigma_0^2 + g\,T\), where \(\sigma_0\) is the additive (read) noise standard deviation, \(g\) the shot-noise coefficient and \(T\) the (clipped, non-negative) clean trace intensity. Samples are clipped at zero, matching a detector that cannot record negative counts (a censored Gaussian, not an exact Poisson model, and not a model of background-subtracted negative values).

Parameters:

• sigma (float or numpy.ndarray) – Additive noise standard deviation. A scalar is interpreted as a fraction of the clean-trace peak (so sigma=0.01 is 1 % of peak, matching the convention used in the test suite); an array is taken as absolute per-pixel standard deviation with the same shape as the trace.

• gain (float, optional) – Shot-noise coefficient \(g\) in trace-intensity units (0 ⇒ purely additive noise). Adds a signal-dependent variance term gain * intensity.

Notes

Pixels are drawn independently, so correlations introduced by filtering or interpolation are not modelled.

std(clean)

Return the per-pixel standard deviation for a given clean trace.

Parameters:

clean (NDArray[float64])

Return type:

NDArray[float64]

draw(clean, rng)

Draw one noisy trace clip(clean + noise, 0, None).

Parameters:

• clean (NDArray[float64])

• rng (Generator)

Return type:

NDArray[float64]

croak.uncertainty.noise_from_residual(measured, result, *, robust=True)

Estimate a NoiseModel from the post-fit residual.

Uses the residual measured - result.trace after a converged retrieval. The residual is noise plus any model misfit, so an imperfect forward model inflates the estimate, while the fit absorbing some noise deflates it slightly. Compare with an independently measured background where one is available. Needs no extra data.

Parameters:

• measured (array_like) – The measured trace (Nomega, Ndelay) (same normalisation as result.trace).

• result (RetrievalResult) – A completed retrieval; result.trace is the scaled clean model trace.

• robust (bool, optional) – If true (default), estimate the standard deviation from the median absolute deviation (MAD × 1.4826), which resists the structured misfit near the signal; otherwise use the plain sample standard deviation.

Returns:

NoiseModel – An additive model whose scalar sigma is the residual std as a fraction of the measured peak.

Return type:

NoiseModel

croak.uncertainty.noise_from_background(raw_trace, region)

Estimate a read-noise NoiseModel from a signal-free region.

Computes the standard deviation of a background patch of the raw trace. Captures only additive read noise (not shot noise on the signal), and the patch must be genuinely signal-free.

Parameters:

• raw_trace (array_like) – The trace to take the background from (ideally the unfiltered measurement).

• region (numpy.ndarray of bool, or tuple of slice) – A boolean mask selecting background pixels, or a (rows, cols) slice tuple.

Returns:

NoiseModel – An additive model whose scalar sigma is the background std as a fraction of the trace peak.

Return type:

NoiseModel

class croak.uncertainty.UncertaintyResult(statistic, point_estimate, samples, interval_68, interval_95, interval_method, method, n_resamples, n_converged, frog_errors, base_error, noise_floor=None, warm_start=True, profiles=None, t_profile=None, thicknesses=None, thickness_central=None, thickness_sigma=None, components=None, propagation_label=None)

Bases: object

Outcome of an uncertainty estimation: the FWHM distribution and its intervals.

Parameters:

• statistic (str)

• point_estimate (float)

• samples (NDArray[float64])

• interval_68 (tuple[float, float])

• interval_95 (tuple[float, float])

• interval_method (str)

• method (str)

• n_resamples (int)

• n_converged (int)

• frog_errors (NDArray[float64])

• base_error (float)

• noise_floor (float | None)

• warm_start (bool)

• profiles (NDArray[float64] | None)

• t_profile (NDArray[float64] | None)

• thicknesses (NDArray[float64] | None)

• thickness_central (float | None)

• thickness_sigma (float | None)

• components (tuple[str, ...] | None)

• propagation_label (str | None)

statistic

The quantity bootstrapped (currently always "fwhm", the temporal intensity FWHM in seconds).

Type:

str

point_estimate

Full-data best-fit FWHM (s) — the value the interval is centred reporting around.

Type:

float

samples

Per-replicate FWHM (s), converged replicates only.

Type:

numpy.ndarray

interval_68, interval_95

Lower/upper bounds (s) of the 68 % and 95 % intervals.

Type:

tuple of float

interval_method

"bias-corrected" or "percentile".

Type:

str

method

"parametric", "delay", "frequency" (statistical bootstrap) or "thickness" (substrate-thickness systematic propagation).

Type:

str

n_resamples

Number of replicates attempted.

Type:

int

n_converged

Number kept after the convergence filter (the effective B).

Type:

int

frog_errors

Per-replicate trace error R (converged replicates).

Type:

numpy.ndarray

base_error

Full-data trace error R (the convergence reference).

Type:

float

noise_floor

Median per-replicate noise floor R0 (Method A only; None otherwise).

Type:

float or None

warm_start

Whether replicates were warm-started from the full-data solution.

Type:

bool

profiles

(n_converged, n_t) peak-aligned, oversampled temporal intensities, if collected.

Type:

numpy.ndarray or None

t_profile

The shared oversampled time axis (s) for profiles.

Type:

numpy.ndarray or None

components

For method="combined", the methods that were combined (e.g. ("parametric", "thickness")); None otherwise.

Type:

tuple of str or None

propagation_label

Human-readable description of the beamline point the FWHM/band was evaluated at (e.g. "+1.0 mm SiO2"), when a propagation transfer function was applied; None for the measurement plane.

Type:

str or None

thicknesses

Per-replicate substrate thickness drawn from the prior (m), converged replicates only (method="thickness" only; None otherwise). Pairs element-wise with samples/frog_errors for the FWHM-vs-thickness and error-vs-thickness diagnostics.

Type:

numpy.ndarray or None

thickness_central

Central (measured) substrate thickness \(L_0\) (m) of the thickness prior (method="thickness" only).

Type:

float or None

thickness_sigma

Standard deviation \(\sigma_L\) (m) of the measured thickness prior (method="thickness" only).

Type:

float or None

property median: float

Median of the bootstrap FWHM samples (s).

property std: float

Standard deviation of the bootstrap FWHM samples (s).

property plus_minus: float

Symmetric half-width of the 68 % interval (s).

summary()

Return a one-line "x ± y fs (…)" summary tagged with the error class.

The caveat differs by method: the statistical bootstraps (parametric/delay/frequency) flag that systematics are excluded, whereas thickness reports a systematic error bar from the measured substrate-thickness prior and names that prior.

Return type:

str

croak.uncertainty.parametric_bootstrap(result, *, noise, n_resamples=300, warm_start=True, rng=None, oversampling=8, convergence_tol=3.0, interval='bias-corrected', collect_profiles=False, propagation=None, propagation_label=None, callback=None, **retrieve_kwargs)

Estimate FWHM uncertainty by the parametric Monte-Carlo bootstrap (Method A).

Simulates n_resamples synthetic traces from the best-fit clean trace (result.trace) plus draws from noise, re-retrieves each (warm-started by default), and returns the spread of the per-replicate FWHM. Conditions on the fitted model as truth, so it yields a statistical interval only.

Parameters:

• result (RetrievalResult) – A completed retrieval (must carry result.trace).

• noise (NoiseModel or float) – The measurement noise model; a bare float is taken as a fraction-of-peak additive NoiseModel.

• n_resamples (int, optional) – Number of bootstrap replicates B (≈200–300 for a 1σ interval).

• warm_start (bool, optional) – Warm-start each replicate from result.spectrum (recommended — isolates measurement variance from optimiser basin-hopping). If false, each replicate starts from a random guess.

• rng (numpy.random.Generator, optional) – Random generator (defaults to a fresh one).

• oversampling (int, optional) – Temporal oversampling for the FWHM (default 8, matching process_result).

• convergence_tol (float, optional) – Drop replicates whose error exceeds this multiple of the median replicate error.

• interval ({"bias-corrected", "percentile"}, optional) – Interval estimator.

• collect_profiles (bool, optional) – Also collect peak-aligned oversampled intensity profiles (for a confidence band).

• propagation (numpy.ndarray of complex, optional) – Spectral transfer function \(H(\omega)\) (same length as result.spectrum) propagating the pulse to a different beamline point; each replicate’s spectrum is multiplied by it before the FWHM/band is taken (see _apply_propagation). None (default) reports at the measurement plane.

• propagation_label (str, optional) – Description of the propagation point recorded on the result (e.g. "+1.0 mm SiO2").

• callback (callable, optional) – callback(done, total, R) after each replicate; raising aborts the run.

• **retrieve_kwargs – Forwarded to croak.retrieve() — must reproduce the original retrieval settings (maxiters, material/thickness, reg_*, phase_basis …).

Returns:

UncertaintyResult

Return type:

UncertaintyResult

croak.uncertainty.resampling_bootstrap(result, measured, *, unit='delay', n_resamples=300, leave_out=None, warm_start=True, rng=None, oversampling=8, convergence_tol=3.0, interval='bias-corrected', collect_profiles=False, propagation=None, propagation_label=None, callback=None, **retrieve_kwargs)

Estimate FWHM uncertainty by resampling the measured trace (Method B).

Needs no noise model. Two resampling units are supported:

• unit="delay" — resample whole delay columns with replacement (a block bootstrap; the natural unit for delay-scanned data, where each column is one laser exposure). The resampled delays may contain duplicates; they are sorted so the delay axis stays monotone.

• unit="frequency" — resample frequency rows via multinomial pixel counts, passed to the weighted retrieval as weights = √counts (weight √k reproduces a multiplicity-k row in the weighted least-squares objective).

Parameters:

• result (RetrievalResult) – A completed retrieval (provides the grid, delays, interaction and warm-start guess).

• measured (array_like) – The measured trace (Nomega, Ndelay) to resample.

• unit ({"delay", "frequency"}, optional) – The resampling unit.

• n_resamples (int, optional) – Number of bootstrap replicates B.

• leave_out (float, optional) – If given, instead of resampling with replacement, drop this fraction of columns (unit="delay") or rows (unit="frequency") each replicate.

• warm_start (bool) – As in parametric_bootstrap().

• rng (Generator | None) – As in parametric_bootstrap().

• oversampling (int) – As in parametric_bootstrap().

• convergence_tol (float) – As in parametric_bootstrap().

• interval (str) – As in parametric_bootstrap().

• collect_profiles (bool) – As in parametric_bootstrap().

• callback – As in parametric_bootstrap().

• propagation (NDArray[complex128] | None) – As in parametric_bootstrap(): a spectral transfer function \(H(\omega)\) propagating each replicate to a different beamline point, and its description.

• propagation_label (str | None) – As in parametric_bootstrap(): a spectral transfer function \(H(\omega)\) propagating each replicate to a different beamline point, and its description.

• **retrieve_kwargs – Forwarded to croak.retrieve() (must reproduce the original settings).

Returns:

UncertaintyResult

Return type:

UncertaintyResult

Notes

Resampled delay columns contain repeats and non-uniform spacings, so "delay" resampling is incompatible with models that need a uniform delay axis (the smearing convolution, marginal-peak centring); use "frequency" or the parametric bootstrap there. Neither method accounts for calibration errors.

croak.uncertainty.thickness_bootstrap(result, measured, *, central_thickness, thickness_sigma, material, n_resamples=300, warm_start=True, rng=None, oversampling=8, convergence_tol=1000000000.0, interval='bias-corrected', collect_profiles=False, clip_min=0.0, propagation=None, propagation_label=None, callback=None, **retrieve_kwargs)

Propagate the measured substrate-thickness uncertainty onto the FWHM.

A dispersive retrieval fits the trace through a forward model that propagates the field through a substrate of an assumed thickness \(L\); the retrieved field — and hence its temporal-intensity FWHM — depends on \(L\). When the substrate thickness is an independently measured input with a Gaussian prior \(L \sim \mathcal N(L_0, \sigma_L)\) (e.g. 9.952 ± 0.539 µm), the FWHM inherits a systematic uncertainty that the statistical bootstraps (parametric_bootstrap(), resampling_bootstrap()) cannot see — every one of their replicates reuses the same fixed thickness.

This estimator marginalises that prior by brute-force Monte-Carlo: it draws n_resamples thicknesses from \(\mathcal N(L_0, \sigma_L)\) (clipped at clip_min), rebuilds the dispersive forward model at each, re-retrieves the same measured trace (warm-started from the central-thickness solution), and returns the spread of the per-draw FWHM. It is the systematic complement to the statistical bootstrap, and the two combine in quadrature when independent (see combine_uncertainties()).

Because thickness is a solver-constructor parameter (the dispersion model is built in ForwardModel), the solver is rebuilt for each draw via _build_runner — cheap next to the retrieval iterations.

Parameters:

• result (RetrievalResult) – A completed dispersive retrieval (provides the grid, delays, interaction, carrier omega0 and the warm-start guess). Its algorithm and the forwarded retrieve_kwargs must reproduce the original settings.

• measured (array_like) – The measured trace (Nomega, Ndelay) to re-fit at each drawn thickness (the data are held fixed; only the assumed thickness varies).

• central_thickness (float) – Central (measured) substrate thickness \(L_0\) in metres.

• thickness_sigma (float) – Standard deviation \(\sigma_L\) of the thickness prior in metres (>= 0; 0 gives a degenerate, zero-width distribution).

• material (str) – Substrate material name (e.g. "SiO2"), as accepted by croak.materials.refractive_index().

• n_resamples (int, optional) – Number of Monte-Carlo draws B (≈200–300 for a 1σ interval).

• warm_start (bool, optional) – Warm-start each draw from the central-thickness solution (recommended: isolates the thickness-induced movement from optimiser basin-hopping). If false, each draw starts from a random guess.

• rng (numpy.random.Generator, optional) – Random generator (defaults to a fresh one).

• oversampling (int, optional) – Temporal oversampling for the FWHM (default 8, matching process_result).

• convergence_tol (float, optional) – Drop draws whose error exceeds this multiple of the median draw error. Defaults to a large value (effectively disabled): unlike the statistical bootstrap, an off-nominal thickness fits the trace genuinely worse, so a raised \(R\) in the tails is mostly physical misfit, not stagnation, and must not be filtered out. Only non-finite (failed) draws are dropped.

• interval ({"bias-corrected", "percentile"}, optional) – Interval estimator.

• collect_profiles (bool, optional) – Also collect peak-aligned oversampled intensity profiles (for a temporal band).

• clip_min (float, optional) – Lower clip (m) for drawn thicknesses (default 0.0; thickness cannot be negative).

• propagation (NDArray[complex128] | None) – As in parametric_bootstrap(): a spectral transfer function \(H(\omega)\) propagating each draw to a different beamline point, and its description. Stacks with the substrate dispersion the retrieval already models (the draw is fit at the measurement plane, then propagated).

• propagation_label (str | None) – As in parametric_bootstrap(): a spectral transfer function \(H(\omega)\) propagating each draw to a different beamline point, and its description. Stacks with the substrate dispersion the retrieval already models (the draw is fit at the measurement plane, then propagated).

• callback (callable, optional) – callback(done, total, R) after each draw; raising aborts the run.

• **retrieve_kwargs – Forwarded to the solver constructor (maxiters, npoints, reg_*, phase_basis …), exactly reproducing the original retrieval. Any material/thickness here are ignored — they are injected per draw.

Returns:

UncertaintyResult – With method="thickness", thicknesses the per-draw thickness (m), and thickness_central/thickness_sigma recording the prior.

Raises:

ValueError – If thickness_sigma < 0, n_resamples < 1 or material is falsy.

Return type:

UncertaintyResult

croak.uncertainty.combine_uncertainties(results, *, interval='bias-corrected', n_samples=2000, rng=None)

Combine independent FWHM uncertainty contributions into one error bar.

Treats each input as an independent source of FWHM deviation (e.g. a statistical bootstrap and the substrate-thickness systematic) and forms the total by convolving their centred sample distributions: a Monte-Carlo draw of the total deviation is the sum of one independently-drawn deviation from each source. This preserves skew and reduces to the familiar quadrature \(\sigma_\text{tot}=\sqrt{\sum_i\sigma_i^2}\) in the Gaussian limit. If every input carries a temporal-intensity band (profiles on a common time axis), the bands are combined the same way, per time sample.

Warning

Only combine sources that estimate different, independent contributions. parametric/delay/frequency all estimate the same statistical precision — combining them double-counts it. A sound budget combines one statistical estimate with each independent systematic (e.g. thickness).

Parameters:

• results (sequence of UncertaintyResult) – Two or more independent contributions (each must carry samples).

• interval ({"bias-corrected", "percentile"}, optional) – Interval estimator applied to the combined samples.

• n_samples (int, optional) – Number of Monte-Carlo draws of the combined distribution.

• rng (numpy.random.Generator, optional) – Random generator (defaults to a fresh one).

Returns:

UncertaintyResult – With method="combined", components listing the input methods, the combined samples and 68 %/95 % intervals, and the combined temporal band in profiles/t_profile when all inputs carried one.

Raises:

ValueError – If fewer than two results are given, or none carries any samples.

Return type:

UncertaintyResult

Notes

The function assumes the contributions are independent; it does not test this. Correlated contributions (e.g. a shared calibration) need joint propagation. Profile intervals are pointwise, not simultaneous bands.

croak.uncertainty.estimate_fwhm_uncertainty(result, measured=None, *, method='parametric', noise=None, **kwargs)

Estimate FWHM uncertainty by statistical resampling or thickness variation.

Parameters:

• result (RetrievalResult) – A completed retrieval.

• measured (array_like, optional) – The measured trace. Required for delay, frequency and thickness. For parametric, it is used to estimate noise if noise is omitted; simulated replicates are generated from result.trace.

• method ({"parametric", "delay", "frequency", "thickness"}, optional) – "parametric" → parametric_bootstrap(); "delay"/"frequency" → resampling_bootstrap(); "thickness" → thickness_bootstrap() (the substrate-thickness systematic, not a statistical bootstrap).

• noise (NoiseModel or float, optional) – Noise model for the parametric method. If omitted there and measured is given, it is estimated from the residual via noise_from_residual(). Ignored by the other methods.

• **kwargs – Forwarded to the chosen function. For method="thickness" this must include central_thickness, thickness_sigma and material. Pass propagation (a spectral transfer function H(ω)) and propagation_label to report the FWHM/band at a different beamline point (see _apply_propagation).

Returns:

UncertaintyResult

Return type:

UncertaintyResult

Examples

>>> import numpy as np, croak
>>> g = croak.Grid(256, dt=0.5e-15)
>>> ew = croak.gaussian_pulse(g, 6e-15)
>>> delays = np.linspace(-40e-15, 40e-15, 40)
>>> trace = croak.maketrace(g.omega, delays, ew, "pg")
>>> res = croak.retrieve(trace, g.omega, delays, "pg", guess=ew, progress=False)
>>> u = croak.estimate_fwhm_uncertainty(
...     res, trace, method="parametric", noise=0.01,
...     n_resamples=20, maxiters=30, rng=np.random.default_rng(0))
>>> bool(u.point_estimate > 0)
True

croak.uncertainty.percentile_interval(samples, level=0.68)

Return the central level percentile interval of samples.

Parameters:

• samples (ArrayLike)

• level (float)

Return type:

tuple[float, float]

croak.uncertainty.bias_corrected_interval(samples, point_estimate, level=0.68)

Return a bias-corrected (BC) bootstrap interval.

The BC interval shifts the percentile cut-offs by a bias correction \(z_0 = \Phi^{-1}(\#\{b: \theta^*_b < \hat\theta\}/B)\) to account for median bias in a skewed bootstrap distribution. The acceleration term of the full BCa interval is omitted (a = 0) because it requires a jackknife over the original data — each point of which is a full retrieval. Falls back to percentile_interval() when the bias correction is degenerate (all replicates on one side of the point estimate).

Parameters:

• samples (array_like) – Bootstrap replicates of the statistic.

• point_estimate (float) – The full-data estimate \(\hat\theta\) (the bias-correction reference).

• level (float, optional) – Central probability of the interval (0.68 ≈ 1σ).

Return type:

tuple[float, float]

croak.uncertainty.synthetic_experiment(*, fwhm, phases=(), shape='gaussian', n=128, dt=4e-16, lambda0=8e-07, ndelay=80, delay_max=6e-14, interaction='pg', oversampling=8)

Build a synthetic experiment with a known true FWHM, from public primitives.

Returns a dict with grid, ew (true spectrum), omega0, delays, trace (clean), interaction and true_fwhm (s). Mirrors the test fixture make_synthetic_experiment but uses only the public API so it can drive coverage_calibration().

Parameters:

• fwhm (float) – Intensity FWHM of the transform-limited pulse (s).

• phases (sequence of float, optional) – Taylor spectral-phase coefficients [GDD, TOD, …] (s², s³, …).

• shape ({"gaussian", "sech"}, optional) – Pulse envelope shape.

• n (int) – Grid, axis and forward-model parameters.

• dt (float) – Grid, axis and forward-model parameters.

• lambda0 (float) – Grid, axis and forward-model parameters.

• ndelay (int) – Grid, axis and forward-model parameters.

• delay_max (float) – Grid, axis and forward-model parameters.

• interaction (str) – Grid, axis and forward-model parameters.

• oversampling (int) – Grid, axis and forward-model parameters.

Return type:

dict

class croak.uncertainty.CoverageResult(level, n_trials, n_resamples, method, coverage, mean_bias, biases, covered, true_fwhms=<factory>)

Bases: object

Empirical coverage of a nominal bootstrap interval on synthetic ground truth.

Parameters:

• level (float)

• n_trials (int)

• n_resamples (int)

• method (str)

• coverage (float)

• mean_bias (float)

• biases (NDArray[float64])

• covered (NDArray[bool])

• true_fwhms (NDArray[float64])

summary()

Return a one-line coverage summary.

Return type:

str

croak.uncertainty.coverage_calibration(*, noise, n_trials=50, n_resamples=200, level=0.68, method='parametric', interval='bias-corrected', fwhm_range=(4e-15, 1e-14), gdd_range=(0.0, 4e-29), shapes=('gaussian',), algorithm='warm-lbfgs', guess_fwhm=None, rng=None, callback=None, **retrieve_kwargs)

Validate interval coverage on synthetic experiments with known FWHM.

For each trial, draws a random true pulse (FWHM, GDD, shape), synthesises and noises its trace, retrieves it, runs the bootstrap, and records whether the nominal-level interval contains the true FWHM. This converts a “spread” into a calibrated interval: if the nominal 68 % interval covers ≈68 % of trials, it is calibrated for the tested pulse population and noise model.

Parameters:

• noise (NoiseModel or float) – Noise injected into each synthetic trace (and used by the parametric method).

• n_trials (int, optional) – Number of synthetic experiments.

• n_resamples (int, optional) – Bootstrap replicates per trial.

• level (float, optional) – Nominal interval probability to test (e.g. 0.68).

• method ({"parametric", "delay", "frequency"}, optional) – Bootstrap method to validate.

• interval ({"bias-corrected", "percentile"}, optional) – Interval estimator (applied to the bootstrap samples at the requested level).

• fwhm_range (tuple of float, optional) – Uniform sampling ranges for the true FWHM (s) and GDD (s²).

• gdd_range (tuple of float, optional) – Uniform sampling ranges for the true FWHM (s) and GDD (s²).

• shapes (tuple of str, optional) – Pulse shapes to sample from.

• algorithm (str, optional) – Retrieval algorithm for each trial.

• guess_fwhm (float, optional) – FWHM of the transform-limited warm-start guess (defaults to the midpoint of fwhm_range).

• rng (numpy.random.Generator, optional) – Random generator.

• callback (callable, optional) – callback(trial, n_trials) after each trial.

• **retrieve_kwargs – Forwarded to croak.retrieve() and the bootstrap (e.g. maxiters).

Returns:

CoverageResult

Return type:

CoverageResult

Linearised (Gauss–Newton/Laplace) covariance of a retrieved pulse.

Pulse retrieval is a nonlinear least-squares fit: a damped Gauss–Newton / Levenberg–Marquardt solver (croak.lm.LM, lm-optx) drives the flattened FROG residual \(r(u)\) to a minimum \(u^\*\), where u is the real parameter vector (spectral amplitude+phase, or B-spline phase coefficients) of croak._jax_pulse. At that optimum the standard nonlinear-least-squares covariance applies:

\[\mathrm{cov}(u) \approx \hat\sigma^2\,(J^\mathsf{T}J)^{+}, \qquad J = \left.\frac{\partial r}{\partial u}\right|_{u^\*},\]

the Gauss–Newton approximation to the inverse Hessian. With the residual whitened by a NoiseModel (each pixel divided by its noise σ) the scale \(\hat\sigma^2 = 1\) and cov(u) = (JᵀJ)⁺ directly; without a noise model it falls back to the internal reduced-χ² estimate \(\hat\sigma^2 = \lVert r(u^\*)\rVert^2/(m-n)\).

This is the fast, analytic complement to the resampling bootstrap in croak.uncertainty: one Jacobian and one eigensolve instead of hundreds of re-retrievals. Two points govern how to read it:

• Gauge degeneracy. FROG’s continuous trivial ambiguities — the absolute phase (CEP) and the absolute timing (a linear spectral phase) — are flat directions of the objective, so \(J^\mathsf{T}J\) is rank-deficient. We use the pseudo-inverse (eigenvalues below rcond × the largest are dropped). The gauge directions carry no uncertainty for any gauge-invariant quantity (FWHM, \(|\tilde E(\omega)|^2\), the aligned temporal intensity), whose gradient is orthogonal to them. Any weakly identified physical direction that also falls below rcond is dropped too, which understates its uncertainty.

• Linearisation. It is a first-order (Laplace) Gaussian approximation about \(u^\*\). We propagate it by drawing parameter samples \(u_k\sim\mathcal N(u^\*,\mathrm{cov}(u))\) and pushing each through the exact (nonlinear) field → FWHM map, which captures the curvature of that map under the Gaussian parameter posterior. It is fast and accurate for well-constrained, high-SNR retrievals; near the transform limit or at low SNR the full croak.uncertainty bootstrap is the reference and this is the cross-check.

Because the covariance is rebuilt from the retrieved spectrum and the JAX forward model, it works for any result whose model it supports, not only one produced by an LM solver. It does not support the focal/collection models, depth weighting or marginal-peak centring, and it omits the regularisation curvature, so for a penalised fit it is the covariance of the unpenalised trace misfit.

class croak.covariance.CovarianceResult(cov, u, n_params, n_residuals, rank, null_dim, reduced_chisq, noise_scaled, sigma_thickness=None, sigma_tau0=None, sigma_smear=None, sigma_smear_delta=None)

Bases: object

Linearised parameter covariance of a least-squares retrieval.

Parameters:

• cov (NDArray[float64])

• u (NDArray[float64])

• n_params (int)

• n_residuals (int)

• rank (int)

• null_dim (int)

• reduced_chisq (float)

• noise_scaled (bool)

• sigma_thickness (float | None)

• sigma_tau0 (float | None)

• sigma_smear (float | None)

• sigma_smear_delta (float | None)

cov

Parameter covariance \(\mathrm{cov}(u)\) (n_params × n_params), the pseudo-inverse Gauss–Newton estimate \(\hat\sigma^2 (J^\mathsf{T}J)^+\).

Type:

numpy.ndarray

u

The parameters \(u^\*\) (the retrieved spectrum in the solver’s parameterisation) the covariance is taken about.

Type:

numpy.ndarray

n_params

Number of free parameters n (length of u).

Type:

int

n_residuals

Number of residual rows m (the FROG block, Nomega × Ndelay).

Type:

int

rank

Numerical rank of \(J^\mathsf{T}J\) (eigenvalues above rcond × max).

Type:

int

null_dim

Dimension of the gauge/flat null space, n_params - rank (≈ 2 for a pointwise full/phase parameterisation: CEP + absolute timing).

Type:

int

reduced_chisq

\(\lVert r(u^\*)\rVert^2/(m-n)\). When a noise model whitened the residual this is the reduced χ² (≈ 1 for a good fit and correct noise); otherwise it is the internal variance scale applied to the covariance.

Type:

float

noise_scaled

True when cov was scaled by the internal reduced-χ² estimate (no noise model given); False when an absolute noise model set \(\hat\sigma^2 = 1\).

Type:

bool

sigma_thickness

Standard error of the fitted dispersive-slab thickness (m), when fit_thickness was set; None otherwise. Already converted from the internal dimensionless offset to physical metres.

Type:

float or None

sigma_tau0

Standard error of the fitted delay-zero offset (s), when fit_tau0 was set; None otherwise.

Type:

float or None

sigma_smear

Standard error of the fitted geometric-smearing multiplier (dimensionless), when fit_smearing was set; None otherwise. For a split fit (fit_smearing_split) this is the gate-shape (p) channel.

Type:

float or None

sigma_smear_delta

Standard error of the delay-channel multiplier of a split smearing fit, when fit_smearing_split was set; None otherwise.

Type:

float or None

property stderr: NDArray[float64]

Per-parameter standard error \(\sqrt{\mathrm{diag}\,\mathrm{cov}}\).

croak.covariance.parameter_covariance(result, measured, *, noise=None, weights=None, material=None, thickness=0.0, npoints=20, omega0=None, quadrature='gausslegendre', phase_only=False, phase_basis='pointwise', n_nodes=20, R_omega=False, fit_thickness=False, fit_tau0=False, smearing=None, fit_smearing=False, fit_smearing_split=False, rcond=1e-10)

Gauss–Newton parameter covariance of a least-squares retrieval.

Rebuilds the JAX forward model and the solver parameterisation at the retrieved spectrum, forms the residual Jacobian \(J = \partial r/\partial u\) by automatic differentiation, and returns \(\mathrm{cov}(u) = \hat\sigma^2 (J^\mathsf{T}J)^{+}\) (see the module docstring). The forward-model and parameterisation keywords must reproduce the original retrieval (the same material/thickness/phase_basis …), exactly as the uncertainty bootstraps require their retrieve_kwargs.

Parameters:

• result (RetrievalResult) – A completed retrieval (provides the spectrum, grid, delays, interaction and carrier omega0).

• measured (array_like) – The measured trace (Nomega, Ndelay) that was fit.

• noise (NoiseModel or float, optional) – Measurement-noise model used to whiten the residual (each pixel divided by its σ), giving an absolutely scaled covariance (JᵀJ)⁺. A bare float is a fraction-of-peak additive NoiseModel. When omitted the covariance is scaled by the internal reduced-χ² estimate instead.

• weights (numpy.ndarray, optional) – Per-frequency weight vector (length Nomega); defaults to ones.

• material (str | None) – Dispersive forward-model parameters, as in croak.lm.LM (omega0 defaults to result.omega0).

• thickness (float) – Dispersive forward-model parameters, as in croak.lm.LM (omega0 defaults to result.omega0).

• npoints (int) – Dispersive forward-model parameters, as in croak.lm.LM (omega0 defaults to result.omega0).

• omega0 (float | None) – Dispersive forward-model parameters, as in croak.lm.LM (omega0 defaults to result.omega0).

• quadrature (str) – Dispersive forward-model parameters, as in croak.lm.LM (omega0 defaults to result.omega0).

• phase_only (bool) – Pulse parameterisation, as in croak._jax_pulse.ls_parameterisation().

• phase_basis (str) – Pulse parameterisation, as in croak._jax_pulse.ls_parameterisation().

• n_nodes (int) – Pulse parameterisation, as in croak._jax_pulse.ls_parameterisation().

• R_omega (bool, optional) – Per-frequency adaptive scaling (must match the retrieval).

• smearing (SmearingKernel or None, optional) – Geometric time-smearing kernel used by the forward model; must match the one the retrieval ran with, or the covariance describes a different model.

• fit_smearing (bool, optional) – Include the smearing multiplier in the augmented parameter vector, giving sigma_smear. Requires smearing.

• fit_smearing_split (bool, optional) – The retrieval fitted the two smearing channels separately (fit_smearing_split); adds the delay-channel multiplier to the vector, giving sigma_smear_delta. Requires fit_smearing.

• fit_thickness (bool, optional) – Include the dispersive-slab thickness and/or delay-zero offset in the covariance (must match the retrieval that fitted them). The fitted values are read from result; their standard errors are returned as sigma_thickness / sigma_tau0.

• fit_tau0 (bool, optional) – Include the dispersive-slab thickness and/or delay-zero offset in the covariance (must match the retrieval that fitted them). The fitted values are read from result; their standard errors are returned as sigma_thickness / sigma_tau0.

• rcond (float, optional) – Relative eigenvalue floor for the pseudo-inverse: eigenvalues of \(J^\mathsf{T}J\) below rcond × the largest are treated as flat directions and dropped (default 1e-10).

Returns:

CovarianceResult

Return type:

CovarianceResult

Notes

With noise supplied, the residual is divided by its per-pixel standard deviation. Currently weights then affect only the fitted intensity gain: they do not mask or weight the whitened residual. Without noise, the weights are included in the normalised residual. The degrees-of-freedom estimate is max(m - n, 1), not the numerical rank.

croak.covariance.covariance_uncertainty(result, measured, *, noise=None, propagation=None, propagation_label=None, n_samples=500, oversampling=8, collect_profiles=False, interval='bias-corrected', rng=None, weights=None, material=None, thickness=0.0, npoints=20, omega0=None, quadrature='gausslegendre', phase_only=False, phase_basis='pointwise', n_nodes=20, R_omega=False, fit_thickness=False, fit_tau0=False, smearing=None, fit_smearing=False, fit_smearing_split=False, rcond=1e-10, callback=None)

FWHM uncertainty from the linearised covariance (the analytic estimator).

Builds the Gauss–Newton parameter covariance (parameter_covariance()), draws n_samples parameter vectors \(u_k \sim \mathcal N(u^\*, \mathrm{cov}(u))\), optionally propagates each to a different beamline point with the transfer function propagation (see croak.session.dispersion.transfer_function()), and returns the spread of the per-sample temporal-intensity FWHM as an UncertaintyResult (method="covariance") — drop-in compatible with the bootstrap output, so the same intervals, temporal band and plotting apply.

Sampling from the local Gaussian parameter approximation and pushing each draw through the exact field → FWHM map is the delta method carried to the nonlinearity of that map: fast (no re-retrieval), and the natural analytic cross-check on the resampling bootstraps. Gauge null directions carry zero variance, so the draws do not jitter the CEP or timing.

result, measured, noise, weights and the forward-model / parameterisation keywords (material, thickness, npoints, omega0, quadrature, phase_only, phase_basis, n_nodes, R_omega, rcond) are as in parameter_covariance() and must reproduce the retrieval.

Parameters:

• propagation (numpy.ndarray of complex, optional) – Spectral transfer function \(H(\omega)\) propagating each draw to a different beamline point; None reports at the measurement plane.

• propagation_label (str, optional) – Description of the propagation point recorded on the result.

• n_samples (int, optional) – Number of Gaussian parameter draws B (default 500; cheap — each is a spectrum evaluation, no re-retrieval).

• oversampling (int, optional) – Temporal oversampling for the FWHM (default 8, matching process_result).

• collect_profiles (bool, optional) – Also collect peak-aligned oversampled intensity profiles (for a band).

• interval ({"bias-corrected", "percentile"}, optional) – Interval estimator.

• rng (numpy.random.Generator, optional) – Random generator (defaults to a fresh one).

• callback (callable, optional) – callback(done, total, R) after the draws are evaluated; for parity with the bootstrap API (called once here, the draws being a single vectorised step).

• result (RetrievalResult)

• measured (ArrayLike)

• noise (NoiseModel | float | None)

• weights (NDArray[float64] | None)

• material (str | None)

• thickness (float)

• npoints (int)

• omega0 (float | None)

• quadrature (str)

• phase_only (bool)

• phase_basis (str)

• n_nodes (int)

• R_omega (bool)

• fit_thickness (bool)

• fit_tau0 (bool)

• smearing (SmearingKernel | None)

• fit_smearing (bool)

• fit_smearing_split (bool)

• rcond (float)

Returns:

UncertaintyResult – With method="covariance".

Return type:

UncertaintyResult

Sessions (headless replay)

A whole retrieval — load, preprocess, retrieve, dispersion compensation and uncertainty — captured as serialisable per-stage parameters and replayed without the GUI. run_session() drives the same stage mappings the wizard uses; retarget() re-points an experimental session at a new dataset folder. SessionOptions also records which stage-1 loader the trace came from (entry, see ENTRIES), so a simulated session replays through assemble_simulated_load_data() rather than the measured loader — see Saving and loading.

Per-stage parameter dataclasses for a retrieval session.

These plain dataclasses describe a complete retrieval — loading, preprocessing, retrieval, dispersion compensation and uncertainty — independently of any user interface. The GUI (The wizard GUI) binds its controls to them and the headless engine (croak.session.engine) replays them; both share this one schema, so a session saved from the GUI can be rerun from a script and vice versa.

The dataclasses are deliberately Qt-free: they carry only data and a little classification logic (which load paths follow the dataset folder vs. which are fixed, used by croak.session.retarget).

class croak.session.params.LoadParams(frog_path='', ifrog_name='', lam_name='', scan_name='', transpose=False, lam_unit='nm', scan_unit='fs', scan_type='delay', reverse_trace=False, interaction='shg', third_order=False, third_order_exp=4.0, energy_j=0.0, use_spec=False, spec_path='', spec_is_csv=True, spec_lam_name='', spec_int_name='', spec_lam_col=0, spec_int_col=1, spec_lam_unit='nm', frog_bg_use=False, frog_bg_path='', frog_bg_lam_name='', frog_bg_int_name='', frog_calib_use=False, frog_calib_path='', frog_calib_lam_unit='nm', spec_bg_use=False, spec_bg_path='', spec_calib_use=False, spec_calib_path='', spec_calib_lam_unit='nm')

Bases: object

Stage 1: which files/datasets/units make up the measured trace.

Holds the FROG trace file and its wavelength/scan datasets and units, the trace orientation (transpose) and delay-sign (reverse_trace) conventions, the optional independent spectrum, the optional backgrounds and spectral-response calibration curves, and the optional measured pulse energy used to rescale the temporal plots to absolute power.

Parameters:

• frog_path (str)

• ifrog_name (str)

• lam_name (str)

• scan_name (str)

• transpose (bool)

• lam_unit (str)

• scan_unit (str)

• scan_type (str)

• reverse_trace (bool)

• interaction (str)

• third_order (bool)

• third_order_exp (float)

• energy_j (float)

• use_spec (bool)

• spec_path (str)

• spec_is_csv (bool)

• spec_lam_name (str)

• spec_int_name (str)

• spec_lam_col (int)

• spec_int_col (int)

• spec_lam_unit (str)

• frog_bg_use (bool)

• frog_bg_path (str)

• frog_bg_lam_name (str)

• frog_bg_int_name (str)

• frog_calib_use (bool)

• frog_calib_path (str)

• frog_calib_lam_unit (str)

• spec_bg_use (bool)

• spec_bg_path (str)

• spec_calib_use (bool)

• spec_calib_path (str)

• spec_calib_lam_unit (str)

reverse_trace: bool = False

negate the scan axis on loading, i.e. mirror the trace about τ = 0. Which end of a stage scan is “gate late” is a wiring convention of the measurement, and a measured file carries no marker recording it (unlike SimulatedLoadParams.reverse_trace, which can auto-detect from /grid/delay_convention), so this is a plain off-by-default override. It matters for PG/SD, whose kernels fix the direction of time: getting it wrong there time-reverses the retrieved pulse and flips the sign of every phase order. A thin-medium SHG trace is symmetric in delay (Interaction.time_reversal_ambiguous), so reversing it changes nothing and cannot settle SHG’s direction-of-time ambiguity — that is croak.processing.resolve_time_direction()’s job. See Conventions.

Type:

Flip the delay-sign convention

static dataset_path_fields()

Path fields that live inside the per-dataset folder.

These move together when an retarget() points an existing session at a new dataset folder: the measured trace, the independent spectrum, and their backgrounds all share that folder and keep their filenames across acquisitions.

Return type:

tuple[str, …]

static fixed_path_fields()

Path fields that stay put when retargeting (shared calibration files).

Return type:

tuple[str, …]

class croak.session.params.SimulatedLoadParams(frog_path='', window_key='Iω_win', z_index=-1, z_thickness_um=None, lam_min_nm=140.0, lam_max_nm=800.0, third_order=True, third_order_exp=2.0, reverse_trace=None, uniform_delay_core=False, interaction='pg', use_spectrum=True, spectrum_source='beamlet', truth_source='beamlet', raw_direct=False)

Bases: object

Stage 1 (alternative): load a numerically simulated FROG scan.

Describes how to turn a Luna scansave HDF5 file (see croak.io.read_simulated_scan()) into the same load_data dict the measured-trace loader produces. Unlike LoadParams the file layout is fixed, so there are no dataset-name/unit pickers — only the physics options: which trace window, the loaded wavelength band, the third-order efficiency correction, the delay-sign convention, the nonlinear interaction, and which independent spectrum (if any) to feed the retrieval.

Notes

The simulated trace is stored as an angular-frequency density \(I_\omega\); the loader divides it by \(\lambda^2\) so the wavelength pipeline’s \(\times\lambda^2\) regrid recovers it undistorted (see croak.preprocess.omega_to_lambda_density()).

The chromatic vignetting imprinted by the physical input mask is a linear filter on the gating beam, not a trace artefact, so it belongs in the reference spectrum: spectrum_source="beamlet" (the post-mask, transverse-integrated grid/Iω_beamlet) constrains the retrieval with the spectrum that actually gates. Choose the independent spectrum and the truth in the same field reference as the model: an on-axis model (the focal mixture) needs the on-axis beamlet, or a frame conversion via RetrieveParams.spectrum_frame_p. Detector response, the radiated source factor and wavelength/frequency density conversion are separate operations. See Loading numerically simulated traces and Conventions.

Parameters:

• frog_path (str)

• window_key (str)

• z_index (int)

• z_thickness_um (float | None)

• lam_min_nm (float)

• lam_max_nm (float)

• third_order (bool)

• third_order_exp (float)

• reverse_trace (bool | None)

• uniform_delay_core (bool)

• interaction (str)

• use_spectrum (bool)

• spectrum_source (str)

• truth_source (str)

• raw_direct (bool)

window_key: str = 'Iω_win'

"Iω_win" (integrated through the selected aperture) or "Iω_win_reimaged" (on-axis).

z_index: int = -1

Propagation slice index (-1 = the last saved slice, i.e. substrate exit); the low-level/back-compat knob, used only when z_thickness_um is None.

z_thickness_um: float | None = None

Material thickness (µm) to load from a multi-thickness scan. When set, the nearest saved /grid/zsave slice is loaded (takes precedence over z_index); None falls back to z_index (the exit slice), which also covers older single-thickness files that lack zsave.

uniform_delay_core: bool = False

Keep only the longest uniformly spaced contiguous stretch of the delay axis (croak.preprocess.uniform_core_indices()). Simulated scans often extend a uniform core with coarser wing points (e.g. 1 fs-step wings to ±40 fs for the Raman wake); the smearing kernel’s delay-axis convolution needs a uniform grid, so retrievals that model smearing must drop the wings, and so must delay_origin="marginal_peak", whose Fourier re-centring also needs uniform delays. Off by default: the wings are real data that the other models can use.

interaction: str = 'pg'

"pg" (transient-grating/PG) or "sd"; both keep the carrier.

spectrum_source: str = 'beamlet'

"beamlet" (post-mask grid/Iω_beamlet, the transverse-integrated beam that gates — the 1-D models’ frame), "beamlet_reimaged" (the ON-AXIS beamlet grid/Iω_beamlet_reimaged, stored by pnps files — the chromatic focal mixture’s own frame, so no spectrum_frame_p is then needed) or "source" (pre-mask grid/Iω, the ideal input). Falls back along reimaged → beamlet → source when the file lacks the requested one.

Type:

which spectrum when use_spectrum

truth_source: str = 'beamlet'

which stored pulse to overlay as the time-domain truth (and to seed truth_init from): "beamlet" (post-mask grid/It_beamlet, the beam that gates — default), "beamlet_reimaged" (its on-axis counterpart, grid/It_beamlet_reimaged — the like-for-like truth for a focal-mixture retrieval) or "source" (pre-mask grid/It, the ideal input). Falls back when the file stores fewer. See croak.io.read_simulated_truth_keys().

raw_direct: bool = False

load the as-simulated ω-density trace on its native grid straight into the retrieval, skipping the marginal-check and preprocess (filtering + regrid) stages. Fully raw: no third-order correction; window/thickness/delay-sign/ known-spectrum selections are still honoured. For sanity-checking the retrieval free of any preprocessing doubt (see croak.session.pipeline.assemble_simulated_tracedata()).

class croak.session.params.PreprocParams(lam_min_nm=700.0, lam_max_nm=900.0, lam_bound_min_nm=350.0, lam_bound_max_nm=1100.0, trange_fs=1000.0, tau_min_fs=-500.0, tau_max_fs=500.0, use_subsample=False, subsample_tau=1, resample_t=False, prefilter=True, lamm_min_nm=300.0, lamm_max_nm=1100.0, lamm_bound_min_nm=300.0, lamm_bound_max_nm=1100.0, taum_min_fs=-500.0, taum_max_fs=500.0, tau_bound_min_fs=-1000.0, tau_bound_max_fs=1000.0, filter_fringes=False, filter_dc=False, dtau_min_fs=0.0, defringe=False, defringe_fraction=0.5, defringe_on_alias='clamp', baseline=False, baseline_smoothness=100.0, baseline_ratio=0.001, baseline_clip=False, baseline_tau_exclude_fs=0.0, use_threshold=False, threshold=0.0001, marginal_correct=False)

Bases: object

Stage 2: the grid, windowing and filtering used to clean and regrid.

Wavelengths are in nm and delays in fs (converted to SI when applied). The grid fields set the retrieval ω/τ grid; the windowing fields the measurement λm/τm extents; the filtering fields the fringe/DC/low-pass and thresholding.

Parameters:

• lam_min_nm (float)

• lam_max_nm (float)

• lam_bound_min_nm (float)

• lam_bound_max_nm (float)

• trange_fs (float)

• tau_min_fs (float)

• tau_max_fs (float)

• use_subsample (bool)

• subsample_tau (int)

• resample_t (bool)

• prefilter (bool)

• lamm_min_nm (float)

• lamm_max_nm (float)

• lamm_bound_min_nm (float)

• lamm_bound_max_nm (float)

• taum_min_fs (float)

• taum_max_fs (float)

• tau_bound_min_fs (float)

• tau_bound_max_fs (float)

• filter_fringes (bool)

• filter_dc (bool)

• dtau_min_fs (float)

• defringe (bool)

• defringe_fraction (float)

• defringe_on_alias (str)

• baseline (bool)

• baseline_smoothness (float)

• baseline_ratio (float)

• baseline_clip (bool)

• baseline_tau_exclude_fs (float)

• use_threshold (bool)

• threshold (float)

• marginal_correct (bool)

class croak.session.params.RetrieveParams(solver='warm-lbfgs', maxiters=300, reltol=1e-05, abstol=1e-08, full=True, R_omega=True, random_init=False, perfect_init=False, truth_init=False, tl_init=False, reuse_result=False, live_preview=True, perfect_init_fs=10.0, dispersive=False, material='SiO2', thickness_um=10.0, npoints=10, delta_z_um=1.0, smearing=False, smear_hole_diameter_mm=1.0, smear_hole_spacing_mm=0.5, smear_sigma_delay_fs=0.357, copra_alpha=0.25, copra_stall_patience=5, smear_nodes=5, smear_scale=1.0, smear_layout='boxcars', focal=False, focal_n_radial=24, focal_n_azimuth=16, focal_rmax_units=4.0, focal_f_mm=100.0, collection='off', collection_mode='integrated', collection_diam_mm=0.5, spectrum_frame_p=0.0, envelope=False, envelope_zr_um=0.0, envelope_lw_um=0.0, envelope_zf_um=0.0, envelope_gouy=1.0, fit_thickness=False, fit_tau0=False, fit_smearing=False, fit_smearing_split=False, polish_two_phase=False, delay_origin='coincidence', tau0_bound_fs=0.0, reg_amp=0.03, reg_phase=0.0, reg_spectrum=0.01, reg_time=0.0, reg_time_lo_fs=-500.0, reg_time_hi_fs=500.0, phase_basis='pointwise', n_nodes=20, cma_strategy='cma', cma_popsize=0, cma_std_init=0.5, post_filter_lam=False, post_filter_tau=False)

Bases: object

Stage 3: the solver, initial guess, dispersive slab and regularisation.

Covers the algorithm and its tolerances, the guess mode, optional dispersive propagation through a slab (the forward model, distinct from the stage-4 dispersion compensation), the regularisation weights, the phase parameterisation, the CMA-ES global-search knobs and the post-filter toggles.

Parameters:

• solver (str)

• maxiters (int)

• reltol (float)

• abstol (float)

• full (bool)

• R_omega (bool)

• random_init (bool)

• perfect_init (bool)

• truth_init (bool)

• tl_init (bool)

• reuse_result (bool)

• live_preview (bool)

• perfect_init_fs (float)

• dispersive (bool)

• material (str)

• thickness_um (float)

• npoints (int)

• delta_z_um (float)

• smearing (bool)

• smear_hole_diameter_mm (float)

• smear_hole_spacing_mm (float)

• smear_sigma_delay_fs (float)

• copra_alpha (float)

• copra_stall_patience (int)

• smear_nodes (int)

• smear_scale (float)

• smear_layout (str)

• focal (bool)

• focal_n_radial (int)

• focal_n_azimuth (int)

• focal_rmax_units (float)

• focal_f_mm (float)

• collection (str)

• collection_mode (str)

• collection_diam_mm (float)

• spectrum_frame_p (float)

• envelope (bool)

• envelope_zr_um (float)

• envelope_lw_um (float)

• envelope_zf_um (float)

• envelope_gouy (float)

• fit_thickness (bool)

• fit_tau0 (bool)

• fit_smearing (bool)

• fit_smearing_split (bool)

• polish_two_phase (bool)

• delay_origin (str)

• tau0_bound_fs (float)

• reg_amp (float)

• reg_phase (float)

• reg_spectrum (float)

• reg_time (float)

• reg_time_lo_fs (float)

• reg_time_hi_fs (float)

• phase_basis (str)

• n_nodes (int)

• cma_strategy (str)

• cma_popsize (int)

• cma_std_init (float)

• post_filter_lam (bool)

• post_filter_tau (bool)

solver: str = 'warm-lbfgs'

“warm-lbfgs” is the extended model seeded from COPRA’s local projection sweep: the same forward model as “lbfgs-ad” with a better starting point. It is the default because it is close to neutral where the landscape is benign and substantially better where it is not (see croak.warm_lbfgs for the measured comparison), so it is the safer choice for a user who does not want to reason about initial guesses. The papers’ validation protocol instead runs “lbfgs-ad” from a transform-limited start (tl_init) for 400 iterations.

maxiters: int = 300

Maximum solver work budget; its unit depends on the algorithm (iterations, residual evaluations or objective evaluations). 300 and a 1e-5 relative tolerance are the settings our smearing and solver comparisons were measured with; 100 and 1e-4 stopped some dispersive retrievals short of convergence.

truth_init: bool = False

Seed from the known complex field carried by a synthetic or ModelPNPS load. The selected truth must contain complex spectral datasets; intensity- only files cannot provide this initial condition.

tl_init: bool = False

its transform limit. Preferred over perfect_init for near-transform-limited pulses: the two perform equivalently, but this one uses the amplitude that was actually measured and assumes no duration. Ignored if random_init or perfect_init is set.

Type:

Start from the measured spectrum with flat phase

copra_alpha: float = 0.25

COPRA global-step size factor (Geib’s alpha, dimensionless). Inert in practice: on our dispersive TG-FROG tests a factor of five in it moved the complex error by a few percent, non-monotonically.

copra_stall_patience: int = 5

Iterations without improvement before COPRA’s local projection sweep hands over to its global descent. Not a smooth knob but a switch between the two stages, and with a smeared model the stages differ in whether they can represent the physics at all: the local projection cannot (a mixture has no single signal to project onto) while the global descent can. A very large value never leaves the local sweep.

smear_scale: float = 1.0

Fixed multiplier on both reduced-kernel widths (dimensionless), applied as a FIXED value (unlike fit_smearing, which optimises it). It exists because the mask geometry can only express d/D = (gap + D)/2D >= 0.5, so a kernel narrower than that is not reachable by adjusting the gap at all.

Use 1.0 for the apertured instrument: in our earlier scans of the kernel width against the 3D instrument simulation at the known truth, the optimum sat at 0.85-1.00 at both 9.5 and 20 um and at both d/D = 1.0 and 1.5. Gaussian-beam traces wanted about 0.5, on top of the 0.748 focal-weight rescale their effective gap already carries. A fitted or chosen value tests sensitivity; it need not identify a single physical mechanism.

smear_layout: str = 'boxcars'

Aperture layout the kernel is built from. "boxcars" is the four-hole BOXCARS mask used for TG/PG (and for a three-hole SD). "sd2" is a genuine TWO-beam self-diffraction: two collinear apertures, the probe entering twice and the gate once with conjugation.

The distinction matters. In a two-beam SD the two unconjugated arms are the SAME beam, so their relative internal delay p (proportional to the difference of their tilts) is identically zero and the smearing collapses to a pure delay-axis convolution. Building such a measurement’s kernel from the boxcars closed form instead attaches a spurious sigma_p of order 0.4 fs that the experiment does not have.

focal: bool = False

Chromatic focal mixture (croak.focal) instead of the reduced smearing kernel — the full model. Each focal-plane node carries the chromatic Airy filter A(r, omega) and its own definite (p, theta), so the r-dependent (bluer-on-axis) local spectra that the reduced kernel averages away are modelled explicitly. Mutually exclusive with smearing (the mixture REPLACES the kernel; asking for both raises). Autodiff solvers only. Roughly 40x the reduced kernel’s cost per iteration at the default 24x16 node grid.

focal_n_azimuth: int = 16

the coherent (apertured) sum needs the finer azimuthal quadrature (8 costs 4e-5 of trace peak). Refine with the radial and collection quadratures to check convergence.

Type:

Number of azimuthal nodes. 16, not 8

focal_f_mm: float = 100.0

Focusing focal length (mm). The reduced kernel never needs it (f cancels between spot size and tilt); the explicit mixture does. For simulated scans it is stored in the file (f_foc when available) and adopted with the rest of the geometry.

collection: str = 'off'

"off" sums incoherently over the focal plane (exact only for full-beam collection), "file" rebuilds the collection aperture from a simulated scan’s own window_def_* record, "manual" builds a centred hole at the phase-matched corner from collection_diam_mm. A tightly apertured instrument sits far from the incoherent limit: the finite aperture changes the spatial interference, and without this the mixture over-blurs.

Type:

Collection model riding on the mixture

spectrum_frame_p: float = 0.0

the measured spectrum used for the initial guess and the reg_spectrum target is reweighted in intensity by ((omega + omega0)/omega0)**(2*p) (absolute over central frequency) before retrieval. The focal mixture’s field is the ON-AXIS focal field, while a measured/beamlet spectrum is spatially integrated; in the declared diffraction convention p = 1 converts the integrated aperture-diffraction beamlet to the on-axis entrance field. Opening the collection reduces the effective value (in our earlier scans, ~1.0 at a 0.5 mm gap and ~0.5 at 1.0 mm). Leave 0 for the reduced kernel, whose implicit frame is the mixture-averaged one. GDD and spectral span are insensitive to it; durations move by a few percent. See Conventions.

Type:

Spectral-frame exponent p

delay_origin: str = 'coincidence'

"coincidence" (gate–probe coincidence, the forward model’s native convention) or "marginal_peak" (every model trace re-centred on its own row-gain weighted delay-marginal peak, the convention the preprocessing gives the data, with no free parameter; needs uniform delays). The latter is the single-cycle choice: a coherently collected trace peaks tens of attoseconds from coincidence, growing with the chirp, so a fixed origin costs duration accuracy (~10 % in our earlier instrument scans) while a free fit_tau0 wanders on chirped pulses. See croak.delay_origin.

Type:

Where the model’s delay zero sits (AD solvers)

tau0_bound_fs: float = 0.0

|τ0| ≤ tau0_bound_fs (fs) for a fitted delay offset; 0 = unbounded.

reg_amp: float = 0.03

Amplitude-smoothness weight (second difference of |E(omega)|, max-normalised). Used together with reg_spectrum: on our dispersive TG-FROG tests, 0.03 alongside a spectrum weight of 0.01 beat every single-penalty setting on every metric at once, across eight arms at two depths: 12 % better in complex error, 2.4 points better in retrieved duration and 2.5x smoother in the worst case than reg_spectrum = 0.1 alone. Adding it never hurt in any configuration measured, noiseless or noisy. The same objective-scale caveat as reg_spectrum applies (the LM family uses a different scale).

reg_spectrum: float = 0.01

Spectrum-divergence weight (penalises mismatch between the retrieved and independently measured spectral amplitude). 1e-2 is the gradient-family default: on our dispersive TG-FROG tests, scanning it over three decades moved the trace error by only 14 % but the pulse error by 44 %, with a clear interior optimum. The penalty does not remove misfit so much as move it out of the spectral amplitude, where it is constrained, and into the phase, where it is not, so both too small a weight (amplitude error dominates) and too large a one (phase error dominates) are worse. Do NOT tune this by watching the trace error: it barely responds.

The papers’ validation protocol uses 0.03 for both this and reg_amp (with “lbfgs-ad” from a transform-limited start, 400 iterations).

This value is for the L-BFGS/gradient family, which minimises the scalar R + lambda P. The LM solvers add the penalty as residual rows (R² + lambda’ P), so the same number is far stronger there; see croak.solver.REG_DEFAULTS for the per-family values and the lambda’ = 2 R lambda mapping between them. croak does not translate automatically, so reset both weights when changing solver to “lm”/”lm-optx” (the GUI does this for you).

class croak.session.params.BeamPathMirror(name='', polarization='s', bounces=0, direction='remove', custom_file='', custom_lam_col=0, custom_val_col=1, custom_unit='nm', custom_datatype='phase')

Bases: object

One mirror in the dispersion stage’s beam-path stack (a serialisable row).

A stable description of a single mirror reflection, so the stack survives a save/reload (unlike the live registry keys the GUI uses at runtime). A built-in mirror is named (name in croak.mirrors.BUILTIN_COATINGS or croak.mirrors.BUILTIN_MIRRORS); a custom mirror leaves name empty and points custom_file at a CSV of wavelength + phase/GDD.

Parameters:

• name (str) – Built-in coating/mirror name; empty for a custom-file mirror.

• polarization ({"s", "p"}) – Polarisation for coatings (ignored by phase-only mirrors).

• bounces (int) – Number of reflections (0 disables the row).

• direction ({"remove", "add"}) – "remove" back-propagates the mirror out of the pulse (negative bounce count); "add" propagates it in.

• custom_file (str) – CSV path for a custom mirror (empty for a built-in).

• custom_lam_col (int) – 0-based wavelength and value column indices in custom_file.

• custom_val_col (int) – 0-based wavelength and value column indices in custom_file.

• custom_unit (str) – Wavelength unit of custom_file (e.g. "nm").

• custom_datatype ({"phase", "gdd"}) – Whether the value column is spectral phase (rad) or GDD (fs²).

class croak.session.params.DispersionParams(gdd_fs2=0.0, tod_fs3=0.0, fod_fs4=0.0, gas_name='', gas_pressure_bar=1.0, gas_path_cm=0.0, material_thickness_mm=<factory>, mirrors=<factory>)

Bases: object

Stage 4: post-retrieval dispersion compensation applied to the pulse.

Taylor coefficients (fs², fs³, fs⁴), an optional active gas (pressure + path), built-in material thicknesses (mm by name), and the beam-path mirror stack. This stage tunes/compresses the retrieved pulse; it does not re-run retrieval.

Notes

Scalar fields are declared before the material_thickness_mm table and the mirrors array-of-tables because TOML requires a table’s scalar keys to precede its sub-tables.

Parameters:

• gdd_fs2 (float)

• tod_fs3 (float)

• fod_fs4 (float)

• gas_name (str)

• gas_pressure_bar (float)

• gas_path_cm (float)

• material_thickness_mm (dict[str, float])

• mirrors (list[BeamPathMirror])

class croak.session.params.UncertaintyParams(method='parametric', n_resamples=200, noise_source='residual', noise_sigma=0.01, warm_start=True, collect_profiles=True, interval='bias-corrected', thickness_sigma_um=0.0, propagate_to_dispersion_point=False)

Bases: object

Stage 5: the FWHM-uncertainty bootstrap configuration.

The method selects a statistical bootstrap — "parametric" (add noise and re-retrieve) or "delay"/"frequency" (resample the trace) — the substrate-thickness systematic "thickness" (Monte-Carlo over the measured thickness prior), or the fast analytic "covariance" (linearised Gauss–Newton covariance; see croak.covariance.covariance_uncertainty()). The remaining fields tune the resampling, noise model and confidence interval.

Parameters:

• method (str)

• n_resamples (int)

• noise_source (str)

• noise_sigma (float)

• warm_start (bool)

• collect_profiles (bool)

• interval (str)

• thickness_sigma_um (float)

• propagate_to_dispersion_point (bool)

croak.session.params.reg_defaults(solver)

Return the default spectrum and amplitude penalty weights.

(reg_spectrum, reg_amp) suited to solver; see REG_DEFAULTS.

Parameters:

solver (str)

Return type:

tuple[float, float]

croak.session.params.reg_family(solver)

Which objective solver minimises: "gradient" (R) or "lm" (R^2).

Regularisation weights are only comparable within a family.

Parameters:

solver (str)

Return type:

str

The SessionOptions bundle: the serialisable session.

SessionOptions groups the per-stage parameter dataclasses, and the name of the stage-1 loader that produced the trace, into one object that round-trips to the options.toml layout (one TOML table per stage). It is the single source of truth for session serialisation — the GUI’s WizardState and the headless engine both go through it — so a session saved in one place reloads identically in the other.

croak.session.options.ENTRIES: tuple[str, ...] = ('experimental', 'simulated', 'synthetic')

Which loader produced the trace — i.e. which of the mutually exclusive stage-1 entries the session started from. "experimental" reads [load], "simulated" reads [simulated]; "synthetic" has no serialised parameters yet (the generator’s settings live only in its GUI controls) and is recorded so a reload can say so rather than silently loading nothing.

class croak.session.options.SessionOptions(entry='experimental', load=<factory>, simulated=<factory>, preproc=<factory>, retrieve=<factory>, dispersion=<factory>, uncertainty=<factory>)

Bases: object

A complete retrieval session: which loader it used and the stage parameters.

Round-trips to the options.toml layout (one TOML table per stage, plus a scalar entry key naming the loader) via to_dict()/from_dict(). The dispersion table carries a nested material_thickness_mm sub-table and a mirrors array-of-tables, which from_dict() rebuilds into BeamPathMirror.

load and simulated are alternatives, selected by entry (see ENTRIES): a measured session fills [load] and leaves [simulated] at its defaults, and vice versa. Both are always written, so the file records which trace the retrieval actually came from — without entry a simulated session would look like an empty experimental one.

Notes

The recipe refers to the external input and calibration files rather than embedding them, and records no retrieval RNG seed; synthetic-generator settings are not serialised for headless replay. See Saving and loading.

Parameters:

• entry (str)

• load (LoadParams)

• simulated (SimulatedLoadParams)

• preproc (PreprocParams)

• retrieve (RetrieveParams)

• dispersion (DispersionParams)

• uncertainty (UncertaintyParams)

to_dict()

Serialise to the {entry, section: {key: value}} options layout.

entry comes first because TOML requires scalar keys to precede the tables they share a level with.

Return type:

dict[str, Any]

classmethod from_dict(data)

Build from a nested options mapping (version header keys are ignored).

A file written before entry existed has no simulated workflow to restore, so it reads back as an experimental session.

Parameters:

data (dict[str, Any])

Return type:

SessionOptions

to_toml(path)

Write the session to path as an options TOML (with version header).

Parameters:

path (str)

Return type:

str

classmethod from_toml(path)

Read a session from an options.toml file.

Parameters:

path (str)

Return type:

SessionOptions

with_(**changes)

Return a copy with the given per-stage params replaced (e.g. load=).

Parameters:

changes (Any)

Return type:

SessionOptions

Pure stage mappings: parameter dataclasses → core croak calls.

Each function turns one stage’s croak.session.params dataclass into the corresponding call on croak’s public API (assemble_load_data → arrays, build_tracedata → load_and_clean(), run_retrieval → retrieve_from_tracedata(), run_uncertainty → estimate_fwhm_uncertainty()). They are Qt-free and side-effect-light (only the loaders read files), so the GUI workers and the headless engine call exactly the same mapping — there is one obvious way to turn a saved session into a retrieval.

croak.session.pipeline.assemble_load_data(p)

Build a load_data dict (SI arrays) from the load parameters.

Returns {lam, scanaxis, input_unit, trace (Nlambda, Ndelay), lam_spec, Ilam_spec, interaction}. The scale curve (calibration × 3rd-order) is applied eagerly to the trace, so everything downstream sees calibrated counts and the calibration never has to be re-applied. With p.reverse_trace the scan axis is negated first, mirroring the trace about τ = 0 (the measured counterpart of SimulatedLoadParams’s reverse_trace).

Parameters:

p (LoadParams)

Return type:

dict

croak.session.pipeline.assemble_simulated_load_data(p)

Build a load_data dict from a numerically simulated FROG scan.

Reads a Luna scansave HDF5 file (croak.io.read_simulated_scan()) and returns the same {lam, scanaxis, input_unit, trace, lam_spec, Ilam_spec, interaction} shape as assemble_load_data(), so the simulated trace flows into build_tracedata() and the retrieval exactly like a measured one.

For a multi-thickness scan (one that stores /grid/zsave), p selects which propagation thickness to load: z_thickness_um (µm) picks the nearest saved slice, and falls back to the substrate-exit slice when None (also the only behaviour for the older single-thickness format).

The simulated trace is an angular-frequency density \(I_\omega\); the measured-data pipeline (regrid) multiplies the loaded trace by \(\lambda^2\). To land on the correct \(I_\omega\) the trace is divided by \(\lambda^2\) here (croak.preprocess.omega_to_lambda_density()), so the round-trip is undistorted. The pipeline is: read → drop ω ≤ 0 and crop to the band → optional ×(λ/µm)^n third-order correction → ÷ λ² → optional delay reversal → normalise. The chosen reference spectrum (post-mask beamlet by default, or the pre-mask source) is offered as the independent spectrum; the mask vignetting is handled by that choice, not by a trace correction. The focal-mixture models need the on-axis reference (beamlet_reimaged) or a frame conversion; see apply_spectrum_frame().

Parameters:

p (SimulatedLoadParams) – The simulated-load options (file, window, band, corrections, …).

Returns:

dict – load_data with SI arrays; trace is (Nlambda, Ndelay) and input_unit is "delay".

Raises:

ValueError – If the requested wavelength band contains no positive-frequency bins.

Return type:

dict

croak.session.pipeline.build_tracedata(p, data, *, fast_baseline=False)

Clean and regrid a loaded trace into a TraceData.

Maps the (nm/fs) preprocess parameters onto croak.preprocess.load_and_clean() SI arguments.

fast_baseline (off by default) solves the arPLS baseline on a coarse delay grid and skips near-zero rows — a ~3× speedup at ~1 % baseline difference, for the GUI’s interactive preview. Headless replay and direct callers get the exact, full-resolution baseline.

Parameters:

• p (PreprocParams)

• data (dict)

• fast_baseline (bool)

Return type:

TraceData

croak.session.pipeline.resolve_guess(p, override=None)

Pick the initial-guess mode (or an explicit spectrum) from the parameters.

Initial-condition flags are mutually exclusive. An explicit override supplies the spectrum for reuse_result; without a compatible previous result that GUI-only choice falls back to "auto".

Parameters:

p (RetrieveParams)

croak.session.pipeline.extra_param_centres(p, seed)

Return the starting centres (thickness, tau0, smear_scale, smear_delta).

Without a seed these are the nominal settings: the slab thickness from the parameters, a zero delay-zero offset and an unscaled smearing kernel.

With a seed (the previous retrieval, when reusing it) each extra that is being fitted again starts from the value that fit found, so re-running is a true warm start: the reused spectrum was optimised jointly with those extras, so restarting them from their priors would re-fit it against, e.g., a differently-referenced delay axis. An extra whose fit_* flag is off keeps the explicit nominal value instead — an unticked box means “hold this at what I set”, not “hold it at whatever the last fit drifted to”.

Parameters:

• p (RetrieveParams) – Retrieval settings (supplies the nominal thickness and the fit flags).

• seed (RetrievalResult or None) – Previous result to warm-start from, or None for a cold start.

Returns:

tuple of float – Thickness (m), delay-zero offset (s), smearing-width multiplier, and the delay-channel multiplier of a split smearing fit (fit_smearing_split; equal to the joint multiplier’s default 1.0 otherwise).

Return type:

tuple[float, float, float, float]

croak.session.pipeline.run_retrieval(p, td, *, guess_override=None, extras_seed=None, truth=None, rng=None, callback=None, focal=None, scan_path=None, scan_window=None)

Run retrieve_from_tracedata() from a RetrieveParams.

Parameters:

• p (RetrieveParams) – Retrieval settings (solver, regularisation, dispersive slab, …).

• td (TraceData) – Cleaned, regridded trace.

• guess_override (array_like or Pulse, optional) – Explicit initial spectrum (reuse-previous-result); overrides the modes.

• extras_seed (RetrievalResult, optional) – Previous result to warm-start the fitted extra parameters from — passed together with guess_override when reusing a result. See extra_param_centres().

• truth (TruthPulse, optional) – Known pulse from the load stage. Required when p.truth_init is set; its complex spectrum is aligned to td and used as the solver seed.

• rng (numpy.random.Generator, optional) – Generator for the random guess / solver; a fresh default is used if None.

• callback (callable, optional) – callback(iteration, R, best_R, snapshot=...) progress hook; the optional snapshot is a zero-argument in-progress-result builder (or None) for a live full-plot preview.

• focal (FocalMixture, optional) – Explicit chromatic focal-plane mixture (optionally with a collection aperture), forwarded to retrieve_from_tracedata(). When None (the default) and p.focal is set, the mixture is built from the params by focal_mixture_from_params() — the route the GUI uses; pass an explicit object to override it from a script. The mixture replaces the reduced kernel, so p.smearing must be off either way.

• scan_path (str, optional) – Path of the loaded scan file, needed only for p.collection == "file" (the aperture is rebuilt from the file’s own window_def_* record).

• scan_window (str or int, optional) – Which collection hole of a multi-aperture scan the trace came from (the loader’s window_key); with p.collection == "file" the aperture is rebuilt from that hole’s record. None = the first window.

Return type:

RetrievalResult

croak.session.pipeline.apply_post_filter(p, result, td)

Window the result to the measurement λ/τ ranges if the toggles are set.

Mirrors the retrieve stage’s post-filter pane: a spectral and/or temporal window to the measurement extents (td.lamm_lims / td.tau_lims), not the full grid. Returns the result unchanged when neither toggle is on.

Parameters:

• p (RetrieveParams)

• result (RetrievalResult)

• td (TraceData)

Return type:

RetrievalResult

croak.session.pipeline.process(result, td, *, energy=None)

Process a result against its measured trace (FWHM, dispersion fit, marginals).

energy is an optional measured pulse energy (J); when given, the result carries absolute peak power (see croak.processing.process_result()).

Parameters:

• result (RetrievalResult)

• td (TraceData)

• energy (float | None)

Return type:

ProcessedResult

croak.session.pipeline.bootstrap_solver_kwargs(p, td)

Solver kwargs for the uncertainty bootstrap (matching the retrieval).

estimate_fwhm_uncertainty() rebuilds the solver internally and filters these by what the chosen algorithm accepts, so it is safe to pass the full candidate set.

Parameters:

• p (RetrieveParams)

• td (TraceData)

Return type:

dict

croak.session.pipeline.smearing_kernel(p, td)

Build the stage-4 geometric-smearing kernel, or None when it is off.

The mask geometry is the stored truth (the displayed delay width is a derived GUI convenience), and the carrier comes from the measurement itself, so the kernel is always evaluated at the wavelength the trace was actually taken at.

Parameters:

• p (RetrieveParams)

• td (TraceData)

Return type:

SmearingKernel | None

croak.session.pipeline.focal_mixture_from_params(p, td, scan_path=None, scan_window=None)

Build the chromatic focal mixture from RetrieveParams, or None.

scan_window names which collection hole of a multi-aperture scan the trace was loaded from (the window_key, e.g. "Iω_win_5"), so that collection="file" rebuilds that hole’s record rather than the first one’s; None keeps the first window, which is right for single-window files.

The focal-model counterpart of smearing_kernel(): the mask geometry comes from the same smear_hole_* fields (it is the same mask), the carrier from the measurement, and the focusing focal length — which the reduced kernel never needs but the explicit mixture does — from focal_f_mm. The collection model rides on the mixture per p.collection:

• "off": incoherent full-beam sum (exact only for full collection);

• "file": the aperture rebuilt from the simulated scan’s own window_def_* record (needs scan_path);

• "manual": a hole of collection_diam_mm at the phase-matched BOXCARS corner (-d, -d), d = (gap + D)/2, hard-edged — the physical pinhole of an experiment. ("file" reproduces a simulated scan’s recorded edge and weighting convention; a soft quadrature edge there is a grid-resolution device, not a physical transmission.)

Raises:

ValueError – If p.smearing is also set (the mixture replaces the reduced kernel), if the interaction is not PG (the mixture’s arm layout is the four-hole BOXCARS TG/PG one), if collection="file" without a scan_path, or on an unknown collection mode.

Parameters:

• p (RetrieveParams)

• td (TraceData)

• scan_path (str | None)

• scan_window (str | int | None)

Return type:

FocalMixture | None

croak.session.pipeline.apply_spectrum_frame(p, td)

Reweight an independent spectrum into a declared retrieval field reference.

Multiply td.Iomega by ((omega + omega0)/omega0)**(2*p) with p=spectrum_frame_p. The initial amplitude and spectral-penalty target therefore move together into the frame the forward model’s field lives in. The focal mixture’s field is the ON-AXIS focal field, bluer than a spatially integrated spectrum; in the stated diffraction convention p=1 converts an integrated beamlet to the on-axis entrance field, and p=0.5 is an empirical effective convention. The reduced kernel’s implicit frame is the integrated one, so leave the exponent at 0 there. The function does not modify a stored complex truth automatically. See Conventions for like-for-like comparisons. Zero exponent or an absent independent spectrum leaves the input unchanged.

Parameters:

• p (RetrieveParams)

• td (TraceData)

Return type:

TraceData

croak.session.pipeline.generation_envelope(p)

Build the generation envelope g(z) at the depth-quadrature nodes, or None.

The Gaussian-crossing form: a Rayleigh amplitude droop and beamlet walk-through, with a net Gouy phase,

\[\begin{split}g(z) = (1+u^2)^{-1/2}\\, e^{-((z-z_f)/L_w)^2/(1+u^2)}\\, e^{i\\gamma \arctan u},\\qquad u = (z - z_f)/z_R,\end{split}\]

evaluated at the Gauss–Legendre nodes of the depth quadrature. All lengths are in-medium values, measured from the slab entrance (envelope_* fields of RetrieveParams). Negligible for slabs much thinner than z_R; see docs/explanation/forward_model.md.

Parameters:

p (RetrieveParams)

croak.session.pipeline.covariance_kwargs(p, td=None)

Build the supported model and parameterisation options for covariance.

This includes the core slab, reduced smearing and phase parameterisation. It does not reconstruct focal/collection, depth-weight or marginal-peak centred models; do not interpret its output as matching those retrievals. td supplies the carrier needed to rebuild a reduced smearing kernel.

Parameters:

• p (RetrieveParams)

• td (TraceData | None)

Return type:

dict

croak.session.pipeline.run_uncertainty(up, rp, result, td, *, dispersion=None, rng=None, callback=None)

Run an FWHM-uncertainty estimate from the saved parameters.

Dispatches on up.method: the resampling bootstraps (parametric / delay / frequency), the substrate-thickness systematic, or the analytic covariance estimator. When up.propagate_to_dispersion_point is set and dispersion applies any dispersion, every replicate is propagated to that beamline point via transfer_function() before the FWHM/band is taken.

Parameters:

• up (UncertaintyParams)

• rp (RetrieveParams)

• result (RetrievalResult)

• td (TraceData)

• dispersion (DispersionParams | None)

• rng (Generator | None)

Return type:

UncertaintyResult

croak.session.pipeline.assemble_simulated_tracedata(p)

Build a TraceData directly from a simulated scan.

Loads the as-simulated angular-frequency-density trace onto the simulation’s native FFT grid and packages it as a retrieval-ready TraceData, bypassing both the assemble_simulated_load_data() ÷λ² and the build_tracedata() regrid/filtering. This is the “skip marginal-check and preprocess” path: it lets the retrieval be sanity-checked with no doubt about those steps.

It is fully raw — no third-order (λ/µm)^n correction — but still honours the window/thickness selection, the delay-sign convention (reverse_trace) and the known-spectrum choice (use_spectrum / spectrum_source). lam_min_nm/lam_max_nm are ignored; the whole native grid is kept.

Parameters:

p (SimulatedLoadParams) – The simulated-load options.

Returns:

TraceData – On the simulation’s native grid (grid.omega + omega0_pulse == scan.omega); trace is the raw ω-density (Nomega, Ndelay) normalised to unit peak.

Raises:

ValueError – If the file’s ω axis is not a uniform FFT grid centred on the carrier (the retrieval needs a uniform baseband grid).

Return type:

TraceData

Apply a saved DispersionParams to a retrieved pulse.

The dispersion stage compensates residual dispersion on the already-retrieved pulse: Taylor terms, propagation through built-in materials / a gas, and the beam-path mirror stack. The GUI does this interactively against live registry entries; here we resolve the persisted parameters into the same croak.dispersion.apply_dispersion() call, registering the gas / mirrors under transient keys and cleaning them up afterwards so no module-global state leaks.

croak.session.dispersion.compose_dispersion(p, result)

Apply the persisted dispersion compensation to result’s spectrum.

Parameters:

• p (DispersionParams) – The stage-4 settings (Taylor terms, materials, gas, RI material, mirrors).

• result (RetrievalResult) – The retrieved pulse to compensate.

Returns:

RetrievalResult – A copy of result with the dispersed spectrum. When p applies no dispersion the input is returned unchanged.

Return type:

RetrievalResult

Notes

The stored trace, error and gain remain those of the original retrieval at the measurement plane.

croak.session.dispersion.has_dispersion(p)

Whether p applies any dispersion (else composition is a no-op).

Parameters:

p (DispersionParams)

Return type:

bool

croak.session.dispersion.transfer_function(p, grid, omega0)

Build the spectral transfer function H(ω) = a(ω)·exp(iφ_prop(ω)) of p.

Propagation between two beamline points is the deterministic linear map \(\tilde E'(\omega) = H(\omega)\,\tilde E(\omega)\). This resolves the persisted dispersion settings — exactly the ones compose_dispersion() applies — into that transfer function by applying them to an all-ones spectrum, so the phase (Taylor + material + gas + mirror) and any mirror reflectivity match the dispersion stage bit-for-bit. The result is the propagation argument of the croak.uncertainty estimators and croak.covariance.covariance_uncertainty(), which propagate the retrieval uncertainty to that point. Returns the identity (all-ones) when p applies no dispersion.

Parameters:

• p (DispersionParams) – The dispersion-stage settings (Taylor terms, materials, gas, RI material, mirrors).

• grid (Grid) – The retrieval grid (its centred omega axis defines H).

• omega0 (float) – Carrier angular frequency (rad/s).

Returns:

numpy.ndarray – Complex transfer function on grid.omega (length grid.n).

Return type:

NDArray[complex128]

Headless session engine: replay a SessionOptions.

run_session drives the same per-stage mappings the GUI uses (croak.session.pipeline) without any Qt dependency, so a saved options.toml can be replayed from a script or the croak CLI. It returns a SessionResult holding every intermediate (loaded arrays, cleaned trace, retrieved + dispersion-compensated pulse, processed summary, uncertainty).

Reads happen at the load stage; the engine itself writes nothing — saving is the caller’s job (see croak.cli), keeping the core pure.

class croak.session.engine.SessionResult(options, data=None, tracedata=None, result=None, dispersed=None, processed=None, uncertainty=None)

Bases: object

Every intermediate produced by run_session().

Fields are populated up to the last stage that ran; later ones stay None. result is the (post-filtered) retrieval; dispersed is it after the stage-4 dispersion compensation; processed summarises the final pulse (dispersed when dispersion ran, else the retrieval).

Parameters:

• options (SessionOptions)

• data (dict | None)

• tracedata (TraceData | None)

• result (RetrievalResult | None)

• dispersed (RetrievalResult | None)

• processed (ProcessedResult | None)

• uncertainty (UncertaintyResult | None)

croak.session.engine.run_session(options, *, stages=('load', 'preproc', 'retrieve', 'dispersion', 'process'), rng=None, retrieve_callback=None, uncertainty_callback=None)

Replay a session from its parameters, returning all intermediates.

Parameters:

• options (SessionOptions) – The session to replay (typically SessionOptions.from_toml(path)).

• stages (tuple of str, optional) – Which stages to run; see DEFAULT_STAGES. Always executed in dependency order. Include "uncertainty" to also bootstrap error bars.

• rng (numpy.random.Generator, optional) – Generator threaded through retrieval and the uncertainty bootstrap for reproducibility; a fresh default is used if None.

• retrieve_callback (callable, optional) – Progress hooks forwarded to the retrieval and bootstrap respectively.

• uncertainty_callback (callable, optional) – Progress hooks forwarded to the retrieval and bootstrap respectively.

Returns:

SessionResult

Return type:

SessionResult

croak.session.engine.DEFAULT_STAGES: tuple[str, ...] = ('load', 'preproc', 'retrieve', 'dispersion', 'process')

Stages run by default. Dispersion compensation is included (it is deterministic and cheap). Uncertainty is opt-in: it is a slow, stochastic bootstrap (many retrievals) whose results are saved separately, so a plain replay reproduces the retrieved + compressed pulse and leaves error bars to an explicit request (stages=(..., "uncertainty") or croak replay --uncertainty).

Point an existing session at a new dataset folder.

Many acquisitions share a layout: the measured trace, independent spectrum and their backgrounds live together in a per-dataset folder with stable filenames, while the spectral-response calibration files sit in a fixed shared location. retarget() rewrites only the dataset-folder paths (dataset_path_fields()) to a new base folder, keeping their filenames, and leaves the calibration paths (fixed_path_fields()) untouched — so a working retrieval can be re-pointed at the next dataset and rerun.

croak.session.retarget.retarget(options, new_base_dir)

Re-point a session’s dataset files at new_base_dir.

Parameters:

• options (SessionOptions) – The session whose dataset folder should change.

• new_base_dir (str) – The new per-dataset folder. Each dataset file keeps its basename and is looked up here; calibration files are left as-is.

Returns:

SessionOptions – A copy with the load paths retargeted (all other stages unchanged).

Raises:

ValueError – If the session did not come from the experimental loader. Retargeting is a property of the measured per-dataset folder layout: a simulated session has one scansave file with no sibling spectrum/background convention, and a synthetic one has no file at all, so there is nothing to re-point and silently returning the session unchanged would mislead.

Return type:

SessionOptions

Notes

No filesystem check is done here; a missing retargeted file raises an error when the session is replayed (the loader raises on the absent path).

Command line and script generation

The croak console script (croak.cli) replays sessions and retargets them at new datasets; croak.scripting.generate_script() writes a standalone, editable Python script that reproduces a retrieval.

The croak command-line interface: replay a saved retrieval session.

croak replay options.toml reruns a session headlessly (load → preprocess → retrieve → dispersion → process), saving result.h5 and the options.toml used. Point it at one or more new dataset folders with --base-dir to rerun the same settings on similarly-laid-out acquisitions (calibration files are kept fixed; see croak.session.retarget.retarget()).

croak.cli.main(argv=None)

Entry point for the croak console script.

Parameters:

argv (list[str] | None)

Return type:

int

Generate a standalone, editable Python script from a saved session.

croak script options.toml emits a small script that reproduces the retrieval by calling the public croak API stage by stage. Every parameter is rendered as an explicit dataclass literal so the script is transparent and hand-editable, while the pipeline steps call the same croak.session.pipeline functions the GUI and the replay command use, so a generated script never drifts from the real mapping. Keep the input files and software version with the script.

croak.scripting.generate_script(options, *, base_dir=None, output='result.h5', seed=0, script_name='replay.py')

Render an editable Python script that replays options.

The loader section follows the session’s entry: an experimental session calls assemble_load_data(), a simulated one assemble_simulated_load_data() (or assemble_simulated_tracedata() for its raw “skip filtering & regrid” option), matching croak.session.engine.run_session().

Parameters:

• options (SessionOptions) – The session to turn into a script.

• base_dir (str, optional) – If given, the script’s load paths are retargeted to this dataset folder (calibration paths are kept fixed), as croak replay --base-dir does.

• output (str, optional) – The result.h5 path the generated script writes to.

• seed (int, optional) – RNG seed embedded in the script for a reproducible run.

• script_name (str, optional) – The script’s own filename, used only in its docstring.

Returns:

str – The script source.

Raises:

ValueError – If the session came from the synthetic generator (nothing to reload).

Return type:

str