chc.uncertainty

Stability tier: evolving.

Calibrated pessimism: a deep-ensemble residual whose disagreement, conformalised for coverage, scores predictive uncertainty for the offline-safety penalty (its U term).

chc.support scores density distance D((x,u),D); this module supplies the complementary calibrated predictive uncertainty U(x,u). Fit K residuals as a deep ensemble; their member disagreement is the epistemic uncertainty of the learned dynamics (large where the members, trained on the same offline data, extrapolate apart), and split conformal turns it into interval widths that cover with probability at least 1 - alpha on data exchangeable with the calibration split, given at least (1 - alpha) / alpha calibration transitions. Both plug into the same penalty channel of chc.support.pessimistic_control via penalty_trajectory(xs, us) -> scalar, so model exploitation is bounded, not merely discouraged.

A third scorer, WassersteinPenalty, targets the deployment-shift failure mode rather than in-distribution epistemic spread: a Wasserstein-1 distributionally-robust margin (radius * Sigma ||d r / d x||) that keeps control where a small shift of the state distribution cannot move the learned dynamics much.

EnsembleResidual

Bases: Module

Deep ensemble of residuals: K members whose mean is the field, whose spread is uncertainty.

Duck-types the Dynamics protocol -- __call__ returns the member mean, so it drops into HybridDynamics.residual. disagreement is the total variance across members (epistemic uncertainty of the learned part), ~0 where members agree and large where they extrapolate apart.

member_outputs(t, x, u)

Stack the K member vector fields at (t, x, u) -> shape (K, out_dim).

disagreement(t, x, u)

Total variance across members at (t, x, u) -- the epistemic uncertainty scalar.

ShardedEnsembleCertificate dataclass

Evidence that the stacked, sharded ensemble fit is the serial per-member recursion.

EnsembleUncertainty

Bases: Module

Predictive-uncertainty penalty U from an ensemble residual -- its member disagreement.

Mirrors SupportModel.penalty_trajectory so it plugs into the same lam_unc channel of chc.support.pessimistic_control. Penalises visiting states/actions where the learned dynamics are uncertain (members disagree) -- exactly where model exploitation happens.

penalty_trajectory(xs, us)

Total member disagreement over the visited (x, u) pairs (xs: (H,n), us: (H,m)).

NestedCVaRPenalty

Bases: Module

Time-consistent (nested) CVaR of the ensemble's disagreement, as a PenaltyModel.

EnsembleUncertainty sums the member variance along the trajectory, which is a risk neutral aggregation: one very bad step averages away against many quiet ones. This replaces the aggregation with a nested conditional risk measure, rho_t = c_t + CVaR_alpha[rho_{t+1}] (Ruszczynski 2010; Shapiro-Dentcheva-Ruszczynski 2021), evaluated on the per-member deviations ||r_k - mean_k r_k||. Same penalty_trajectory(xs, us) -> scalar shape, so it drops into the lam_unc channel of chc.support.pessimistic_control unchanged.

The distinction that matters is which adversary the number prices, and both are available:

• penalty_trajectory (nested) -- the adversary re-picks the worst alpha fraction of members at every step, so a model may be wrong here and a different one wrong there.
• static_penalty_trajectory -- the adversary must commit to one member for the whole horizon, then the tail is taken once over those horizon-summed scenarios.

Nested is therefore never smaller (subadditivity of CVaR, Sum_t CVaR[c_t] >= CVaR[Sum_t c_t] -- a standard coherence property, not a result of this repo), and the gap is the price of time consistency. Time consistency is the reason to pay it: a plan optimal for the static measure at t = 0 need not be the plan you would re-choose at t = 1 after seeing the first step, so a receding-horizon controller that re-solves -- which is what chc.mpc.mpc_control does -- can chase its own tail. The nested measure has no such gap by construction.

HONEST SCOPE: with the members re-evaluated independently at each visited pair, the recursion collapses to Sum_t CVaR_alpha[c_t]; this is the aggregation rule made risk-averse, not a dynamic-programming solve of a nested-risk MDP, and it does not certify a bound on the cost.

member_costs(xs, us)

(H, K) per-member deviation from the ensemble mean field at each visited pair.

penalty_trajectory(xs, us)

Nested: the tail is taken per step, so the worst members may differ across steps.

static_penalty_trajectory(xs, us)

Static: one tail over horizon-summed per-member scenarios (time-INCONSISTENT).

NestedRiskCertificate dataclass

Numeric evidence that the nested measure dominates the static one, and by how much.

SplitConformal

Bases: Module

Split-conformal calibration of the ensemble's next-state uncertainty (coverage guarantee).

Calibrated on a held-out split: the normalised nonconformity score is s = ||x_next - mean|| / (sigma + eps) with sigma the ensemble predictive std; q_hat is the k-th smallest of the n scores, k = ceil((n + 1)(1 - alpha)). interval_width = q_hat*(sigma + eps) has marginal coverage >= 1 - alpha on data exchangeable with the calibration split (measured by coverage). That needs k <= n, so n >= (1 - alpha) / alpha: with fewer scores no finite width carries the guarantee, and calibrate refuses rather than read the largest score, which covers n / (n + 1).

