mixle.analysis.health_risk module

Dose-response and population health-risk models (K3, work-plan Workstream K).

Given an exposure/dose – a bare number, an array of dose realisations, or an IC-1 Posterior over a receptor field (K1/K2 transport output) – these push it through a named dose-response curve into an outcome-probability distribution, so downstream liability/constraint code (K6) always has a distribution, never a bare point estimate:

• DoseResponse – a named model (loglinear / logit / hill / threshold_linear) with caller-supplied params; DoseResponse.probability() maps a dose into an IC-1-shaped DerivedQuantity (samples + credible interval + the prior_dominated honesty flag). When dose is a Posterior, the pushforward runs through the posterior’s own derived_quantity so the flag propagates correctly from the exposure uncertainty.

• cumulative_exposure() – trapezoidal time-integration of an exposure series, with an optional first-order biological-decay discount (older exposure counts less toward the current body burden) – the intake feeding a chronic dose-response evaluation.

• population_risk() – aggregates per-receptor dose-response probabilities (one draw per posterior sample, or a single point evaluation for a bare array) into an expected-case-count DerivedQuantity.

This module supplies the dose-response machinery; it ships no regulatory/clinical dose-response table – params are always supplied by the caller or a knowledge lookup (see Non-goals).

K6 (this file’s second addition, appended after K3) turns a K3/K4 risk distribution into the two things J’s objective/optimizer actually need – a priced cost term and a hard feasibility screen:

• health_liability() – prices a case-count/exceedance-probability DerivedQuantity (e.g. population_risk’s output, or K4’s safety_risk_surface) into a dollar-liability DerivedQuantity at cost_per_case per case, discounted – the health_cost term J6’s priced_liabilities/risk_adjusted_plan (analysis/valuation.py) sums alongside carbon (L6) and remediation (G-side) liabilities.

• exposure_constraints() – screens a list of candidate operating options against named occupational/community exposure limits, marking each feasible (or not, naming the binding limit(s)) so H4’s two_stage_stochastic_plan (stochastic_opt.py) only ever optimizes over the surviving, feasible candidate set.

K4 (this file’s third addition) turns a G4 ground-deformation posterior into a spatial safety-risk surface and a people-weighted incident probability:

• safety_risk_surface() – the spatial gradient (tilt) of a deformation field, pushed through the field’s own IC-1 Posterior.derived_quantity pushforward so P(tilt > gradient_limit) per cell comes back as a DerivedQuantity (samples + credible interval + the prior_dominated honesty flag), not a point estimate. An optional static terrain slope adds to the deformation-induced tilt before the exceedance test – a location that is already steep needs less differential settlement to become unsafe.

• incident_probability() – combines a per-cell hazard probability surface with a people- exposure_map (who is where) through a logistic link into a per-cell probability that a hazard cell becomes an actual incident (nobody is at risk on a steep-but-empty cell).

safety_risk_surface accepts either a raw np.ndarray deformation field (a deterministic, already-inverted grid – treated as a single degenerate draw with prior_dominated=False) or an IC-1 Posterior over deformation (the intended G4 case). The frozen Posterior.derived_quantity method does the sampling and the honesty-flag bookkeeping; this module only supplies the gradient/exceedance pushforward function, so the uncertainty accounting is always the posterior implementation’s own (A2), never re-derived here.

Non-goals: no deformation physics (mixle_pde.poroelastic owns the InSAR inversion itself), no economic liability (health_liability, K6).

K5 (this file’s fourth addition) turns a scalar monitoring series into a calibrated real-time exceedance alert – the sibling of mixle.analysis.coverage / mixle.analysis.extreme for the health & safety pillar. It answers a monitoring-shift question – “is this reading trending toward, or past, a regulatory/occupational exposure limit, and can I trust the alert?” – without pretending a single noisy reading settles it:

• exposure_exceedance_monitor() – per-timestep P(exposure > limit) treating the OBSERVED reading itself as a noisy measurement of the true exposure level, scaled by a causal (strictly past-only) local predictive fit around it (the same distribution-over-a-threshold idea as IC-8’s mixle_pde.decision_quantities.prob_exceed, here applied to a scalar monitoring series rather than a spatial posterior field, since a live sensor stream is not itself an IC-1 Posterior) – a sudden extreme reading enters its own score and can trigger at the timestep it happens, not only once later readings drag a lagging rolling mean toward it. The raw probability is then run through mixle.inference.conformal.split_conformal() against a held-out calib reference (known-safe, sub-limit history, drawn from a stable reference period and scored by the same past-only rule) so the alert threshold is distribution-free calibrated: under exchangeability with an adequately-sized calib, the empirical false-alarm rate is bounded by alpha, not just “probably fine” from an untested normal-theory cutoff. When calib is omitted or too small for alpha to have real resolution, the monitored series scores itself as a graceful degradation, but the report comes back explicitly calibrated=False rather than silently implying the bound holds.

An alert firing (ExceedanceReport.alerts.any()) is the hook the mlops drift/retrain half (G7, mixle_mlops.drift_retrain) watches for a re-check/retrain trigger – that wiring is a signal (read this array), not a code dependency: this module never imports anything from mixle_mlops.

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]

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.