cellstate_core/

decay.rs

//! Unified parametric decay module.
//!
//! All temporal and spatial decay in CELLSTATE is unified under the **Weibull
//! survival function**, a single parametric family that subsumes exponential
//! decay, heavy-tail decay, and sharp-cutoff decay:
//!
//! ```text
//! S(x; τ, β) = exp( -ln(2) · (x / τ)^β )
//! ```
//!
//! where:
//! - `x ≥ 0` — independent variable (time in days, graph depth in hops, …)
//! - `τ > 0` — **half-life parameter**: `S(τ; τ, 1) = 0.5`
//! - `β > 0` — **shape parameter**:
//!   - `β = 1.0` → standard memoryless exponential (current recency default)
//!   - `β > 1.0` → super-exponential with sharp sigmoid-like cutoff
//!   - `β < 1.0` → sub-exponential with heavy tail (graph proximity, etc.)
//!
//! # Mathematical Properties (all verified by property tests)
//!
//! - **P1 (Identity)**: `S(0; τ, β) = 1` for all `τ, β > 0`
//! - **P2 (Half-life)**: `S(τ; τ, 1) = 0.5`
//! - **P3 (Monotonicity)**: `x₁ ≤ x₂ ⟹ S(x₁) ≥ S(x₂)`
//! - **P4 (Bounds)**: `0 ≤ S(x) ≤ 1` for all `x ≥ 0`
//! - **P5 (Asymptotic)**: `lim_{x→∞} S(x) = 0`
//! - **P6 (Smoothness)**: `S` is C^∞ for `x > 0`
//!
//! # Subsumption of Legacy Decay Functions
//!
//! | Legacy function               | Equivalent call                           |
//! |-------------------------------|-------------------------------------------|
//! | `0.5^(t/7)` (note recency)   | `parametric_decay(t, &{7.0, 1.0})`       |
//! | `0.5^(t/30)` (artifact rec.) | `parametric_decay(t, &{30.0, 1.0})`      |
//! | `0.5^(t/7)` (warmth decay)   | `parametric_decay(t, &{7.0, 1.0})`       |
//! | `1 - d/D` (graph linear)     | Replaced by `parametric_decay(d, &{…})`  |
//!
//! Re-export path: `cellstate_core::decay::*`

use serde::{Deserialize, Serialize};

/// Configuration for a parametric decay function.
///
/// Encapsulates the two parameters of the Weibull survival function.
/// Builders configure these per scoring factor via `ScoringDecayConfig`.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct DecayParams {
    /// Half-life parameter τ: the value of `x` at which decay reaches 0.5
    /// when `shape = 1.0`. Must be positive.
    pub half_life: f32,

    /// Shape parameter β controlling the decay curve:
    /// - `1.0` = standard exponential (memoryless)
    /// - `> 1.0` = super-exponential (sharp cutoff near τ)
    /// - `< 1.0` = sub-exponential (heavy tail past τ)
    ///
    /// Must be positive.
    pub shape: f32,
}

impl DecayParams {
    /// Create params for standard exponential decay with the given half-life.
    #[inline]
    pub fn exponential(half_life: f32) -> Self {
        Self {
            half_life,
            shape: 1.0,
        }
    }

    /// Create params for heavy-tail decay (useful for graph proximity).
    #[inline]
    pub fn heavy_tail(half_life: f32) -> Self {
        Self {
            half_life,
            shape: 0.65,
        }
    }

    /// Create params for sharp cutoff decay (useful for time-sensitive windows).
    #[inline]
    pub fn sharp_cutoff(half_life: f32) -> Self {
        Self {
            half_life,
            shape: 2.0,
        }
    }

    /// Validate that parameters are in the domain.
    pub fn validate(&self) -> bool {
        self.half_life > 0.0
            && self.shape > 0.0
            && self.half_life.is_finite()
            && self.shape.is_finite()
    }
}

