API reference#

Warning

The public API is not frozen and may change before v0.1. The supported Finder/Result path is M0 versus M1. An opt-in M0/M1/M2 ladder is workable in a production workflow, but its statistical estimator remains experimental; searches above two absorbers are not implemented.

Finder and results#

The high-level entry point: give it a spectrum and get a typed Result back.

The finder runs a Spectrum through the quality policy, evidence calculation, and absorber prior, then returns a typed result that the catalog writer can consume.

It covers:

  • the null versus one-absorber comparison, in exact or FILTER mode;

  • all four statuses — completed, insufficient_data, quality_rejected and failed — as results, not exceptions;

  • the legacy probability fields, computed as the reference defines them rather than derived from one another;

  • full provenance: preset, model, prior, grid, backend, compatibility profile, quality policy, evidence mode and the per-spectrum evaluated-sample count.

With max_absorbers=2 and experimental_multi_absorber=True, it also computes the two-absorber evidence and returns an M0/M1/M2 ladder with aligned priors, posteriors and a selected model. The model is workable in a production workflow, while its current statistical estimator remains experimental. It reproduces the reference evidences bitwise under a controlled seed, but close pairs and low-signal pairs are known weaknesses. An M2 result writes to a catalogue as two flat rows sharing a TARGETID; the ladder itself travels in gp_dla_finder.io.structured. See gp_dla_finder.multi.

It does not yet cover:

  • more than two absorbers. max_absorbers > 2 is refused rather than truncated, so a configuration cannot claim a multiplicity that was never evaluated;

  • FILTER together with two absorbers, refused for the same reason: the hybrid is neither full-grid M2 evidence nor the reference’s multi-absorber FILTER path;

  • validated point estimates. The current redshift and column density come from the best evaluated grid point. They are usable for inspection, but not yet reliable enough to quote as science measurements; uncertainty fields stay NaN;

  • multi-absorber early stops beyond the NaN rung.

Per-spectrum empirical-Bayes mean-flux fitting is implemented – see gp_dla_finder.mean_flux.

Unavailable quantities remain absent from the result as NaN or an empty list. They are not replaced with values that could be read as measurements.

class gp_dla_finder.finder.AbsorberCandidate(grid_z_abs, grid_log_nhi, model=1, z_abs_err=nan, log_nhi_err=nan)[source]#

One absorber candidate: an evaluated grid point, not a detection.

A completed result carries a candidate whenever the evidence path evaluated the absorber model, including spectra whose absorber posterior is far below a useful threshold. It is where that model fits best, which is still a well-defined place in a spectrum with no absorber.

In other words, len(result.absorber_candidates) == 1 means “the search ran”, not “a DLA was found”. Whether a candidate is a detection is a policy question answered by Result.detected() at an explicit threshold, and applied at the catalogue boundary – never silently inside the scientific result.

grid_z_abs and grid_log_nhi are the best evaluated grid point. They provide a usable preliminary location, but they are not a validated MAP estimate. The field names keep grid explicit so their status is clear.

Parameters:
model: int = 1#

1 for the one-absorber model, 2 for a member of the best two-absorber pair. Without it, two candidates in a list say nothing about which fit produced them.

Type:

Which model this candidate belongs to

z_abs_err: float = nan#

Uncertainties, if a validated estimator ever supplies them. NaN for now.

class gp_dla_finder.finder.Finder(config=None, *, model=None, prior=None, grid=None, warn_about_threads=True)[source]#

Run the current inference path for one spectrum at a time.

Construct once and reuse: the model, prior and sample grid are loaded here, not per spectrum, and loading them per call would dominate the runtime.

Parameters:
run(spectrum, *, targetid=0)[source]#

Score one spectrum. Never raises for an ordinary bad input.

A processing failure is a result with a status and a reason code, not an exception, because a batch layer needs to record it and continue — and because a failure silently reported as “no absorber” would corrupt any population statistic built from the output.

Parameters:
Return type:

Result

class gp_dla_finder.finder.Result(targetid, z_qso, status, reason='', ladder=None, mean_flux=None, absorber_candidates=<factory>, log_evidence_null=nan, log_evidence_absorber=nan, p_absorber=nan, p_null=nan, legacy_two_model_p_absorber=nan, logp_absorber=nan, logp_null=nan, evidence_mode='', n_evaluated=0, screening_score=nan, screening_n_evaluated=0, screening_log_evidence_null=nan, screening_log_evidence_absorber=nan, n_usable_pixels=0, quality_fraction=nan, provenance=<factory>)[source]#

The inference record for one spectrum, including an incomplete run.

Immutable, and immutable through its containers. The dataclass is frozen, the candidate list becomes a tuple, and provenance is recursively frozen: nested mappings become read-only proxies, nested sequences become tuples, and arrays are copied and marked non-writable. Mutating whatever the caller passed in cannot reach the result afterwards.

This keeps the provenance tied to the run that produced it, even if the caller later changes the objects originally passed to the package.

Parameters:
ladder: ModelLadder | None = None#

Log evidences for M0, M1 and (when max_absorbers >= 2) M2.

None when only the null-versus-one comparison ran. Check ModelLadder.complete before reading the posteriors: a stopped rung is a model that was never measured, not one measured as impossible.

mean_flux: MeanFluxFit | None = None#

The per-spectrum mean-flux scan, when one ran. None otherwise.

Carried whole rather than summarised: the winner alone cannot tell a reader whether the grid was decisive or whether the maximum sat against an edge, and both change how much the fitted value should be trusted.

absorber_candidates: Sequence[AbsorberCandidate]#

Evaluated absorber candidates – NOT detections. See AbsorberCandidate; use detected() for the policy question.

log_evidence_null: float = nan#

Model evidences. NaN unless the status is completed.

p_absorber: float = nan#

Two-model posterior of at least one absorber, under the named prior. Posterior probability of AT LEAST ONE absorber.

On an M2 run this is the multi-model value, summed over the completed absorber models. On a null-versus-one run it is the two-model value. One number, one meaning – previously the ladder and this field could disagree, and detected() silently used the two-model one.

legacy_two_model_p_absorber: float = nan#

The two-model M0/M1 posterior, kept for legacy comparability.

Equal to p_absorber when only M0 and M1 were evaluated. On an M2 run it is the value the OLD two-model calculation would have given, and it is never used for selection or detection.

logp_absorber: float = nan#

log posterior of the one-absorber model. Carried independently of p_absorber: the reference treats them as different quantities.

evidence_mode: str = ''#

"exact" or "filter". Empty when the status is not completed. This describes the evidence fields above, which a later refinement stage may replace.

screening_score: float = nan#

FILTER-prefix log Bayes factor, absorber over null. A ranking statistic: not a probability, not an evidence. NaN means no screening stage ran, never “screened and scored zero”.

screening_n_evaluated: int = 0#

Samples the screening stage evaluated. 0 when it did not run.

screening_log_evidence_null: float = nan#

The two evidences the screening stage produced, kept so the screening decision can be re-derived after a refinement stage replaces the finals.

property log_bayes_factor: float#

log evidence ratio, absorber over null, in this result’s mode.

property was_screened: bool#

Whether a FILTER screening stage ran for this spectrum.

refined(*, log_evidence_null, log_evidence_absorber, p_absorber, p_null, logp_absorber, logp_null, n_evaluated, provenance=None)[source]#

A copy with full-grid evidences replacing the screened ones.

The forward-compatibility path for the two-stage workflow: screen with FILTER, then re-run candidates near a decision boundary on the full grid. The screening fields are carried through unchanged, so the product records both what screening said and what refinement concluded.

Nothing calls this yet – the refinement stage is not implemented. It exists so the representation is testable now rather than discovered to be lossy later.

Parameters:
Return type:

Result

detected(threshold)[source]#

Whether this counts as a detection at threshold.

threshold is required and has no default. A detection depends on the analysis threshold, so callers must state that policy explicitly.

Parameters:

threshold (float)

Return type:

bool

gp_dla_finder.finder.results_to_catalogue(results, *, detection_threshold)[source]#

Turn results into the catalogue model the FITS writers consume.

detection_threshold is required and has no default. It decides which completed results contribute an absorber row, so it determines what the catalog contains. This is a policy decision that belongs to the caller, not to a library default. It is recorded in the run record as GPDLF_DETECTION_THRESHOLD so a consumer can recover the selection that produced the file.

All results must agree on RUN_DEFINING_PROVENANCE; see _check_run_provenance().

Parameters:
Return type:

Catalogue

gp_dla_finder.finder.screening_score(result)[source]#

The stored FILTER screening statistic, or NaN if no screening ran.

Reads Result.screening_score. It used to derive the value from the result’s current evidence fields, which broke in two ways the two-stage workflow would have hit immediately: a refinement stage rewrites those fields and the mode, so the derivation returned NaN for a spectrum that certainly had been screened; and a screened non-detection has no absorber row to carry it.

Definition. The log Bayes factor – one-absorber over null – computed from the FILTER prefix estimate. It is a ranking statistic: larger means the one-absorber model fits better on the samples that were evaluated. It is deliberately not a probability and not an evidence, which is why the column is named GPDLF_SCREENING_SCORE rather than anything resembling P_ or LOGZ_.