calibrate(model, data, dt, alpha=0.1, eps=1e-06) classmethod

Calibrate on data's transitions, held out from the ensemble's fit.

Raises:

Type	Description
TypeError	when model's residual is not an EnsembleResidual.
ValueError	when alpha is not in (0, 1), eps is not positive, or data holds fewer than (1 - alpha) / alpha transitions.

interval_width(x, u)

Calibrated prediction-interval half-width q_hat * (sigma(x,u) + eps) (scalar).

coverage(data)

Empirical coverage on data: fraction of true next-states within the interval.

penalty_trajectory(xs, us)

Total calibrated interval width over the visited (x, u) pairs (a penalty).

HorizonConformal dataclass

One conformal scale for every step of a horizon: a simultaneous band scale * r_k.

Experimental: it may change or be withdrawn in any release. Each calibration window is scored by max_k ||x_{t+k} - x_hat_{t+k}|| / r_k over the horizon, r_k a fixed profile such as the plan's tube (CausalPlan.uncertainty_tube[1:]), and scale is the ceil((n + 1)(1 - delta))-th smallest of the n scores. Cleaveland, Lee, Pappas and Lindemann (AAAI 2024, Theorem 3 and Corollary 1): on windows drawn i.i.d., a fresh window's every step lies within scale * r_k with probability at least 1 - delta, and, its scores tied with probability 0, at most 1 - delta + 1 / (n + 1). The profile and the model must not be fitted on the calibration windows. Per step at 1 - delta each, the steps together hold less; split delta / H per step (Bonferroni), they hold more than needed.

On one series rather than i.i.d. windows, calibrate with mixing states 1 - delta - eta, eta from mixing_penalty, and refuses where eta >= delta.

What it covers: windows of the policy that logged the calibration windows, under that policy's actions. Not the plan's own actions, which move the plant elsewhere; and it moves no number of the tube, which it reads as a profile only.

band property

scale * r_k at each step: every step's radius at the joint level.

level property

The joint level the band holds: 1 - delta, less eta on one series.

calibrate(errors, profile, delta=0.05, *, mixing=None, train_lag=None) classmethod

Calibrate on errors, (n, H, states), one row per window as horizon_errors returns them, against profile, one positive radius per step.

Parameters:

Name	Type	Description	Default
mixing	Callable[[int], float] | None	on windows cut from one series, beta(j) at a lag of j windows (mixing_penalty); None where the windows are i.i.d.	None
train_lag	int | None	with mixing, the lag in windows from the data the model and the profile were fitted on to the test window (mixing_penalty).	None

Raises:

Type	Description
ValueError	when delta is not in (0, 1), the errors and the profile do not pair, too few windows certify 1 - delta, train_lag is given without mixing, or the series' eta is at least delta.

covers(errors)

Whether each window of errors, (W, H, states), lies within the band at every step.

WassersteinPenalty

Bases: Module

Wasserstein-1 distributionally-robust penalty on the learned residual.

Where SupportModel scores in-distribution density distance D and the ensemble scores epistemic spread U, this scores robustness to a distribution shift of the states. By Kantorovich-Rubinstein W1 duality the worst-case cost over a Wasserstein-1 ball of radius radius equals the empirical cost plus radius times the cost's Lipschitz constant (Gao-Kleywegt; Blanchet-Kang-Murthy; Shafieezadeh-Abadeh-Esfahani-Kuhn). The fragile quantity is the learned residual r (the known physics f_known is trusted, hence excluded), so the penalty is radius * Sigma_t ||d r / d x|| -- a gradient-norm margin that steers control away from states where a small deployment shift moves the learned dynamics a lot. Satisfies the PenaltyModel protocol, so it drops into the lam_unc channel of chc.support.pessimistic_control unchanged; radius is the distribution-shift budget.

from_model(model, radius) classmethod

Penalise the learned residual of a hybrid model (f_known is trusted, so excluded).

local_lipschitz(x, u)

Local Lipschitz constant of the residual in x: Frobenius norm of d r / d x.

penalty_trajectory(xs, us)

W1-DRO margin radius * Sigma_t ||d r/d x|| over the visited (x, u) pairs.

ConfoundingRobustPenalty

Bases: Module

Pessimism against a HIDDEN-CONFOUNDING effect error, for the general OC stack (§34/§38).

The §35 controller robustifies a scalar tracker; this lifts the same sensitivity radius into the multivariate hybrid solver. Under hidden confounding the control-effect matrix B is only partially identified: the offline estimate B_hat sits within a sensitivity half-width radius (the §32 Delta(Gamma)) of the truth. By the §34 dimensional insight a per-step control u_t then carries a per-step transition error ||Delta_B @ u_t|| <= radius * ||u_t||, and radius * Sigma_t ||u_t|| accumulates that over a trajectory. Penalising it steers control away from exploiting a partially-identified effect with large actions -- the confounding analogue of WassersteinPenalty's deployment-shift margin. Satisfies the PenaltyModel protocol (penalty_trajectory(xs, us) -> scalar), so it drops into the lam_unc channel of chc.support.pessimistic_control unchanged.

