mixle.analysis package

Applied statistical methods that are not probability-distribution families.

These are estimators / analysis routines (extreme-value, kernel density, species coverage, geostatistics, rank aggregation, spatial mixtures, max-stable processes, covariance shrinkage, mine economics, health/exposure monitoring) that operate on data but are not SequenceEncodableProbabilityDistribution families and are not used by any distribution. They were previously scattered under mixle.stats; collected here so the distributions package stays focused on distribution families.

This package’s __all__ is the union of every submodule’s own declared __all__ (MXR-080-1593). The façade previously omitted seven names – including result types returned by exported functions (TransitionRiskResult, SmithMaxStableFit, VoiEstimate) and argument protocols needed to call them (WaterBudget, GaussianObservationModel) – so from mixle.analysis import ... could reach a function but not the type it hands back or the type it demands. A submodule’s __all__ is that module’s declaration of its intentional stable surface; re-exporting anything less here makes the package a lossy view of it rather than the entry point it is documented to be. mixle/tests/api_consistency_test.py asserts the two stay in agreement.

class GPDFit(shape, scale, threshold, n_exceedances, n_total, method, n_dropped_nonpositive=0)

Bases: object

Fitted Generalized Pareto Distribution for exceedances over a threshold.

Parameters:

• shape (float)

• scale (float)

• threshold (float)

• n_exceedances (int)

• n_total (int)

• method (str)

• n_dropped_nonpositive (int)

shape

tail index xi (> 0 heavy/Pareto tail, = 0 exponential, < 0 bounded).

Type:

float

scale

scale beta (> 0).

Type:

float

threshold

the threshold u the exceedances were measured over.

Type:

float

n_exceedances / n_total

exceedance and full-sample sizes (for return levels).

method

"mle" or "pwm".

Type:

str

n_dropped_nonpositive

count of finite, non-positive (<= 0) input values excluded because they are not actually above the threshold (a receipt for gpd_fit()’s filtering; NaN and non-finite input are rejected outright rather than silently dropped, so they never contribute to this count).

Type:

int

property endpoint: float

Right endpoint of the support (finite when shape < 0, else inf).

gpd_fit(exceedances, *, threshold=0.0, method='mle', n_total=None)

Fit a Generalized Pareto Distribution to threshold exceedances.

Parameters:

• exceedances (ndarray) – the excesses x - u for observations above the threshold (all positive). NaN or non-finite entries are invalid data, not sub-threshold observations, and raise. Finite entries <= 0 are legitimately not above the threshold, so they are dropped and the count is receipted on the returned GPDFit (n_dropped_nonpositive) rather than silently vanishing. If you have raw data, use peaks_over_threshold() instead.

• threshold (float) – the threshold u (stored for return-level computation).

• method (str) – "mle" (Nelder-Mead on the GPD likelihood) or "pwm" (probability-weighted moments, closed form, robust for xi < 0.5).

• n_total (int | None) – full sample size before thresholding (defaults to the number of valid, positive exceedances; if given explicitly it must be at least that many).

Returns:

A GPDFit.

Raises:

ValueError – exceedances contains NaN/non-finite values; fewer than two positive exceedances remain after filtering; n_total is smaller than the exceedance count; method is not "mle"/"pwm"; the MLE optimizer fails to converge; or the fitted (shape, scale) is not a valid GPD (see _validate_gpd_params()).

Return type:

GPDFit

peaks_over_threshold(data, threshold, *, method='mle')

Peaks-over-threshold: select exceedances above threshold and fit a GPD to the excesses.

Parameters:

• data (ndarray)

• threshold (float)

• method (str)

Return type:

GPDFit

return_level(fit, period)

POT return level: the level exceeded on average once per period observations.

x_m = u + (beta/xi) [ (m zeta_u)^xi - 1 ] with zeta_u = n_exceed/n_total the exceedance rate (period = m). For xi = 0 it reduces to u + beta log(m zeta_u).

Raises:

ValueError – period is not positive (a non-positive return period is not meaningful – m = 0 sends the level to +/-infinity and negative m is undefined for non-integer xi).

Parameters:

• fit (GPDFit)

• period (float)

Return type:

float

hill_estimator(data, k)

Hill estimator of the tail index xi = 1/alpha from the top k order statistics.

xi_hat = (1/k) sum_{i=1}^{k} log(X_(n-i+1) / X_(n-k)) – consistent for heavy (Pareto-type, xi > 0) tails. For a Pareto tail with exponent alpha this estimates 1/alpha.

Parameters:

• data (ndarray)

• k (int)

Return type:

float

moment_estimator(data, k)

Dekkers–Einmahl–de Haan moment estimator of the extreme-value index (any tail sign).

Generalises Hill to xi of either sign by combining the first two log-moments of the top k exceedances; works for heavy, light, and bounded (xi < 0) tails. Needs k >= 2: with a single log-spacing the second moment is identically the square of the first (zero variance), which sends the estimator’s denominator to zero.

Parameters:

• data (ndarray)

• k (int)

Return type:

float

mean_residual_life(data, thresholds)

Mean-excess (mean-residual-life) function for POT threshold selection.

e(u) = mean(X - u | X > u). Over a range where the GPD fits, e(u) is approximately linear in u (slope xi/(1-xi)); the lowest threshold from which the plot is linear is the choice.

Returns:

{'threshold', 'mean_excess', 'n_exceed'}.

Parameters:

• data (ndarray)

• thresholds (ndarray)

Return type:

dict[str, ndarray]

endpoint_estimator(data, k, *, method='gpd')

Estimate the finite right endpoint of a bounded support (frontier / boundary estimation).

Fits a GPD to the top k exceedances over X_(n-k); when the tail index xi is negative the support is bounded and the endpoint is x+ = X_(n-k) - beta/xi (which, by the GPD support constraint, necessarily exceeds the observed maximum). Generic to econometric frontier analysis, reliability limits, and image-edge localisation. Returns inf if the estimated tail is unbounded (xi >= 0).

Parameters:

• data (ndarray) – the sample.

• k (int) – number of upper order statistics (exceedances) used.

• method (str) – "gpd" – GPD-MLE endpoint. No other method is implemented.

Returns:

The estimated right endpoint (inf if unbounded).

Raises:

ValueError – method is not "gpd", k is out of range, or the upper order statistics used are non-finite.

Return type:

float

n_records(data)

Number of upper records. For an i.i.d. sequence of length n the expectation is H_n.

Parameters:

data (ndarray)

Return type:

int

record_times(data)

Indices at which a new running maximum (upper record) occurs, including the first observation.

data must be a finite one-dimensional series, like the tail estimators above. A single NaN used to silently truncate the answer rather than raise: np.maximum.accumulate propagates the NaN forward, every later x[i] > running[i-1] comparison against it is false, and every record after the NaN disappears – record_times([1, nan, 2, 3]) returned just [0], which n_records() then reported as an authoritative count of one.

An empty data has no observations and so, vacuously, no records: returns an empty index array rather than raising.

Parameters:

data (ndarray)

Return type:

ndarray

class BMDResult(bmd, bmdl, bmr, model, dof, status='ok', bmd_se=nan, _coef=<factory>)

Bases: object

A fitted benchmark-dose analysis: the BMD, its lower confidence bound (BMDL), and fit metadata.

status is one of:

• "ok": the curve fit converged, the BMD was found and bracketed, and the BMDL lower bound was established. bmd and bmdl are real doses.

• "unidentifiable": the curve fit itself did not converge, or the fitted curve never reaches the benchmark target within a bounded search (flat, wrong-signed, or requiring an implausibly large dose). bmd and bmdl are nan – never a search boundary or other placeholder dressed up as a real dose.

• "bmdl_unavailable": the BMD itself was identified (bmd is real), but the BMDL computation did not converge to a valid bound. bmdl is nan.

bmd_se is the delta-method standard error of the BMD (nan unless status == "ok"); bmdl = bmd - z * bmd_se clipped at 0, z the one-sided normal quantile for the requested confidence level.

Parameters:

• bmd (float)

• bmdl (float)

• bmr (float)

• model (str)

• dof (int)

• status (str)

• bmd_se (float)

• _coef (ndarray)

property converged: bool

True only when both the BMD and the BMDL were genuinely established.

benchmark_dose(dose, n_affected, n_total, *, bmr=0.10, model='loglogistic', risk='extra', ci_level=0.95)

Fit a quantal dose-response and report the benchmark dose (BMD) and its lower bound (BMDL).

dose/n_affected/n_total are per-dose-group arrays (n_affected <= n_total). The curve is fit by maximum likelihood (DR-ALG K8); the BMD solves for the dose giving bmr extra (or added) risk over the fitted background rate. The BMDL is the one-sided ci_level lower confidence bound on the BMD by the delta method: the gradient of the BMD with respect to the fitted coefficients is obtained by implicit differentiation of the BMD-defining equation, the coefficient covariance by inverting the observed Fisher information at the MLE, and Var(BMD) by propagating one through the other (grad @ Cov @ grad); BMDL = BMD - z * SE(BMD) with z the one-sided normal quantile for ci_level, clipped at 0 (dose cannot be negative). See BMDResult.status for what happens when any step of that chain fails to converge.

Parameters:

• dose (ndarray)

• n_affected (ndarray)

• n_total (ndarray)

• bmr (float)

• model (str)

• risk (str)

• ci_level (float)

Return type:

BMDResult

rfd_exceedance(exposure, bmd, *, uf=100.0, n=2000, rng=None)

P(exposure > RfD) as an IC-8 DerivedQuantity, where RfD = BMDL / uf (EPA convention).

exposure may be an IC-1 Posterior (the pushforward runs through its own derived_quantity so prior_dominated propagates), an array of exposure draws (resampled to n if its length differs), or a bare scalar (a degenerate point mass).

uf must be finite and strictly positive, n a positive exact sample count, and exposure draws finite and nonnegative – all rejected up front, before any division or resampling. bmd.bmdl must also be finite: a BMDResult with status != "ok" (see BMDResult) has no real BMDL to divide by an uncertainty factor, and a NaN silently compared with draws > nan would quietly evaluate to all-False – a fabricated-looking “0% exceedance” result – rather than surfacing the underlying unidentifiability. The finite/nonnegative exposure check applies uniformly to a Posterior’s own draws too, not only to a plain array/scalar: it lives inside fn below (the pushforward handed to Posterior.derived_quantity), not in a separate top-level check on the raw exposure – otherwise a mis-specified exposure posterior could emit a negative or NaN draw that flowed straight through draws > rfd (silently False for NaN) into a confident-looking “not exceeding” result with no warning.

Parameters:

• exposure (Posterior | np.ndarray)

• bmd (BMDResult)

• uf (float)

• n (int)

• rng (np.random.Generator | None)

Return type:

DerivedQuantity

class DoseResponse(model, params=<factory>)

Bases: object

A named dose-response model: model selects the functional form, params its coefficients.

model in {"loglinear", "logit", "hill", "threshold_linear"}:

• loglinear: P = 1 - exp(-beta * dose) (params: {"beta"}) – the EPA-style linear low-dose cancer/chronic form.

• logit: P = sigmoid(a * dose + b) (params: {"a", "b"}, both optional).

• hill: P = emax * dose^n / (ec50^n + dose^n) (params: {"ec50"}, "emax"/"n" optional) – saturating receptor-occupancy form.

• threshold_linear: P = clip(slope * (dose - threshold), 0, 1) (params: {"slope"}, "threshold" optional) – no response below threshold.

No regulatory dose-response table ships here – params are supplied by the caller (see the module Non-goals).

Frozen, holding its own copy of params. The eager domain validation below is only meaningful if the validated coefficients are the ones actually used: with a mutable record and a shared dict, a caller could construct a valid model and then swap model or edit params – and every later response_fn() would rebuild the closure from coefficients that never passed the construction-time gate.

Parameters:

• model (str)

• params (dict[str, Any])

response_fn()

The elementwise dose -> outcome-probability function for this model + params.

Return type:

Callable[[ndarray], ndarray]

probability(dose, *, n=2000, rng)

Push dose through the dose-response model into an outcome-probability DerivedQuantity.

dose may be an IC-1 Posterior over exposure (the pushforward runs through the posterior’s own derived_quantity, so prior_dominated propagates from the exposure uncertainty), an array of dose draws/ensemble members (“array-with-UQ”, resampled to n if its length differs), or a bare scalar (a degenerate point mass – the returned quantity still carries a (trivial) credible interval). Bare ensembles are one-dimensional. Posterior draws explicitly use shape (n,) for scalar dose or (n, n_receptors) for a receptor vector. Every dose must be finite and nonnegative, and n is an exact non-Boolean positive integer.

Parameters:

• dose (Any)

• n (int)

• rng (np.random.Generator)

Return type:

DerivedQuantity

cumulative_exposure(series, dt, *, decay=0.0)

Time-integrated exposure (trapezoidal rule), with optional first-order biological decay.

decay=0 is the plain trapezoidal integral of series over its dt-spaced timesteps (area under the exposure-rate curve). decay > 0 discounts each sample toward the final timestep by exp(-decay * (t_end - t)) before integrating – the way a biological half-life would – so a spike long ago contributes less to the current cumulative body burden than an equally large spike near the end of the series. Feeds a chronic dose-response evaluation (e.g. via DoseResponse.probability()).

series must be a one-dimensional finite nonnegative time series; multidimensional receptor/time data require an explicit aggregation before this scalar integration. dt must be finite and positive (MXR-080-0098: a zero or negative timestep makes the integral either trivially zero or sign-flipped, silently), and decay finite and non-negative (a negative “decay” would amplify, not discount, older readings without bound).

Parameters:

• series (ndarray)

• dt (float)

• decay (float)

Return type:

float

population_risk(exposure, dr, *, n, rng)

Aggregate per-receptor dose-response probabilities into an expected-case-count DerivedQuantity.

exposure is a per-receptor dose: an IC-1 Posterior whose draws are (n, n_receptors) dose vectors (K1/K2 transport output propagated through UQ), or a plain (n_receptors,) array of point doses. Each posterior draw is pushed through dr’s response function and summed over receptors, so the returned quantity carries the expected-case-count distribution (not just its mean); a bare array has no exposure uncertainty and yields a degenerate (constant) distribution.

Parameters:

• exposure (Posterior | np.ndarray)

• dr (DoseResponse)

• n (int)

• rng (np.random.Generator)

Return type:

DerivedQuantity

health_liability(risk, *, cost_per_case, discount=0.0)

Price a K3/K4 risk DerivedQuantity into an expected-liability DerivedQuantity (K6, work-plan §7-K).

risk is any IC-1-shaped DerivedQuantity over an expected case count or exceedance probability (population_risk()’s output, or K4’s safety_risk_surface); every draw in risk.samples is multiplied by cost_per_case (dollars per case/incident) and divided by (1 + discount) – a single-period present-value factor (discount=0 is undiscounted; multi-period accounting applies its own per-period factor before summing across periods, this function prices one period/one risk term at a time). Pricing does not change how uncertain the underlying risk is: the returned quantity keeps risk’s sample count/shape and its prior_dominated flag unchanged, so a liability that is prior-dominated upstream is still honestly flagged as such downstream. Handed to J6’s priced_liabilities/risk_adjusted_plan (analysis/valuation.py) as the health_cost callable’s output.

risk.samples must be a nonempty finite nonnegative (n,) case-count series or (n, n_receptors) risk surface; negative harm is never interpreted as a benefit. cost_per_case must be finite and non-negative (a dollar price cannot be negative), and discount finite and strictly greater than -1.0 (MXR-080-0098: the 1 / (1 + discount) present-value factor is only well-defined for a denominator that is positive; discount == -1 divides by zero, and anything below it silently flips the liability’s sign).

Parameters:

• risk (DerivedQuantity)

• cost_per_case (float)

• discount (float)

Return type:

DerivedQuantity

exposure_constraints(options, limits, *, treat_unmodeled_as_safe=False)

Screen candidate operating options against named exposure/exceedance limits (K6).

limits maps an occupational/community exposure metric name (e.g. "silica_pm4", an 8-hour TWA, or an exceedance probability) to the regulatory/policy limit for that metric. options is a list of plain dicts – each one candidate operating configuration, carrying (among whatever other plan data the caller needs, e.g. block cost or grade) a value for zero or more of those metric keys.

Returns a new list (the input dicts are never mutated), one entry per option, each the original option’s key/value pairs plus:

• "status": one of "safe", "violating", or "unknown". "violating" means at least one limited metric the option DOES carry a finite value for breaches its (finite) limit. "unknown" means no limit was confirmed breached, but at least one limited metric could not be evaluated for this option – missing from the option entirely, a non-finite (NaN/Inf) measurement, or a non-finite/missing limit – “not modeled” and “safe” are never the same outcome here, because this is a hard occupational/community exposure screen, not a soft optimization preference (MXR-080-0095). A confirmed violation always dominates an unknown: one metric’s missing data never erases a different metric’s real breach.

• "feasible": True iff status == "safe" – unless treat_unmodeled_as_safe=True, in which case status == "unknown" also counts as feasible (see below). Kept as a plain boolean for the existing contract H4’s two_stage_stochastic_plan filtering relies on.

• "binding": the sorted list of limit names actually confirmed breached (empty unless status == "violating") – naming exactly which limit made the option infeasible.

• "unmodeled": the sorted list of limited metric names that could not be evaluated for this option (empty unless status == "unknown") – naming exactly which evidence was missing.