Why it is not redundant with the row’s log Bayes factor. The two are currently equal because FILTER evidence columns contain the prefix values. A future two-stage workflow can refine candidates near a decision boundary on the full grid while retaining this value from the screening stage. The catalog can then record both stages.

NaN in full-grid mode, because there was no screening stage. A NaN here means “not screened”, never “screened and scored zero”.

Populated only for completed results; a failed one has no evidences to combine.

Parameters:

result (Result)

Return type:

float

Configuration#

Inference configuration and the named operating points we support.

A Config contains the knobs that actually affect inference. It is smaller than the reference pipeline’s Parameters, which also mixes in training, file-loading, and catalog-filtering settings. Keeping those out makes it easier to see what changed the numerical result.

Model-coupled quantities are not part of the configuration. The GP rank, rest-frame grid, and flux normalization band belong to the trained model. Config.validate_against() checks that the model and configuration agree.

Presets#

In practice, start from a preset instead of assembling a config by hand. The preset records the operating point you started from.

class gp_dla_finder.config.Config(min_lambda=911.75, max_lambda=1250.0, num_forest_lines=31, num_lines=3, prev_tau_0=0.00246, prev_beta=3.62, max_z_cut_kms=3000.0, min_z_cut_kms=3000.0, min_z_separation_kms=3000.0, prior_z_qso_increase_kms=30000.0, max_absorbers=4, num_samples=50000, sample_grid='pw14_172_225_50000', log_nhi_range=(17.2, 22.5), log_nhi_prior_alpha=0.97, filter_low_likelihood=False, filter_n_initial_floor=5000, filter_empty_mask_fallthrough=False, early_stop_mode='baseline', voigt_backend='numpy', lsf_kernel='desi-r3000-7tap', broadening=True, pixel_spacing_dex=0.0001, enable_tau_eb=True, tau_eb_factors=(0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0), tau_eb_objective='null', tau_eb_apply_hcd_mask=False, tau_eb_mask_threshold_sigma=1.5, seed=0, quality_policy=None, compatibility='reference-d5b306e6', preset='\x00unset', experimental_multi_absorber=False, base_preset='')[source]#

Everything the inference path reads, apart from the trained model.

Distances in velocity are given in km/s and converted internally, matching the reference implementation’s kms_to_z.

Parameters:

  • min_lambda (float)

  • max_lambda (float)

  • num_forest_lines (int)

  • num_lines (int)

  • prev_tau_0 (float)

  • prev_beta (float)

  • max_z_cut_kms (float)

  • min_z_cut_kms (float)

  • min_z_separation_kms (float)

  • prior_z_qso_increase_kms (float)

  • max_absorbers (int)

  • num_samples (int)

  • sample_grid (str)

  • log_nhi_range (tuple[float, float])

  • log_nhi_prior_alpha (float)

  • filter_low_likelihood (bool)

  • filter_n_initial_floor (int)

  • filter_empty_mask_fallthrough (bool)

  • early_stop_mode (Literal['baseline', 'no_null_stop', 'pre_occam'])

  • voigt_backend (str)

  • lsf_kernel (str)

  • broadening (bool)

  • pixel_spacing_dex (float)

  • enable_tau_eb (bool)

  • tau_eb_factors (tuple[float, ...])

  • tau_eb_objective (Literal['null', 'absorber'])

  • tau_eb_apply_hcd_mask (bool)

  • tau_eb_mask_threshold_sigma (float)

  • seed (int | None)

  • quality_policy (str | None)

  • compatibility (str)

  • preset (str)

  • experimental_multi_absorber (bool)

  • base_preset (str)

min_lambda: float = 911.75#

Rest-frame window the GP is evaluated on, angstroms.

num_forest_lines: int = 31#

Lyman-series members used for the mean-flux suppression / forest noise.

num_lines: int = 3#

Lyman-series members used for the absorber’s own Voigt profile.

prev_tau_0: float = 0.00246#

Effective-optical-depth normalisation, tau_eff = tau_0 (1+z)^beta. Default is Turner et al. (2024), arXiv:2405.06743.

max_z_cut_kms: float = 3000.0#

z_abs search window is inset from the quasar redshift and the Lyman limit by these velocities.

min_z_separation_kms: float = 3000.0#

Minimum velocity separation between two absorbers in a multi-absorber model.

prior_z_qso_increase_kms: float = 30000.0#

The absorber-existence prior counts catalogued quasars with z < z_qso + this, expressed as a velocity.

max_absorbers: int = 4#

Maximum number of absorbers modelled per spectrum.

num_samples: int = 50000#

Number of quasi-Monte-Carlo samples in the evidence integral.

sample_grid: str = 'pw14_172_225_50000'#

Name of the QMC sample grid this operating point uses. Kept explicit so a preset never inherits a grid by accident.

log_nhi_range: tuple[float, float] = (17.2, 22.5)#

Column-density prior range, log10(N_HI / cm^-2).

The deployed range reaches well below the DLA threshold. The low end regularises the inference and is required to reproduce the production catalogs, but performance there is not independently validated. Treat log10 N_HI > 20 as the trusted regime.

log_nhi_prior_alpha: float = 0.97#

Mixture weight of the Prochaska et al. (2014) component in the N_HI prior.

filter_low_likelihood: bool = False#

Evaluate the one-absorber evidence on a truncated prefix of the sample grid instead of the whole configured grid.

Off by default for v0.1. The full-grid path is the conservative one and is what a production preset selects; FILTER is a screening approximation that a caller must ask for by name. It changed classification in 3 of 15 constructed cases against the adopted 100k full-grid reference, so it is not interchangeable with the full-grid evidence. See docs/filter.md.

filter_n_initial_floor: int = 5000#