HONEST SCOPE, and it now depends on which constructor was used. Plain ConfoundingRobustPenalty(radius=...) is an identification-radius regulariser: the §34 inequality bounds the state-transition error, and converting that into a bound on the objective needs the sensitivity multiplier Delta J <= Sigma_t L_{V,t+1} * radius * ||u_t||, which this path does not supply -- so lam_unc absorbs it as an unidentified scale. certified supplies it -- the adjoint norm ||lambda_{t+1}||, the RK4 injection gain and a second-order deviation tube -- and on that path lam_unc = 1 is the bound rather than a knob. confounding_cost_bound_certificate is the gate that can fail it.

Either way the COEFFICIENT is derived from the §32 sensitivity rather than being an arbitrary actuation budget; Gamma and the CVaR-gap calibration remain the analyst's inputs, and it does NOT test for confounding.

from_sensitivity(cvar_gap, gamma) classmethod

Radius = the §32 bounded-density-ratio inflation (Gamma-1)/(Gamma+1) * cvar_gap.

certified(radius, dyn, x0, us_reference, dt, cost) classmethod

Supply the missing multiplier, so lam_unc = 1 is the bound rather than a free knob.

The §34 inequality bounds the per-step transition error; turning it into a bound on the objective needs Delta J <= Sigma_t L_{V,t+1} * radius * ||u_t||, and L_{V,t+1} is locally the adjoint norm ||lambda_{t+1}||. That is a quantity the planner already computes, so leaving it to be absorbed by lam_unc was throwing away information, not avoiding an assumption. chc.adjoint.costate_norms returns it along us_reference.

The first-order term alone is not an upper bound, and measuring rather than assuming is what showed it: at the optimum the Cauchy-Schwarz step is nearly tight, so the dropped curvature term -- positive, O(radius^2) -- pushes the realised worst case above it at every radius: an adversary reaches 1.014 of it at radius 0.005 and 1.69 at 0.2 on a two-lever oscillator, while the weights below hold at 0.995 and 0.953. chc.adjoint.perturbation_cost_weights therefore carries the deviation tube and the cost curvature as well, which is why radius is an argument here: a second-order term cannot be folded into a radius-free weight.

Three things this is and is not:

• Second order, exact for LQ. For a linear plant with a quadratic cost the objective is exactly quadratic along a perturbation direction, so the expansion closes and the only slack left is Cauchy-Schwarz. Nonlinear plants keep an O(radius^3) remainder; confounding_cost_bound_certificate measures whether it matters.
• At a reference. The weights are frozen at us_reference, which breaks the circularity of weights that depend on the plan that depends on the weights. Iterating (re-solve, re-weight) is available to the caller and is not done here, because a fixed point is a different claim from a bound at a named trajectory.
• Not a confounding test. radius still comes from an assumed Gamma.

penalty_trajectory(xs, us)

Confounding pessimism radius * Sigma_t L_t ||u_t|| (xs unused).

L_t is 1 when no weights were supplied -- the identification-radius regulariser -- and chc.adjoint.perturbation_cost_weights when they were, on which path the sum is a certified upper bound on the cost gap rather than a scale-free direction.

LipschitzRolloutCertificate dataclass

Numeric evidence: the measured rollout deviation stays under the certified Gronwall bound.

ContractiveRolloutCertificate dataclass

Evidence a contractive residual has log-norm mu<0 and a flat (bounded) rollout radius.

TimeVaryingRolloutCertificate dataclass

Evidence the time-varying tube is a valid, tighter-than-constant envelope with a cutoff.

ConfoundingRobustCertificate dataclass

Evidence the closed-form MSM worst-case equals the brute-force sharp bound and behaves.

ConfoundingRobustClosedLoopCertificate dataclass

Evidence confounding widens the closed-loop tube and shortens the certified-safe horizon.

ConfoundingBoundCurve dataclass

Predicted vs adversarially realised cost gap, radius by radius. ok iff none exceeded.

GammaBenchmark dataclass

Gamma expressed in units of the confounding the OBSERVED covariates actually carry.

GammaBenchmarkCertificate dataclass

Falsifiable gates for benchmark_gamma and negative_control_gamma.

fit_ensemble(model, data, dt, n_members=5, width=16, depth=2, bootstrap=True, seed=0, steps=2000, lr=0.01)

Fit K residuals into an EnsembleResidual (split keys + optional bootstrap), in parallel.

Each member is a fresh MLPResidual (dims inferred from data = {x, u, x_next}) trained by the same Adam recursion chc.train.fit_residual runs. With bootstrap each member sees a row-resample of the data, so disagreement reflects both init randomness and data variation. Returns HybridDynamics(model.known, EnsembleResidual) and the per-member loss histories.

The members are trained together: one stacked parameter pytree, jax.vmap over the member axis, jax.lax.scan over the Adam steps, and the stack sharded over the member axis of _member_mesh. This is one XLA program rather than K x steps dispatches, and it is the same recursion -- sharded_ensemble_certificate pins the agreement in ULP.

sharded_ensemble_certificate(seed=0, n_members=8, n_samples=96, parity_steps=60, dt=0.05)