This is a HARD safety screen, so it FAILS CLOSED by default: an option that breaches no limit it was evaluated against, but leaves one or more limits unevaluated, is "unknown" and therefore feasible=False – unmodeled evidence is never silently treated as compliant, and a NaN measurement (which would otherwise compare False, and so silently pass, against > limit) is treated the same as a missing one. Pass treat_unmodeled_as_safe=True to explicitly opt into the old permissive behavior; this is a deliberate policy override the caller must request, and it defaults to off. A caller filters the returned list down to the feasible options before handing the survivors’ blocks to H4’s two_stage_stochastic_plan (stochastic_opt.py) – an infeasible (violating OR, by default, unknown) option is dropped from the candidate set entirely, so the optimizer never has the chance to select it (see the K6 DoD).

Parameters:

• options (list[dict])

• limits (dict[str, float])

• treat_unmodeled_as_safe (bool)

Return type:

list[dict]

class KDE(data, *, bandwidth='silverman', bounds=None, adaptive=False)

Bases: object

A fitted kernel density estimate, 1-D or multivariate.

Use kde() to construct. Evaluate with evaluate() (or call the instance). Supports a Gaussian kernel – an axis-aligned product kernel with a per-dimension bandwidth for multivariate (n, d) data, rows treated as observations and never reshaped – adaptive bandwidths (any dimension), and, for 1-D data only, normalized boundary correction (bounds; MXR-080-0099).

Parameters:

• data (np.ndarray)

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

• adaptive (bool)

data

the fitted sample as (n, d) observations (d == 1 for a 1-D sample).

n

number of observations.

d

dimensionality (1 for a 1-D sample).

bandwidth

the resolved bandwidth – a scalar float when d == 1, else a length-d numpy.ndarray (one bandwidth per dimension).

evaluate(x)

Density at points x – (m,) for 1-D, (m, d) for multivariate – with reflection boundary correction if bounds was set (1-D only). Always returns a length-m array.

Parameters:

x (ndarray)

Return type:

ndarray

kde(data, *, bandwidth='silverman', bounds=None, adaptive=False)

Construct a kernel density estimate (Gaussian kernel; product kernel for multivariate data).

Parameters:

• data (ndarray) – (n,) sample, or (n, d) for a d-dimensional sample – rows are observations and must stay paired (never flatten an (n, d) array before calling this; MXR-080-0099).

• bandwidth – "silverman", "scott", a single positive finite number (used in every dimension), or – for (n, d) data – a length-d sequence of positive finite numbers (one bandwidth per dimension). Automatic selection needs >= 2 observations with nonzero variation in every dimension.

• bounds – (lo, hi) support limits for normalized boundary correction; either may be None for an unbounded side (e.g. (0.0, None) for a positive variable). Half-lines use an exactly normalized reflected image; finite intervals normalize every truncated kernel. Only supported for 1-D (d == 1) data.

• adaptive (bool) – use Abramson variable bandwidths (wider where the pilot density is low); supported for any dimension.

Returns:

A KDE.

Raises:

ValueError – data is empty, non-finite, or more than 2-D; bandwidth is an unrecognized method name, the wrong length, or not strictly positive and finite; automatic bandwidth selection is requested on fewer than 2 observations or a constant sample in some dimension (MXR-080-0100); or bounds is given together with d > 1, out of order, or non-finite (MXR-080-0099 / MXR-080-0100).

Return type:

KDE

kde_mode(data, *, bandwidth='silverman', bounds=None, grid=None, ci=False, n_boot=500, ci_level=0.95, seed=0)

Estimate the mode (peak location) of a 1-D density, optionally with a bootstrap CI.

Univariate only: a dense grid search over d dimensions is exponential in d and is not implemented here (unlike KDE, this never accepted multivariate data as a documented feature, so an (n, d) input is now rejected explicitly instead of being silently flattened with .ravel() – the same silent-flattening defect KDE itself had; MXR-080-0099).

Parameters:

• data (ndarray) – (n,) sample.

• bandwidth – passed to kde().

• bounds – passed to kde().

• grid (ndarray | None) – evaluation grid; defaults to 512 points spanning the data range.

• ci (bool) – if True return a percentile bootstrap interval for the mode.

• n_boot (int) – bootstrap controls.

• ci_level (float) – bootstrap controls.

• seed (int | RandomState | None) – bootstrap controls.

Return type:

float | dict

Bootstrap policy for degenerate resamples (MXR-080-1587): a resample of a perfectly valid, nonconstant sample can itself come out constant – [0, 0, 1] resamples to [0, 0, 0] roughly one draw in four – and an automatic bandwidth selector has no spread to estimate from, so it raises. That made the interval fail at random, depending on the seed, for a fit that had already succeeded. A degenerate resample is now scored at the ORIGINAL fit’s bandwidth (the sample the bandwidth was selected from is the full one, and a bootstrap replicate is a resample of it, not a new dataset), rather than skipped – skipping would quietly narrow the interval by dropping exactly the least-dispersed replicates. The count of replicates that needed the fallback is reported as n_degenerate_resamples so a caller can see when the interval rests largely on them.

Returns:

The mode (float), or {'mode', 'ci_low', 'ci_high', 'n_boot', 'n_degenerate_resamples'} when ci is True.

Raises:

ValueError – data is empty or not 1-D; or, when ci is True, n_boot is not a positive integer or ci_level is not in the open interval (0, 1) (MXR-080-0100). Also propagates any KDE construction error from kde() on the ORIGINAL sample (e.g. automatic bandwidth selection on a constant input) – that is a genuine property of the data, unlike a degenerate resample, which is handled by the policy above.

Parameters:

• data (ndarray)

• grid (ndarray | None)

• ci (bool)

• n_boot (int)

• ci_level (float)

• seed (int | RandomState | None)

Return type:

float | dict

intensity(events, grid, *, bandwidth='silverman', domain=None, edge_correct=True)

Kernel intensity lambda(t) of a 1-D inhomogeneous Poisson / point process.

Unlike a density (which integrates to 1), the intensity integrates to the expected number of events: lambda_hat(t) = sum_i K_h(t - t_i). With edge_correct the estimate is divided by the fraction of the kernel falling inside domain, removing the downward bias near the boundary.

Univariate only: like kde_mode(), this never accepted multivariate data as a documented feature, so an (m, d) input is now rejected explicitly instead of being silently flattened with .ravel() (MXR-080-0099).

Parameters:

• events (ndarray) – (m,) event locations.

• grid (ndarray) – nonempty finite one-dimensional points t at which to evaluate the intensity.

• bandwidth – "silverman", "scott", or a float.

• domain (tuple[float, float] | None) – (lo, hi) observation window (defaults to the event range); used for edge correction.

• edge_correct (bool) – divide by the in-window kernel mass at each t.

Returns:

The intensity evaluated on grid.

Raises:

ValueError – events is empty, non-finite, or not 1-D; bandwidth is an unrecognized method name or a non-positive/non-finite number; automatic bandwidth selection is requested on fewer than 2 events or constant events; or (when edge_correct is True) the effective domain – explicitly passed, or defaulted to the event range – is non-finite or has lo >= hi (MXR-080-0100; a collapsed-to-a-point domain, e.g. from constant events with an explicit numeric bandwidth, previously produced a silent divide-by-near-zero blowup).

Return type:

ndarray

silverman_bandwidth(data, *, d=1)

Silverman’s rule-of-thumb bandwidth 0.9 min(sd, IQR/1.34) n^{-1/(d+4)} for a 1-D sample.

data is always treated as one 1-D sample (e.g. a single column of a multivariate dataset). d only sets the exponent: it is the dimensionality of the joint sample this bandwidth will be used in a product kernel for – pass it explicitly when selecting one dimension’s bandwidth for a multivariate KDE, so the n^{-1/(d+4)} rate reflects the true joint dimension rather than the univariate default d=1 (n^{-1/5}) rate (MXR-080-0099).

Raises:

ValueError – data has fewer than 2 observations, non-finite entries, or zero variation (see _require_variation(); MXR-080-0100).

Parameters:

• data (ndarray)

• d (int)

Return type:

float

scott_bandwidth(data, *, d=1)

Scott’s rule-of-thumb bandwidth sd * n^{-1/(d+4)} for a 1-D sample.

data is always treated as one 1-D sample; d only sets the exponent – see silverman_bandwidth() (whose d parameter has exactly the same meaning here).

Raises:

ValueError – data has fewer than 2 observations, non-finite entries, or zero variation (see _require_variation(); MXR-080-0100).

Parameters:

• data (ndarray)

• d (int)

Return type:

float

exception CoverageInsufficientDataError

Bases: ValueError

The requested coverage/diversity estimand is not identified by an empty sample.

turing_coverage(counts)

Turing’s sample-coverage estimate and the complementary unseen probability mass.

C = 1 - f1 / n where f1 is the number of singletons and n the total count: the estimated probability that the next observation is a previously seen species. 1 - C = f1/n is the Good–Turing estimate of the total probability of all unseen species.

Returns:

{'coverage', 'unseen_mass', 'n', 'f1'}.

Parameters:

counts (ndarray)

Return type:

dict[str, float | int]

good_turing(counts)

Simple Good–Turing smoothed probabilities (Gale & Sampson 1995).

Reallocates probability from seen to unseen items using the frequencies of frequencies. Empirical Turing estimates r* = (r+1) N_{r+1}/N_r are used for small r and a smoothed log-linear fit S(r) takes over once the two diverge (the Gale switch), giving stable discounts in the sparse tail.

The log-linear fit needs at least two distinct abundance classes to be identifiable. With exactly one – e.g. an all-singleton sample like [1, 1, 1], a canonical Good–Turing input, not a malformed one – there is nothing to fit a slope against, so smoothing is skipped and the seen mass is reallocated by raw frequency instead (the Simple Good–Turing small-support fallback; previously this raised a LinAlgError out of np.polyfit). With no observed counts at all, coverage is not identifiable and the common CoverageInsufficientDataError is raised.

Parameters:

counts (ndarray) – per-species abundances (zeros ignored).

Returns:

{'p0', 'proba', 'r_star', 'r', 'insufficient_evidence', 'reason'} – p0 is the total probability assigned to unseen species; proba are the smoothed probabilities of the input species (aligned to the positive entries of counts, summing to 1 - p0); r_star / r are the discounted and raw frequencies for the distinct abundance classes. insufficient_evidence is always False for a returned estimate; empty samples raise CoverageInsufficientDataError.

Return type:

dict[str, ndarray | float | bool | str]

chao1(counts, *, ci_level=0.95)

Chao1 nonparametric richness estimator from abundance data (bias-corrected).

S_chao1 = S_obs + f1 (f1 - 1) / (2 (f2 + 1)) (Chao 1984, bias-corrected form), a lower bound on total richness driven by the singleton (f1) and doubleton (f2) counts. Returns a standard error and a log-normal confidence interval for the number of undetected species (Chao 1987), so the interval respects S_chao1 >= S_obs.

Returns:

{'estimate', 'observed', 'f1', 'f2', 'se', 'ci_low', 'ci_high'}.

Parameters:

• counts (ndarray)

• ci_level (float)

Return type:

dict[str, float]

chao2(incidence, *, ci_level=0.95)

Chao2 richness estimator from replicated incidence (presence/absence) data.

Parameters:

• incidence (ndarray) – finite binary (n_species, n_sites) 0/1 matrix (validated by _incidence(); every entry must be exactly 0 or 1). At least _MIN_REPLICATED_SITES sampling units are required – Chao2 estimates unseen richness from replication across units, so one unit does not identify it (CoverageInsufficientDataError).

• ci_level (float) – confidence level for the log-normal interval.

Returns:

{'estimate', 'observed', 'q1', 'q2', 'se', 'ci_low', 'ci_high', 'sites'} where q1/q2 are the numbers of species found in exactly one / two sites.

Return type:

dict[str, float]

ace(counts, *, rare_threshold=10)

ACE: Abundance-based Coverage Estimator of richness (Chao & Lee 1992).

Splits species into abundant (> rare_threshold) and rare, estimates sample coverage from the rare group’s singletons, and corrects for the coefficient of variation of the rare abundances.

Returns:

{'estimate', 'observed', 's_rare', 's_abund', 'c_ace'}.

Raises:

CoverageInsufficientDataError – when the estimated sample coverage c_ace is zero (every rare individual is a singleton) – the 1 / c_ace correction is undefined there and unseen richness is unbounded, so no finite estimate is reported.

Parameters:

• counts (ndarray)

• rare_threshold (int)

Return type:

dict[str, float]

ice(incidence, *, rare_threshold=10)

ICE: Incidence-based Coverage Estimator of richness (the Chao–Lee estimator for incidence data).

Parameters:

• incidence (ndarray) – finite binary (n_species, n_sites) 0/1 matrix (validated by _incidence(); every entry must be exactly 0 or 1). At least _MIN_REPLICATED_SITES sampling units are required, for the same identifiability reason as chao2().

• rare_threshold (int) – species found in <= rare_threshold sites are treated as infrequent.

Returns:

{'estimate', 'observed', 's_infreq', 's_freq', 'c_ice'}.

Raises:

CoverageInsufficientDataError – for a single-unit design, or when the estimated incidence coverage c_ice is zero (every infrequent species seen in exactly one unit) – the correction term s_infreq / c_ice is undefined there.

Return type:

dict[str, float]

hill_numbers(counts, q=(0.0, 1.0, 2.0))

Hill numbers (effective number of species) of order q.

The unified diversity profile: q=0 is observed richness, q=1 is the exponential of Shannon entropy, q=2 is the inverse Simpson concentration. Larger q weights common species more, so the profile D(q) summarises evenness as well as richness.

Parameters:

• counts (ndarray) – per-species abundances.

• q (float | ndarray) – a scalar order or an array of orders.

Returns:

Array of Hill numbers, one per requested order (scalar input still returns a length-1 array).

Return type:

ndarray

rarefaction_curve(counts, sizes=None)

Individual-based rarefaction: expected richness when subsampling m individuals (Hurlbert).

E[S(m)] = sum_i (1 - C(n - x_i, m) / C(n, m)) – the expected number of species seen in a random subsample of m of the n individuals. Used to compare richness between samples at a common sample size (or coverage).

Parameters:

• counts (ndarray) – per-species abundances.

• sizes (ndarray | None) – subsample sizes m to evaluate; defaults to 1 .. n. Must be exact integers in [0, n] (validated by _rarefaction_sizes()); the result is returned in the same order as sizes (not sorted or deduplicated).

Returns:

{'sizes', 'expected_richness'}.

Return type:

dict[str, ndarray]

class CohortAttribution(hazard_ratio, hr_ci, attributable_fraction, af_ci, cif, provenance)

Bases: object

Cox-PH hazard ratio + attributable fraction for one exposure, with intervals and provenance.

Parameters:

• hazard_ratio (float)

• hr_ci (tuple[float, float])

• attributable_fraction (float)

• af_ci (tuple[float, float])

• cif (dict[int, ndarray])

• provenance (dict)

hazard_ratio

exp(coef[exposure_col]) from the fitted Cox model.

Type:

float

hr_ci

Wald (lo, hi) confidence interval on hazard_ratio (log scale, from CoxResult.se).

Type:

tuple[float, float]

attributable_fraction

fraction of hazard among the exposed attributable to exposure, (hazard_ratio - 1) / hazard_ratio.

Type:

float

af_ci

bootstrap (lo, hi) confidence interval on attributable_fraction.

Type:

tuple[float, float]

cif

cause-specific Aalen-Johansen cumulative incidence, {cause: array} (empty unless competing=True).

Type:

dict[int, numpy.ndarray]

provenance

fit diagnostics, the RNG’s full reproducible state (rng_state), and af_distribution – an _AFDistribution (IC-1 DerivedQuantity-shaped) when the bootstrap produced adequate evidence, or an _InsufficientBootstrapEvidence when it did not (see _MIN_BOOTSTRAP_DRAWS).

Type:

dict

The record is frozen and owns its cumulative-incidence curves (read-only copies), so the point estimates, the intervals and the provenance that documents how they were fitted cannot drift apart after the fact: rebinding hazard_ratio without hr_ci – or editing a cif array in place through the caller’s own handle – left a record whose parts described different fits.

Construction now also enforces the record’s own domain (MXR-080-1900). It previously enforced none: a hazard_ratio of -2.0 (impossible – it is exp(coef)), an attributable_fraction of 4.2 (impossible – (HR - 1) / HR approaches but never exceeds 1), and an inverted hr_ci of (9.0, 1.0) all constructed silently and were indistinguishable downstream from a real fit. The invariants, each one the producer already guarantees:

• hazard_ratio finite and strictly positive.

• hr_ci finite, strictly positive, and ordered lo <= hi.

• attributable_fraction finite and <= 1. It is NOT bounded below by 0: a protective exposure has HR < 1 and a legitimately negative attributable fraction.

• af_ci either two finite values ordered lo <= hi, or (nan, nan) – the documented marker for the insufficient-bootstrap-evidence case (see _InsufficientBootstrapEvidence), which is a real state the producer emits. A half-NaN interval is not.

Deliberately NOT checked: that attributable_fraction equals (hazard_ratio - 1) / hazard_ratio exactly, or that either interval brackets its point estimate. Both are true of cohort_attribution()’s own output, but a bootstrap percentile interval is not required to contain the point estimate in general, and pinning the AF identity would turn one documented parameterization into the only constructible one.

provenance is copied DEEPLY: the previous dict(...) detached only the top level, so nested entries – coef/se lists, the rng_state mapping that exists precisely so a fit can be replayed – stayed aliased to the caller’s own objects and could be rewritten after the attribution was reported.

cohort_attribution(covariates, time, event, *, exposure_col=0, competing=False, latency=0.0, n_boot=1000, rng=None)

Attribute a cohort’s hazard to one exposure via Cox PH, with an Aalen-Johansen competing-risks CIF.

Parameters:

• covariates (ndarray) – (n, p) covariates (finite); column exposure_col is the exposure of interest.

• time (ndarray) – (n,) event/censoring times (finite, non-negative).

