1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
/// Typed errors emitted by the survival location-scale family pipeline.
///
/// Each variant carries a pre-formatted `reason` string so `Display` is
/// byte-equivalent to the original `format!(...)` outputs the module used
/// before the typed-error migration. The category split lets callers
/// pattern-match on the failure kind without dragging the string apart.
impl_reason_error_boilerplate!
// ---------------------------------------------------------------------------
// Overflow-safe arithmetic for the survival exact-Newton chain
// ---------------------------------------------------------------------------
//
// The survival location-scale model computes inv_sigma = exp(-eta_ls) and
// multiplies it through many intermediate quantities (q0, qdot, g, ...).
// When eta_ls is very negative (sigma → 0, distribution very concentrated),
// exp(-eta_ls) can overflow to inf, poisoning downstream sums with NaN via
// inf * 0 or inf - inf patterns.
//
// The protection strategy is layered:
//
// Layer 1 – `exp_neg_stable`: exact exp(-eta_ls) for every value that is
// representable in binary64, saturating near f64::MAX only past the
// representability boundary (`EXP_SATURATION_MAX_ARG` ≈ ln(f64::MAX))
// instead of overflowing to +inf. Underflow (exp(-x) → 0 for large
// positive x) is allowed because it is the mathematically correct
// limit. The saturation never rewrites a finite model: it engages
// only where the true value has no f64 representation.
//
// Layer 2 – `survival_q0_from_eta`: uses exact log-space arithmetic to
// detect when |eta_t * inv_sigma| genuinely exceeds f64::MAX and
// saturates to ±MAX instead of overflowing; when exp(-eta_ls) alone is
// unrepresentable but the product is finite it evaluates the product in
// the log domain, so the value channel matches the mathematical
// function on the whole representable range.
//
// Layer 3 – factorized time-derivative algebra and compensated subtraction:
// the base dq/dt chain is evaluated as exp(-eta_ls) * (eta_t*eta_ls' - eta_t')
// so the shared exp(-eta_ls) factor is applied only once, and
// d_eta/dt = d_raw + qdot is formed with a compensated sum that
// carries an explicit roundoff bound into the monotonicity gate.
//
// Layer 4 – `safe_product` / `safe_sum2` plus `exact_row_kernel`: the generic
// arithmetic guards still clamp inf products to MAX/MIN and map
// inf + (-inf) → 0 as defense in depth, and the row kernel splits the old
// `!g.is_finite()` hard error
// into NaN (hard error for genuinely bad data) and ±inf (clamped to MAX
// so the monotonicity guard can apply).
//
// The invariant: no NaN ever reaches the solver; all overflow paths saturate
// to large finite values that the monotonicity floor and penalty then control.
// ---------------------------------------------------------------------------
// Layer 1 (one-sided overflow guard on the inverse-sigma link), its
// helper `exp_neg_stable`, and `exp_sigma_inverse_from_eta_scalar` now
// live in `crate::sigma_link` so every consumer — solver
// internals here, `main.rs` callers, and any Rust↔Python boundary
// code — picks up the same clamp. Keeping a local copy here previously
// allowed silent semantic divergence between the canonical sigma_link
// version (unclamped) and the survival-local clamped version.