impl Default for DecayParams {
    fn default() -> Self {
        Self::exponential(7.0)
    }
}

/// Weibull survival function: `S(x; τ, β) = exp(-ln(2) · (x/τ)^β)`
///
/// Returns a value in `[0.0, 1.0]` representing the decay of a signal
/// at distance/time `x` from origin.
///
/// # Panics
///
/// Debug-asserts that `params.half_life > 0` and `params.shape > 0`.
/// In release builds, invalid params produce `NaN` or `0.0`.
///
/// # Properties
///
/// All properties below are verified by the property test suite:
///
/// - **P1**: `parametric_decay(0, _) == 1.0`
/// - **P2**: `parametric_decay(τ, {τ, 1.0}) ≈ 0.5`
/// - **P3**: Monotone decreasing in `x`
/// - **P4**: Output ∈ `[0.0, 1.0]`
/// - **P5**: Approaches `0` as `x → ∞`
#[inline]
pub fn parametric_decay(x: f32, params: &DecayParams) -> f32 {
    // Runtime guards: return 0.0 (fully decayed) for invalid / non-finite params
    // so that invalid configurations fail safe instead of silently producing NaN.
    if params.half_life <= 0.0
        || params.shape <= 0.0
        || !params.half_life.is_finite()
        || !params.shape.is_finite()
        || !x.is_finite()
    {
        return 0.0;
    }

    if x <= 0.0 {
        return 1.0;
    }

    // S(x; τ, β) = exp(-ln(2) · (x/τ)^β)
    let ratio = x / params.half_life;
    let exponent = -(core::f32::consts::LN_2) * ratio.powf(params.shape);

    // Clamp to [0, 1] to guard against floating-point edge cases
    exponent.exp().clamp(0.0, 1.0)
}

/// Convenience: standard exponential decay with the given half-life.
///
/// Equivalent to `parametric_decay(x, &DecayParams::exponential(half_life))`.
///
/// This is the direct replacement for the legacy `0.5_f32.powf(x / half_life)`.
#[inline]
pub fn exponential_decay(x: f32, half_life: f32) -> f32 {
    parametric_decay(x, &DecayParams::exponential(half_life))
}

/// Per-factor decay configuration.
///
/// Each scoring factor can use the unified parametric decay with its own
/// half-life and shape. Builders override these via `ScoringDecayConfig`
/// in the pack bundle manifest or environment variables.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct ScoringDecayConfig {
    /// Recency decay for notes (default: τ=7 days, β=1.0 exponential)
    pub recency_notes: DecayParams,

    /// Recency decay for artifacts (default: τ=30 days, β=1.0 exponential)
    pub recency_artifacts: DecayParams,

    /// Warmth (access frequency) decay (default: τ=7 days, β=1.0 exponential)
    pub warmth: DecayParams,

    /// Graph proximity decay over hop distance (default: τ=5 hops, β=1.0)
    ///
    /// Replaces the legacy linear decay `1 - d/D` with smooth exponential.
    pub graph_proximity: DecayParams,

    /// Causal proximity decay over event recency (default: τ=10 events, β=0.8 heavy-tail)
    pub causal_proximity: DecayParams,
}

impl Default for ScoringDecayConfig {
    fn default() -> Self {
        Self {
            recency_notes: DecayParams::exponential(7.0),
            recency_artifacts: DecayParams::exponential(30.0),
            warmth: DecayParams::exponential(7.0),
            graph_proximity: DecayParams::exponential(5.0),
            causal_proximity: DecayParams {
                half_life: 10.0,
                shape: 0.8,
            },
        }
    }
}