Check fit_ensemble against a serial oracle built from the untouched public trainer.

Two claims, each pinned by the measurement that can falsify it. Numerics: every member of the stacked fit matches the member chc.train.fit_residual produces from the same key and the same bootstrap rows, reported in ULP of the working dtype rather than as an absolute -- the threshold is a precision claim, so it is stated in the unit precision actually has. Layout: the member axis is genuinely distributed, so a stacked array committed to _member_mesh and consumed by a jitted vmap comes back spanning the whole mesh; on a one-device host that is trivially 1, and the subprocess test forces eight.

parity_steps is deliberately short. Adam on this loss is chaotic: past a few hundred steps a one-ULP difference in reduction order is amplified without bound (measured 6.6 ULP at 200 float32 steps, 7557 at 400), so a long-horizon parity threshold would be either vacuous or flaky. The equivalence certified is of the recursion, which a short horizon tests exactly.

cvar_upper(values, alpha)

Upper CVaR (superquantile) of equally-weighted values at level alpha.

Mean of the worst alpha fraction, with the exact fractional-tail weight rather than a rounded count, so it agrees with _top_tail_mean's convention and is continuous in alpha. Differentiable: jnp.sort is a permutation, so this is a fixed linear functional of the order statistics almost everywhere. alpha = 1 gives the mean, alpha -> 0 the max.

nested_risk_certificate(xs, us, ensemble, alpha=0.2)

Check nested >= static >= risk-neutral, and that alpha = 1 collapses the three.

horizon_errors(model, states, actions, dt, horizon, *, stride=None)

The (W, horizon, n) errors of the model's open loop over one unit's logged windows.

Window i starts at t = i * stride (stride defaults to horizon, so the windows tile the log and each transition is in one) and runs the model's RK4 steps from the logged x_t under the logged u_t, ..., u_{t + horizon - 1}. Its row k - 1 is x_{t + k} - x_hat_{t + k}, the plant's state less the open loop's at step k. HorizonConformal.calibrate scores them.

Parameters:

Name	Type	Description	Default
model	Dynamics	the field the plan rolls out, read as chc.integrate.rollout reads it.	required
states	ArrayLike	(T + 1, n), one unit's consecutive logged states.	required
actions	ArrayLike	(T, m), the actions logged between them.	required
dt	float	the step.	required
horizon	int	the window's number of steps, at least 1.	required
stride	int | None	the steps between window starts, at least 1.	None

Raises:

Type	Description
ValueError	on states and actions of lengths that do not pair, a horizon or stride below 1, entries that are not finite, or a log too short for one window.

mixing_penalty(mixing, calibration, train_lag=None)

eta, the coverage that split conformal pays on one beta-mixing series.

Oliveira, Orenstein, Ramos and Romano (JMLR 2024, Theorem 4): on a stationary beta-mixing sequence, the score at a test index i falls within the empirical 1 - alpha quantile of calibration consecutive scores with probability at least 1 - alpha - eta, eta = eps_cal + eps_train + delta_cal, where eps_train = beta(i - n_train) and

eps_cal = inf over (a, m, r) of sigma(a) sqrt(4 L / (n - r + 1)) + L / (3 m) + (r - 1) / n

with L = log(4 / (delta_cal - 4 (m - 1) beta(a) - beta(r))), 2 m a = n - r + 1, delta_cal > 4 (m - 1) beta(a) + beta(r) and sigma(a) = sqrt(1/4 + (2 / a) sum_{j < a} (a - j) beta(j)). The theorem holds at every delta_cal > 0, so this takes the least eta over every block (a, m, r) and a grid of delta_cal: any point of it is a valid penalty, the least one the tightest found.

Parameters:

Name	Type	Description	Default
mixing	Callable[[int], float]	beta(j), the beta-mixing coefficient of the score sequence at a lag of j scores, in [0, 1] (the supremum over events of the gap between the joint law of the past and the future and their product, Bradley's convention). A value outside [0, 1] is clipped into it. gaussian_mixing gives a bound for windows of a stable Gaussian VAR(1).	required
calibration	int	n, the number of consecutive calibration scores.	required
train_lag	int | None	the lag, in scores, from the last score the model or the profile was fitted on to the test window, which adds beta(train_lag); None where they were fitted on data independent of the series, which adds nothing.	None

Raises:

Type	Description
ValueError	on fewer than 2 calibration scores or a train_lag below 1.

gaussian_mixing(state_matrix, noise_covariance, *, stride=1, span=0)

A bound on the beta-mixing coefficient of a stable Gaussian VAR(1) s' = F s + w, w ~ N(0, Q), at a lag of j windows that start stride steps apart and each span span steps: beta(j * stride - span), or 1 where the windows meet.

For a stationary Markov chain beta(k) is the mean over the stationary law of the total variation from P^k(s, .) to it. By Pinsker's inequality and Jensen's, it is at most sqrt(E KL / 2), and for this chain E KL = -log det(I - S^-1 F^k S F^k') / 2, S the stationary covariance: beta(k) <= sqrt(-log det(I - S^-1 F^k S F^k')) / 2 (validation/tube_levels.mac). On one state of correlation rho it is sqrt(-log(1 - rho^(2k))) / 2.