• event (ndarray) – (n,) event indicator. Binary (0/1) when competing=False; an integer cause label (0 = censored, 1 = the outcome of interest, 2..K = competing causes) when competing=True. Either way, the Cox fit is cause-specific for cause 1: competing-cause events are censored at their time for the hazard-ratio fit (only the CIF, when requested, reports the competing causes’ own cumulative incidence). Values must be exact integers – fractional codes are rejected outright, not truncated (a truncated 1.5 would silently relabel a subject from “the event of interest” to “censored”).

• exposure_col (int) – column of covariates treated as the exposure (must be a valid column index).

• competing (bool) – if True, also run aalen_johansen for cause-specific cumulative incidence.

• latency (float) – exposure-to-effect lag (left-truncation, via cox_ph’s start); 0 disables it. Must be non-negative – negative latency is rejected, not silently treated as 0.

• n_boot (int) – cohort-resample bootstrap draws for the attributable-fraction interval. If fewer than _MIN_BOOTSTRAP_DRAWS converge (or too small a fraction of n_boot, see _MIN_BOOTSTRAP_FRACTION), af_ci is (nan, nan) and provenance["af_distribution"] is an _InsufficientBootstrapEvidence rather than a real interval.

• rng – seed, numpy.random.Generator, or None.

Returns:

A CohortAttribution.

Raises:

ValueError – if the cohort is malformed – see _validate_cohort() for the exact checks (mismatched array shapes, non-finite values, negative time, an out-of-vocabulary or fractional event code, negative latency, or an out-of-bounds exposure_col).

Return type:

CohortAttribution

class Variogram(model, nugget, psill, rng, nu=1.5, anisotropy=None)

Bases: object

A fitted variogram model gamma(h) = nugget + psill * shape(h).

model

"spherical", "exponential", "gaussian" (aka "squared_exponential" / "rbf" – covariance psill * exp(-(h/rng)**2)), or "matern".

Type:

str

nugget

discontinuity at h=0 (measurement error / micro-scale variance).

Type:

float

psill

partial sill (correlated variance); nugget + psill is the total sill.

Type:

float

rng

range (correlation length).

Type:

float

nu

Matern smoothness (ignored by other models).

Type:

float

anisotropy

optional (angle_rad, ratio) geometric anisotropy – coordinates are rotated by angle and the minor axis scaled by 1/ratio before distances are taken.

Type:

tuple[float, float] | None

Raises:

ValueError – if model is not implemented; nugget/psill are not finite and >= 0; rng/nu are not finite and > 0; or anisotropy is set with a non-finite angle or a ratio that is not finite and > 0 (a zero or negative ratio collapses or flips the minor axis, which previously produced NaN predictions and variances downstream instead of an error).

Parameters:

• model (str)

• nugget (float)

• psill (float)

• rng (float)

• nu (float)

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

Frozen. Every kriging solver reads these fields directly and trusts the construction-time validation below without re-checking; while the record was mutable, setting rng = 0, a negative sill, an unimplemented model, or a collapsed anisotropy ratio after construction bypassed the whole check and put the invalid parameter straight into the solve. Freezing keeps the validated state true for the variogram’s lifetime, which is what those solvers assume. anisotropy is normalized to a (angle, ratio) tuple so it cannot be a mutable pair either.

gamma(h)

Evaluate the semivariogram at lag distances.

Parameters:

h (ndarray)

Return type:

ndarray

cov_field(h)

Covariance of the correlated field part (excludes the nugget): psill (1 - shape).

Parameters:

h (ndarray)

Return type:

ndarray

exception VarianceCalibrationError

Bases: ValueError

Base error for a variance-calibration request that has no publishable finite factor.

exception VarianceCalibrationUnidentifiable

Bases: VarianceCalibrationError

Raised when the requested coverage holds without any identifiable positive scale.

exception VarianceCalibrationOutOfRange

Bases: VarianceCalibrationError

Raised when the required scale is outside the finite floating-point domain.

empirical_variogram(coords, values, *, n_bins=15, max_dist=None)

Binned empirical (semi-)variogram: mean 0.5 (z_i - z_j)^2 by separation distance.

max_dist defaults to half the largest pairwise distance. Bins are n_bins equal-width intervals, left-closed and right-open ([edges[b], edges[b+1])) except the last, which is closed on both ends – so a pair separated by exactly max_dist (routine, not an edge case: e.g. the endpoints of an evenly spaced line always sit exactly on the default cutoff) lands in the last bin instead of being classified as “beyond the last bin” and silently dropped.

Returns:

{'lag', 'semivariance', 'count', 'insufficient_evidence', 'reason'} – lag / semivariance / count cover each non-empty distance bin. insufficient_evidence is True (with reason set and the other fields empty) only when no valid distinct pair falls within [0, max_dist] – e.g. a max_dist narrower than every pairwise distance – so no bin is identifiable; it is False (with reason empty) whenever at least one bin is populated. Invalid geometry, fewer than two observations, or no distinct locations raise before distance computation.

Parameters:

• coords (ndarray)

• values (ndarray)

• n_bins (int)

• max_dist (float | None)

Return type:

dict[str, ndarray | bool | str]

fit_variogram(coords, values, *, model='spherical', n_bins=15, nu=1.5)

Fit a variogram model to data by least squares on the empirical variogram.

Returns:

A fitted Variogram (nugget, partial sill, range).

Raises:

ValueError – if the empirical variogram has fewer than _MIN_POPULATED_BINS_FOR_FIT populated bins, so nugget/partial sill/range are not identifiable from the data (see empirical_variogram()).

Parameters:

• coords (ndarray)

• values (ndarray)

• model (str)

• n_bins (int)

• nu (float)

Return type:

Variogram

ordinary_kriging(coords, values, variogram, query, *, noise=None)

Ordinary kriging: BLUP of the field at query under an unknown constant mean.

Parameters:

• coords (ndarray) – (n, d) data locations.

• values (ndarray) – (n,) measured field.

• variogram (Variogram) – a fitted Variogram.

• query (ndarray) – (q, d) prediction locations.

• noise (ndarray | None) – optional (n,) per-observation measurement variance (heteroscedastic nugget); if None the homoscedastic variogram.nugget is used on the diagonal.

Returns:

{'prediction', 'variance'} arrays of length q.

Raises:

ValueError – if coords/values/query are empty, mismatched in shape, or contain non-finite values, or if noise is provided and is mismatched in shape, non-finite, or negative.

Return type:

dict[str, ndarray]

universal_kriging(coords, values, variogram, query, *, degree=1, noise=None)

Universal kriging: kriging with a polynomial spatial trend (drift) of the given degree.

degree=0 uses a constant drift, degree=1 removes a linear trend, and degree=2 a quadratic one. No other degree is implemented. Use when the field has a large-scale drift on top of the stationary residual the variogram describes. The selected drift basis must be full-rank with at least one observation per coefficient.

Returns:

{'prediction', 'variance'}.

Raises:

ValueError – if coords/values/query are empty, mismatched in shape, or contain non-finite values, or if noise is provided and is mismatched in shape, non-finite, or negative.

Parameters:

• coords (ndarray)

• values (ndarray)

• variogram (Variogram)

• query (ndarray)

• degree (int)

• noise (ndarray | None)

Return type:

dict[str, ndarray]

calibrate_variance(predicted_var, residuals, *, target=0.9)

Scale factor that makes kriging predictive intervals hit a target coverage.

Finds the smallest observed order-statistic scale c such that standardised residuals residual / sqrt(c * predicted_var) achieve at least the target empirical central coverage under a Gaussian predictive. Returns the variance multiplier c; multiply predicted_var by it to recalibrate (generic GP/kriging variance recalibration). No fixed search caps are used.

Parameters:

• predicted_var (ndarray) – (m,) held-out kriging variances (must be strictly positive).

• residuals (ndarray) – (m,) held-out actual - predicted.

• target (float) – desired central coverage (e.g. 0.9); must be finite and strictly in (0, 1).

Returns:

The variance multiplier c (> 0).

Raises:

• ValueError – if target is not finite and strictly in (0, 1) – a target <= 0, >= 1, or NaN previously converged to a silent boundary scale factor (1e-6 or 1e6) instead of raising; if predicted_var and residuals do not have the same nonempty one-dimensional shape (a valid paired held-out sample is required); or if either contains non-finite values, or predicted_var is not strictly positive everywhere.

• VarianceCalibrationUnidentifiable – if zero residuals already meet the target for every positive multiplier.

• VarianceCalibrationOutOfRange – if the required factor is outside finite representable range.

Return type:

float

borda_count(rankings)

Borda positional aggregation: each item scores (m - 1 - position) summed over voters.

Returns:

{'consensus', 'scores'} – the consensus ordering (best first) and per-item Borda scores.

Parameters:

rankings (ndarray)

Return type:

dict[str, ndarray]

copeland(rankings)

Copeland pairwise-majority aggregation.

Each ordered pair contributes a win/loss by majority across voters; an item’s score is wins minus losses. Closely tracks the Condorcet winner when one exists.

Returns:

{'consensus', 'scores', 'wins'}.

Parameters:

rankings (ndarray)

Return type:

dict[str, ndarray]

kemeny_consensus(rankings, *, exact_max_items=8)

Kemeny median ranking: minimise the total Kendall-tau distance to all input orderings.

The Kemeny consensus is the maximum-likelihood aggregation under the Mallows–Kendall model and the Condorcet-consistent choice. Exact by enumeration when m <= exact_max_items; otherwise a local search (adjacent transpositions) from the Borda ordering. Exact enumeration evaluates m! permutations, so exact_max_items is an explicit factorial resource boundary; it must be an exact non-Boolean nonnegative integer.

Returns:

{'consensus', 'distance', 'exact', 'search_mode', 'exact_max_items'} – the consensus, total Kendall distance, whether the result is exact, and the requested/achieved search mode.

Parameters:

• rankings (ndarray)

• exact_max_items (int)

Return type:

dict

mallows_fit(rankings, *, exact_max_items=8)

Fit a Mallows model (Kendall): central ranking + dispersion theta.

The central ranking is the kemeny_consensus(); theta is the MLE concentration, found by matching the observed mean Kendall distance to its expectation under the model. Larger theta means tighter agreement (theta -> 0 is uniform/no-consensus).

theta is nan when fewer than two items are ranked (MXR-080-1594): with m < 2 there is exactly one possible ordering, so every ranking has Kendall distance zero under every theta and the likelihood is flat – the dispersion is unidentifiable, not infinite. Reporting inf there (the perfect-agreement branch below) would claim maximal evidence for concentration from data that carries no evidence at all. The trivial center is still returned, since it is correct. inf remains the honest answer for genuine perfect agreement among two or more items, where the observed zero distance really was one of many attainable outcomes.

Returns:

{'center', 'theta', 'mean_distance', 'consensus_distance', 'exact', 'search_mode', 'exact_max_items'}, preserving the Kemeny search guarantee used to obtain the center.

Parameters:

• rankings (ndarray)

• exact_max_items (int)

Return type:

dict

kendall_distance(a, b)

Kendall-tau distance: the number of item pairs ordered oppositely by a and b.

Parameters:

• a (ndarray)

• b (ndarray)

Return type:

int

spearman_footrule(a, b)

Spearman footrule distance: the sum of absolute position differences across items.

Parameters:

• a (ndarray)

• b (ndarray)

Return type:

int

cayley_distance(a, b)