Floor on the coarse-scan budget; the scan uses max(num_samples // 20, this).

filter_empty_mask_fallthrough: bool = False#

False stops early with the coarse 1-absorber evidence (deployed), True falls through to the full sample set.

Type:

When the coarse scan finds no viable region

early_stop_mode: Literal['baseline', 'no_null_stop', 'pre_occam'] = 'baseline'#

Multi-absorber early-stop policy. "baseline" is deployed.

voigt_backend: str = 'numpy'#

Voigt backend name. "numpy" is the official backend and is always available; "libcerf" exists only where the optional compiled extension was built, and reproduces the Faddeeva implementation behind the deployed catalogues. Naming a backend that was not built raises at construction time rather than silently falling back to a different forward model.

lsf_kernel: str = 'desi-r3000-7tap'#

Named line-spread function; see gp_dla_finder.voigt.

broadening: bool = True#

Apply instrumental broadening to absorber profiles.

pixel_spacing_dex: float = 0.0001#

Pixel spacing, dex, used to pad the spectrum for the LSF convolution.

enable_tau_eb: bool = True#

Fit tau_0 per spectrum before inference instead of using prev_tau_0.

tau_eb_factors: tuple[float, ...] = (0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0)#

Candidate multipliers of prev_tau_0 scanned by the fit.

tau_eb_objective: Literal['null', 'absorber'] = 'null'#

the null-model evidence (cheap) or the absorber model.

Type:

Objective for the fit

tau_eb_apply_hcd_mask: bool = False#

at population scale it over-corrects.

Type:

Mask strongly negative residuals during the fit. Off in production

seed: int | None = 0#

Seed for the multi-absorber resampler. An integer makes a run reproducible; None opts into nondeterminism explicitly.

quality_policy: str | None = None#

Named catalog-selection policy, or None for no policy. The production presets name the deployed DESI requirement; a custom configuration must choose one deliberately, because a quality cut decides which spectra reach a catalogue and must never be inferred from an instrument label. None leaves selection entirely to the caller and is not a validated catalogue rule – see gp_dla_finder.quality.

compatibility: str = 'reference-d5b306e6'#

Named compatibility profile; see gp_dla_finder.compat. The default reproduces the pinned reference bitwise, including two floating-point no-ops it performs. Selecting "clean" drops them and changes results, which is why the profile name travels in result provenance.

preset: str = '\x00unset'#

Name recording what this configuration reproduces.

There is no default. Constructing a Config without naming a preset raises, so no workflow can silently inherit the DESI production operating point. Use desi_y3(), or pass preset="custom" to state deliberately that this configuration reproduces nothing standard.

experimental_multi_absorber: bool = False#

Opt in to the EXPERIMENTAL two-absorber path.

max_absorbers=2 alone is not enough: this must also be True. Two deliberate choices, because the M2 estimator matches the legacy implementation on the tested surface but is not a validated close-pair method – it is weak on close pairs and low signal, and its bounded benchmark is not a survey calibration. A user should not reach it by nudging a number.

base_preset: str = ''#

The named preset this configuration STARTED from, kept even after scientifically consequential overrides.

preset is the effective identity and becomes "<base>+modified" as soon as anything consequential is overridden; base_preset answers the other question a reader needs – which operating point was the starting point. Config.desi_y3(filter_low_likelihood=True) is not the canonical desi_y3 and must not claim that name, but it did not come from nowhere either.

property compatibility_profile: CompatibilityProfile#

The arithmetic-fidelity profile this configuration runs under.

property selected_quality_policy: QualityPolicy | None#

The named quality policy, or None when selection is the caller’s.

property backend_provenance: Mapping[str, object]#

What the selected Voigt backend is, for a result’s provenance block.

For a compiled backend this includes the libcerf version, the sha256 of the shared library it linked, the compiler and the redacted optimisation flags – because a compiled backend’s numbers depend on all of them.

property evidence_mode: str#

How the one-absorber evidence integral is evaluated.

"exact" evaluates the whole configured sample grid and is the v0.1 default for every production preset.

The name "exact" is an API compatibility label, not a mathematical claim. It means “the package’s full configured QMC-grid estimator” – not an analytic integral, not zero numerical error, and not a formal convergence proof. The grid itself is finite and its own residual error is unquantified (see docs/filter.md).

"filter" evaluates only a prefix of that grid and is a screening approximation. It must be requested explicitly; no production preset selects it. The chosen mode is a scientific choice, not an optimisation detail, so it travels in result provenance and labels every value it produces.

static kms_to_z(kms)[source]#

Velocity in km/s as a redshift difference.

Parameters:

kms (float)

Return type:

float

property convolution_half_width: int#

Pixels of padding each side needed for the LSF convolution.

property n_models: int#

null plus 1..max_absorbers.

Type:

Number of models compared

property model_labels: tuple[str, ...]#

Labels for the model-posterior vector, in order.

INERT_FIELDS: ClassVar[frozenset[str]] = frozenset({'base_preset', 'preset'})#

Fields whose value cannot change what the inference computes, so overriding them does not make a configuration “modified”.

replace(**changes)[source]#

A modified copy, relabelled if anything consequential changed.

The effective preset becomes "<base>+modified"; base_preset is preserved, so provenance answers both “which named operating point was the starting point” and “what exactly ran”.

Overriding only inert bookkeeping does not relabel: a name is a claim about the numerics, and re-stating a value that cannot move a number is not a change to them.

Return type:

Config

property is_modified: bool#

Whether this configuration diverges from its named base preset.

normalized()[source]#

Every scientifically consequential field, in a stable order.

The complete answer to “what exactly produced this result”, as opposed to the preset name, which is only a label.

Return type:

dict[str, object]

property digest: str#

Stable short hash of normalized().

Two configurations with the same digest compute the same thing, whatever they are called. Comparing digests is how results are checked for compatibility before being combined into one catalogue.

validate_against(model)[source]#

Check a trained model can serve this configuration.

Raises:

ValueError – If the model’s rest-frame grid does not span the search window, or if the model does not record a normalisation band. Both would otherwise fail deep inside the likelihood, or worse, silently produce a wrong prediction.

Parameters:

model (GPModel)

Return type:

None

classmethod desi_y3(**overrides)[source]#

The deployed DESI Y3 production operating point.

This is the configuration behind the production catalogs: a two-way model comparison (null versus k absorbers, no separate sub-DLA channel), a column-density prior spanning [17.2, 22.5], and the full configured evidence grid.

Return type:

Config

classmethod desi_y3_refined(**overrides)[source]#

As desi_y3() but with the 100,000-sample QMC grid.

A denser numerical integration grid, selected explicitly. At least one archived mock-catalog production used 100k rather than 50k, so reproducing that run means choosing this preset deliberately and recording it.

“Refined” describes the integration grid only. It does not establish identity with any deployed production array, and does not by itself make a result more accurate scientifically.

Return type:

Config

classmethod desi_y3_fast(**overrides)[source]#

As desi_y3() but with a smaller sample budget.

For exploration and tutorials. The evidence integral is noisier, so multi-absorber results in particular will differ from production.

Return type:

Config

Forward model#

Voigt absorption profiles for the Lyman series.

The absorption profile of an absorber at redshift z_dla with column density nhi is

\[\tau(\lambda) = N_\mathrm{HI} \sum_j a_j\, V(v_j(\lambda); \sigma, \gamma_j)\]

evaluated as exp(-tau) and convolved with a discrete instrumental line-spread function (LSF). V is the Voigt function, computed from the Faddeeva function scipy.special.wofz — the same function the reference C implementation obtains from libcerf.

Backends#

The numerics live behind VoigtBackend. Two backends exist:

numpy

NumpyVoigtBackend, using scipy.special.wofz. Always present, and the official v0.1 backend.

libcerf

LibcerfVoigtBackend, a compiled extension using libcerf’s Faddeeva function. Present only where the optional extension was built. It exists for fidelity, not speed: it is the Faddeeva implementation behind the deployed DESI catalogs, and it is measurably not the same function as SciPy’s. End to end it is about 3% slower than the NumPy backend (measured; the Voigt evaluation is roughly a third of the per-sample cost and NumPy’s wofz ufunc is already a tight C loop).

Backend selection is explicit. Requesting a backend that was not built raises an error. A compiled backend that disagrees with the NumPy backend beyond BACKEND_AGREEMENT_ATOL is not registered at all.

Instrumental line-spread function#

The LSF is a named, tabulated kernel, not a global default, because the kernel is a property of the instrument/model configuration and choosing the wrong one changes the inferred profile shape. The production kernel is "desi-r3000-7tap".

Warning

The historical pure-Python module in the reference implementation (gpy_dla_detection/voigt.py) carries a BOSS R=2000 kernel while the compiled extension it stands in for uses the DESI R≈3000 kernel above. Those are different forward models: the profile shapes differ by up to ~4e-2 at log10 N_HI = 19 on a DESI grid, largest in the LLS/sub-DLA regime. Users must select a kernel that describes their instrument.

References

Garnett et al. (2017), arXiv:1605.04460 – note this is the correct identifier; some upstream docstrings cite 1605.04538, which is an unrelated paper. Ho, Bird & Garnett (2020), arXiv:2003.11036. Ho, Bird & Garnett (2021), arXiv:2103.10964.

gp_dla_finder.voigt.available_backends()[source]#

Names of Voigt backends available in this installation.

Always contains "numpy". Contains "libcerf" only where the optional compiled extension was built and passed its agreement check at import.

Return type:

tuple[str, …]

gp_dla_finder.voigt.backend_provenance(name)[source]#

What a backend is, for the record a result carries.

Includes the Faddeeva implementation, whether the backend is compiled, and – for a compiled backend – its measured agreement with the NumPy backend.

Parameters:

name (str)

Return type:

Mapping[str, object]

gp_dla_finder.voigt.backend_rejections()[source]#

Backends that were built but failed their agreement checks, and why.

Return type:

Mapping[str, str]

gp_dla_finder.voigt.gaussian_lsf_kernel(resolving_power, pixel_scale, wavelength, *, truncate_sigma=4.0)[source]#

Build a Gaussian LSF kernel from a resolving power.

Conventions, all of which are testable and none of which are implicit:

  • FWHM = wavelength / resolving_power in angstroms, so sigma = FWHM / (2 sqrt(2 ln 2));

  • the kernel is sampled on the pixel grid, sigma_pixels = sigma / pixel_scale, at integer offsets from the centre;

  • it is truncated at +/- ceil(truncate_sigma * sigma_pixels) pixels, so it always has odd length and a well-defined centre;

  • it is normalised to sum to exactly 1 after truncation.

Warning

This is an approximation: a Gaussian of constant resolving power. Real instruments have a line-spread function that varies with wavelength and across spectrograph arms, and is not exactly Gaussian. Use it for exploration and for instruments this package has no named kernel for. It is not a reproduction path – for that use a named historical kernel (PRODUCTION_KERNEL, BOSS_KERNEL), whose values are pinned and parity-tested.

Parameters:

  • resolving_power (float) – R = lambda / FWHM, dimensionless. Must be finite and positive.

  • pixel_scale (float) – Wavelength step per pixel, angstroms. Must be finite and positive.

  • wavelength (float) – Wavelength at which to evaluate R, angstroms. For a constant-R instrument use a representative wavelength of the search window.

  • truncate_sigma (float) – Half-width of the kernel in standard deviations.

Returns:

Normalised, symmetric, odd-length, read-only kernel.

Return type:

numpy.ndarray

Raises:

ValueError – On non-finite or non-positive inputs, or if the resulting kernel would be narrower than one pixel – an under-sampled LSF is silently wrong, so it is refused rather than rounded up.

gp_dla_finder.voigt.voigt_absorption(wavelengths, nhi, z_dla, num_lines=3, *, kernel='desi-r3000-7tap', broadening=True, backend='numpy')[source]#

Instrument-convolved Lyman-series absorption profile exp(-tau).

Parameters:

  • wavelengths (ndarray) – Observed-frame wavelengths, angstroms. When broadening is true these must be padded by kernel_half_width(kernel) pixels at each end.

  • nhi (float) – Neutral hydrogen column density, cm^-2 (linear, not log10).

  • z_dla (float) – Absorber redshift.

  • num_lines (int) – Number of Lyman-series members to include, starting at Lyman-alpha.

  • kernel (str) – Named LSF kernel. Defaults to the production DESI kernel.

  • broadening (bool) – Convolve with the LSF. When false the bare profile is returned at full length.

  • backend (str) – Voigt backend name.

Returns:

exp(-tau) in [0, 1]; length len(wavelengths) - 2 * half_width when broadening, else len(wavelengths).

Return type:

numpy.ndarray

gp_dla_finder.voigt.lsf_kernel(name)[source]#

Return a named LSF kernel.

The returned array is read-only and shared; copy it if you need to modify one (lsf_kernel(name).copy()).

Raises:

KeyError – If name is not a known kernel. There is deliberately no default and no nearest-match behaviour: silently substituting a kernel changes the forward model.

Parameters:

name (str)

Return type:

ndarray

Spectra and evidence#

Validate and prepare one spectrum for the GP likelihood.

Turns a raw quasar spectrum into the arrays the Gaussian-process likelihood needs:

  1. structural validation — is this a spectrum at all;

  2. normalisation and masking — divide by the median flux over the model’s own training band, then keep only unmasked pixels inside the search window;

  3. padded grid — extend the window by the line-spread function’s half-width so absorber profiles can be convolved without edge effects.

The arithmetic follows the reference implementation’s NullGP.set_data. Two details affect numerical fidelity and are preserved exactly:

  • the padding is built with np.logspace at a fixed spacing in dex. On a linear-Å grid that is not the local pixel scale, so the padded pixels are not evenly spaced with their neighbours. It affects only the outermost half_width pixels of the convolution, but it is part of the reference arithmetic;

  • the median is nanmedian over pixels that are inside the normalisation band and unmasked, computed before the search-window cut.

exception gp_dla_finder.gp.spectrum.InsufficientData(reason, detail='')[source]#

A valid spectrum that cannot support inference.

Not an error in the input and not a non-detection. Carries a stable reason code so a batch layer can aggregate causes.

Parameters:
Return type:

None

class gp_dla_finder.gp.spectrum.PreparedSpectrum(rest_wavelength, wavelength, flux, noise_variance, window_wavelength, mask_in_window, padded_wavelength, normalization_median, z_qso, diagnostics=<factory>)[source]#

A spectrum reduced to the pixels the likelihood will actually use.

Parameters:
rest_wavelength: ndarray#

Rest-frame wavelengths of the kept pixels.

wavelength: ndarray#

Observed-frame wavelengths of the kept pixels.

flux: ndarray#

Normalised flux at the kept pixels.

noise_variance: ndarray#

Normalised noise variance at the kept pixels.

window_wavelength: ndarray#

Observed wavelengths inside the search window before masking. The absorber profile is evaluated here, then reduced by mask_in_window.

mask_in_window: ndarray#

Which pixels of window_wavelength survive the mask.

padded_wavelength: ndarray#

window_wavelength extended by the LSF half-width at each end.

normalization_median: float#

The median used to normalise, in the input’s flux units.

class gp_dla_finder.gp.spectrum.Spectrum(wavelength, flux, ivar, z_qso, mask=None)[source]#

One quasar spectrum, in observed-frame arrays.

Parameters:

  • wavelength (numpy.ndarray) – Observed-frame wavelengths, angstroms, strictly increasing.

  • flux (numpy.ndarray) – Calibrated flux. Units are irrelevant: the model self-normalises.

  • ivar (numpy.ndarray) – Inverse variance. Zero marks a pixel with no information; such pixels are folded into the mask, matching the reference pipeline.

  • mask (numpy.ndarray | None) – Optional additional bad-pixel mask, True meaning bad.

  • z_qso (float) – Quasar emission redshift.

property rest_wavelength: ndarray#

Rest-frame wavelengths, lambda / (1 + z_qso).

property noise_variance: ndarray#

1 / ivar, with NaN where ivar == 0, as the reference expects.

gp_dla_finder.gp.spectrum.prepare_spectrum(spectrum, model, config, *, min_usable_pixels=1)[source]#

Validate, normalise, mask, window and pad a spectrum.

Parameters:

  • min_usable_pixels (int) – Below this many unmasked pixels in the search window the spectrum cannot support inference. The default of 1 is a structural floor only: it is deliberately not a science quality cut, because the reference pipeline’s selection (>20 % unmasked in a fixed rest-frame window) lives in its DESI I/O layer, not in the inference, and inventing a threshold here would silently change which spectra are searched.

  • spectrum (Spectrum)

  • model (GPModel)

  • config (Config)

Raises:

InsufficientData – The spectrum is valid but unusable: no normalisation coverage, nothing left after masking, or too few usable pixels.

Return type:

PreparedSpectrum

Assemble the GP model and calculate null and one-absorber evidences.

The workflow first interpolates the trained model onto the spectrum, applies the Lyman-series mean-flux suppression and absorption noise, and then evaluates the null evidence. The one-absorber evidence is a quasi-Monte-Carlo average over the absorber sample grid.

The arithmetic follows NullGP.get_interp, NullGP.log_model_evidence and the reference DLAGP sample-likelihood path. A few details that look incidental affect reference fidelity:

  • mu, log_omega and each of the k columns of M are interpolated separately with 1-D linear interpolation, matching the reference’s per-eigenvector loop;

  • the absorption-noise term is multiplied by the squared mean-flux suppression, because the whole model is re-levelled to mu * a_lya rather than mu;

  • every per-sample log-likelihood carries - log(N), and the evidence adds + log(N) back. The pair cancels analytically, but the order changes the floating-point result and is therefore preserved.

gp_dla_finder.gp.evidence.coarse_scan_size(config)[source]#

How many samples the FILTER coarse scan evaluates.

max(num_samples // 20, filter_n_initial_floor), verbatim from the reference’s parallel_log_model_evidences. The floor of 5000 reproduces the historical hard-coded value, so at the 10k operating point the “20x reduction” is actually a factor of two, and only at 100k does it become 20x.

Parameters:

config (Config)

Return type:

int

gp_dla_finder.gp.evidence.assemble_model(prepared, model, config, *, tau_0=None, beta=None)[source]#

Interpolate the model and apply the forest terms.

Parameters:

  • tau_0 (float | None) – Mean-flux prior. Defaults to the configuration’s values; the per-spectrum empirical-Bayes fit will pass its own.

  • beta (float | None) – Mean-flux prior. Defaults to the configuration’s values; the per-spectrum empirical-Bayes fit will pass its own.

  • prepared (PreparedSpectrum)

  • model (GPModel)

  • config (Config)

Return type:

AssembledModel

gp_dla_finder.gp.evidence.null_log_evidence(prepared, assembled)[source]#

Stage 6: the null-model log evidence.

Parameters:
Return type:

float

gp_dla_finder.gp.evidence.one_absorber_log_evidence(prepared, assembled, grid, config, *, return_samples=False, mode=None)[source]#

Stage 7: the one-absorber log evidence, by QMC average over the grid.

\[p(D \mid 1) \approx \frac{1}{N} \sum_i p(D \mid z_i, N_{\mathrm{HI},i})\]

computed in log space with the reference’s exact bookkeeping: each sample likelihood carries - log N, and the estimator adds + log N back after the log-mean-exp.

Parameters:

  • mode (str | None) – None takes the mode from Config.evidence_mode. Production presets use the full configured grid. Passing a mode explicitly overrides the configuration, as in the reference-parity tests. “exact” evaluates all num_samples grid points. “filter” evaluates only the first coarse_scan_size() of them. For the one-absorber evidence that is the whole of what the reference’s FILTER=1 path does: its adaptive region-A machinery selects a high-likelihood subset for refinement, but the reference’s “FILTER fix #5” then discards those refined samples when forming the 1-absorber evidence and uses the coarse scan alone. The refinement only reaches the multi-absorber evidences, so at k = 1 FILTER is not an adaptive approximation at all: it is the same estimator on a fixed prefix of the QMC sequence. Note the normalisation: the log(N) in the round trip stays the full sample count even when only a prefix is evaluated, exactly as the reference does it.

  • prepared (PreparedSpectrum)

  • assembled (AssembledModel)

  • grid (AbsorberSampleGrid)

  • config (Config)

  • return_samples (bool)

Returns:

The log evidence, or a pair of the log evidence and the per-sample log likelihoods when return_samples is set. In filter mode the un-evaluated tail of that array is NaN, so a caller can see what was and was not computed. The per-sample array is not retained by default.

Return type:

float or tuple

gp_dla_finder.gp.evidence.absorber_search_window(prepared, config)[source]#

The z_abs interval searched for this spectrum.

Verbatim from the reference’s Parameters.min_z_dla / max_z_dla: the window is the part of the modelled region where a Lyman-alpha absorber could fall, inset from the quasar and from the Lyman limit by the configured velocities.

Parameters:
Return type:

tuple[float, float]

Gaussian-process likelihood for a quasar spectrum.

The null (no-absorber) model for the observed flux is

\[p(y \mid \lambda, \sigma^2) = \mathcal{N}(y;\; \mu \cdot a_{\mathrm{Ly}\alpha},\; M M^{\mathsf{T}} + \Omega^2 + V)\]

where mu and M come from the trained model, a_lya is the Lyman-series mean-flux suppression, Omega^2 is the Lyman-forest absorption noise, and V is the pipeline’s own per-pixel noise variance.

The covariance is a rank-k update to a diagonal, so the Woodbury identity evaluates the log-likelihood in O(n k^2) rather than O(n^3) — the single hot operation in the whole finder.

The numerics follow the reference implementation (gpy_dla_detection/effective_optical_depth.py and gpy_dla_detection/null_gp.py). Expression order, summation order and the LAPACK routine all affect bitwise reference fidelity, even when an alternative would be mathematically equivalent.

gp_dla_finder.gp.likelihood.effective_optical_depth(wavelengths, beta, tau_0, z_qso, num_forest_lines)[source]#

Per-pixel, per-line Lyman-series effective optical depth.

\[\tau_i(\lambda) = \tau_0 \frac{f_i \lambda_i}{f_1 \lambda_1} (1 + z_i)^{\beta}, \qquad 1 + z_i = \lambda / \lambda_i\]

Absorbers beyond the quasar are switched off by an indicator z_i <= z_qso. The reference multiplies by the indicator rather than masking, so pixels outside a line’s forest contribute exactly zero rather than NaN; that choice is preserved because exp(-sum) then stays finite everywhere.

Parameters:

  • wavelengths (ndarray) – Observed-frame wavelengths, angstroms.

  • beta (float) – Effective-optical-depth power law, tau_eff = tau_0 (1 + z)^beta.

  • tau_0 (float) – Effective-optical-depth power law, tau_eff = tau_0 (1 + z)^beta.

  • z_qso (float) – Quasar redshift.

  • num_forest_lines (int) – Number of Lyman-series members to include.

Returns:

Shape (n_pixels, num_forest_lines). Sum over axis 1 for the total.

Return type:

numpy.ndarray

gp_dla_finder.gp.likelihood.log_mvnpdf_low_rank(y, mu, M, d)[source]#

log N(y; mu, M M^T + diag(d)) via the Woodbury identity.

With D = diag(d) and B = I + M^T D^-1 M,

\[K^{-1} = D^{-1} - D^{-1} M B^{-1} M^{\mathsf{T}} D^{-1}, \qquad \log\det K = \log\det D + \log\det B\]

so the cost is O(n k^2) instead of O(n^3).

Ported verbatim, including two choices that look incidental and are not: B gets its identity added by strided in-place addition on the raveled array, and B^-1 is applied via two lapack.dtrtri triangular inversions rather than triangular solves. Both affect the floating-point result at the last bits, which matters because the reference-equivalence tests are bitwise.

Parameters:

  • y (ndarray) – Observed values, shape (n,).

  • mu (ndarray) – Mean, shape (n,).

  • M (ndarray) – Low-rank covariance factor, shape (n, k).

  • d (ndarray) – Diagonal of the noise covariance, shape (n,). Must be positive.

Returns:

The log density.

Return type:

float

Priors, grids and models#

Absorber-existence prior: P(k absorbers | z_qso).

The prior is empirical. It comes from counting, among catalogued sightlines below a quasar’s redshift, what fraction host a damped absorber:

\[P(\ge 1 \text{ absorber} \mid z_\mathrm{QSO}) = M / N\]

with M the number of catalogued absorbers and N the number of quasars with z < z_QSO + delta. Multi-absorber priors follow by assuming independence, P(>= k) = (M/N)^k, and differencing.

Representation#

The reference pipeline recomputes this from ~115 MB of SDSS catalogues on every run, but touches them through a single counting call. Since that call is a monotone step function of one scalar, it is stored here as a sorted redshift array plus a cumulative absorber count – an exact representation, about 0.1 MB, with no interpolation and no tolerance. tools/build_prior_table.py builds it and proves the equivalence at every breakpoint and on a dense grid.

class gp_dla_finder.prior.AbsorberPrior(name, z_qsos, cumulative_absorbers, *, provenance=<factory>)[source]#

An empirical absorber-existence prior as an exact step table.

Parameters:
name#

Asset name, or "<external>".

Type:

str

z_qsos#

Sorted quasar redshifts of the selected catalogue sample, shape (N,).

Type:

numpy.ndarray

cumulative_absorbers#

cumulative_absorbers[i] is the number of sightlines among z_qsos[:i+1] that host an absorber.

Type:

numpy.ndarray

provenance#

Read-only record of sources, checksums, selection, and the equivalence proof.

Type:

Any

supports(z_qso)[source]#

Whether z_qso lies within the catalogue’s redshift support.

Parameters:

z_qso (float)

Return type:

bool

counts(z_qso, z_increase)[source]#

Number of catalogued absorbers and quasars below z_qso + z_increase.

This is the low-level reference-parity operation. Below the catalogue floor it reproduces the reference’s historical clamp, which the reference source describes as temporary, and emits a warning. Above the top of the catalog, the counts saturate, which is equally an extrapolation.

The public finder rejects unsupported redshifts by default; use supports() to test before calling.

Parameters:
Return type:

tuple[int, int]

absorber_fraction(z_qso, z_increase)[source]#

M / N: P(at least one absorber | z_qso).

Parameters:
Return type:

float

log_prior_no_absorber(z_qso, z_increase)[source]#

log P(no absorber | z_qso).

Parameters:
Return type:

float

log_priors(z_qso, max_absorbers, z_increase)[source]#

log P(exactly k absorbers | z_qso) for k = 1 .. max_absorbers.

P(>= k) = (M/N)^k under independence, differenced to give exactly-k. The last entry is left as P(>= max_absorbers), matching the reference: the top model absorbs the tail.

Parameters:
Return type:

ndarray

gp_dla_finder.prior.available_priors()[source]#

Names of prior tables bundled with this installation.

Return type:

tuple[str, …]

gp_dla_finder.prior.load_prior(name='dr9q_concordance', *, path=None)[source]#

Load an absorber-existence prior table.

Parameters:
Return type:

AbsorberPrior

Quasi-Monte-Carlo absorber sample grids.

The evidence for a k-absorber model is an integral over absorber parameters, approximated by a QMC average:

\[p(D \mid k) \approx \frac{1}{N} \sum_i p(D \mid z_i, N_{\mathrm{HI},i})\]

A AbsorberSampleGrid holds those samples. Column densities are drawn from the Prochaska et al. (2014) CDDF prior mixed with a uniform component; absorber redshifts are stored as offsets in [0, 1) and stretched onto each spectrum’s own search window at inference time, so one grid serves every spectrum.

Provenance status#

Grids carry an explicit identity status. A grid marked "regenerated, production-array identity unverified" reproduces the reference generator bitwise. It has not been shown byte-identical to the deployed arrays, whose QMC state was not recorded reproducibly. You can use such a grid for development and tutorials, but it does not support a production-equivalence claim. AbsorberSampleGrid.is_verified exposes that distinction to callers.

class gp_dla_finder.samples.AbsorberSampleGrid(name, offset_samples, log_nhi_samples, *, nhi_samples=None, integrity_problems=None, provenance=<factory>)[source]#

QMC samples of (redshift offset, log10 N_HI).

Parameters:
name#

Asset name, or "<external>".

Type:

str

offset_samples#

Absorber-redshift offsets in [0, 1), shape (N,). Mapped onto a spectrum’s search window by sample_redshifts().

Type:

numpy.ndarray

log_nhi_samples#

log10(N_HI / cm^-2) samples, shape (N,).

Type:

numpy.ndarray

provenance#

Read-only record of the prior, the QMC construction, the generating environment, per-array checksums, and the identity status.

Type:

Any

integrity_problems: tuple[str, ...] | None = None#

Integrity problems found when this grid was loaded, or None if the check never ran.

None and () are deliberately different. () means the sidecar was checked against the file and the arrays and matched; None means nothing was verified, which is what a hand-constructed grid gets. Only the first is inference-ready – see usable_for_inference.

property log_nhi_range: tuple[float, float]#

Support actually spanned by the samples.

property declared_support: tuple[float, float] | None#

The prior support the grid was generated over, from provenance.

None when the grid carries no provenance, or when what it carries is not the expected shape. A malformed record is not a support range, and reporting one would be worse than reporting nothing.

property inference_metadata: tuple[str, ...]#

Which required provenance fields this grid is missing.

Empty means the grid can be checked against a configuration, and so can be used for inference. Anything else names what a caller has to supply before it can.

The fields are the ones a consistency check needs, not everything the builder records: a stable name, the sample count, the declared support, the prior family and mixture weight, how it was generated, and a hash to pin the arrays.

property unusable_because: tuple[str, ...]#

Everything standing between this grid and inference, in one list.

Missing metadata and failed integrity checks are different faults with the same consequence, and a caller fixing one wants to see the other in the same message rather than after another attempt.

property usable_for_inference: bool#

Whether this grid may be used to evaluate a spectrum.

Two conditions, and both are required:

  • every field in REQUIRED_GRID_METADATA is present; and

  • the recorded digests were checked against this file and these arrays and matched.

Metadata presence alone is not enough. A sidecar can be copied from a different grid, or the arrays edited after it was written, and the field names would look correct either way. A grid loaded from a bare .npz, or one built by hand, fails the second condition because nothing was ever verified.

Such a grid still loads and can be inspected – the arrays are all there. It is inference that Finder refuses.

property declared_prior_alpha: float | None#

The PW14 mixture weight the grid was generated with, from provenance.

None when the grid carries no provenance, which is the case for an .npz loaded without its sidecar. A caller that needs the check must keep the sidecar; see load_sample_grid().

property is_verified: bool#

Whether this grid’s identity with a deployed production grid is proven.

Fails closed. Only a status on VERIFIED_IDENTITY_STATUSES counts as verified; missing, unknown, misspelled, regenerated or "unverified" all return False. An earlier implementation asked whether the status failed to start with "regenerated", which would have treated a typo — or any new status string — as verified.

sample_redshifts(z_min, z_max)[source]#

Stretch the stored offsets onto a spectrum’s absorber search window.

Mirrors the reference: z_i = z_min + (z_max - z_min) * offset_i.

Parameters:
Return type:

ndarray

gp_dla_finder.samples.load_sample_grid(name='pw14_172_225_50000', *, path=None)[source]#

Load a QMC absorber sample grid.

With name, one of the packaged grids. With path, an .npz built by tools/build_sample_grid.py anywhere on disk — you do not need to put a custom grid inside the installed package.

An external grid keeps its provenance. The builder writes <name>.json beside the .npz; if that sidecar is present it is loaded too, and the grid takes its recorded name.

Without the sidecar the arrays still load, but the grid cannot be used for inference. It has no declared support, prior mixture or stable name, so nothing can check it against a configuration, and a run using it would record a configuration describing a different grid. Inspect such a grid freely; Finder refuses it. Check AbsorberSampleGrid.usable_for_inference if you need to know before constructing one.

Parameters:
Return type:

AbsorberSampleGrid

Trained Gaussian-process quasar-emission models.

A GPModel is the learned prior over unabsorbed quasar spectra: a mean flux mu on a rest-frame wavelength grid, a rank-k factor M giving the low-rank covariance K = M M^T, a per-pixel amplitude log_omega for the Lyman-forest absorption noise, and three learned scalars (c_0, tau_0, beta) parameterising that noise.

Model-coupled quantities, including the rank, rest-frame grid, and flux normalization band, are properties of the model. They are not global defaults; the package checks them against the configuration before inference.

Assets#

Packaged models live in gp_dla_finder/data/models as an .npz of inference arrays plus a .json of provenance. Arrays may be stored as float32 where that is bitwise lossless (see tools/convert_model.py); they are always loaded as float64 so the arithmetic matches the reference implementation exactly.

class gp_dla_finder.model.GPModel(name, rest_wavelengths, mu, M, log_omega, log_c_0, log_tau_0, log_beta, normalization_min_lambda=None, normalization_max_lambda=None, provenance=<factory>)[source]#

A trained GP prior over quasar emission.

Parameters:
name#

Asset name, or "<external>" for a model loaded from an arbitrary path.

Type:

str

rest_wavelengths#

Rest-frame grid the model is defined on, angstroms, shape (W,).

Type:

numpy.ndarray

mu#

Mean flux on that grid, shape (W,).

Type:

numpy.ndarray

M#

Low-rank covariance factor, shape (W, k); K = M M^T.

Type:

numpy.ndarray

log_omega#

Log amplitude of the absorption-noise term, shape (W,).

Type:

numpy.ndarray

log_c_0, log_tau_0, log_beta

Learned scalars of the Lyman-forest noise model.

normalization_min_lambda, normalization_max_lambda

Rest-frame band the training spectra were normalised over. None for older models that do not record it; nan if trained unnormalised.

provenance#

Read-only provenance mapping; empty for externally loaded models.

Type:

collections.abc.Mapping[str, Any]

property rank: int#

Rank k of the low-rank covariance factor.

covers(min_lambda, max_lambda)[source]#

Whether the model’s grid spans a rest-frame window, inclusive.

Parameters:
Return type:

bool

gp_dla_finder.model.load_model(name='phase2_2lpt_loa124_nohcd_nobal_wide_m', *, path=None)[source]#

Load a trained GP model.

Parameters:

  • name (str) – Bundled asset name. Ignored when path is given.

  • path (str | Path | None) – Load from disk instead: either a packaged .npz or an HDF5 .h5 / MATLAB v7.3 .mat file from the reference pipeline. HDF5 input needs the optional legacy extra (pip install 'gp_dla_finder[legacy]').

Return type:

GPModel

Examples

>>> model = load_model()
>>> model.rank
30
gp_dla_finder.model.model_provenance(name='phase2_2lpt_loa124_nohcd_nobal_wide_m')[source]#

Provenance record for a bundled model: source checksum, training run, grid.

Parameters:

name (str)

Return type:

Mapping[str, Any]

Catalogue output#

Catalog schemas for both DESI compatibility and complete run accounting.

The deployed DESI DLA catalog is a flat, one-row-per-DLA FITS table, and downstream analyses already read it. We preserve its column names, data types, row granularity, and identifier conventions as the compatibility surface.

A DLA-only table has one important blind spot: a spectrum with no selected DLA produces no row. That is enough for a DLA list, but it cannot distinguish a null result from a quality rejection, an inference failure, or a spectrum that was never processed. The three schemas are:

ABSORBER_SCHEMA

one row per absorber. The legacy compatibility surface.

SPECTRUM_SCHEMA

one row per attempted spectrum, carrying its status and reason code. Null detections and rejections live here.

RUN_SCHEMA

run-level configuration and provenance, recorded once rather than repeated on every absorber row.

Two FITS products are built from them: a strict legacy export containing only the absorber table with exactly the historical columns, for readers that require that schema, and an extended product carrying all three. Both are flat: one row per absorber, one row per spectrum, no nesting and no variable-length arrays.

FITS is the compact catalogue, not the inference record#

The FITS product is the compact summary used by the DESI workflow. The full M0/M1/M2 ladder – priors, evidences, posteriors, the selected model, and which rungs were evaluated – travels in the structured JSON output, gp_dla_finder.io.structured.write_structured_results(). The in-memory Result keeps the same information either way.

A spectrum for which two absorbers were selected contributes two ordinary rows sharing a TARGETID, with DLAID values <targetid>000 and <targetid>001. That is the whole multi-absorber representation in FITS.

Legacy semantics#

The reference implementation defines the probability columns independently (dlasearch.py, in the per-absorber loop):

  • P_DLA and P_NULL are spectrum-level scalars, repeated unchanged on every absorber row of that spectrum;

  • LOGP_DLA is log_posteriors_dla[n] — a per-absorber-index quantity, the log posterior of the n-absorber model, not the log of P_DLA;

  • LOGP_NULL is a spectrum-level scalar;

  • MODEL_P is model_posteriors[1 + num_subdla + n], the posterior of the specific absorber-count model.

These quantities coincide only in special cases and diverge when more than one absorber is modeled. We therefore preserve all five independently.

Z_DLA_ERR and NHI_ERR#

Present in both products, because the reference produces them and downstream readers expect the columns. The row model accepts and preserves them when a validated estimator supplies a value. Until this package has such an estimator they are written as documented NaN. The optional emcee validation samples do not by themselves define a production uncertainty.

Other invariants#

  • No sample arrays, ever. Catalogue output is summary-only, enforced by an explicit field allowlist rather than by guessing from array shapes.

  • A FILTER-derived value may occupy an evidence field, but only when it is labelled: every row carries GPDLF_EVIDENCE_MODE, and constructing a row without a valid label raises. Run-level labelling alone would be insufficient, because one file can hold both modes.

class gp_dla_finder.catalogue.ModelRow(targetid, model_index, log_evidence, log_prior, posterior, evaluated, selected)[source]#

One model of one spectrum’s ladder.

Parameters:
class gp_dla_finder.catalogue.Column(name, dtype, unit, description, legacy=False, provisional=False)[source]#

One catalogue column, with the intent behind it recorded.

Parameters:
dtype: str#

FITS/numpy datatype code, as astropy’s Table understands it.

legacy: bool = False#

True when the column exists to satisfy an existing downstream reader and its meaning is fixed by history rather than by this package.

provisional: bool = False#

True when the column may be absent from the extended product because the quantity has no validated estimator yet.

gp_dla_finder.catalogue.schema_for(table)[source]#

Return a named schema.

Raises:

KeyError – If table is unknown.

Parameters:

table (str)

Return type:

tuple[Column, …]

class gp_dla_finder.catalogue.AbsorberRow(targetid, dlaid, z_qso, z_dla, nhi, p_dla, p_null, logp_dla, logp_null, model_p, log_evidence_absorber, log_evidence_null, evidence_mode, z_dla_err=nan, nhi_err=nan, ra=nan, dec=nan, snr_forest=nan, snr_redside=nan, dlaflag=0, screening_score=nan, model_index=1)[source]#

One absorber. Summary quantities only, by construction.

The five legacy probability fields are carried independently. None is derived from another, because in the reference they are different quantities: p_dla/p_null are spectrum-level, logp_dla is per-absorber-index, and model_p is the posterior of that absorber-count model.

Parameters:
p_dla: float#

Spectrum-level posterior of the absorber model, repeated on each row.

p_null: float#

Spectrum-level posterior of the null model.

logp_dla: float#

log posterior of the n-absorber model. NOT log(p_dla).

logp_null: float#

Spectrum-level log posterior of the null model.

model_p: float#

Posterior of this absorber-count model.

evidence_mode: str#

see EVIDENCE_MODES.

Type:

“exact” or “filter”. Required

z_dla_err: float = nan#

Uncertainties, when a validated estimator supplies them. NaN otherwise.

model_index: int = 1#

1 for the one-absorber model, 2 for a member of the best two-absorber pair.

Not a FITS column, deliberately. The DESI product is flat and stays flat; membership is implicit there in how many rows share a TARGETID. This field exists so the structured JSON can state it, and so it does not have to be re-derived from the row index – which gives the wrong answer, because both members of a selected M2 pair belong to M2 while their row indices are 0 and 1.

Type:

Which absorber-count model this candidate belongs to

property log_bayes_factor: float#

log evidence ratio, absorber over null, in this row’s mode.

class gp_dla_finder.catalogue.Catalogue(absorbers=<factory>, spectra=<factory>, models=<factory>, run=<factory>)[source]#

Absorber rows, spectrum rows, the model ladder, and the run record.

Parameters:
models: Sequence[ModelRow]#

One row per model per spectrum, when a ladder was evaluated.

class gp_dla_finder.catalogue.SpectrumRow(targetid, z_qso, status, reason='', n_absorbers=0, p_absorber=nan, log_evidence_null=nan, log_evidence_absorber=nan, evidence_mode='', quality_fraction=nan, n_usable_pixels=0, ra=nan, dec=nan, n_evaluated=0, screening_score=nan, screening_n_evaluated=0)[source]#

One attempted spectrum, whatever happened to it.

Parameters:
n_evaluated: int = 0#

Samples actually evaluated for THIS spectrum. Recorded per spectrum because FILTER can stop after different counts for different spectra, and a run-level number would falsely claim one count for all of them.

screening_score: float = nan#

FILTER-prefix log Bayes factor. Carried on the SPECTRUM row, not only on absorber rows, so a screened spectrum that did not pass the detection threshold still records what screening said about it – it has no absorber row to carry it. NaN means no screening stage ran.

screening_n_evaluated: int = 0#

Samples the screening stage evaluated; 0 when it did not run.

FITS writers for the compact legacy and extended catalog products.

Use the legacy product when old downstream code needs the historical table. Use the extended product when you also need to know what happened to every spectrum.

write_legacy_catalogue()

exactly the historical DESI columns, in the historical order, one row per absorber, in a single binary table. This product supports downstream software that expects the historical schema. It cannot represent a null or rejected spectrum, which is why the extended product also exists.

write_catalogue()

the same absorber table plus a per-spectrum status table and a run/provenance extension.

The no-samples guarantee#

Catalog output contains no QMC or posterior sample arrays. An explicit schema allowlist enforces this restriction.

astropy is an optional dependency. It is imported inside the functions so the inference core keeps needing nothing but NumPy and SciPy.

gp_dla_finder.io.fits.read_catalogue_metadata(path)[source]#

Read schema version and run provenance back out of a written catalogue.

Raises:

ValueError – If the file’s schema major version is not one this package understands. A reader that silently accepts an unknown major version is how a column whose meaning changed gets read as though it had not.

Parameters:

path (str | Path)

Return type:

Mapping[str, object]

gp_dla_finder.io.fits.write_catalogue(path, catalogue)[source]#

Write the extended product: absorbers, per-spectrum status, and run info.

The absorber HDU keeps the legacy columns and order, so a reader that only knows the historical schema still works on it, and adds the new, unambiguously named evidence fields after them. Three flat tables, and no model-ladder HDU: the ladder is gp_dla_finder.io.structured.write_structured_results().

Parameters:
Return type:

Path

gp_dla_finder.io.fits.write_legacy_catalogue(path, catalogue)[source]#

Write the strict legacy view: historical columns only, one HDU.

The DESI-compatible catalogue product. Flat, one row per absorber: a spectrum with two selected absorbers contributes two ordinary rows sharing a TARGETID, with contiguous DLAID values. There is no nesting, no variable-length column and no model-ladder table.

Z_DLA_ERR and NHI_ERR are written as NaN, documented as such in the header. No estimator for them has been chosen or validated, and emitting a plausible-looking number would be worse than emitting nothing.

This product cannot represent a spectrum with no absorber. Use write_catalogue() when that matters — which is whenever the question is about a population rather than a list. Neither product carries the model ladder: that is gp_dla_finder.io.structured.write_structured_results().

Parameters:
Return type:

Path

The structured result: everything the flat FITS catalogue does not carry.

The FITS product follows the DESI catalogue: flat, fixed columns, and one row per absorber. We keep that layout because downstream analyses already read it.

The M0/M1/M2 ladder does not fit a flat table. Its length varies with the run, and its posteriors are model-level rather than absorber-level. So it travels here instead, in JSON:

  • every model’s log evidence, log prior and posterior;

  • whether each rung was evaluated at all, and which was selected;

  • which absorber candidate belongs to which model – taken from the candidate itself, not inferred from its position in the list;

  • the run provenance, unchanged from the FITS run record.

JSON works here because the structured result is small – one object per spectrum, not per pixel – and it needs no dependency beyond the standard library. If you later retain posterior samples or full inference state, use a separate HDF5 product instead.

Strict JSON, and how a missing value stays distinguishable#

The output is strict RFC 8259: every non-finite number is written as null, and serialisation runs with allow_nan=False so a future non-finite value cannot silently reintroduce a bare NaN or Infinity literal. Any JSON parser reads these files, not only Python’s.

null on its own would lose something, so it does not carry the meaning alone. Each model states whether it ran:

evaluated = false, posterior = null   # the rung was never evaluated
evaluated = true,  posterior = 0.0    # evaluated, and its posterior is zero

That is the distinction the ladder exists to record, and evaluated is what records it. A null says only “no number here”.

gp_dla_finder.io.structured.read_structured_results(path)[source]#

Read a structured result back.

Raises:

ValueError – If the document is not a structured result, or its major format version is one this package does not understand. Reading an unknown major version silently is how a field whose meaning changed gets misread.

Parameters:

path (str | Path)

Return type:

Mapping[str, object]

gp_dla_finder.io.structured.selected_models(payload)[source]#

{targetid: 'M0' | 'M1' | ...} for every spectrum that has a ladder.

A spectrum with no ladder is absent rather than defaulted: a run that never evaluated the models has no selection to report, and reporting M0 for it would claim a measurement that was not made.

Parameters:

payload (Mapping[str, object])

Return type:

dict[int, str]

gp_dla_finder.io.structured.structured_payload(catalogue)[source]#

The JSON document, as a dict, without writing it.

Separate from the writer so it can be embedded in another document, sent over a wire, or asserted against in a test without a temporary file.

Parameters:

catalogue (Catalogue)

Return type:

dict

gp_dla_finder.io.structured.write_structured_results(path, catalogue)[source]#

Write the complete structured result, ladder included.

The companion to the FITS catalogue rather than a replacement for it: write both when the ladder matters, and point analyses that only need a DLA list at the FITS file.

Parameters:
Return type:

Path

Policies and diagnostics#

Named data-quality policies, kept separate from basic spectrum validation.

There are two different questions here: can we run on this spectrum, and should it enter a selected catalog? This module keeps those questions separate.

Structural validation asks whether the input is a usable spectrum: correct shapes, increasing wavelengths, non-negative inverse variance, and whether enough of it survives masking to compute anything. That lives in gp_dla_finder.gp.spectrum and is not negotiable.

Quality selection asks whether a well-formed spectrum is suitable for a catalog. This is a survey decision rather than an inference decision. The deployed DESI pipeline requires at least 20% of pixels in rest-frame 900–1230 Å to have non-zero inverse variance. This module provides that named policy.

We enforce three properties:

  • the low-level inference API never applies a survey cut on its own. A configuration selects a named policy, or none at all;

  • a quality rejection is not an inference failure. It has a separate status and reason code;

  • whether a policy ran, which one, its threshold and window, the measured usable fraction and the verdict all go into result provenance. A catalog built with a cut and one built without must not be indistinguishable after the fact.

Disabling the policy transfers quality selection to the caller. Structural validation alone is not a scientifically validated catalog-selection rule and must not be presented as one.

class gp_dla_finder.quality.QualityAssessment(policy, policy_version, passed, usable_fraction, n_usable, n_in_window, threshold, rest_lambda_min, rest_lambda_max)[source]#

What a policy measured on one spectrum.

Parameters:
provenance()[source]#

Flat, JSON-friendly record for a result’s provenance block.

Return type:

Mapping[str, object]

class gp_dla_finder.quality.QualityPolicy(name, version, summary, min_usable_fraction, rest_lambda_min, rest_lambda_max)[source]#

A named, versioned catalogue-selection rule.

Parameters:
version: str#

Bumped whenever the meaning changes. A policy must never decide differently under a fixed name and version.

min_usable_fraction: float#

Minimum fraction of pixels in the window that must be usable.

rest_lambda_min: float#

Rest-frame window the fraction is measured over, angstroms.

assess(spectrum)[source]#

Measure spectrum against this policy.

Measures rather than decides: the caller reads QualityAssessment.passed. A policy never raises on a well-formed spectrum, because “this does not belong in the catalogue” is a result, not an error.

Return type:

QualityAssessment

gp_dla_finder.quality.quality_policy(name)[source]#

Return a named quality policy.

Raises:

KeyError – If name is unknown.

Parameters:

name (str)

Return type:

QualityPolicy

Named compatibility profiles for reproducing the reference arithmetic.

Reproducing the reference implementation bitwise required reproducing two pieces of arithmetic that are mathematically no-ops:

rest_frame_round_trip

The reference derives observed wavelengths as wave / (1 + z) * (1 + z) rather than reusing the input array. In floating point that round trip is not the identity: it moves wavelengths by ~1e-13 A, which propagates into the effective optical depth and then into every evidence.

log_norm_round_trip

Every per-sample log-likelihood carries - log(N) and the quasi-Monte-Carlo estimator adds + log(N) back. The pair cancels analytically. It does not cancel in floating point, because the subtraction happens before a log-mean-exp and the addition after it.

Neither operation changes the analytic model. We retain them because this port must reproduce a pinned reference bitwise; comparisons at the $10^{-12}$ level identified both effects.

We isolate these operations in a named, versioned compatibility profile and record the profile in result provenance. Two consequences follow:

  • the package reproduces the arithmetic of a specific commit. If the reference is ever corrected, REFERENCE_D5B306E6 must not be edited – add a new profile, so old results stay reproducible;

  • a clean-arithmetic mode exists (CLEAN) for comparison. Selecting it changes numerical results, and the provenance records that choice.

class gp_dla_finder.compat.CompatibilityProfile(name, version, summary, rest_frame_round_trip, log_norm_round_trip, reference_repo=None, reference_commit=None)[source]#

A named set of reference-fidelity arithmetic choices.

Parameters:

  • name (str)

  • version (str)

  • summary (str)

  • rest_frame_round_trip (bool)

  • log_norm_round_trip (bool)

  • reference_repo (str | None)

  • reference_commit (str | None)

name: str#

Stable identifier, recorded in result provenance.

version: str#

Bumped whenever the meaning of the profile changes. A profile’s numbers must never change under a fixed name and version.

rest_frame_round_trip: bool#

Reproduce observed = rest * (1 + z) instead of reusing the input grid.

log_norm_round_trip: bool#

Carry - log(N) per sample and add + log(N) back in the estimator.

reference_repo: str | None = None#

The reference this profile was measured against, if any.

provenance()[source]#

Flat, JSON-friendly record for a result’s provenance block.

Return type:

Mapping[str, object]

gp_dla_finder.compat.compatibility_profile(name)[source]#

Return a named compatibility profile.

Raises:

KeyError – If name is unknown. Substituting a profile would silently change the arithmetic, so there is no fallback.

Parameters:

name (str)

Return type:

CompatibilityProfile

BLAS thread diagnostics for a workload where more threads can be slower.

The GP likelihood multiplies matrices of shape (1105, 30) — about 2 Mflop per product. Whether extra BLAS threads help or hurt at that size is not obvious and is not monotonic, so it was measured rather than assumed. On the audited host, freshly started 2-, 4- and 8-thread processes all beat one thread.

Two sweeps on the same 10-core machine (OpenBLAS 0.3.23, order-balanced), and they do not agree:

  • 1 thread: live pool 0.2438 ms, fresh process 0.2341 ms

  • 2 threads: live pool 0.2196 (-9.9 %), fresh process 0.1882 (-19.6 %)

  • 4 threads: live pool 0.1693 (-30.6 %), fresh process 0.1547 (-33.9 %)

  • 8 threads: live pool 0.4272 (+75 %), fresh process 0.2040 (-12.9 %)

  • 10 threads: live pool 1.7508 (+618 %), fresh process 1.6164 (+591 %)

The separate-process measurements best represent normal use because thread limits are set before a process starts. Resizing a live pool makes eight threads appear 75% slower, while a process started with eight threads remains 13% faster than one thread. We therefore keep the live-pool sweep as a separate diagnostic.

On this host, a few threads improve performance, eight threads still help, and performance degrades when the pool uses all ten cores. This result describes one machine and should not be treated as a universal optimum.

What this does and does not show#

It shows single-process scaling as a function of pool size. It does not show the multi-worker regime (worker processes x BLAS threads > cores), which has not been measured here. One BLAS thread per worker remains sound guidance for catalogue production, but that is reasoning about how independent processes contend, not a measurement this project has made.

A shared 2-core CI runner measured 20-23 % overhead at 2 threads on one occasion and 1.6 % on another, with different host CPUs — and on a 2-core machine, 2 threads is the whole machine. Small-pool behaviour is sensitive to the host.

Runtime integration#

Finder calls warn_once_about_blas_threads() at the inference boundary. The diagnostic runs once per process, not at import time or inside the likelihood loop.

What this module does, and deliberately does not#

It detects and warns once. It does not set environment variables, impose a process-wide limit, or change the caller’s thread policy.

Detection needs threadpoolctl, an optional dependency (pip install 'gp_dla_finder[performance]'). Without it the package works normally and simply cannot offer the diagnostic: no warning, no behavior change, no noise.

The caveat that belongs with the number#

Neither measurement settles what a production Linux machine will do. The Linux number above comes from a 2-core shared CI runner, which cannot spawn enough threads to produce the failure mode at all — it bounds the small-pool case and says nothing about a 64-core node. That is why the warning says “measured performance risk” rather than “your run is four times slower”.

exception gp_dla_finder.performance.BLASPerformanceWarning[source]#

A BLAS thread configuration likely to be slow for this workload.

Its own class so it can be silenced without also silencing scientific or numerical warnings:

warnings.filterwarnings("ignore", category=BLASPerformanceWarning)

It is a performance advisory. It never indicates an inference error, and it is not proof that a particular workload is slower — only that the configuration carries a measured risk for small-matrix work.

gp_dla_finder.performance.blas_thread_report()[source]#

What BLAS pools this process has, as far as they can be detected.

Returns a mapping with available false when threadpoolctl is not installed, rather than guessing at the active thread pools. NumPy and SciPy frequently ship separate BLAS libraries, so pools may hold more than one entry with different thread counts, and the advisory keys off the largest.

Return type:

Mapping[str, object]

gp_dla_finder.performance.warn_once_about_blas_threads(stacklevel=3)[source]#

Warn at most once if a threaded BLAS looks badly configured for this work.

Call this from the high-level inference boundary on first use — not at import, and never per likelihood evaluation. Returns whether a warning was issued, which is what makes it testable.

Silent, and returns False, when threadpoolctl is absent, when the pools look fine, or when it has already fired.

Parameters:

stacklevel (int)

Return type:

bool

Exception taxonomy for package and spectrum-level failures.

Callers need to distinguish four situations reliably:

  • the run is misconfigured or its assets are wrong — nothing will work, fail the batch (ConfigurationError, AssetError);

  • this particular spectrum is structurally invalid — the caller passed something that is not a usable spectrum (SpectrumError);

  • this spectrum is well formed but cannot support inference — fully masked, no normalisation coverage, too few usable pixels. That is not an error and not a non-detection; see InsufficientData;

  • the numerics failed on otherwise valid input (NumericalError).

A processing failure must never be reported as a no-absorber result. Those are different states and downstream population statistics depend on telling them apart.

exception gp_dla_finder.errors.AssetError[source]#

A packaged or supplied asset is missing, malformed, or inconsistent.

Global, like ConfigurationError.

exception gp_dla_finder.errors.ConfigurationError[source]#

The configuration is invalid, or incompatible with the chosen assets.

Global: it invalidates the whole run, not one spectrum.

exception gp_dla_finder.errors.GPDLAError[source]#

Base class for every error this package raises deliberately.

exception gp_dla_finder.errors.NumericalError[source]#

The numerics failed on structurally valid input.

Spectrum-local. A non-positive-definite covariance, an all-NaN likelihood slice. Raised rather than returned as a result, because a numerical failure is not a scientific conclusion.

exception gp_dla_finder.errors.SpectrumError[source]#

The input is not a structurally valid spectrum.

Spectrum-local: a batch layer may record the failure and continue. Wrong shapes, non-monotonic wavelengths, negative inverse variance, a non-finite redshift. Distinct from a valid spectrum that merely cannot support inference.

On this page
Show Source