s must be the whole state the scores are read off: on a log under a policy with its own persistence, the policy's state belongs in it.

Raises:

Type	Description
ValueError	on a matrix that is not square, a covariance not of its size, not symmetric or not positive definite, a spectral radius of 1 or more, or stride below 1 or span below 0.

lipschitz_rollout_bound(lipschitz, model_error, dt, horizon)

Certified H-step trajectory-error bound from a Lipschitz field via discrete Gronwall.

Two Euler rollouts of an L-Lipschitz field, one perturbed per step by <= model_error, deviate by e_k obeying e_{k+1} <= (1 + L*dt)*e_k + dt*model_error, e_0 = 0, so

``e_H <= model_error * ((1 + L*dt)^H - 1) / L``   (``L>0``; the ``L->0`` limit is linear,
``model_error * dt * H``, no blow-up).

Turns the CERTIFIED constant L of LipschitzResidual (plus the trusted known-field norm) into a certified pessimism radius that propagates a per-step error BUDGET through time. This COMPLEMENTS (does not replace) chc.support.SupportModel: the support density-distance scores WHERE the one-step error model_error is large (off-support), and the bound propagates it (see support_calibrated_error). model_error should be a one-step bound (a conformal upper quantile or a bounded-disturbance envelope), not a validation average. Machine-checked in proofs/lipschitz_rollout.v (rollout_error_bound); derived in validation/lipschitz_rollout.mac. HONEST SCOPE: the bound is exp(L*T) (T = H*dt) -- tight for small L*T (bounded-gain residual, short horizon, safety-critical), loose else; a contraction metric (ContractiveResidual, log-norm mu<0) removes the exponential; the open loop assumes a FIXED action sequence -- for re-planning use closed_loop_rollout_bound; for per-step L_k use time_varying_rollout_bound.

Raises:

Type	Description
ValueError	if lipschitz or model_error is negative, infinite or nan, if dt is not a finite positive step, or if horizon is not a whole number of steps. A negative budget, step or horizon reads a negative radius, which certifies every step. A bound past the float range is inf.

lipschitz_rollout_certificate(seed=0, state_dim=2, horizon=8, dt=0.05, model_error=0.1)

Roll two fields differing by a bounded per-step error; confirm the deviation obeys the bound.

The field is f_known + r with f_known a linear map of known spectral norm and r a LipschitzResidual with certified constant L_r; L = ||A||_2 + L_r upper-bounds its x-Lipschitz constant. A fixed error of norm model_error is added each step (the learned part is the fragile one); the Euler deviation from a shared start stays under lipschitz_rollout_bound(L, model_error, dt, horizon).

contractive_rollout_bound(contraction_rate, lipschitz, model_error, dt, horizon)

UNIFORMLY BOUNDED explicit-Euler rollout radius for a CONTRACTING field (log-norm mu=-c<0).

The discrete Euler contraction factor is q = sqrt(1+2*mu*dt+L^2*dt^2) (NOT 1+mu*dt; L = the full Lipschitz constant), since two Euler steps deviate by |d+dt*(f(x)-f(y))|^2 <= (1+2*mu*dt+L^2*dt^2)|d|^2 (Maxima contractive_euler.mac; Rocq contractive_euler.v). Contraction (q<1) needs the step dt<2c/L^2 (sufficient, from the worst-case (mu,L) bounds); under it the discrete Gronwall bound e_H <= model_error*dt*(1 - q^H)/(1 - q) is capped at model_error*dt/(1-q) for ALL horizons -- no e^{L*T} blow-up -- tends to the continuous radius model_error/c as dt -> 0. The payoff of a certified negative log-norm (ContractiveResidual) over the non-negative ||.||-Lipschitz of LipschitzResidual; the q < 1 cap uses gronwall_bounded. Returns +inf if not contracting or the step is too large (dt >= 2c/L^2), where Euler overshoots.

Raises:

Type	Description
ValueError	if contraction_rate is not finite, if lipschitz or model_error is negative, infinite or nan, if dt is not a finite positive step, or if horizon is not a whole number of steps.

contractive_rollout_certificate(seed=0, state_dim=3, horizon=60, dt=0.05, model_error=0.1)

Confirm a ContractiveResidual contracts and its radius stays flat.

Checks (i) the empirical one-sided Lipschitz <r(x)-r(y),x-y>/||x-y||^2 <= -rate (certified contraction), and (ii) the perturbed-rollout deviation stays under the BOUNDED explicit-Euler radius eps*dt/(1-q) (step dt<2c/L^2) even as the horizon grows -- contrasted with the e^{L*T} that the same field's norm-Lipschitz L would incur (lipschitz_blowup).

time_varying_rollout_bound(lipschitz, model_error, dt, *, integrator='rk4')

Per-step certified error tube e_0..e_H for time-varying L_k and BUDGET eps_k.

Two rollouts by integrator from one start, of an L_k-Lipschitz field and of the field off by at most eps_k wherever step k reads it, stay within e_k of each other, e_0 = 0:

• "rk4", the integrator of chc.integrate.rollout and so of every plan: e_{k+1} = e_k + dt (L_k e_k + eps_k) phi(L_k dt), phi(z) = 1 + z/2 + z^2/6 + z^3/24, from bounding RK4's four stages one after another (Rocq rk4_rollout_error_bound, validation/rk4_rollout.mac). Its growth 1 + z phi(z) is RK4's stability polynomial, and on a linear scalar field with a constant error the tube is the two rollouts' gap exactly;
• "euler": e_{k+1} = (1 + L_k dt) e_k + dt eps_k (Rocq gronwall_var_comparison), a bound on explicit-Euler rollouts only. On RK4 rollouts it is not one: x' = x against x' = x + 0.1 at dt = 1 reads 0.1, 0.3, 0.7 where the RK4 gaps are 0.1708, 0.6335, 1.8866.

Unlike the constant-L lipschitz_rollout_bound, this exposes WHICH step / channel drives the growth -- feed certified_horizon for the honest certified_until_step. eps_k is a budget: it should be a CERTIFIED per-step bound (a conformal quantile, a bounded-disturbance envelope, or support_calibrated_error combining a base error with chc.support.SupportModel), not a validation-set average.

Raises:

Type

Description

ValueError

if an L_k is negative, infinite or nan. L_k bounds ||f(x) - f(y)|| / ||x - y||, which no field takes below zero; a contracting field's negative log-norm shrinks a tube only through linear_rollout_bound's matrices or contractive_rollout_bound. In either recursion it can turn the radii negative, which would certify every step, and an infinite one reads inf * 0, a nan. Also if an eps_k is negative, infinite or nan, which no budget is, or if dt is not a finite positive step.

certified_horizon(lipschitz, model_error, dt, tolerance, *, integrator='rk4')

The largest step H whose certified error e_H stays within tolerance.

Past it the plan is flagged uncertain; e is monotone, so this is the first crossing. A radius that overflows to inf is within no tolerance, inf included. The tube is time_varying_rollout_bound's for integrator, which raises as it does; a tolerance that is negative or nan raises too.

linear_rollout_bound(state_matrix, model_error, dt)

The RK4 error tube e_0..e_H of a field affine in the state, x' = A x + g(t, u).

Two RK4 rollouts from one start, of the field and of the field off by at most eps_k wherever step k reads it, under the same actions, stay within

``e_k = sum_{j<k} eps_j g_{k-1-j}``,  ``g_m = dt sum_i ||R(dt A)^m c_i(dt A)||_2``,

R(z) = 1 + z + z^2/2 + z^3/6 + z^4/24 being RK4's stability polynomial and c_1..c_4 the weights its stages put on the field's errors, (1 + z + z^2/2 + z^3/4)/6, (2 + z + z^2/2)/6, (2 + z)/6 and 1/6, which sum to phi: one step's gap is exactly R(dt A) e + dt sum_i c_i(dt A) delta_i (validation/rk4_rollout.mac). The matrices, not a norm of A, carry the growth, so a contracting A shrinks the tube where a norm-Lipschitz one could only widen it; on a scalar A = L >= 0 it is time_varying_rollout_bound's RK4 tube.

Raises:

Type	Description
ValueError	if state_matrix is not a finite square matrix, an eps_k is negative, infinite or nan, or dt is not a finite positive step.

closed_loop_rollout_bound(state_lipschitz, control_lipschitz, policy_lipschitz, model_error, dt, horizon)

Closed-loop rollout radius when the plan is RE-PLANNED each step by an L_pi-Lipschitz pi.

A state error perturbs the action, which perturbs the next state, so the growth rate is L_x + L_u*L_pi (not L_x alone): e_{k+1} <= (1 + dt(L_x + L_u*L_pi)) e_k + dt*eps. Reduces to lipschitz_rollout_bound with the combined constant. CAVEAT: for a clipping / active-set / threshold MPC L_pi may be unbounded (not globally Lipschitz); use this only where the controller is Lipschitz (a fixed feedback gain, or a smooth relaxation).

Raises:

Type	Description
ValueError	if a Lipschitz constant is negative, infinite or nan: a negative one lowered the combined rate below the state's own. Otherwise as lipschitz_rollout_bound.

support_calibrated_error(base_error, scale, support_distances)

Per-step error BUDGET eps_k = base_error + scale * D(x_k, u_k) from a support distance.

This COMBINES (does not replace) chc.support.SupportModel: the certified rollout bound propagates a per-step budget through time, and the SupportModel says WHERE that budget is large (off-support, where the one-step model error grows). SupportModel scores the local risk, the Gronwall tube propagates it.

constraint_margin(nominal_g, lipschitz_g, error_radii)

Tightened constraint margin g(x_hat_k) + L_g * e_k for an L_g-Lipschitz constraint.

By Rocq constraint_tightening, if this margin is <= 0 at step k then the TRUE trajectory (within the certified radius e_k of the nominal x_hat_k) satisfies g(x_k) <= 0 -- the tube becomes a robust feasibility check. Admissible while all <= 0.

time_varying_rollout_certificate(seed=0, horizon=20, dt=0.05, tolerance=0.3)