Cayley distance: the minimum number of transpositions turning a into b (m - #cycles).

Parameters:

• a (ndarray)

• b (ndarray)

Return type:

int

class SpatialMixture(shape, n_components, emission, beta=1.0)

Bases: object

A grid-structured mixture with a Potts prior on the latent labels and pluggable mixle emissions.

Parameters:

• shape – grid shape, e.g. (nx, ny) or (nx, ny, nz) – defines the neighbour structure.

• n_components (int) – number of mixture components (latent classes).

• emission – a mixle ParameterEstimator for the per-component family, e.g. MultivariateGaussianEstimator() – this is what makes the class domain-agnostic.

• beta (float) – Potts coupling (>= 0); larger smooths the labels more. 0 is an ordinary mixture.

Raises:

ValueError – if any shape dimension is not an exact positive integer, if the in-memory grid exceeds _MAX_GRID_CELLS, if n_components is not in [1, prod(shape)] (MXR-080-0115: more components than grid cells guarantees at least one permanently empty component by pigeonhole – rejected here rather than left to fail unpredictably during fitting), or if beta is not finite and non-negative.

fit(observations, *, max_iter=40, mf_iter=3, seed=0)

Fit by mean-field variational EM. observations is a length-prod(shape) sequence of per-cell observations (row order matches shape.ravel()); each is a single emission datum.

Robustness: components are initialized by a short hard-assignment pass and the Potts coupling is annealed from 0 to beta over the first iterations, so components form before the smoothness prior is applied (a strong prior on a degenerate init otherwise collapses every cell into one). Every partition is repaired to be nonempty (see _repair_empty_components()) BEFORE it is re-estimated (MXR-080-0115).

Raises:

ValueError – if max_iter or mf_iter is not a positive integer, or if observations does not have exactly prod(shape) entries (MXR-080-0116: a mismatched count would otherwise let responsibilities, neighbors, and encoded data describe different numbers of cells, surfacing later as an unrelated shape-mismatch crash deep inside the M-step instead of a clear error here).

Parameters:

• max_iter (int)

• mf_iter (int)

• seed (int)

Return type:

SpatialMixture

responsibilities()

The posterior label probabilities, (prod(shape), n_components) – a simplex per cell.

Return type:

ndarray

labels()

The MAP label field, reshaped to shape.

Return type:

ndarray

entropy()

Per-cell posterior entropy (label uncertainty), reshaped to shape.

Return type:

ndarray

component(j)

The fitted mixle emission distribution of component j.

Parameters:

j (int)

Return type:

Any

class SmithMaxStable(sigma)

Bases: object

The Smith max-stable process with Gaussian storm-profile covariance sigma (d x d, SPD).

A spatial process, not an i.i.d. leaf distribution: its full likelihood is intractable, so it is not a SequenceEncodableProbabilityDistribution – it exposes the things that do have closed forms (extremal_coefficient, bivariate_cdf) plus a sampler, and is fitted by the module-level fit_smith_maxstable() (composite/madogram estimation), mirroring the functional fit style of the other non-leaf spatial models. Margins are unit Frechet; spatial dependence grows with sigma.

Parameters:

sigma (np.ndarray)

property sigma: ndarray

Owned read-only covariance used consistently by dependence formulas and simulation.

extremal_coefficient(h)

theta(h) = 2 * Phi(a/2) with a the Mahalanobis lag length – 1 at h=0 (full dependence) rising to 2 as the lag grows (independence).

Parameters:

h (ndarray)

Return type:

float

bivariate_cdf(z1, z2, h)

P(Z(s) <= z1, Z(s+h) <= z2) = exp(-V(z1, z2)) – the Smith bivariate distribution.

z1/z2 are thresholds on unit-Frechet margins and so must be strictly positive (the unit-Frechet support is (0, inf)); the closed form below divides by z1/z2 and takes log(z2/z1), which silently produces nan/inf rather than a real probability at or below 0.

Parameters:

• z1 (float)

• z2 (float)

• h (ndarray)

Return type:

float

sampler(locations, seed=None)

Return a sampler over the requested spatial locations.

Parameters:

• locations (ndarray)

• seed (int | None)

Return type:

SmithMaxStableSampler

class SmithMaxStableFit(model, n_locations, n_replicates, n_pairs, residual, status='ok')

Bases: object

A fitted isotropic Smith max-stable process, plus a diagnostic of whether the pairwise dependence-scale match actually identified an interior optimum.

status is one of:

• "ok": the fitted scale is an interior point of the search bracket – a genuine stationary point of the F-madogram least-squares objective.

• "boundary": the fit landed at (or numerically indistinguishable from) the edge of the search bracket – the signature of an objective with no interior minimum. This happens, for example, when replicates carry almost no independent information (every location’s rank collapses together), which pushes the empirical extremal coefficient to the same value for every pair and drives the least-squares match monotonically toward one edge instead of a genuine optimum. model is still populated with the boundary value, but it is a search artifact, not a dependence-scale estimate, and should not be trusted.

Parameters:

• model (SmithMaxStable)

• n_locations (int)

• n_replicates (int)

• n_pairs (int)

• residual (float)

• status (str)

model

the fitted SmithMaxStable.

Type:

mixle.analysis.max_stable.SmithMaxStable

n_locations

distinct spatial locations used.

Type:

int

n_replicates

replicated fields used.

Type:

int

n_pairs

location pairs the least-squares match was computed over (n_locations choose 2).

Type:

int

residual

the F-madogram objective’s value at the fitted scale (mean squared error between the empirical and model extremal coefficients across pairs; exactly 0 for a 2-location fit, which has one pair and one parameter).

Type:

float

status

see above.

Type:

str

Construction validates the record against its own documented domain (MXR-080-1900). It previously validated nothing at all, so SmithMaxStableFit(model=None, n_locations=-1, n_replicates=0, n_pairs=-5, residual=nan) constructed cleanly and reported converged is True – a record claiming a converged fit that had no model, no data and no objective value. The invariants are exactly what fit_smith_maxstable() produces: a real SmithMaxStable; at least two locations (a dependence scale needs a pair) and one replicate; n_pairs the n_locations choose 2 the docstring above states; a finite, non-negative residual (it is a mean of squares); and a status this class actually defines.

Deliberately NOT changed: the record stays a MUTABLE dataclass. Its frozen siblings in this package were converted under a separate finding, and freezing this one is a broader API change than the shape invariants this validates – the fields can still be rebound after construction and escape the checks below.

property converged: bool

True only when the fit is a genuine interior optimum, not a search-bound artifact.

class SmithMaxStableSampler(dist, locations, seed=None)

Bases: object

Sampler for a fitted Smith max-stable process at fixed locations.

sample() truncates the (in principle unbounded, infinite-storm) Smith spectral representation to a finite spatial box and a tolerance-controlled number of storms; see its docstring for exactly what that truncation does and does not rigorously control, and approximation_diagnostic() for an empirical check of the result.

Parameters:

• dist (SmithMaxStable)

• locations (np.ndarray)

• seed (int | None)

sample(size=None, *, n_storms=20_000, tol=1e-3, box_sigma=5.0)

Draw max-stable field(s) at the locations via a boxed, tolerance-truncated storm sum.

The true Smith representation is a max over an infinite Poisson process of storms (xi_i, U_i) spread over the whole of R^d; only then are margins exactly unit Frechet. This sampler approximates that in two distinct ways, one rigorously controlled and one not:

1. Storm-count truncation (rigorously controlled by ``tol``). Storms are drawn in decreasing order of their Poisson intensity mark xi (xi_i = 1/Gamma_i, the natural arrival order of this construction). Because the Gaussian storm kernel is bounded by its peak density phi_Sigma(0), a storm with mark xi can contribute at most xi * vol * phi_Sigma(0) to any site, no matter where its center U_i lands (vol the sampling box’s volume). Since xi is strictly decreasing, once that bound drops below tol for one storm it has also dropped below tol for every later one – so it is exact, not heuristic, to stop there: every dropped storm’s true contribution is provably below tol at every site, which bounds the total storm-count truncation error (relative to the boxed field’s own true value) by tol. n_storms is a safety cap on iterations, not the primary control; a UserWarning is raised if that cap binds before tol is reached for some replicate, since then the tol guarantee was not actually achieved. Every storm drawn is used before its own tol check is applied, so at least one storm always contributes to the field regardless of how loose tol is – an all-zero field (like the previous n_storms=0 behavior) is just as invalid a “draw” as a negative one, since 0 is outside unit Frechet’s (0, inf) support.

2. Spatial truncation (NOT rigorously controlled). Storm centers are drawn only from a finite box around the query locations, expanded by box_sigma storm-kernel standard deviations in each dimension, rather than the true unbounded domain. This introduces a bias that empirically shrinks as box_sigma grows, but is not exactly bounded here – use approximation_diagnostic() to check it for your particular locations/sigma.

Net effect: simulated margins are only approximately, not exactly, unit Frechet. The storm-count half of that approximation is controlled to within tol (when n_storms doesn’t bind first); the box half is not formally bounded and should be checked empirically.

Parameters:

• size (int | None) – number of independent replicate fields (None for a single field).

• n_storms (int) – safety cap on storms drawn per replicate; must be a positive integer. Storms stop earlier, as soon as tol is satisfied. Must be a positive integer – with n_storms=0 no storms are drawn at all, which cannot produce a valid unit-Frechet draw (0 is outside its (0, inf) support), so that case is rejected outright rather than silently returned as if it were one.

• tol (float) – stop once every remaining storm’s max-possible contribution to every site is below this (absolute, field-scale) tolerance. Must be finite and strictly positive.

• box_sigma (float) – half-width, in storm-kernel standard deviations, of the box storm centers are drawn from beyond the query locations’ own extent. Must be finite and strictly positive.

Returns:

Field(s) at the locations. Shape (len(locations),) if size is None, else (size, len(locations)).

Return type:

ndarray

size (when supplied) and n_storms are exact non-Boolean positive integers; zero does not represent a valid request for a field or a storm process.

approximation_diagnostic(n=2000, probability_levels=(0.25, 0.5, 0.75, 0.9), **sample_kwargs)

Empirically estimate how close this sampler’s truncated simulation is to true unit-Frechet margins.

Draws n replicate fields (with sample_kwargs forwarded to sample(), so a specific n_storms/tol/box_sigma setting can be checked) and compares their empirical CDF, at each site and each requested probability level, against the true unit-Frechet CDF – see FrechetApproximationDiagnostic. n is an exact non-Boolean positive integer and probability_levels is a nonempty finite one-dimensional sequence strictly inside (0, 1).

Parameters:

• n (int)

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

Return type:

FrechetApproximationDiagnostic

class FrechetApproximationDiagnostic(n_replicates, probability_levels, max_abs_error, mean_abs_error, per_site_max_abs_error)

Bases: object

Empirical check of how close a truncated SmithMaxStableSampler draw is to true unit Frechet margins.

At each site and each requested probability level p, compares the empirical CDF of simulated draws (evaluated at the true unit-Frechet quantile for p) against p itself – a QQ/PP-plot-style discrepancy. Both summary statistics should shrink toward 0 as the sampler’s truncation is loosened (smaller tol, larger box_sigma/n_storms) and should NOT shrink just from raising n alone (n only reduces Monte Carlo noise in the estimate of a possibly nonzero bias).

Parameters:

• n_replicates (int)

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

• max_abs_error (float)

• mean_abs_error (float)

• per_site_max_abs_error (ndarray)

n_replicates

number of simulated draws the diagnostic was computed from.

Type:

int

probability_levels

the probability levels checked.

Type:

tuple[float, …]

max_abs_error

worst-case |empirical_cdf(q_p) - p| over every (site, level) pair.

Type:

float

mean_abs_error

mean |empirical_cdf(q_p) - p| over every (site, level) pair.

Type:

float

per_site_max_abs_error

worst-case error at each site, shape (n_locations,).

Type:

numpy.ndarray

fit_smith_maxstable(locations, fields)

Fit an isotropic Smith max-stable process (sigma = s^2 I) to replicated spatial extremes.

locations is (n_locations, d) and fields is (n_replicates, n_locations) of block maxima. Estimation matches the binned empirical extremal coefficient (from the F-madogram) to the model 2 Phi(|h| / (2 s)).

This is a single-parameter (s) fit estimated from pairwise statistics, so it is only identifiable given enough structure to compute those statistics from:

• at least 2 locations, so there is at least one pair (with at least one nonzero lag between some pair) to compute an empirical extremal coefficient from at all;

• at least 2 replicates, so the per-location rank transform is not degenerate – with exactly 1 replicate every location’s rank trivially collapses to a single point, making the empirical extremal coefficient identically 1 (full dependence) for every pair regardless of distance, which is not a real dependence signal.

Returns a SmithMaxStableFit (the fitted model plus fit-quality diagnostics); raises ValueError if the inputs are malformed or too small to attempt a fit at all.

Parameters:

• locations (ndarray)

• fields (ndarray)

Return type:

SmithMaxStableFit

class LedoitWolfEstimator(dim=None, name=None, keys=None)

Bases: ParameterEstimator

Estimate a multivariate Gaussian with a Ledoit-Wolf-shrunk covariance.

Returns a MultivariateGaussianDistribution whose mean is the sample mean and whose covariance is shrunk toward a scaled identity by the data-driven Ledoit-Wolf intensity. The chosen intensity is exposed on the returned distribution as dist.shrinkage for inspection. When the accumulated sufficient statistics carry too little effective weight or no identifiable dispersion, raises LedoitWolfInsufficientData.

Parameters:

• dim (int | None)

• name (str | None)

• keys (str | None)

accumulator_factory()

Return an accumulator factory for streaming Ledoit-Wolf sufficient statistics.

Return type:

LedoitWolfAccumulatorFactory

estimate(nobs, suff_stat)

Estimate a Gaussian distribution from accumulated Ledoit-Wolf sufficient statistics.

Raises:

• LedoitWolfInsufficientData – If the accumulated weight is below _MIN_EFFECTIVE_COUNT or the observations contain no numerically identifiable dispersion.

• ValueError – If the shrinkage computation itself produces a non-finite mean or covariance despite adequate, already-validated sufficient statistics (a numerical-overflow defect, not an insufficient-data condition). A substantively non-positive-semidefinite result is rejected by MultivariateGaussianDistribution’s own construction, which this relies on as the final PSD gate.

Return type:

MultivariateGaussianDistribution

exception LedoitWolfInsufficientData(reason, effective_count, dim)

Bases: ValueError

Specific estimation error raised when the data cannot identify covariance.

effective_count and dim let a caller decide whether to accumulate more data and retry. Raising, rather than returning a foreign non-distribution object, preserves the common ParameterEstimator result protocol.

Parameters:

• reason (str)

• effective_count (float)

• dim (int | None)

cost_curve(depth, grade, throughput, *, params)

Parametric mining + processing cost in $/t, as a function of depth, grade, and throughput.

depth, grade, and throughput are broadcastable array-likes (one entry per block or per scheduling period; scalars broadcast against the others). params recognizes (all optional, defaulting to zero/no-effect):

• base_cost: $/t floor cost.

• haul_cost_per_m: $/t per metre of depth – haulage and pumping/dewatering cost, modeled as linear in depth, so the curve is strictly increasing in depth whenever this is positive.

• grade_complexity_coef: $/t scale of a 1 / grade metallurgical-complexity penalty – lower-grade ore needs proportionally more material handled and processed per unit of recovered metal, so this term falls as grade rises.

• throughput_scale_coef / design_capacity: the plant has one throughput, design_capacity (Q*), at which fixed costs are spread most efficiently; cost rises quadratically away from it in either direction – throughput_scale_coef * ((Q - Q*) / Q*) ** 2 – capturing both under-utilized fixed-cost drag below Q* and overtime/expediting/accelerated-wear cost above it (the classic “decreasing then rising past design capacity” U-shaped average-cost curve).

depth must be finite and non-negative (metres below surface/datum; there is no such thing as negative haul depth); grade and throughput must be finite and strictly positive (both are also divisors elsewhere in the mine-planning stack, and a zero or negative ore grade or throughput is not a physical plan). Every numeric params entry must be finite, and the cost-coefficient entries (base_cost, haul_cost_per_m, grade_complexity_coef, throughput_scale_coef) must be non-negative – a negative coefficient would let one additive term of the curve go negative and silently discount an otherwise non-negative total cost into a subsidy (MXR-080-0119). Raises ValueError if any of these are violated.

Returns the elementwise $/t cost, broadcast to the common shape of the three inputs.

Parameters:

• depth (Any)

• grade (Any)

• throughput (Any)

• params (dict)

Return type:

ndarray

capex_opex(plan, *, params)

Roll a mine plan’s tonnage/depth/grade/throughput profile up into (total capex, total opex).

plan is a mapping or attribute-bearing object exposing, per scheduling period:

• tonnage: array-like, tonnes mined/processed each period (required).

• depth, grade, throughput: array-likes (or scalars, broadcast against tonnage) fed to cost_curve() to get each period’s $/t.

• capex_schedule (optional): array-like of lumpy capital spend per period (e.g. pre-strip, plant construction, fleet purchases); summed into total capex on top of the params below.

params is passed through to cost_curve() for the opex side, plus two capex-only keys: capex_fixed (one-off, tonnage-independent capital) and capex_per_tonne (sustaining capital that scales with total tonnage mined over the plan).

Total opex is sum(tonnage * cost_curve(depth, grade, throughput, params=params)); total capex is capex_fixed + capex_per_tonne * sum(tonnage) + sum(capex_schedule). Returns (capex, opex), both plain floats – the totals monte_carlo_npv() (J2) discounts into a DCF, and the same $/t curve this function calls is what feeds block_cost for H’s optimizers.

tonnage must be finite and non-negative and have at least one period; depth/grade/ throughput must each be a scalar or share tonnage’s exact shape (see _require_scalar_or_period_shape()); capex_schedule, if given, must be finite and non-negative. A negative tonnage or a negative cost/capital term would otherwise turn modeled spending into modeled revenue with no error (MXR-080-0119); raises ValueError if violated.

Parameters:

• plan (Any)

• params (dict)

Return type:

tuple[float, float]

class NPVDistribution(samples, mean, p10, p50, p90, sensitivity)

Bases: _NPVDistributionFields

A Monte-Carlo DCF outcome: the full NPV sample, not just a point estimate.

samples is the length-n array of per-draw NPVs (real dollars, one entry per Monte-Carlo trial) that value_at_risk() / conditional_value_at_risk() (J5) consume directly. mean/p10/p50/p90 are the usual project-finance summary of that distribution. sensitivity decomposes the NPV variance into the grade- and price-uncertainty contributions, as Sobol first-order variance-based sensitivity indices (see monte_carlo_npv() and _sobol_first_order_share()): keyed "grade" / "price" (each factor’s fraction of total variance explained on its own, mathematically in [0, 1] and not required to sum to 1 – NPV is multiplicative in grade and price, so some variance is a grade/price interaction effect attributed to neither) plus the raw "grade_variance" / "price_variance" / "total_variance".

Construction validates samples: non-empty and finite (no NaN/Inf) – defense-in-depth so invalid state can never flow downstream to a caller even if some upstream computation (here, monte_carlo_npv()) fails to validate its own inputs, the same guard applied to carcinogenic_risk.RiskQuantity and the _SampleDerivedQuantity carriers elsewhere in mixle.analysis.

Parameters:

• samples (ndarray)

• mean (float)

• p10 (float)

• p50 (float)

• p90 (float)

• sensitivity (Mapping[str, float])

Return type:

NPVDistribution

monte_carlo_npv(posterior, price_paths, cost_model, schedule, *, discount_rate, n=10000, rng)

Monte-Carlo discounted-cash-flow NPV distribution (DR-ALG J2).

Draws n grade realizations off the IC-1 posterior and pairs them, draw for draw, with n (resampled/aligned as needed) price scenarios from price_paths; for each draw and each of the schedule’s periods, cashflow_t = tonnage_t * grade_t * price_t - opex_t - capex_t, discounted at discount_rate (period 0 undiscounted, period t divided by (1 + discount_rate) ** t) and summed into one NPV per draw:

• posterior: an IC-1 Posterior (frozen mixle.reason.posterior_protocol.Posterior); its .samples(n, rng) are the grade draws. A single project-life grade (d == 1) broadcasts across every period; a per-period posterior (d == len(schedule)) is used period by period.

• price_paths: a scenario-major (m, n_periods) array-like (one row per scenario), OR a J1 PriceForecast.paths passed directly – (n_periods, m), time-major – detected and transposed automatically (see _align_price_paths()). Resampled with replacement to n rows when m != n.

• cost_model: declares exactly one supported positional signature, either cost_model(t, tonnage_t) or cost_model(t), to get that period’s deterministic opex_t. Signature binding selects the form before execution; an interface accepting both is rejected as ambiguous rather than guessed.

• schedule: per-period tonnage (required) and capex (optional, default zero); see _unpack_schedule(). len(schedule)’s tonnage vector fixes n_periods.

• discount_rate: the DCF discount rate per period. Must be finite and strictly greater than -1 – discount_rate <= -1 makes the per-period discount factor 1 / (1 + discount_rate) ** t divide by zero (at exactly -1) or alternate sign every other period (below -1), neither of which is a coherent discount rate.

• n / rng: Monte-Carlo draw count and the shared numpy.random.Generator. n must be a positive int (not merely a positive-valued float); posterior.samples(n, rng) is required to return exactly n draws, checked explicitly rather than trusted – a posterior that returns fewer (e.g. a single row) used to silently numpy-broadcast across every price scenario, fabricating n “independent” draws that all shared one real grade sample (MXR-080-0118).

Every numeric input – the grade draws, price_paths, schedule’s tonnage/capex, and cost_model’s returned opex – must be finite and non-negative; raises ValueError if any of these, n, or discount_rate are violated, before any DCF arithmetic runs.

Returns an NPVDistribution with the raw samples, mean/p10/p50/p90, and a sensitivity dict decomposing NPV variance into grade vs. price contributions as Sobol first-order variance-based sensitivity indices (Saltelli/Jansen two-independent-sample estimator; see _sobol_first_order_share()) – each factor’s own share of the total NPV variance, estimated from a SECOND joint sample independent of the one samples/mean/p10/p50/ p90 are built from, so the indices stay meaningful (population-bounded in [0, 1]) even when a caller’s price_paths happen to be paired against the posterior’s draws in a way that makes the actual joint sample’s own variance small (MXR-080-0117). Not a full ANOVA decomposition (grade and price’s interaction/second-order effect, if any, is not itself reported), but each reported share is individually a real, correctly-estimated first-order Sobol index.

Parameters:

• posterior (Posterior)

• price_paths (Any)

• cost_model (Callable)

• schedule (Any)

• discount_rate (float)

• n (int)

• rng (Generator)

Return type:

NPVDistribution

class EmissionFactors(scope1, scope2, scope3, sigma=None)

Bases: object

Per-activity-key CO2e emission factors, one dict per GHG-Protocol scope.

scope1/scope2/scope3 map an activity key (e.g. "diesel_L", "grid_kWh", "explosives_kg", "transport_t_km") to a CO2e factor expressed per unit of that activity (e.g. kg CO2e per litre of diesel). A key absent from a scope’s dict simply does not contribute to that scope. sigma optionally gives the standard deviation of each factor (not the activity quantity itself), keyed the same way across all three scopes, for Monte-Carlo uncertainty propagation in emissions_footprint(); a key with no entry in sigma is treated as exactly known (std 0).

Parameters:

• scope1 (dict[str, float])

• scope2 (dict[str, float])

• scope3 (dict[str, float])

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

class Footprint(scope1, scope2, scope3, total, ci=None, provenance=<factory>)

Bases: object

A Scope 1/2/3 CO2e footprint with an optional 90% credible interval and full provenance.

scope1/scope2/scope3/total are in the same physical CO2e units as the emission factors (typically kg CO2e / tCO2e). ci is the (lo, hi) 90% Monte-Carlo interval on total when factor uncertainties were propagated, else None (the default – a caller that only needs climate_terms()’ point costs can construct one without ci/provenance). provenance identifies both inputs to the activity x factors product, plus how the interval was produced: factor_source; the 64-hex activity_content_hash and factor_content_hash (sha256 of the canonical activity / factor-model encodings, the same hashing convention IC-2 uses for field artifacts); factor_keys (the union of activity keys the factor model prices, for a readable diff when the two hashes disagree); the scopes actually included in total; and draws / uncertainty_propagated describing the Monte-Carlo configuration behind ci. Hashing the activity alone was not enough to distinguish two runs: the same schedule costed against a supplier-specific factor set and against a national average is a different footprint with different provenance.

The record is frozen: a reported footprint is evidence, and rebinding total or scope2 after construction would leave the numbers, the interval and the provenance describing three different computations. ci is normalized to a (lo, hi) tuple and provenance to a dict this record owns, so a later edit to the caller’s own dict cannot retroactively rewrite what a footprint claims it was computed from.

That copy is DEEP (MXR-080-1900). dict(provenance) detaches only the top level, so a nested container – provenance["scopes"] as a list, a nested factor-model dict – stayed aliased to the caller’s own object and could be rewritten after the footprint was reported. Deliberately NOT checked here: that total equals scope1 + scope2 + scope3. The producer always makes it so, but the record is documented as constructible directly for callers that only need climate_terms()’ point costs, and the test-suite already builds partially-specified footprints; a sum invariant would refuse states the library legitimately accepts.

Parameters:

• scope1 (float)

• scope2 (float)

• scope3 (float)

• total (float)

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

• provenance (dict)

emissions_footprint(activity, factors, *, scopes=(1, 2, 3), n=0, rng=None)

Scope 1/2/3 CO2e footprint of a production activity schedule.

Each activity key (e.g. diesel_L, grid_kWh, explosives_kg, transport_t_km) is priced by the corresponding factor in whichever of factors.scope1/scope2/scope3 includes that key: Scope 1 is direct combustion/blasting, Scope 2 is purchased electricity, Scope 3 is upstream reagents plus downstream transport. scopes selects which of the three scopes are actually included in the returned footprint (a scope not requested reports 0.0 and does not contribute to total – e.g. scopes=(1, 2) for a Scope-1/2-only disclosure).

If n > 0 and factors.sigma is supplied, each priced factor is treated as Normal(mean=factor, std=sigma.get(key, 0.0)) and resampled n times (factors with no sigma entry stay fixed); the resulting distribution of total yields a 90% credible interval in ci. Without both a positive n and a non-empty sigma, ci is None – the point total is still returned, just without an uncertainty band.

provenance always fingerprints both factors of the activity x factors product – activity_content_hash AND factor_content_hash (64-hex sha256 over the canonical encodings, IC-2’s hashing convention) – alongside factor_keys, scopes, draws and uncertainty_propagated. A downstream carbon-cost or transition-risk term (L3/L6) can therefore be traced back to the exact activity schedule and the exact factor table it was priced with; recording only the activity meant a footprint restated against a revised or supplier-specific factor set was indistinguishable from the original by its provenance alone.

scopes must not contain duplicate entries (each GHG-Protocol scope may appear at most once) – a duplicate raises rather than being silently summed once on the point-total path and once per occurrence on the Monte-Carlo path, which previously left the point estimate and the credible interval describing two different footprints. activity values must be finite and nonnegative; factors (scope1/scope2/scope3) values must be finite; factors.sigma values, when supplied, must be finite and nonnegative – a negative sigma is rejected rather than silently treated as an exactly-known (zero-uncertainty) factor. n must be an exact non-Boolean nonnegative integer; zero explicitly requests the point estimate without an interval. All of this is validated before either the point-total or the Monte-Carlo path runs.

Parameters:

• activity (dict[str, float])

• factors (EmissionFactors)

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

• n (int)

• rng (Generator | None)

Return type:

Footprint

class TransitionRiskResult(samples, prior_dominated, scenario_mean, ranking, carbon_cost, provenance)

Bases: object

The carbon-adjusted NPV distribution across carbon-price/policy scenarios (L3).

Satisfies the frozen mixle.reason.posterior_protocol.DerivedQuantity structural protocol – samples, prior_dominated, credible_interval – so a carbon-adjusted value can flow anywhere a DerivedQuantity is expected (J5 tail risk, J2 re-valuation). samples is shaped (n, k): the n baseline npv_samples draws, each re-priced under every one of the k carbon-price scenarios (one column per scenario) – the re-ranking below stays uncertainty-aware rather than collapsing straight to a point estimate. prior_dominated is always False: there is no prior/regulariser here, the distribution’s width is set entirely by npv_samples.

Beyond the protocol, scenario_mean (per-scenario mean carbon-adjusted NPV), ranking (scenario indices sorted best -> worst by scenario_mean), and carbon_cost (the priced-and-discounted carbon cost subtracted from each scenario) carry the scenario-level comparison transition_risk() exists to produce.

Construction validates samples: non-empty and finite (no NaN/Inf) – defense-in-depth so invalid state can never flow downstream to a caller even if some upstream pushforward (here, transition_risk()) fails to validate its own inputs, the same guard applied to carcinogenic_risk.RiskQuantity and the _SampleDerivedQuantity carriers elsewhere in mixle.analysis.

The record is frozen and owns its arrays (validated copies, write-locked), so ranking and scenario_mean cannot drift out of agreement with the samples they were derived from: an unfrozen record let a caller reorder the ranking or overwrite a scenario mean while every summary method kept reporting off the original draws.

Freezing only stopped the drift AFTER construction, though. Construction itself accepted any finite scenario_mean and any permutation ranking, so a record could be born already contradicting its own draws – scenario_mean reporting one ordering of the scenarios while credible_interval(), computed off samples, reported another (MXR-080-1900). Two agreement invariants now close that, both holding exactly for what transition_risk() produces:

• scenario_mean must equal the per-scenario (column) mean of samples, to a tolerance scaled by the magnitude of the draws – summation order, not disagreement, is all that is allowed to separate them.

• scenario_mean[ranking] must be non-increasing, i.e. the ranking really does run best -> worst. Ties may be broken in any order, since np.argsort gives no stable tie order and demanding a particular one would refuse rankings the producer itself emits.

Additionally, when provenance carries the module’s own "carbon_cost" key (the doubled representation transition_risk() writes), it must agree with the carbon_cost field – the two used to be free to disagree about the very cost the summaries were priced from.

Deliberately NOT checked: that samples[:, j] + carbon_cost[j] recovers one common baseline NPV across scenarios. It is true of transition_risk()’s output, but the record is also constructed directly with independent samples/carbon_cost in this package’s own tests, so requiring it would reject states the library legitimately produces.

Parameters:

• samples (ndarray)

• prior_dominated (bool)

• scenario_mean (ndarray)

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

• carbon_cost (ndarray)

• provenance (dict)

credible_interval(level)

Per-scenario central level interval of the carbon-adjusted NPV, each shape (k,).

level goes through the shared mixle.analysis interval contract (MXR-080-1580).

Parameters:

level (float)

Return type:

tuple[ndarray, ndarray]

transition_risk(footprint, carbon_price_paths, *, npv_samples, discount=None)

Carbon-adjusted NPV distribution + scenario ranking under a set of carbon-price paths (L3).

For each of the k scenarios in carbon_price_paths (a (k,) flat price or a (k, t) per-period path – see _coerce_price_paths()), the priced carbon cost is footprint.total * sum_t(price[t] * discount[t]) (discount defaults to all-ones, i.e. no discounting, when omitted – pass period discount factors, e.g. 1 / (1 + r) ** t, to match a J2 monte_carlo_npv DCF). That per-scenario carbon cost is subtracted from every draw of the baseline npv_samples (a J2 NPVDistribution.samples-shaped (n,) array), yielding an (n, k) carbon-adjusted value distribution: one re-priced NPV distribution per scenario, still carrying the original valuation uncertainty.

Scenarios are ranked by mean carbon-adjusted NPV (scenario_mean, descending: best scenario first) so a high-carbon-price/policy scenario reliably re-ranks below a low-price one, with the gap between any two scenarios scaling linearly in footprint.total – a bigger footprint pays proportionally more carbon cost under the same price paths. The returned TransitionRiskResult satisfies IC-1’s DerivedQuantity protocol so the re-ranking is always inspectable with a credible interval, not just a point estimate; it feeds J5 tail risk and J2 re-valuation directly.

Validated before any ranking happens: carbon_price_paths must have at least one scenario and at least one period, with every price finite and nonnegative; discount, when supplied, must be finite and nonnegative; footprint.total must be finite; and npv_samples must be a non-empty, finite, one-dimensional (n,) array (a multi-dimensional array is rejected, not silently flattened). Letting any of these through used to allow NaN scenario means and a meaningless ranking to construct successfully, failing only later when a caller requested a credible interval.

Parameters:

• footprint (Footprint)

• carbon_price_paths (ndarray)

• npv_samples (ndarray)

• discount (ndarray | None)

Return type:

DerivedQuantity

class WaterBudget(*args, **kwargs)

Bases: Protocol

Structural stand-in for L2’s WaterBudget (mixle-pde) – see the module’s repo-boundary note.

Never constructed here; a real WaterBudget (or anything shaped like one) is passed in and duck-typed via getattr(water, ..., default) throughout this module, so a caller-supplied object satisfying only shortfall_m3 (the one required field) still works – the optional fields default to absent, matching the “or None if absent” fallbacks already documented on _water_demand_m3() and _shortfall_time_fraction().

climate_terms(footprint, water, *, carbon_price, water_limit_m3=None)

Fold an emissions footprint and a water budget into J6’s carbon cost + H4’s water constraint.

Returns {"carbon_cost": float, "water_feasible": bool | None, "water_binding": bool | None, "shortfall_time_fraction": float | None}:

• carbon_cost = footprint.total * carbon_price – the priced carbon term J6’s objective (analysis/valuation.py) subtracts alongside capex_opex’s cost roll-up. Always a definite float; it does not depend on water evidence.

• water_binding is set when the L2 WaterBudget.shortfall_m3 > 0 (the routed water balance already ran dry at some point over the plan) OR when a caller-supplied water_limit_m3 is given and the budget’s demand exceeds it; water_feasible is its negation – a hard constraint H4 (stochastic_opt.py) folds into a per-option cost (derating or dropping the option) rather than into the mean-grade objective. Either can come back None (unknown) instead of a definite bool – see below.

• shortfall_time_fraction is the fraction of the water budget’s own per-step storage trajectory sitting at a binding zero (_shortfall_time_fraction()) – a deterministic duration/occupancy statistic over ONE trajectory, reported independently of water_limit_m3 even for options that are feasible on average but fragile. This is NOT a probability: it carries no uncertainty information. A genuine cross-realization probability requires an explicit ensemble of trajectories (_shortfall_probability_over_ensemble()), which is not available from a single water object.

Missing or non-finite water evidence is an explicit UNKNOWN state, never a permissive pass:

• When water is None (no water budget available for this option – e.g. an inland option with no catchment exposure), water_feasible, water_binding, and shortfall_time_fraction are all None: climate risk for that option is carbon-only, but the water terms are reported as “not evaluated”, not silently “known feasible”.

• When shortfall_m3, the water-limit demand, or water_limit_m3 itself is non-finite (NaN/inf) or negative, that check contributes None (unknown) rather than silently comparing and passing as non-binding. All three are physical volumes in m3 and cannot be negative under the documented contract, so a negative one is malformed evidence rather than evidence of comfortable headroom (or, for a negative demand/limit, of a violation). The two binding checks combine via three-valued (Kleene) OR: a confirmed binding signal from either check always wins, an unknown never gets rounded down to “feasible”.

• When a water object is supplied but carries no readable shortfall – the attribute is absent, None, a Boolean, or non-numeric – the shortfall check likewise contributes None. This is the case the promise above used to skip (MXR-080-1900): the water is None branch abstained correctly, but a present-but-empty water object took getattr(..., shortfall_m3, 0.0)’s numeric default and reported water_feasible=True, a definite all-clear on a hard constraint from an object holding no water evidence at all. shortfall_time_fraction reports None in the same case rather than a confident 0.0 (“never short over the whole trajectory”) – unless a real storage trajectory is present, which is genuine evidence and is still used.

Parameters:

• footprint (Footprint)

• water (WaterBudget | None)

• carbon_price (float)

• water_limit_m3 (float | None)

Return type:

dict

resistance_raster(habitat, *, floor=1e-3)

Per-cell movement cost from a fitted suitability field: cost_c = 1 / max(suitability_c, floor).

Reads only habitat.mean (N1’s fitted intensity field lambda_c); floor keeps near-zero suitability cells at a large-but-finite cost (near-impermeable) rather than blowing up to inf, so the resulting raster is always safe to feed straight into least_cost_corridor() / habitat_connectivity() without any further cleanup.

Raises:

ValueError – if floor is not finite and strictly positive, or if habitat.mean contains a NaN/Inf/negative entry. An earlier revision silently clamped negative suitability into the same near-impermeable cost as zero suitability (MXR-080-0073), hiding invalid upstream model output; a negative suitability now always raises instead of being clamped. Zero suitability remains a legitimate input (clamped to floor, as documented above).

Parameters:

• habitat (HabitatModel)

• floor (float)

Return type:

np.ndarray

least_cost_corridor(resistance, patch_a, patch_b)

Cheapest movement path between two cells of a resistance raster (the least-cost corridor).

resistance is an n-dimensional cost raster (e.g. from resistance_raster()); patch_a/ patch_b are flat (row-major/np.ravel) cell indices into it. Delegates to mixle.relations.ShortestPath over the rook-adjacency grid graph with each edge weighted by _edge_cost() (mean resistance of its two endpoint cells); a cell resistance of inf (a mined footprint, fragmentation_impact()) makes every edge touching it non-finite and so unusable, effectively removing that cell from the graph. Lower resistance = better connected.

ShortestPath’s underlying engine (mixle.relations.best_first_paths) is a pure best-first search with no closed/visited set – intentional there, since it is also used to enumerate k best (not just the single best) paths, where revisiting a state via a costlier path is legitimate output. On a cyclic grid graph asked for only its single best path (k=1), that same lack of a closed set is combinatorially explosive (every back-and-forth detour is a distinct, never-pruned path). Since successors(state) is called exactly once per pop and pops come out in non-decreasing cost order, the first call for a given state is provably its cheapest arrival (Dijkstra’s invariant, valid here because every edge cost is non-negative) – so a small closed set local to this call, never touching relations.py, turns the search back into ordinary Dijkstra: each state’s successors are computed once and [] (already settled) on every later call.

Returns:

the total path resistance and the visited flat cell indices from patch_a to patch_b inclusive; (inf, []) if no finite-cost path exists.

Return type:

(corridor_resistance, path)

Parameters:

• resistance (ndarray)

• patch_a (int)

• patch_b (int)

effective_conductance(resistance, sources, sinks)

Circuit-theory effective conductance between sources and sinks over a resistance raster.

Solves the standard Dirichlet/Laplacian boundary-value problem for landscape connectivity (McRae et al. 2008; the “Circuitscape” formulation): every sources cell is held at potential 1 and every sinks cell at potential 0 – equivalent to shorting each terminal’s cells together with an ideal conductor, with no arbitrary coupling constant needed (unlike a max-flow super-source/-sink). The remaining free cells’ potentials solve the reduced linear system L_free @ v_free = -L_free,fixed @ v_fixed for the graph Laplacian L (edge weight = conductance = 1 / edge_cost, _edge_cost()); the total current flowing out of the source terminal at this unit potential difference is the effective conductance (I = G @ V, V = 1), read off directly from the Laplacian: current out of node i is (L @ v)[i], which is exactly 0 at every free node (that is what the linear solve enforces) and the true injected/withdrawn current at every fixed node.

This is a genuinely different quantity from maximum flow (max_flow_connectivity()): conductances in series add reciprocally (two unit conductances in series give an effective conductance of 0.5, not 1), so a max-flow solver – which only sees bottleneck capacity and ignores series resistance entirely – systematically overstates connectivity whenever current has to cross more than one edge to reach the sink. On a unit-conductance square between opposite corners this function returns the textbook value 1; a max-flow solver over the same graph returns 2 (the sum of the two parallel corner-to-corner paths’ raw capacities, MXR-080-0071).

Cells with no path at all to either terminal (an isolated island cut off by +inf-resistance neighbours on every side) carry no current and are excluded from the linear solve rather than left to make it singular; a raster whose sources have no path to its sinks at all correctly yields 0.0 (the linear system still solves – each side simply settles to its own terminal’s potential everywhere reachable, so no potential gradient, hence no current, exists anywhere).

Raises:

ValueError – on a NaN/negative/zero resistance entry (+inf, the impassable-cell sentinel, is fine), an empty raster, an out-of-bounds source/sink index, or overlapping source/sink sets.

Parameters:

• resistance (ndarray)

• sources (Sequence[int])

• sinks (Sequence[int])

Return type:

float

habitat_connectivity(resistance, sources, sinks)

Domain-friendly alias for effective_conductance(), the circuit-theory connectivity metric.

Kept under this name because it is the pre-existing public entry point N4/N6/J callers and tests already use; see effective_conductance() for the full Dirichlet/Laplacian derivation and its McRae et al. (2008) citation. Prior to MXR-080-0071 this function actually computed maximum flow (now max_flow_connectivity()) despite the circuit-theory claim in its docstring – a different, systematically-inflated quantity for any network with more than one edge between source and sink. It now computes the real thing.

Parameters:

• resistance (ndarray)

• sources (Sequence[int])

• sinks (Sequence[int])

Return type:

float

max_flow_connectivity(resistance, sources, sinks)

Maximum-flow connectivity between sources and sinks over a resistance raster’s conductance graph.

Builds a conductance graph (arc capacity 1 / edge_cost between rook-adjacent cells) with a super-source over sources and a super-sink over sinks, then returns mixle.relations.max_flow()’s value between them – the maximum flow (bottleneck-limited throughput) the habitat network can carry from the source patch(es) to the sink patch(es).

This is NOT the circuit-theory effective-conductance quantity (effective_conductance() / habitat_connectivity()): max flow only sees the tightest cut and entirely ignores series resistance, so it systematically overstates connectivity whenever current must cross more than one edge to reach the sink (MXR-080-0071). This function exists as a distinctly-named, independently useful bottleneck-throughput metric in its own right – never use it as a stand-in for the other.

Parameters:

• resistance (ndarray)

• sources (Sequence[int])

• sinks (Sequence[int])

Return type:

float

fragmentation_impact(resistance, footprint_mask, sources, sinks)

Habitat-connectivity impact of a development footprint: baseline vs. disturbed corridor/connectivity.

Sets every footprint_mask cell’s resistance to inf (impassable / zero conductance – an H3 mine plan, in this module’s worked instantiation) and recomputes both least_cost_corridor() (between the first sources/sinks cell, taken as the representative patch-to-patch corridor endpoints – a judgment call: the public API only takes cell sets here, not a single designated pair) and habitat_connectivity() (over the full sources/sinks sets) before and after. This is the population-viability proxy N6/J feed into the biodiversity-offset objective: a footprint that severs the only corridor raises corridor_resistance and drops connectivity sharply; one that misses every real corridor leaves both essentially unchanged.

Returns:

A dict with corridor_resistance_baseline/corridor_resistance_mined, connectivity_baseline/connectivity_mined, delta (connectivity_baseline - connectivity_mined, i.e. connectivity lost), and mincut_edges – the baseline network’s minimum-cut arcs (mixle.relations.min_cut()), i.e. the single weakest link a footprint would need to sever to maximally damage connectivity.

Parameters:

• resistance (ndarray)

• footprint_mask (ndarray)

• sources (Sequence[int])

• sinks (Sequence[int])

Return type:

dict

habitat_offset_liability(plan_footprint, habitat, *, offset_ratio, unit_offset_cost)

Priced biodiversity-offset liability of disturbing plan_footprint (a J6 priced-objective term).

lost_equivalents = sum_{c in footprint} suitability_c * area_c (suitability = N1’s fitted HabitatModel.mean); the liability is offset_ratio * lost_equivalents * unit_offset_cost – an additive dollar term the same shape J6’s priced_liabilities already sums for carbon/health/ remediation (workstream-J.md J6). offset_ratio=0 or unit_offset_cost=0 reduces this to zero, i.e. no biodiversity-offset requirement – both remain valid, meaningful inputs and are not rejected.

Because lost_equivalents is linear in the boolean footprint, the per-cell rate offset_ratio * unit_offset_cost * suitability_c * area_c is itself a valid per-block deduction a MILP-based optimizer (H4/J6’s risk_adjusted_plan) can net directly out of expected per-block profit – this function is the scalar evaluator for a given (candidate or solved) footprint.

Raises:

ValueError – if offset_ratio/unit_offset_cost is negative or non-finite, or if _lost_equivalents rejects plan_footprint/habitat. A negative ratio or cost used to silently turn habitat damage into a “profit” (MXR-080-0073); both are now rejected. Also if the priced product itself is not representable in float64 – finite inputs near the float64 ceiling used to publish an infinite liability (MXR-080-1572).

Parameters:

• plan_footprint (np.ndarray)

• habitat (HabitatModel)

• offset_ratio (float)

• unit_offset_cost (float)

Return type:

float

no_net_loss_constraint(plan_footprint, habitat, *, offset_ratio)

Hard no-net-loss constraint payload: created/purchased offsets >= offset_ratio * lost_equivalents.

Returns a dict carrying both the raw quantities (lost_equivalents, per_cell_lost_equivalents, required_offset) and a solver-agnostic linear-constraint row in this repo’s standard coeffs @ x <= bound convention (mixle.relations/mixle.stochastic_opt’s a_ub rows), expressed over a single offsets_created decision variable: coeffs=[-1.0], bound=-required_offset encodes -offsets_created <= -required_offset, i.e. offsets_created >= required_offset. Placing that row (and the offsets_created column it references) into the wider extraction/offset-purchase decision space is H4/J6’s job – this module never edits their MILP variable indexing, only hands them the row.

Because that row is emitted ALREADY normalized into the <= convention, its "sense" label is "<=" – the label describes the row that is actually here, not the inequality the row expresses about offsets_created. It used to read ">=" (MXR-080-1570), which made the payload internally contradictory: mixle.analysis.objective.hard_constraints() and mixle.stochastic_opt.risk_adjusted_plan() both trust the label and negate a ">=" row into <= form, so a second negation turned the no-net-loss FLOOR into the CEILING offsets_created <= required_offset – permitting zero offsets and forbidding any offset purchase above the bare minimum, i.e. exactly the opposite of a no-net-loss guarantee.

Raises:

ValueError – if offset_ratio is negative or non-finite, or if _lost_equivalents rejects plan_footprint/habitat (NaN/negative suitability or area). A NaN or negative lost_equivalents used to be able to propagate into required_offset, silently weakening or inverting a constraint that is supposed to be a hard floor (MXR-080-0073); both inputs are now validated before this payload is built. Also if required_offset itself is not representable in float64 – finite inputs near the float64 ceiling used to publish an infinite bound, which is not a hard constraint but an unsatisfiable row (MXR-080-1572).

Parameters:

• plan_footprint (np.ndarray)

• habitat (HabitatModel)

• offset_ratio (float)

Return type:

dict

class SpeciesObservation(species_id, detection, location, crs=None, covariates=<factory>, modality='occurrence', provenance=<factory>)

Bases: object

One presence/absence record for a species (an IC-4 Observation specialisation).

location is the (already discretised) study-area coordinate: for a presence record consumed by fit_sdm(), the first component of location is the fractional cell index ([0, K)) on the same grid as the covariates/cell_area passed to fit_sdm – resolving a real-world crs-referenced location onto that grid is covariate/CRS ingest (B-series), out of scope here.

Parameters:

• species_id (str)

• detection (bool)

• location (ndarray)

• crs (str | None)

• covariates (dict[str, Any])

• modality (str)

• provenance (dict[str, Any])

class HabitatModel(beta, beta_cov, design, cell_area, *, species_id, var_scale=1.0, prior_dominated=False, provenance=None)

Bases: object

Fitted IPP habitat-suitability field; satisfies IC-1 Posterior over the suitability field.

beta/beta_cov are the Laplace-approximate posterior over the log-linear intensity coefficients; design is the (K, p) covariate design matrix (intercept + environmental covariates); cell_area is the per-cell area used as the Poisson offset during fitting. The suitability field is the fitted intensity lambda_c = exp(design_c @ beta) (mean); samples()/cov/credible_interval() push the beta-posterior forward through the same log-link (a delta-method / lognormal approximation), scaled by a held-out-calibrated variance multiplier, so every downstream consumer (N2’s no-mine mask, N4’s resistance raster) sees one calibrated field posterior.

species_id is required (MXR-080-1590): a suitability field is a statement about ONE species, and a model with no identity on it cannot be checked against the species whose legal exclusion it is being used to impose (MXR-080-1591) or traced back to the dataset it was fitted from. provenance carries immutable input receipts (fit_sdm() fills in the occurrence, background, covariate and cell-area digests) so a published field can be tied to its source data.

Parameters:

• beta (np.ndarray)

• beta_cov (np.ndarray)

• design (np.ndarray)

• cell_area (np.ndarray)

• species_id (str)

• var_scale (float)

• prior_dominated (bool)

• provenance (dict | None)

samples(n, rng)

Draw n intensity-field realisations by pushing beta draws through exp(X @ beta).

Returns:

(n, K) array of intensity-field draws.

Raises:

ValueError – n is not a positive exact integer (MXR-080-0114).

Parameters:

• n (int)

• rng (Generator)

Return type:

ndarray

property mean: ndarray

Fitted intensity field lambda_c = exp(design_c @ beta) – the suitability surface.

Uses the same overflow-safe exp as fitting (MXR-080-0114): a linear predictor that never overflowed under the training-time clip can still overflow plain exp() at a different design row (e.g. an extrapolated covariate) or for a directly-constructed model.

property cov: ndarray | LinearOperator

Delta-method covariance of the intensity field, recalibrated.

Domains with at most 2,048 cells return a dense (K, K) array. Larger domains return a symmetric LinearOperator backed by the (K, p) intensity Jacobian and (p, p) coefficient covariance, avoiding quadratic field storage.

Raises:

ValueError – the Jacobian or requested covariance operation is not finite. In particular, finite model inputs can still imply a covariance outside float64’s representable range; such a matrix is not published as posterior evidence.

credible_interval(level)

Per-cell central credible interval of the suitability field (lognormal delta-method).

Raises:

ValueError – level is not finite and strictly in (0, 1) (MXR-080-0114).

Parameters:

level (float)

Return type:

tuple[ndarray, ndarray]

derived_quantity(fn, n, rng)

Pushforward fn over n intensity-field draws into a DerivedQuantity (IC-1).

Parameters:

• fn (Callable[[ndarray], ndarray])

• n (int)

• rng (Generator)

Return type:

_PushforwardQuantity

critical_habitat_mask(threshold)

Boolean mask, True where fitted suitability lambda_c >= threshold.

The hard no-mine constraint N2 hands to H (same shape/role as a G9 seepage polygon).

Parameters:

threshold (float)

Return type:

ndarray

fit_sdm(occurrences, covariates, cell_area, *, species_id, background=None, ridge=1e-3)

Fit a presence-only habitat-suitability model as an inhomogeneous-Poisson point process.

Discretizes the study area into K = covariates.shape[0] cells, bins occurrences (only detection=True records) into per-cell presence counts, and fits a log-linear intensity log lambda_c = [1, covariates_c] @ beta by maximizing the frozen IPP log-likelihood with a ridge * ||beta||^2 penalty and a log(area_c) offset (Poisson GLM equivalent to MaxEnt). If background quadrature points are given, they are binned the same way and folded into the offset as extra effective area (an effort/thinning correction for presence-only sampling bias). The returned HabitatModel’s field variance is recalibrated on an internal held-out cell fold via mixle.analysis.kriging.calibrate_variance() so its credible intervals hit their nominal (90%) coverage.

species_id names the ONE species being modelled and is required (MXR-080-1590). Every occurrences record must carry that same id: the fit used to bin every truthy detection regardless of which species it belonged to, so a single fitted intensity could silently pool records from unrelated species into one field that claimed to describe all of them. The id is carried onto the returned model, together with immutable input digests in provenance, so a published scientific result can be traced back to the exact dataset it came from and so a downstream legal exclusion can check the field it is applying actually belongs to the species it is applying it for (MXR-080-1591).

Parameters:

• occurrences (list[SpeciesObservation]) – presence records for species_id; each location’s first component is the fractional cell index. A record for any other species is rejected, not silently pooled.

• covariates (ndarray) – (K, p - 1) environmental covariates per cell (an intercept column is prepended).

• cell_area (ndarray) – (K,) per-cell area (the Poisson offset).

• species_id (str) – the species this field is being fitted for.

• background (ndarray | None) – optional quadrature/background point locations (same cell-index convention as occurrences), used to correct for uneven survey effort.

• ridge (float) – L2 penalty strength on beta (also regularizes the Laplace covariance).

Returns:

A fitted HabitatModel carrying species_id and input receipts in provenance.

Raises:

ValueError – ridge is not finite and non-negative; covariates is not finite; cell_area does not have exactly one finite, strictly positive entry per covariate row; any presence or background location is non-finite or outside the cell-index domain [0, K) (MXR-080-0111); species_id is not a non-empty string or any occurrence record names a different species (MXR-080-1590); or the internal beta fit does not converge to finite coefficients (MXR-080-0113).

Return type:

HabitatModel

critical_habitat_exclusion(habitat, listed, *, suitability_cut, buffer_cells=0)

Boolean no-mine mask over blocks (work-plan algorithm steps 1-2).

For every listed record whose critical_habitat flag is set, THAT SPECIES’ OWN fitted suitability field contributes its lambda_c >= suitability_cut cells (via IC-12’s critical_habitat_mask) to the exclusion; a record with critical_habitat=False (a species that is tracked but has no critical-habitat designation) contributes nothing. If no listed record qualifies, nothing is excluded. The per-species contributions are unioned (OR).

habitat is either a single HabitatModel or an explicit {species_id: HabitatModel} mapping, and every qualifying listed species must have a model whose species_id matches it exactly (MXR-080-1591). Previously any qualifying record applied whatever single field happened to be passed in, and no comparison was made at all: a model fitted for species A could impose species B’s legal exclusion, and several listed species all collapsed onto species A’s same field. A qualifying species with no matching model raises rather than borrowing another species’ range.

When a matched fit is prior-dominated (habitat’s own honesty flag: not enough presence data has yet informed the field), there isn’t enough evidence to clear any block, so the mask is excluded conservatively in full (every block is treated as potential critical habitat) rather than optimistically reporting whatever the under-determined field happens to say.

buffer_cells conservatively dilates the resulting mask outward (_dilate_conservatively()) to approximate a regulatory buffer around a designation boundary.

A real CriticalHabitatDesignation polygon (a regulator’s own mapped boundary, independent of the fitted suitability field) is not rasterized here: that requires resolving its SpatialBounds (real crs coordinates) onto this function’s abstract block-index grid, which is exactly the covariate/CRS-ingest machinery N1 scoped out (B-series, not yet landed) – see the PR notes for this documented gap. Once a grid-registration utility lands, folding a designation in is an additional OR into mask before the buffer dilation, using the same conservative-inclusion rule.

Parameters:

• habitat (HabitatModel | dict[str, HabitatModel]) – N1’s fitted habitat-suitability field(s) (IC-12 HabitatModel, satisfies IC-1 Posterior) – one model, or a {species_id: HabitatModel} mapping.

• listed (Sequence[ListedSpecies]) – ListedSpecies records (mixle_knowledge.contracts); critical_habitat, species_id (to match the record to its own fitted field) and, when informativeness must be checked, the matched model’s own honesty flag are consulted here – every other field (citation, jurisdiction, listing status) is provenance the caller and downstream audit trail carry, not logic this function branches on.

• suitability_cut (float) – the fitted-intensity threshold passed to critical_habitat_mask.

• buffer_cells (int) – conservative dilation radius, in blocks, applied after the union.

Returns:

A (K,) boolean array, True where the block is excluded (no-mine).

Return type:

np.ndarray

apply_habitat_constraints(network, exclusion_mask)

Fold a critical-habitat exclusion mask into an IC-9-shaped network payload (work-plan algorithm step 3); H1/H4 read this, this module never calls the solver.

network is a plain mapping over the reference block network, in exactly mixle.relations.min_cost_flow()’s frozen (cap, cost, supply) shape: "cap"/"cost" are (n, n) arc matrices, "supply" is the optional length-n node supply vector, and "block_nodes" (optional, defaults to arange(len(exclusion_mask))) maps each block index to its node index in that network – the same convention G9’s apply_env_constraints uses.

Excluded blocks become forbidden nodes: every arc touching one has its capacity zeroed (no flow can originate from, or land back on, a no-mine block), and, if that node carried supply, the supply is zeroed too (there is nothing left to extract there). Any nodes/arcs/fixed_costs/ demands the caller supplied (the shape mixle.relations.network_design() itself takes) are passed through unchanged alongside forbidden_nodes, so a fixed-charge caller can exclude the same nodes from its own arc set.

Every input is validated before any exclusion is applied, so a rejected call never leaves a partially-mutated result (MXR-080-0092): cap/cost must be equal-shape, square, and finite, and cap must additionally be nonnegative (a negative arc capacity has no meaning); supply, when supplied, must be a finite length-n vector; and block_nodes must be exact, unique node ids in [0, n) (see _validate_block_nodes() for the fractional/negative/out-of-range failure modes this replaces).

Parameters:

• network (dict[str, Any]) – the reference network payload (cap/cost/supply/block_nodes, plus any network_design-shaped pass-through fields).

• exclusion_mask (ndarray) – the (K,) boolean no-mine mask from critical_habitat_exclusion().

Returns:

cap/cost (habitat-adjusted), forbidden_nodes (sorted node-id list), and supply when the caller provided one, plus any pass-through fields.

Return type:

A new dict

Raises:

ValueError – if cap/cost are not equal-shape finite square arrays, if cap contains a negative entry, if supply is not a finite length-n vector, or if block_nodes is not a length-matched array of exact, unique, in-range node ids.

priced_liabilities(plan, *, carbon_price, health_cost, remediation_cost)

Net a plan’s per-block environmental/health/carbon liabilities into one additive dollar array.

plan is a mapping or attribute-bearing object exposing three per-block, length-n_blocks array-likes:

• grade: the block grade (or other geology/complexity proxy) fed to remediation_cost – a G9-style environmental remediation-cost model keyed on what is actually mined.

• exposure: the block’s public-health exposure proxy fed to health_cost – a K6-style health-cost model.

• emissions: the block’s carbon-equivalent emissions, priced directly at carbon_price (an L6-style market/regulatory carbon price is a multiplier, not a fitted cost model, so no callable is needed for it).

Returns a dict with per-block arrays "remediation", "health", "carbon", their elementwise sum "total" (the single per-block liability mixle.stochastic_opt.risk_adjusted_plan() nets out of revenue), and the scalar roll-ups "remediation_total"/"health_total"/ "carbon_total"/"grand_total" for reporting.

Every plan field is validated BEFORE it is handed to a cost callable (MXR-080-1573). The callables are opaque user code: a model that ignores its argument, saturates, or returns a constant turns a NaN grade or an infinite exposure into a perfectly ordinary finite cost, so the invalid source evidence disappears and the resulting liability record carries no trace that the plan it was priced from was malformed. grade/exposure are checked for finiteness only, not sign – both are declared as arbitrary geology/complexity and exposure proxies, so a negative value is a legitimate reading; it is only the prices that carry the non-negativity convention.

Every assembled output is likewise validated after the fact: three individually finite per-block costs near the float64 ceiling can overflow total, and finite per-block totals can overflow grand_total’s reduction, publishing an infinite liability from a record that claims to be validated (MXR-080-1573).

Raises:

ValueError – if carbon_price, plan’s emissions, or either cost callable’s return value is non-finite or negative – a negative liability would net out of revenue as a profit increase rather than a cost, and a NaN would propagate unguarded into the objective – if plan’s grade/exposure is non-finite, if any of plan’s grade/exposure/emissions or a cost callable’s return value has the wrong shape, or if the netted per-block total or any scalar roll-up is not representable in float64.

Parameters:

• plan (Any)

• carbon_price (float)

• health_cost (Callable[[ndarray], ndarray])

• remediation_cost (Callable[[ndarray], ndarray])

Return type:

dict

hard_constraints(*, no_mine_mask=None, caps=None)

Assemble the constraints dict mixle.stochastic_opt.risk_adjusted_plan() consumes.

• no_mine_mask: an optional boolean array-like, length n_blocks. True entries are a G9-style no-mine/buffer-zone polygon’s enclosed blocks – hard-fixed to x_b = 0 (stronger than an inequality row: those blocks can never be selected, regardless of the rest of the objective).

• caps: an optional list of linear activity caps – a K6 exposure budget, an L6 water budget, or any other coeffs @ x <sense> bound row over the block decision vector. Each entry is a dict with "coeffs" (length-n_blocks array-like) and "bound" (float), plus an optional "sense" ("<=" by default; ">=" is negated into the solver’s standard <= form so risk_adjusted_plan never has to special-case it). This is the same coeffs/bound/ sense field naming mixle.analysis.biodiversity.no_net_loss_constraint() already uses for its own linear-constraint payload.

Returns a dict with whichever of "no_mine_mask"/"caps" were supplied (both are omitted, i.e. an empty dict, when neither argument is given – meaning no hard constraint at all).

Raises:

ValueError – if no_mine_mask is not a genuine 1-D boolean array (data that is merely truthy under an implicit bool cast, e.g. a float array or a NaN entry, is rejected rather than silently coerced), if any cap’s "coeffs" is not a finite 1-D array, if any cap’s "bound" is not finite, or if a cap’s "sense" is not "<="/">=".

Parameters:

• no_mine_mask (Any | None)

• caps (list[dict] | None)

Return type:

dict

class OptionValue(value, exercise_boundary, premium_over_npv)

Bases: NamedTuple

The priced option, alongside the exercise policy that generated it.

value is the total value of holding the option (project value + optionality). exercise_boundary is the per-time-step critical underlying value at which immediate exercise first becomes optimal (nan at any step where it is never optimal to exercise). premium_over_npv is value minus the naive npv_dist.mean – the dollar amount the flexibility to defer/expand/abandon is worth over just committing to the point-estimate NPV today.

Parameters:

• value (float)

• exercise_boundary (ndarray)

• premium_over_npv (float)

value: float

Alias for field number 0

exercise_boundary: ndarray

Alias for field number 1

premium_over_npv: float

Alias for field number 2

real_option_value(npv_dist, *, volatility, horizon, kind='defer', rate, n_steps=None, expand_fraction=0.3, expansion_cost=0.0, salvage_value=0.0)

Price the defer/expand/abandon option on a project via an additive binomial lattice on NPV.

The NPV process is modeled as a driftless arithmetic random walk started at npv_dist.mean: each of n_steps (default: max(horizon, 1)) steps up to horizon moves the value by +-h with h = volatility * scale * sqrt(dt), scale = abs(npv_dist.mean) (or 1.0 if the mean is exactly zero, so volatility is never degenerate) – i.e. volatility is a fractional per-sqrt-period dispersion relative to the project’s own scale, matching how the parameter is usually quoted (e.g. 0.3 ~ “30% dispersion”), while still letting the walk go negative the way a signed NPV must be able to. The risk-neutral measure is taken driftless (an NPV is already the discounted expectation of a self-financing project, so there is no further risk-neutral drift to add); rate only discounts the continuation value, which is what makes immediate exercise optimal once volatility vanishes. At every node the holder may exercise (payoff via _intrinsic(), depending on kind) or continue holding the option; the American value is the larger of the two, by backward induction from the horizon.

Parameters:

• npv_dist (NPVDistribution) – anything exposing .mean (a float) – typically J2’s NPVDistribution. Only the mean is used; real_option_value prices the option on top of the point estimate, not a re-derivation of the distribution itself. Must be finite.

• volatility (float) – fractional (per-sqrt-period) dispersion of NPV around its mean. Must be finite and >= 0; 0 means no dispersion and the option collapses to immediate exercise of kind’s own payoff on npv_dist.mean (_intrinsic()) – max(npv_dist.mean, 0) for "defer", max(npv_dist.mean, salvage_value) for "abandon", npv_dist.mean + max(expand_fraction * max(npv_dist.mean, 0) - expansion_cost, 0) for "expand".

• horizon (int) – number of periods over which the option may be exercised. 0 means “decide now”. Must be a non-Boolean integral scalar >= 0.

• kind (str) – one of "defer", "expand", "abandon".

• rate (float) – per-period discount rate applied to the continuation value. Must be finite.

• n_steps (int | None) – lattice steps (defaults to max(horizon, 1), i.e. one step per period). Must be a non-Boolean positive integral scalar when given explicitly.

• expand_fraction (float) – for kind="expand", the fractional capacity bonus applied to a positive underlying value when the expansion is exercised. Must be finite and nonnegative.

• expansion_cost (float) – for kind="expand", the investment cost of actually exercising the expansion, netted against expand_fraction’s capacity bonus – without it, a positive-NPV project always “expands” for free and there is no real exercise decision being priced (see _intrinsic()). Ignored for other kinds. Must be finite and nonnegative; default 0.0 (no cost, matching this function’s behavior before expansion_cost existed).

• salvage_value (float) – for kind="abandon", the value recovered by abandoning the project instead of riding out its NPV, replacing the previous hardcoded assumption of a total write-off. Ignored for other kinds. Must be finite; default 0.0 (total write-off, matching this function’s behavior before salvage_value existed).

Returns:

An OptionValue.

Return type:

OptionValue

class GaussianObservationModel(obs_matrix, obs_cov)

Bases: NamedTuple

A declared linear-Gaussian observation model: y = obs_matrix @ theta + noise, noise ~ N(0, obs_cov), for a latent theta distributed as the current posterior’s (mean, cov).

MXR-080-0109: the opt-in heuristic invents future-posterior centers/spreads by rescaling draws from the CURRENT posterior; it never samples an observation from a likelihood/forward model and conditions on it, so it is not a real Bayesian preposterior for non-Gaussian, multimodal, or constrained beliefs. This is the simplest well-defined experiment likelihood for which the preposterior – the distribution of what TOMORROW’s posterior could look like, before the observation is actually taken – has a closed form: conjugate-Gaussian updating. Passing this to voi_dollars() / voi_estimate() (as observation_model=) is what makes them compute a genuine expected value of sample information (EVSI): simulate a hypothetical y from theta’s actual prior-predictive marginal, analytically condition the CURRENT posterior on that y to get a real future posterior, decide against it, and average over many simulated y – rather than the opt-in heuristic’s rescaling of today’s draws (see voi_estimate()’s docstring for why that heuristic is not a real preposterior in general).

obs_matrix (H, shape (k, d)) and obs_cov (R, shape (k, k), symmetric positive definite) together say: a real experiment would return a k-dimensional noisy linear readout of the d-dimensional latent quantity the posterior is over (d = posterior.mean’s length). This path requires posterior.mean / posterior.cov to fully characterize the belief – i.e. a genuinely Gaussian posterior – since that is what makes the conjugate update exact; a dense posterior.cov is required (a matrix-free LinearOperator covariance is not supported here). obs_matrix = 0 (or obs_cov large) declares a literally uninformative experiment – VOI comes back exactly zero for that case, not just approximately (see voi_estimate()’s docstring).

Parameters:

• obs_matrix (ndarray)

• obs_cov (ndarray)

obs_matrix: ndarray

Alias for field number 0

obs_cov: ndarray

Alias for field number 1

class VoiEstimate(value, standard_error, ci_low, ci_high, method)

Bases: NamedTuple

A Monte Carlo value-of-information estimate, reported honestly alongside its own sampling noise (MXR-080-0108: independent no-info/with-info draws made a should-be-zero difference noisy, and flooring that noise at zero upward-biased it; see voi_estimate()’s docstring for the fix).

value is the point estimate, E[value | info] - E[value | no info] – NOT floored at zero (see voi_estimate()’s docstring for why flooring a noisy estimate at zero systematically biases it upward, rather than merely rounding it). standard_error is the paired-difference Monte Carlo standard error of value across the n_outer hypothetical-information replicates (paired because each replicate’s with-info and no-info evaluations share their underlying random draw – see voi_estimate()). ci_low/ci_high are an approximate two-sided 95% confidence interval (value +- 1.96 * standard_error; a normal approximation, not exact for small n_outer). When ci_low <= 0 <= ci_high, value is statistically indistinguishable from “no additional value” – a more honest read than a bare point estimate that happens to be (noisily) positive or negative.

method names which computation produced the estimate: "variance_rescaling_heuristic" (an explicit opt-in Gaussian-dispersion approximation, exact only when the posterior is close to Gaussian, unimodal, and unconstrained) or "gaussian_conjugate_evsi" (a genuine sample-then-condition preposterior computation under a declared linear-Gaussian observation_model).

Parameters:

• value (float)

• standard_error (float)

• ci_low (float)

• ci_high (float)

• method (str)

value: float

Alias for field number 0

standard_error: float

Alias for field number 1

ci_low: float

Alias for field number 2

ci_high: float

Alias for field number 3

method: str

Alias for field number 4

voi_estimate(posterior, decision_fn, drill_info, *, rng, n_outer=64, n_inner=256, observation_model=None)

Value of information, honestly: the point estimate together with its own Monte Carlo uncertainty.

decision_fn maps a set of posterior draws (an (n, d) array, physical units) to the dollar value of the single best decision made using that belief state – e.g. a risk-neutral go/no-go choice is max(samples.mean(), 0), not an average of a per-draw payoff: the latter implicitly assumes the realization is already known, which is exactly the perfect-information case this function is pricing the gap to, not the belief itself.

Two methods are available:

• observation_model=None together with explicit drill_info["method"] == "variance_rescaling_heuristic": a pre-posterior Monte Carlo built on the law of total variance, NOT real Bayesian conditioning on a simulated observation. Today’s posterior variance splits into “where the post-drill posterior will be centered” (unknown until the drillhole is actually put in) and “how much spread remains once it lands there.” If the drillhole is expected to remove a fraction r of variance (_variance_reduction()), then for n_outer hypothetical drill outcomes a center is drawn with standard deviation scaled by sqrt(r) around today’s posterior mean (how much the belief could plausibly shift), and a shared base draw (see “Common random numbers” below) is re-centered there and rescaled by sqrt(1 - r) for the remaining spread. This INVENTS future-belief “centers” by rescaling today’s draws rather than drawing them from an actual forward/likelihood model, so it is accurate only when the posterior is (close to) Gaussian, unimodal, and unconstrained – outside that regime it is flagged with a UserWarning (see _warn_if_not_gaussian_like()), not silently trusted.

• observation_model= a GaussianObservationModel (method="gaussian_conjugate_evsi"): a genuine sample-then-condition preposterior. A hypothetical observation y is simulated from the latent’s real prior-predictive marginal (y = H @ theta + noise for theta ~ N(posterior.mean, posterior.cov)), the CURRENT posterior is analytically conditioned on that y via closed-form conjugate-Gaussian updating to get a genuine future posterior, and the decision is re-made against it – averaged over many simulated y for E[value | info]. This is the textbook EVSI construction for the declared-likelihood case, but it still requires the PRIOR (posterior.mean / posterior.cov) to be genuinely Gaussian for the conjugate update to be exact – also screened by the same warning.

Common random numbers (MXR-080-0108): each of the n_outer replicates draws ONE base/innovation sample and uses it for BOTH the no-info decision and the with-info (refined) decision for that replicate, differing only by a deterministic transform (the heuristic’s re-centering/rescaling, or the conjugate update’s affine map) – not independent draws for the two sides. This matters when there is truly no information to gain (r = 0 for the heuristic; obs_matrix = 0 or obs_cov arbitrarily large for the declared-likelihood path): with independent draws, the no-info and with-info values would be two separately noisy estimates of the same quantity, so their difference would be noise centered at zero; flooring that noise at zero (the previous behavior) kept only its positive half, upward-biasing the reported VOI (a standard-normal fake posterior with a go/no-go decision spuriously returned positive VOI at zero information for most tested seeds). With a shared draw, the with-info value is bit-for-bit identical to the no-info value whenever there is truly no information (checked directly for the heuristic via an analytic short-circuit below, not just left to fall out of the Monte Carlo – and provably exact by construction for the declared-likelihood path too), and the two stay correlated – hence lower-variance in their difference – away from that boundary as well.

Parameters:

• posterior (Posterior) – the current (pre-drill) belief, satisfying IC-1’s Posterior protocol.

• decision_fn (Callable[[ndarray], float]) – belief state ((n, d) draws) -> expected dollar value of the best decision.

• drill_info (dict[str, Any]) – describes the hypothetical drillhole. Recognized keys: method (must explicitly be "variance_rescaling_heuristic" when no observation_model is supplied), variance_reduction (direct fractional variance reduction, default 0.5), or candidate_geometry + forward_op (+ optional region / cell_volumes) to route through the C8 VOI hook when available; n_outer_samples / n_inner_samples override the Monte Carlo sample counts. Ignored (aside from the sample-count overrides) when observation_model is given.

• rng (Generator) – seeded random generator for reproducibility.

• n_outer (int) – number of hypothetical-information replicates (overridden by drill_info).

• n_inner (int) – posterior draws per replicate (overridden by drill_info).

• observation_model (GaussianObservationModel | None) – when given, a declared GaussianObservationModel routes the computation through the genuine "gaussian_conjugate_evsi" path instead of the opt-in heuristic; see above. Requires a dense (non-LinearOperator) posterior.cov.

Returns:

the (unfloored) point estimate plus its Monte Carlo uncertainty.

Return type:

A VoiEstimate

voi_dollars(posterior, decision_fn, drill_info, *, rng, n_outer=64, n_inner=256, observation_model=None)

The value-of-information point estimate, in dollars: E[value | info] - E[value | no info].

A thin convenience wrapper around voi_estimate() for callers that only want the point estimate; see its docstring for the full method (the opt-in heuristic vs. the genuine observation_model EVSI path, common random numbers, the zero-information analytic equality, and what method the estimate used) and for how to also get the Monte Carlo standard error / CI.

Parameters:

• posterior (Posterior) – forwarded to voi_estimate() unchanged.

• decision_fn (Callable[[ndarray], float]) – forwarded to voi_estimate() unchanged.

• drill_info (dict[str, Any]) – forwarded to voi_estimate() unchanged.

• rng (Generator) – forwarded to voi_estimate() unchanged.

• n_outer (int) – forwarded to voi_estimate() unchanged.

• n_inner (int) – forwarded to voi_estimate() unchanged.

• observation_model (GaussianObservationModel | None) – forwarded to voi_estimate() unchanged.

Returns:

The value-of-information point estimate, in dollars. NOT floored at zero – see voi_estimate()’s docstring for why a floor is dishonest, not merely conservative; this can be a small negative number when the true VOI is at or near zero, which is the honest reflection of Monte Carlo noise, not a defect. Use voi_estimate() if you need to tell that case apart from a real, larger difference.

Return type:

float

class VoiStoppingDecision(voi_dollars, sample_cost, net_value, keep_sampling, standard_error)

Bases: NamedTuple

A real decision-theoretic answer to “should we sample again?”: compare the value of one more sample against its cost, rather than an arbitrary uncertainty threshold picked by hand.

standard_error is the underlying VoiEstimate’s Monte Carlo standard error, carried through so a caller can tell a close-to-the-cost net_value (statistical noise) apart from a real, larger gap – the same honesty voi_estimate() reports at the VOI level (MXR-080-0108/0110).

Parameters:

• voi_dollars (float)

• sample_cost (float)

• net_value (float)

• keep_sampling (bool)

• standard_error (float)

voi_dollars: float

Alias for field number 0

sample_cost: float

Alias for field number 1

net_value: float

Alias for field number 2

keep_sampling: bool

Alias for field number 3

standard_error: float

Alias for field number 4

voi_stopping_decision(posterior, decision_fn, drill_info, *, sample_cost, rng, n_outer=64, n_inner=256, observation_model=None)

Should the next sample (drillhole, monitoring well, survey station, …) actually be taken?

A real, principled alternative to picking an uncertainty threshold by hand and telling an LLM loop-controller to stop when it’s cleared (the pattern used, and flagged as a real gap, in experiments/adaptive-groundwater-monitoring and experiments/adaptive-gravity-survey-design): sample again iff voi_dollars() – the actual expected dollar value the next sample would add to the decision – exceeds what that sample costs. As uncertainty tightens, voi_dollars shrinks toward zero (there’s less left to learn that would change the decision), so this converges to a stopping rule on its own without any separately-chosen threshold; the only free parameter is sample_cost, which is an actual real-world number (drilling/survey cost) rather than an arbitrary uncertainty width.

Parameters:

• posterior (Posterior) – forwarded to voi_estimate() unchanged – see its docstring for what each means.

• decision_fn (Callable[[ndarray], float]) – forwarded to voi_estimate() unchanged – see its docstring for what each means.

• drill_info (dict[str, Any]) – forwarded to voi_estimate() unchanged – see its docstring for what each means.

• rng (Generator) – forwarded to voi_estimate() unchanged – see its docstring for what each means.

• n_outer (int) – forwarded to voi_estimate() unchanged – see its docstring for what each means.

• n_inner (int) – forwarded to voi_estimate() unchanged – see its docstring for what each means.

• observation_model (GaussianObservationModel | None) – forwarded to voi_estimate() unchanged – see its docstring for what each means.

• sample_cost (float) – the real dollar cost of taking the next sample. Must be finite (MXR-080-0110: an unvalidated NaN cost previously compared as False against everything, silently forcing keep_sampling=False regardless of the actual VOI).

Returns:

the computed VOI (with its standard error), the cost it was compared against, their difference, and whether sampling should continue (voi_dollars > sample_cost).

Return type:

A VoiStoppingDecision

safety_risk_surface(deformation, *, gradient_limit, slope=None)

Map a deformation field into a per-cell P(tilt > gradient_limit) safety-risk surface.

Parameters:

• deformation (Posterior | np.ndarray) – an IC-1 Posterior over a flattened (d,) subsidence/deformation field (the G4 poroelastic InSAR-inversion case), or a plain np.ndarray spatial grid (any shape) of already-point-estimated deformation values.

• gradient_limit (float) – the tilt/gradient magnitude above which a cell is considered geotechnically unsafe (e.g. an angular-distortion or differential-settlement limit).

• slope (np.ndarray | None) – an optional static terrain-slope field, same spatial shape as the deformation grid. When given, it is added to the deformation-induced gradient magnitude before the exceedance test: a cell that is already steep needs less additional differential movement to cross gradient_limit.

Returns:

A DerivedQuantity whose samples are the per-cell exceedance indicator (0.0/1.0) drawn over the posterior’s Monte-Carlo replicates (or a single deterministic replicate for an ndarray input), flattened in the grid’s row-major (C) order, together with a credible interval and the prior_dominated flag. The per-cell risk probability is samples.mean(axis=0).

Return type:

DerivedQuantity

gradient_limit must be finite and non-negative (a tilt magnitude threshold cannot itself be negative), and slope/an ndarray deformation must be finite (MXR-080-0098): a NaN either way would otherwise compare False against > gradient_limit and silently mark an unevaluable cell as “not exceeding” instead of raising. Posterior evidence additionally requires a finite flattened mean and exactly (_MC_SAMPLES, n_cells) finite draws; the returned binary risk samples are validated independently before publication.

incident_probability(hazard, exposure_map, *, model='logit')

Combine a per-cell hazard probability surface with a people-exposure_map into incident risk.

A hazard exceedance is only a safety incident if someone is exposed to it: an empty, unstable cell and a busy, unstable cell are not the same risk. hazard is expected in [0, 1] (e.g. the per-cell mean of safety_risk_surface(…).samples); exposure_map is a non-negative people- density/occupancy weight of the same shape.

exposure_map is on the expected people present scale: a cell’s value is the mean occupancy of that cell, so 0 means nobody, 1 means about one person on average, and larger values are denser. Both models take that scale literally and return 0 for an unoccupied cell.

Parameters:

• hazard (ndarray) – per-cell hazard probability, shape (*grid_shape,), values in [0, 1].

• exposure_map (ndarray) – per-cell non-negative expected-occupancy weight, same shape as hazard.

• model (str) – “logit” (default) is a logistic link that multiplies the hazard’s own odds by the probability that anybody is actually present, 1 - exp(-exposure_map) (the Poisson occupancy probability for that mean): an unoccupied cell yields 0, denser occupancy monotonically raises the result, and it saturates at the hazard probability itself. “linear” is the simple product hazard * exposure_map clipped to [0, 1].

Returns:

Per-cell incident probability, same shape as hazard, values in [0, 1].

Return type:

ndarray

hazard outside [0, 1] is REJECTED, not clipped into a falsely-confident boundary probability (MXR-080-0098): the pre-fix "logit" path clipped any hazard value – 5.0, -3.0, whatever – straight into [eps, 1 - eps] before ever checking it was a valid probability to begin with. The narrow eps-clip that remains below is purely for the logit transform’s own numerical stability at the exact 0/1 boundary of an ALREADY-valid probability, not a substitute for validating the input range.

The "logit" occupancy factor used to be log1p(exposure_map) added to the log-odds (MXR-080-1589). That has no zero element: log1p(0) is 0, so an empty cell added no log-odds at all and hazard 0.8 at zero occupancy returned incident probability 0.8 – flatly contradicting this function’s own stated contract that an exceedance is only an incident when someone is exposed. It was also unbounded above, so a dense-enough cell reported MORE incidents than hazard exceedances. Multiplying the odds by an occupancy PROBABILITY fixes both: zero occupancy is the zero element, and the incident probability can never exceed the hazard probability it is conditioned on.

class ExceedanceReport(alerts, prob_exceed, false_alarm_target, calibrated=True, warmed_up=<factory>)

Bases: object

Per-timestep exceedance call: which points alert, their raw probability, and the target rate.

alerts is the boolean array the caller (or the mlops drift/retrain wiring) acts on; when calibrated is True it is already conformal-calibrated – do not re-threshold prob_exceed again downstream.

calibrated (MXR-080-0097) is False whenever the returned alerts do NOT carry the distribution-free false-alarm-rate guarantee: calib was omitted (the monitored series – including any anomalies it contains – scored itself, which is neither independent of nor exchangeable with what it is calibrating), or the effective calibration set was too small to give false_alarm_target real resolution (see _min_adequate_calib_size). The alerts are still computed on a best-effort basis in that case, but the caller must treat them as a heuristic, not a proven bound.

warmed_up (MXR-080-0096) marks which timesteps had enough strictly-prior history to be scored at all; an unwarmed timestep always reports prob_exceed == 0.0 and never alerts – an explicit “not yet evaluated”, never a confident “safe”.

The record is frozen and owns read-only copies of its three arrays. calibrated and warmed_up are the honesty flags that qualify alerts; while the record was mutable, flipping calibrated to True or overwriting alerts in place upgraded a best-effort heuristic into an apparently-proven false-alarm bound with nothing else changing.

Construction additionally enforces the record’s own domain (MXR-080-1900), which it previously did not check at all – prob_exceed of 5.0, -3.0 or NaN, a false_alarm_target of 7.5, and a prob_exceed/warmed_up array of a different length than alerts all constructed silently, so the record could report an alert whose “probability” was not one, at a “target rate” that was not a rate, for timesteps that did not line up with the alerts:

• prob_exceed finite and in [0, 1]. It is a survival-function value by construction, so this refuses nothing the producer emits.

• false_alarm_target finite and strictly in (0, 1) – the same range exposure_exceedance_monitor() already requires of its alpha, now also true of a directly constructed report.

• prob_exceed the same length as alerts, one probability per alert.

• warmed_up either the same length as alerts or EMPTY. Empty is the field’s own default and means “warm-up was not recorded”; refusing it would reject reports the package’s own tests construct, and it is not the same claim as a per-timestep flag that disagrees with the alerts.

Parameters:

• alerts (ndarray)

• prob_exceed (ndarray)

• false_alarm_target (float)

• calibrated (bool)

• warmed_up (ndarray)

exposure_exceedance_monitor(series, limit, *, alpha=0.05, calib=None)

Flag exceedance excursions in series against limit at a calibrated false-alarm rate.

Parameters:

• series (ndarray) – (n,) monitoring readings (e.g. silica PM4 concentration over time). Must be finite.

• limit (float) – the occupational/community exposure limit being monitored against. Must be finite.

• alpha (float) – target false-alarm rate (ExceedanceReport.false_alarm_target); the empirical alert rate on exchangeable, non-exceeding data is bounded by this via conformal calibration – but only when the returned calibrated is True (see calib below). Must be in (0, 1).

• calib (ndarray | None) – (m,) held-out reference readings known to be exposure-compliant (sub-limit), drawn from a stable, anomaly-free reference period and scored by the same past-only rule as series – the temporal-exchangeability scheme this monitor’s false-alarm-rate guarantee depends on. When omitted, series scores itself (a graceful degradation for callers with no separate holdout) – but the monitored series, including any anomalies it contains, is neither independent of nor exchangeable with what it is calibrating, so the returned report is explicitly marked calibrated=False and the alpha bound is NOT guaranteed (MXR-080-0097). An explicitly-supplied calib that is too small to give alpha real resolution (fewer than roughly _min_adequate_calib_size(alpha) usable, warmed-up points) is likewise reported uncalibrated rather than silently accepted.

Returns:

An ExceedanceReport. When calibrated is False the alerts are still a best-effort heuristic, not a proven distribution-free bound.

Return type:

ExceedanceReport

Algorithm:

1. At every timestep t of series (and, separately, of calib), fit a causal (strictly past-only) local Gaussian SCALE std[t] from series[max(0, t-window):t] – never series[t] itself, never anything at or after t – via _causal_local_scale(). Timesteps with too little prior history are not “warmed up” (ExceedanceReport.warmed_up) and are reported as an explicit not-yet-evaluated state (prob_exceed=0.0, never alerting) instead of reaching past the window into a global – and, for an early index, necessarily partly-future – statistic (MXR-080-0096).

2. Score prob_exceed[t] as the probability that the TRUE exposure level exceeds limit, treating the OBSERVED series[t] itself as a noisy measurement of that level with scale std[t]: P(X > limit) for X ~ Normal(series[t], std[t]) (scipy.stats.norm.sf()). Centering on the observation itself – rather than on the local mean of prior readings alone, which is blind to what was actually just measured – is what lets a single sudden extreme reading enter its own score and trigger at the exact timestep it happens, instead of only once later readings drag a lagging rolling mean toward it (MXR-080-0096).

3. Calibrate the alert threshold: treat each WARMED-UP calibration timestep’s prob_exceed (scored by the identical observed-value rule) as a one-sided conformal nonconformity score (mixle.inference.conformal.split_conformal(), side="upper", against a constant zero “prediction”) on the known-safe calib set. The returned upper bound is the smallest cutoff such that, under exchangeability with calib, at most an alpha fraction of non-exceeding timesteps would clear it – a distribution-free false-alarm-rate guarantee, not a normal-theory approximation, PROVIDED calib genuinely satisfies that exchangeability and is large enough (see calib above and calibrated).

4. Alert wherever a warmed-up series timestep’s prob_exceed clears that threshold.

class SlopeFactor(oral_csf=None, inhalation_iur=None, sigma_log=0.0, source='EPA-IRIS')

Bases: object

A chemical’s cancer potency, EPA-IRIS style.

Parameters:

• oral_csf (float | None)

• inhalation_iur (float | None)

• sigma_log (float)

• source (str)

oral_csf

oral cancer slope factor, (mg/kg-day)^-1. None if the chemical has no established oral potency.

Type:

float | None

inhalation_iur

inhalation unit risk, (ug/m3)^-1. None if no established inhalation potency.

Type:

float | None

sigma_log

log-scale standard deviation of a multiplicative log-normal uncertainty band around the point-estimate slope factor (0.0 = treat the slope factor as fixed). Must be finite and non-negative – validated at construction (a negative or NaN value has no valid meaning as a standard deviation, so it is rejected here rather than silently treated as zero/fixed uncertainty downstream).

Type:

float

source

provenance tag for the potency values (default the EPA IRIS database, the regulatory-toxicology standard; caller-supplied values should override this).

Type:

str

class RiskQuantity(samples, prior_dominated=False)

Bases: object

A pushforward risk distribution: draws + credible interval + the prior_dominated flag.

Structurally satisfies the IC-1 DerivedQuantity protocol (samples, prior_dominated, credible_interval) and additionally exposes mean – the point estimate callers actually read off first.

Construction validates samples: non-empty, finite, and a valid probability-like value in [0, 1] (the range the excess_lifetime_cancer_risk/radon_wlm_risk LNT pushforward always produces from validated non-negative inputs – see _lnt_risk()). This is defense-in-depth: invalid state cannot silently flow downstream to a caller even if some upstream pushforward fails to validate its own inputs.

Parameters:

• samples (ndarray)

• prior_dominated (bool)

property mean: float

the sample-mean excess risk.

Type:

Point estimate

credible_interval(level=0.9)

Central level credible interval of the risk samples (e.g. level=0.9 -> 5%/95%).

level goes through the shared mixle.analysis interval contract (validated_level()), so every IC-1 carrier in the package rejects the same out-of-range levels the same way (MXR-080-1580).

Parameters:

level (float)

Return type:

tuple[float, float]

excess_lifetime_cancer_risk(exposure, sf, *, route='oral', n=2000, rng=None)

Excess lifetime cancer risk under the linear no-threshold model (DR-ALG K7).

route="oral" treats exposure as an LADD (lifetime average daily dose, mg/kg-day) and multiplies through SlopeFactor.oral_csf; route="inhalation" treats exposure as an air concentration (ug/m3) multiplied through SlopeFactor.inhalation_iur.

Parameters:

• exposure (Posterior | ndarray | float) – the lifetime-average dose/concentration. An IC-1 Posterior (its derived_quantity pushforward is used, so prior_dominated propagates from the exposure posterior), a plain array of exposure samples (already representing exposure uncertainty), or a single deterministic scalar. Every exposure draw must reduce to a single value per sample – one route/chemical/time-point/receptor already selected or aggregated by the caller – and must be finite and non-negative; this is enforced uniformly for all three forms, including a Posterior’s own draws (there is no honest way to guess, on the caller’s behalf, which route/chemical/time point/receptor a genuinely multi-dimensional draw was meant to represent, so that case is rejected with a clear error rather than silently collapsed to its first column).

• sf (SlopeFactor) – the chemical’s SlopeFactor. Its potency coefficient for route, and its sigma_log, must be finite and non-negative (sigma_log is also validated at SlopeFactor construction).

• route (str) – "oral" or "inhalation".

• n (int) – number of posterior draws to take when exposure is a Posterior, or the number of slope-factor draws to take when exposure is a bare scalar and sf.sigma_log > 0. Must be a positive exact integer. On the Posterior path this is now a hard delivery receipt (MXR-080-1900): a posterior that returns fewer than n draws is an error, not a quietly smaller risk distribution.

• rng (Generator | None) – numpy random Generator (a fresh default one is created if omitted).

Returns:

A DerivedQuantity of excess lifetime cancer risk (samples + CI + prior_dominated).

Return type:

DerivedQuantity

A Boolean exposure (scalar or dtype=bool array) is rejected rather than coerced: True is a flag, not a dose, and float(True) == 1.0 made one indistinguishable from a measured 1 mg/kg-day LADD once it was inside the LNT product (MXR-080-1900).

radon_wlm_risk(wlm, *, risk_per_wlm=5.38e-4, n=2000, rng=None)

Radon lung-cancer risk from cumulative working-level-months (BEIR-VI linear coefficient).

risk = wlm * risk_per_wlm for small products (the BEIR-VI committee’s linear excess-relative-risk coefficient per WLM of cumulative radon-progeny exposure, default 5.38e-4 per WLM), saturating via the same 1 - exp(-x) low-dose-linear-no-threshold form _lnt_risk() uses once that product would exceed ~0.01 – without it, cumulative exposures large enough to be realistic over a working lifetime push the bare linear form above 1, which is not a valid probability. Both wlm and risk_per_wlm must be finite and non-negative (exposure and potency are never negative).

Parameters:

• wlm (Posterior | ndarray | float) – cumulative working-level-months. An IC-1 Posterior over exposure (its own derived_quantity pushforward is used, so prior_dominated propagates from the exposure posterior), a plain array of exposure samples (already representing exposure uncertainty), or a single deterministic scalar. Each draw must reduce to one value per sample, exactly as excess_lifetime_cancer_risk() requires – a (n, d>1) matrix mixes route/receptor/time axes this function has no honest way to collapse.

• risk_per_wlm (float) – BEIR-VI risk coefficient per WLM.

• n (int) – number of posterior draws to take when wlm is a Posterior. Must be a positive exact integer. Ignored for the array/scalar paths, which carry their own sample count.

• rng (Generator | None) – numpy random Generator used to draw from a Posterior (a fresh default one is created if omitted). Unused for the array/scalar paths, which are deterministic.

Returns:

A DerivedQuantity of radon-attributable lung-cancer risk.

Return type:

DerivedQuantity

n/rng were previously accepted and never used at all, on any path – a negative, fractional or Boolean draw count was silently taken and discarded (MXR-080-1574). They now drive the Posterior path they were reserved for, and n is validated on every path so an impossible request is rejected rather than ignored. On that path n is additionally a delivery receipt: a posterior returning fewer than n draws raises instead of yielding a thinner risk distribution that reads like the full one (MXR-080-1900). A Boolean wlm is rejected for the same reason a Boolean dose is – float(True) == 1.0 is not one working-level-month.

Submodules

• mixle.analysis.biodiversity module
• mixle.analysis.carcinogenic_risk module
• mixle.analysis.covariance_shrinkage module
• mixle.analysis.coverage module
• mixle.analysis.developmental_risk module
• mixle.analysis.emissions module
• mixle.analysis.epidemiology module
• mixle.analysis.extreme module
• mixle.analysis.habitat_constraints module
• mixle.analysis.health_risk module
• mixle.analysis.kde module
• mixle.analysis.kriging module
• mixle.analysis.max_stable module
• mixle.analysis.objective module
• mixle.analysis.rank_aggregation module
• mixle.analysis.real_options module
• mixle.analysis.sdm module
• mixle.analysis.spatial_mixture module
• mixle.analysis.valuation module