impl ScoringDecayConfig {
    /// Validate all decay parameters.
    pub fn validate(&self) -> bool {
        self.recency_notes.validate()
            && self.recency_artifacts.validate()
            && self.warmth.validate()
            && self.graph_proximity.validate()
            && self.causal_proximity.validate()
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    // ── Deterministic unit tests ──────────────────────────────────────

    #[test]
    fn test_identity_at_zero() {
        let params = DecayParams::exponential(7.0);
        assert_eq!(parametric_decay(0.0, &params), 1.0);
    }

    #[test]
    fn test_half_life_exponential() {
        let params = DecayParams::exponential(7.0);
        let result = parametric_decay(7.0, &params);
        assert!((result - 0.5).abs() < 1e-6, "Expected ~0.5, got {}", result);
    }

    #[test]
    fn test_half_life_large() {
        let params = DecayParams::exponential(365.0);
        let result = parametric_decay(365.0, &params);
        assert!((result - 0.5).abs() < 1e-5, "Expected ~0.5, got {}", result);
    }

    #[test]
    fn test_double_half_life() {
        let params = DecayParams::exponential(7.0);
        let result = parametric_decay(14.0, &params);
        assert!(
            (result - 0.25).abs() < 1e-5,
            "Expected ~0.25, got {}",
            result
        );
    }

    #[test]
    fn test_triple_half_life() {
        let params = DecayParams::exponential(7.0);
        let result = parametric_decay(21.0, &params);
        assert!(
            (result - 0.125).abs() < 1e-5,
            "Expected ~0.125, got {}",
            result
        );
    }

    #[test]
    fn test_asymptotic_very_large_x() {
        let params = DecayParams::exponential(7.0);
        let result = parametric_decay(10000.0, &params);
        assert!(result < 1e-10, "Expected ~0, got {}", result);
    }

    #[test]
    fn test_heavy_tail_decays_slower_than_exponential() {
        let exp = DecayParams::exponential(5.0);
        let heavy = DecayParams::heavy_tail(5.0);
        // At x > τ, heavy tail should be higher (decays slower)
        let x = 10.0;
        assert!(
            parametric_decay(x, &heavy) > parametric_decay(x, &exp),
            "Heavy tail should decay slower at x={}",
            x
        );
    }

    #[test]
    fn test_sharp_cutoff_decays_faster_past_half_life() {
        let exp = DecayParams::exponential(5.0);
        let sharp = DecayParams::sharp_cutoff(5.0);
        // At x > τ, sharp cutoff should be lower (decays faster)
        let x = 8.0;
        assert!(
            parametric_decay(x, &sharp) < parametric_decay(x, &exp),
            "Sharp cutoff should decay faster at x={}",
            x
        );
    }

    #[test]
    fn test_backward_compat_with_legacy_exponential() {
        // Legacy: 0.5_f32.powf(days / 7.0)
        let params = DecayParams::exponential(7.0);
        for days in [0.0, 1.0, 3.5, 7.0, 14.0, 30.0, 100.0] {
            let new_val = parametric_decay(days, &params);
            let old_val = 0.5_f32.powf(days / 7.0);
            assert!(
                (new_val - old_val).abs() < 1e-5,
                "Mismatch at days={}: new={}, old={}",
                days,
                new_val,
                old_val
            );
        }
    }

    #[test]
    fn test_backward_compat_artifact_recency() {
        // Legacy: 0.5_f32.powf(days / 30.0)
        let params = DecayParams::exponential(30.0);
        for days in [0.0, 7.0, 15.0, 30.0, 60.0, 90.0] {
            let new_val = parametric_decay(days, &params);
            let old_val = 0.5_f32.powf(days / 30.0);
            assert!(
                (new_val - old_val).abs() < 1e-5,
                "Mismatch at days={}: new={}, old={}",
                days,
                new_val,
                old_val
            );
        }
    }

    #[test]
    fn test_negative_x_returns_one() {
        let params = DecayParams::exponential(7.0);
        assert_eq!(parametric_decay(-5.0, &params), 1.0);
    }

    #[test]
    fn test_very_small_half_life() {
        let params = DecayParams::exponential(0.001);
        // At x=1.0, this should be essentially zero
        let result = parametric_decay(1.0, &params);
        assert!(
            result < 1e-10,
            "Expected ~0 for tiny half-life, got {}",
            result
        );
    }

    #[test]
    fn test_very_large_half_life() {
        let params = DecayParams::exponential(1_000_000.0);
        // At x=1.0, this should be essentially 1.0
        let result = parametric_decay(1.0, &params);
        assert!(
            (result - 1.0).abs() < 0.001,
            "Expected ~1.0 for huge half-life, got {}",
            result
        );
    }

    #[test]
    fn test_convenience_exponential_decay() {
        let direct = parametric_decay(5.0, &DecayParams::exponential(7.0));
        let convenience = exponential_decay(5.0, 7.0);
        assert_eq!(direct, convenience);
    }

    #[test]
    fn test_decay_params_validate() {
        assert!(DecayParams::exponential(7.0).validate());
        assert!(DecayParams::heavy_tail(5.0).validate());
        assert!(DecayParams::sharp_cutoff(3.0).validate());
        assert!(!DecayParams {
            half_life: 0.0,
            shape: 1.0
        }
        .validate());
        assert!(!DecayParams {
            half_life: -1.0,
            shape: 1.0
        }
        .validate());
        assert!(!DecayParams {
            half_life: 1.0,
            shape: 0.0
        }
        .validate());
        assert!(!DecayParams {
            half_life: 1.0,
            shape: -1.0
        }
        .validate());
        assert!(!DecayParams {
            half_life: f32::NAN,
            shape: 1.0
        }
        .validate());
        assert!(!DecayParams {
            half_life: 1.0,
            shape: f32::INFINITY
        }
        .validate());
    }

    #[test]
    fn test_scoring_decay_config_defaults_valid() {
        let config = ScoringDecayConfig::default();
        assert!(config.validate());
    }

    // ── NaN / Infinity production guard tests ────────────────────────

    #[test]
    fn test_nan_half_life_returns_zero() {
        let params = DecayParams {
            half_life: f32::NAN,
            shape: 1.0,
        };
        assert_eq!(parametric_decay(5.0, &params), 0.0);
    }

    #[test]
    fn test_nan_shape_returns_zero() {
        let params = DecayParams {
            half_life: 7.0,
            shape: f32::NAN,
        };
        assert_eq!(parametric_decay(5.0, &params), 0.0);
    }

    #[test]
    fn test_nan_age_returns_zero() {
        let params = DecayParams::exponential(7.0);
        assert_eq!(parametric_decay(f32::NAN, &params), 0.0);
    }

    #[test]
    fn test_infinity_half_life_returns_zero() {
        let params = DecayParams {
            half_life: f32::INFINITY,
            shape: 1.0,
        };
        assert_eq!(parametric_decay(5.0, &params), 0.0);
    }

    #[test]
    fn test_neg_infinity_half_life_returns_zero() {
        let params = DecayParams {
            half_life: f32::NEG_INFINITY,
            shape: 1.0,
        };
        assert_eq!(parametric_decay(5.0, &params), 0.0);
    }

    #[test]
    fn test_infinity_shape_returns_zero() {
        let params = DecayParams {
            half_life: 7.0,
            shape: f32::INFINITY,
        };
        assert_eq!(parametric_decay(5.0, &params), 0.0);
    }

    #[test]
    fn test_infinity_age_returns_zero() {
        let params = DecayParams::exponential(7.0);
        assert_eq!(parametric_decay(f32::INFINITY, &params), 0.0);
    }

    #[test]
    fn test_negative_shape_returns_zero() {
        let params = DecayParams {
            half_life: 7.0,
            shape: -1.0,
        };
        assert_eq!(parametric_decay(5.0, &params), 0.0);
    }

    #[test]
    fn test_zero_half_life_returns_zero() {
        let params = DecayParams {
            half_life: 0.0,
            shape: 1.0,
        };
        assert_eq!(parametric_decay(5.0, &params), 0.0);
    }

    // ── Property tests (proptest) ─────────────────────────────────────

    mod prop {
        use super::*;
        use proptest::prelude::*;

        proptest! {
            #![proptest_config(ProptestConfig::with_cases(2000))]

            /// P1: S(0; τ, β) = 1.0 for all valid τ, β
            #[test]
            fn prop_identity_at_zero(
                half_life in 0.001f32..10000.0,
                shape in 0.01f32..10.0,
            ) {
                let params = DecayParams { half_life, shape };
                let result = parametric_decay(0.0, &params);
                prop_assert!(
                    (result - 1.0).abs() < 1e-6,
                    "P1 violated: decay(0, τ={}, β={}) = {} ≠ 1.0",
                    half_life, shape, result
                );
            }

            /// P2: S(τ; τ, 1) = 0.5 (half-life semantics for exponential shape)
            #[test]
            fn prop_half_life_semantics(
                half_life in 0.01f32..10000.0,
            ) {
                let params = DecayParams { half_life, shape: 1.0 };
                let result = parametric_decay(half_life, &params);
                prop_assert!(
                    (result - 0.5).abs() < 1e-4,
                    "P2 violated: decay(τ={}, τ={}, β=1) = {} ≠ 0.5",
                    half_life, half_life, result
                );
            }

            /// P3: S(x₁) ≥ S(x₂) when x₁ ≤ x₂ (monotone decreasing)
            #[test]
            fn prop_monotone_decreasing(
                half_life in 0.01f32..10000.0,
                shape in 0.01f32..10.0,
                x1 in 0.0f32..10000.0,
                delta in 0.001f32..1000.0,
            ) {
                let params = DecayParams { half_life, shape };
                let x2 = x1 + delta;
                let s1 = parametric_decay(x1, &params);
                let s2 = parametric_decay(x2, &params);
                prop_assert!(
                    s1 >= s2,
                    "P3 violated: S({}) = {} < S({}) = {} with τ={}, β={}",
                    x1, s1, x2, s2, half_life, shape
                );
            }

            /// P4: 0 ≤ S(x) ≤ 1 for all x ≥ 0
            #[test]
            fn prop_bounded_unit_interval(
                half_life in 0.001f32..10000.0,
                shape in 0.01f32..10.0,
                x in 0.0f32..100000.0,
            ) {
                let params = DecayParams { half_life, shape };
                let result = parametric_decay(x, &params);
                prop_assert!(
                    (0.0..=1.0).contains(&result),
                    "P4 violated: S({}, τ={}, β={}) = {} not in [0, 1]",
                    x, half_life, shape, result
                );
            }

            /// P5: S(k·τ) ≈ 0 for sufficiently large k (asymptotic zero)
            ///
            /// The required multiplier depends on shape β. For heavy tails (β < 1),
            /// convergence to zero is slow: we need k^β > ln(1/ε)/ln(2).
            /// We restrict β ≥ 0.3 here so k=1000 suffices for ε=0.01.
            #[test]
            fn prop_asymptotic_zero(
                half_life in 0.01f32..100.0,
                shape in 0.3f32..10.0,
            ) {
                let params = DecayParams { half_life, shape };
                let result = parametric_decay(1000.0 * half_life, &params);
                prop_assert!(
                    result < 0.01,
                    "P5 violated: S(1000τ, τ={}, β={}) = {} ≥ 0.01",
                    half_life, shape, result
                );
            }

            /// Subsumption: parametric_decay(x, {τ, 1.0}) ≡ 0.5^(x/τ)
            #[test]
            fn prop_backward_compat_exponential(
                half_life in 0.01f32..10000.0,
                x in 0.0f32..10000.0,
            ) {
                let params = DecayParams { half_life, shape: 1.0 };
                let new_result = parametric_decay(x, &params);
                let old_result = 0.5_f32.powf(x / half_life);
                prop_assert!(
                    (new_result - old_result).abs() < 1e-4,
                    "Subsumption violated at x={}, τ={}: new={} vs old={}",
                    x, half_life, new_result, old_result
                );
            }

            /// Shape ordering: at x > τ, lower β (heavier tail) yields higher score
            #[test]
            fn prop_shape_ordering_past_half_life(
                half_life in 1.0f32..100.0,
                beta_low in 0.1f32..0.9,
                beta_high in 1.1f32..5.0,
                x_ratio in 1.5f32..10.0,  // x = x_ratio * τ, so x > τ
            ) {
                let x = x_ratio * half_life;
                let heavy = DecayParams { half_life, shape: beta_low };
                let sharp = DecayParams { half_life, shape: beta_high };
                let s_heavy = parametric_decay(x, &heavy);
                let s_sharp = parametric_decay(x, &sharp);
                prop_assert!(
                    s_heavy >= s_sharp,
                    "Shape ordering violated at x={}: heavy(β={})={} < sharp(β={})={}",
                    x, beta_low, s_heavy, beta_high, s_sharp
                );
            }

            /// DecayParams::validate rejects all invalid configurations
            #[test]
            fn prop_validate_rejects_nonpositive(
                half_life in -1000.0f32..0.0,
                shape in 0.01f32..10.0,
            ) {
                let params = DecayParams { half_life, shape };
                prop_assert!(!params.validate());
            }
        }
    }

    // ── Moved from integration tests (scoring_property_tests.rs) ─────

    #[test]
    fn test_decay_exponential_at_half_life_integration() {
        let params = DecayParams::exponential(7.0);
        let d = parametric_decay(7.0, &params);
        assert!(
            (d - 0.5).abs() < 1e-4,
            "exponential decay at half_life should be 0.5, got {}",
            d
        );
    }

    #[test]
    fn test_decay_weibull_shape_effect() {
        // Higher shape = sharper cliff
        let params_gentle = DecayParams {
            half_life: 10.0,
            shape: 0.5,
        };
        let params_sharp = DecayParams {
            half_life: 10.0,
            shape: 3.0,
        };

        // At x < half_life, gentle decays faster (concave up vs sigmoid)
        let gentle_at_3 = parametric_decay(3.0, &params_gentle);
        let sharp_at_3 = parametric_decay(3.0, &params_sharp);
        // Both should be > 0.5 since x < half_life
        assert!(gentle_at_3 > 0.5);
        assert!(sharp_at_3 > 0.5);
    }

    mod prop_integration {
        use super::super::*;
        use proptest::prelude::*;

        proptest! {
            #![proptest_config(ProptestConfig::with_cases(500))]

            /// Property: graph proximity -- deeper depth produces lower or equal score.
            #[test]
            fn prop_graph_proximity_monotone(d1 in 0u32..10, d2 in 0u32..10) {
                let params = DecayParams::exponential(5.0);
                let s1 = parametric_decay(d1 as f32, &params);
                let s2 = parametric_decay(d2 as f32, &params);
                if d1 <= d2 {
                    prop_assert!(s1 >= s2 - 1e-6,
                        "depth {} score {} should be >= depth {} score {}",
                        d1, s1, d2, s2);
                }
            }

            /// Property: causal proximity -- events_ago=0 yields max score.
            #[test]
            fn prop_causal_entry_decay(events_ago in 0u32..200) {
                let params = DecayParams { half_life: 10.0, shape: 0.8 };
                let score_at_0 = parametric_decay(0.0, &params);
                let score_at_n = parametric_decay(events_ago as f32, &params);
                prop_assert!(score_at_0 >= score_at_n - 1e-6,
                    "events_ago=0 score {} should be >= events_ago={} score {}",
                    score_at_0, events_ago, score_at_n);
            }
        }
    }
}