Rising per-step L_k + a bumpy budget: the tube stays under the constant-max-L bound.

Identifies certified_until_step and the safety margin gives a safe_until_step cutoff.

confounding_robust_inflation(cvar_upper, cvar_lower, gamma)

MSM confounding inflation (Gamma-1)/(Gamma+1) * (CVaR_up - CVaR_lo) over a point effect.

Under a BOUNDED DENSITY-RATIO sensitivity model (the weight lies in [1/Gamma, Gamma] with mean 1 -- a marginal, covariate-unconditional special case of Tan's MSM, whose full form constrains the treatment-assignment odds with propensity/covariate-dependent bounds) the sharp worst-case E[wY] puts w=Gamma on the top-tau tail (tau=1/(Gamma+1), mean-preserving) and w=1/Gamma elsewhere. The gap over the nominal mean is this closed form (Maxima confounding_robust_cvar.mac; Rocq msm_inflation_*): 0 at Gamma=1 (point ID), nonnegative, monotone in Gamma. Dorn-Guo (2023); Oprescu et al., B-Learner (2023); Tan (2006).

For a continuous lever the same closed form is the sharp bound on the mean outcome at one dose level under the continuous MSM (Jesson et al. 2022; Frauen, Melnychuk and Feuerriegel 2023, Thm 1). It bounds a level, not a slope: a contrast of two levels is bounded by the upper bound at one less the lower bound at the other (Frauen et al., Lemma 3). A channel radius spent as radius * ||u|| (ConfoundingRobustPenalty) is a structural sensitivity model of its own, not the MSM's bound.

msm_worst_case_mean(outcomes, gamma)

Sharp MSM worst-case (upper) mean of outcomes over the [1/Gamma, Gamma] weight box.

The confounding-robust pessimistic value a controller should plan against when the treated outcome sample may be confounded up to sensitivity Gamma. Equals mean + inflation with the CVaR tails of confounding_robust_inflation; reduces to the sample mean at Gamma=1.

confounding_robust_radius(nominal_radius, outcomes, gamma)

Inflate a pessimism radius by the MSM confounding gap: rho0 + (worst_case - mean).

Rocq robust_radius_ge_nominal / robust_radius_monotone: the returned radius is never below nominal_radius (pessimism only grows under assumed confounding) and is monotone in Gamma. Feeds chc.support.pessimistic_control as a widened uncertainty budget.

confounding_robust_certificate(seed=0, n=8, gamma=3.0)

Confirm the CVaR closed form matches the sharp LP worst-case on a mean-preserving grid.

Chooses (n, Gamma) so tau*n = n/(Gamma+1) is an integer, making the sharp box-LP optimum a pure top-tau*n assignment (weight Gamma on the largest points, 1/Gamma on the rest) that a brute-force sort computes exactly -- a real check, no LP solver needed.

confounding_robust_closed_loop_bound(state_lip, control_lip, policy_lip, base_error, control_magnitude, cvar_gap, gamma, dt, horizon)

Closed-loop rollout radius (Result 31) with the per-step budget inflated by MSM confounding.

Composes the confounding radius (§32) with the closed-loop replan tube (§31): the growth rate is L_x + L_u*L_pi (re-planning feeds state error through an L_pi-Lipschitz policy). The §32 half-width Delta_B = confounding_robust_inflation(cvar_gap, 0, gamma) is an effect (B) error; the per-step state-transition error it induces is ||Delta_B * u_k|| <= control_magnitude*Delta_B, so the budget is eps_k = base_error + control_magnitude*Delta_B (the action multiplier makes the units match -- effect error times control is a state error). gamma=1 recovers the plain closed-loop tube. Rocq closed_loop_confounding_monotone: the tube is monotone in the confounding at every horizon.

confounding_robust_closed_loop_certificate(gamma=3.0, horizon=25, dt=0.05)

Show the §32-inflated §31 closed-loop tube widens and its safe horizon shrinks with Gamma.

A re-planning controller (L_pi-Lipschitz policy) on an L_x/L_u plant, a nominal margin to a constraint g <= 0: as the assumed MSM sensitivity Gamma grows the confounding-inflated tube eats the margin sooner, so the safe planning horizon contracts.

confounding_cost_bound_certificate(radii=(0.005, 0.01, 0.02, 0.05, 0.1, 0.2), horizon=20, dt=0.05, restarts=6, ascent_steps=250, seed=0)

Try to break ConfoundingRobustPenalty.certified, and report where it bends.

A two-lever damped oscillator, planned to optimality, then attacked: the adversary picks the control-channel error M with ||M||_2 <= radius that moves the objective most, by projected gradient ascent from several starts on both signs of Delta J. Random sampling was the first version and it is not good enough -- in four parameters it underestimates the worst case by enough to hide a violated bound.

Two ratios are reported because they answer different questions. ratio is the certificate: realised over the full second-order weight, and ok demands every entry stay at or below one. first_order_ratio is the diagnostic that motivated the second-order term: it exceeds one at every radius, which is why ||lambda|| alone was a calibrated estimate, not a bound.

benchmark_gamma(treated, covariates, assumed_gamma, names=None, quantile=0.95, ridge=1e-06)

Price an assumed MSM Gamma against the observed covariates, in MSM units.

Gamma is the analyst's unfalsifiable input, and Result 32 says so every time it is used. It is unfalsifiable; it is not uncalibrated. The MSM bounds the assignment odds ratio between the true propensity and the modelled one uniformly over units, and dropping an observed covariate produces exactly such a pair -- so

``Gamma_j = exp( quantile_i | logit e(x_i) - logit e_{-j}(x_i) | )``

is the sensitivity a confounder as strong as covariate j would generate, measured rather than assumed.

quantile = 1.0 is the MSM's own quantity and it is not the default, because it grows with the sample. The model bounds the odds ratio uniformly over units, so the sup is the faithful statistic -- but under an unbounded covariate the sup is an extreme order statistic. Measured on a standard-normal design, the sup runs 309 -> 382 -> 397 -> 734 as n goes 500 -> 4000 -> 32000 -> 128000, while the 95th percentile sits at 22, 24, 19, 19 and the median at 2.7, 3.0, 2.7, 2.8. A benchmark that quadruples because more data arrived is not a benchmark, so the default reports a quantile and says which one; pass quantile=1.0 when the uniform bound is the point and the covariates are bounded.

multiples_of_strongest then reads Gamma as "an unobserved confounder k times as strong as the strongest thing we did observe", on the log scale odds ratios compose in. A k far below 1 is an assumption nobody should be impressed by; a k far above 1 is one the analyst has to defend.

SCOPE. Per-covariate and one-at-a-time: this is the marginal strength of each covariate, not a policy-level aggregation over a set of them, and the two differ whenever covariates are correlated. Dropping a covariate from the propensity is a benchmark, not a claim that the unobserved confounder resembles it. The propensity is logistic-linear in the columns as passed; a benchmark is only as good as that model, which is why the covariates should already carry the basis expansion the analyst believes.

ridge penalises each slope of the propensity in units of its column's variance, so an implied Gamma reads the same in any units and from any origin of the covariates.

negative_control_gamma(outcomes, tol=1e-09, gamma_max=1000000.0)

Smallest Gamma whose MSM interval covers zero on a KNOWN-NULL outcome.

The other half of calibration, and the half that can refute. On an outcome whose true effect is zero, any nonzero estimate is confounding, so the smallest Gamma that reconciles the two is a lower bound on the confounding actually present: assuming less than it is refuted by the data rather than merely unappealing. That transfers to the outcome of interest only when the control is U-comparable, confounded by the same unobserved causes (Lipsitch, Tchetgen Tchetgen and Cohen 2010). A control confounded by more than the outcome overstates the outcome's Gamma, and one confounded by less understates it (Baitairian et al., arXiv 2510.16560). Monotonicity of the MSM inflation in Gamma makes the search a bisection, exact to tol.

Returns inf when no Gamma reconciles the null -- the endpoint saturates at the extreme order statistic as Gamma -> inf, so a sample whose values all share the mean's sign cannot be pulled across zero by a bounded density ratio at all, and the negative control has refuted the model class instead of calibrating it.

The interval is NOT symmetric about the mean, and which endpoint has to travel depends on the sign of the estimate: a positive estimate is reconciled by the lower endpoint reaching zero, a negative one by the upper. Reflecting the sample by -sign(mean) maps both onto the single upper-endpoint routine msm_worst_case_mean computes, which is why the sharp CVaR tails are used in the direction that actually binds rather than the convenient one.

gamma_benchmark_certificate(strengths=(0.25, 0.5, 1.0, 2.0), sizes=(500, 4000, 32000), biases=(0.1, 0.2, 0.4), quantile=0.95, samples=4000, seed=0)

Check that the benchmark measures confounding strength and that the calibration binds.

Five gates, each able to fail:

1. Monotone in strength. A covariate whose logit coefficient grows induces a larger implied Gamma. A benchmark that does not order confounders by strength orders nothing.
2. The null covariate's score is a root-n noise floor. Dropping a coefficient that is truly zero must not manufacture sensitivity beyond sampling noise -- and "beyond sampling noise" is a rate, not a constant, so the gate is that log(Gamma_null) * sqrt(n) stays bounded as n grows rather than that a single draw sits below a fixed number. Measured on six seeds per size, that product runs 2.1, 4.0, 5.0, 4.4 at n = 1000 .. 64000 while Gamma_null itself falls 1.069 -> 1.017; a single draw of 1.27 at n = 4000 is inside that distribution, which is exactly why one draw cannot be the gate.
3. Ranking. On one design carrying strong / weak / null covariates, the implied Gamma ranks them in that order.
4. The reported quantile is stabler in n than the sup. This is the measurement the default rests on, so it is run rather than remembered: the sup is an extreme order statistic under an unbounded covariate and grows with the sample; the quantile should not.
5. The negative control lands on the endpoint, and grows with the planted bias. At the returned Gamma the binding interval endpoint sits at zero to solver tolerance, and a larger planted bias needs a larger Gamma -- plus inf when the sample cannot be reconciled at all.