feat(information,save): KL/JS divergence + save overhaul

- Fix digamma Stirling expansion; KLDivergence with JS, TV, cross-entropy + tests
- Discrete + continuous estimator traits; KnnEntropy/MI (default k=3)
- PersistentId(u64), seeded counter, atomic write fix, VecDeque events
- Chained version migration; remove dead helpers; basic_save.rs updated

Author: M1thieu

crates/information/src/measures/divergence.rs

@@ -1,73 +1,185 @@
-/// Kullback-Leibler divergence and related measures
-/// Core information-theoretic distance metrics
+/// Kullback-Leibler divergence and related measures.
+/// Core information-theoretic distance metrics for comparing probability distributions.
+///
+/// LP usage: measure evolutionary distance between species, compare creature behaviour
+/// distributions, quantify trait drift between populations. In multiplayer, these
+/// metrics let the simulation decide when two players' lineages count as separate species.
 pub struct KLDivergence;
 
 impl KLDivergence {
-    /// Calculate KL divergence D(P||Q) = Σ P(i) log₂(P(i)/Q(i))
-    /// Measures how distribution P differs from reference Q
-    /// Returns divergence in bits - NOT symmetric
+    /// D(P||Q) = Σ P(i) log₂(P(i)/Q(i))
+    ///
+    /// How many extra bits are needed to encode P-distributed events using a code
+    /// optimised for Q. NOT symmetric. Returns +∞ when P has mass where Q has none.
     pub fn divergence(p_probs: &[f64], q_probs: &[f64]) -> f64 {
         assert_eq!(
             p_probs.len(),
             q_probs.len(),
             "P and Q must have same length"
         );
 
-        let mut kl_div = 0.0;
+        let mut kl = 0.0;
         for (&p, &q) in p_probs.iter().zip(q_probs) {
             if p > 0.0 && q > 0.0 {
-                kl_div += p * (p / q).log2();
-            } else if p > 0.0 && q == 0.0 {
-                // P has probability where Q doesn't - infinite divergence
+                kl += p * (p / q).log2();
+            } else if p > 0.0 {
+                // P has support where Q does not — infinite divergence by definition
                 return f64::INFINITY;
             }
-            // p == 0.0 contributes nothing to KL divergence
+            // p == 0 contributes 0 regardless of q  (0 · log(0/q) := 0)
         }
-
-        kl_div
+        kl
     }
 
-    /// Calculate Jensen-Shannon divergence - symmetric version of KL
-    /// JS(P,Q) = 0.5 * [D(P||M) + D(Q||M)] where M = 0.5*(P+Q)
-    /// Always finite and bounded [0, 1] bits
+    /// JS(P,Q) = 0.5 · [D(P||M) + D(Q||M)]  where M = 0.5·(P+Q)
+    ///
+    /// Symmetric, always finite, bounded [0, 1] bits (log₂ base).
+    /// Reaches 1 bit only when P and Q have fully disjoint support.
+    /// Preferred over raw KL for comparing creature trait distributions
+    /// because it never blows up and is a proper metric (square root is a distance).
     pub fn jensen_shannon(p_probs: &[f64], q_probs: &[f64]) -> f64 {
         assert_eq!(
             p_probs.len(),
             q_probs.len(),
             "P and Q must have same length"
         );
 
-        // Calculate mixture distribution M = 0.5*(P+Q)
-        let m_probs: Vec<f64> = p_probs
+        let m: Vec<f64> = p_probs
             .iter()
             .zip(q_probs)
             .map(|(&p, &q)| 0.5 * (p + q))
             .collect();
 
-        let kl_pm = Self::divergence(p_probs, &m_probs);
-        let kl_qm = Self::divergence(q_probs, &m_probs);
-
-        0.5 * (kl_pm + kl_qm)
+        // M[i] >= 0.5 · max(P[i], Q[i]), so KL(P||M) and KL(Q||M) are always finite.
+        0.5 * (Self::divergence(p_probs, &m) + Self::divergence(q_probs, &m))
     }
 
-    /// Calculate cross-entropy H(P,Q) = -Σ P(i) log₂(Q(i))  
-    /// Useful for ML loss functions and distribution comparison
+    /// H(P,Q) = -Σ P(i) log₂(Q(i))
+    ///
+    /// Cross-entropy: bits to encode P-distributed events with a Q-optimal code.
+    /// H(P,Q) = H(P) + D(P||Q). Returns +∞ when Q has no mass where P has mass.
     pub fn cross_entropy(p_probs: &[f64], q_probs: &[f64]) -> f64 {
         assert_eq!(
             p_probs.len(),
             q_probs.len(),
             "P and Q must have same length"
         );
 
-        let mut cross_ent = 0.0;
+        let mut ce = 0.0;
         for (&p, &q) in p_probs.iter().zip(q_probs) {
             if p > 0.0 && q > 0.0 {
-                cross_ent -= p * q.log2();
-            } else if p > 0.0 && q == 0.0 {
+                ce -= p * q.log2();
+            } else if p > 0.0 {
                 return f64::INFINITY;
             }
         }
+        ce
+    }
+
+    /// TV(P,Q) = 0.5 · Σ|P(i) - Q(i)|
+    ///
+    /// Total variation distance: the maximum probability gap any single event
+    /// can have between P and Q. Symmetric, bounded [0, 1], no log required.
+    /// Fastest way to ask "how different are these two distributions?" when you
+    /// do not need the information-theoretic interpretation of KL/JS.
+    pub fn total_variation(p_probs: &[f64], q_probs: &[f64]) -> f64 {
+        assert_eq!(
+            p_probs.len(),
+            q_probs.len(),
+            "P and Q must have same length"
+        );
+        0.5 * p_probs
+            .iter()
+            .zip(q_probs)
+            .map(|(&p, &q)| (p - q).abs())
+            .sum::<f64>()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn kl_identical_is_zero() {
+        let p = [0.5, 0.5];
+        assert!(KLDivergence::divergence(&p, &p).abs() < 1e-10);
+    }
+
+    #[test]
+    fn kl_disjoint_support_is_infinity() {
+        let p = [1.0, 0.0];
+        let q = [0.0, 1.0];
+        assert_eq!(KLDivergence::divergence(&p, &q), f64::INFINITY);
+    }
+
+    #[test]
+    fn kl_known_value() {
+        // D([0.5, 0.5] || [0.25, 0.75]) = 0.5*log2(2) + 0.5*log2(2/3) ≈ 0.2075
+        let p = [0.5, 0.5];
+        let q = [0.25, 0.75];
+        let kl = KLDivergence::divergence(&p, &q);
+        let expected = 0.5 * (0.5_f64 / 0.25).log2() + 0.5 * (0.5_f64 / 0.75).log2();
+        assert!((kl - expected).abs() < 1e-10, "got {}", kl);
+    }
+
+    #[test]
+    fn js_identical_is_zero() {
+        let p = [0.25, 0.25, 0.25, 0.25];
+        assert!(KLDivergence::jensen_shannon(&p, &p).abs() < 1e-10);
+    }
+
+    #[test]
+    fn js_disjoint_support_is_one_bit() {
+        let p = [1.0, 0.0];
+        let q = [0.0, 1.0];
+        let js = KLDivergence::jensen_shannon(&p, &q);
+        assert!((js - 1.0).abs() < 1e-10, "expected 1 bit, got {}", js);
+    }
+
+    #[test]
+    fn js_is_symmetric() {
+        let p = [0.7, 0.2, 0.1];
+        let q = [0.1, 0.5, 0.4];
+        let diff =
+            (KLDivergence::jensen_shannon(&p, &q) - KLDivergence::jensen_shannon(&q, &p)).abs();
+        assert!(diff < 1e-12, "JS must be symmetric, diff = {}", diff);
+    }
+
+    #[test]
+    fn js_bounded_between_zero_and_one() {
+        let p = [0.6, 0.3, 0.1];
+        let q = [0.1, 0.2, 0.7];
+        let js = KLDivergence::jensen_shannon(&p, &q);
+        assert!(js >= 0.0 && js <= 1.0, "JS out of [0,1]: {}", js);
+    }
+
+    #[test]
+    fn cross_entropy_of_self_equals_entropy() {
+        // H(P, P) = H(P). For fair coin P=[0.5,0.5], H = 1 bit.
+        let p = [0.5, 0.5];
+        let h = KLDivergence::cross_entropy(&p, &p);
+        assert!((h - 1.0).abs() < 1e-10, "expected 1 bit, got {}", h);
+    }
+
+    #[test]
+    fn total_variation_identical_is_zero() {
+        let p = [0.3, 0.3, 0.4];
+        assert!(KLDivergence::total_variation(&p, &p).abs() < 1e-12);
+    }
+
+    #[test]
+    fn total_variation_disjoint_is_one() {
+        let p = [1.0, 0.0];
+        let q = [0.0, 1.0];
+        assert!((KLDivergence::total_variation(&p, &q) - 1.0).abs() < 1e-12);
+    }
 
-        cross_ent
+    #[test]
+    fn total_variation_bounded() {
+        let p = [0.6, 0.4];
+        let q = [0.2, 0.8];
+        let tv = KLDivergence::total_variation(&p, &q);
+        assert!(tv >= 0.0 && tv <= 1.0, "TV must be in [0, 1], got {}", tv);
     }
 }

crates/information/src/measures/estimators.rs

@@ -1,17 +1,42 @@
 use super::mutual::MutualInfo;
 use super::shannon::Shannon;
 
-/// Trait for discrete entropy estimators.
+// ── Discrete estimator traits ─────────────────────────────────────────────────
+
+/// Compute entropy from a sequence of discrete integer labels (in bits).
 pub trait DiscreteEntropyEstimator {
     fn estimate_entropy_bits(&self, values: &[i32]) -> f64;
 }
 
-/// Trait for discrete mutual-information estimators.
+/// Compute mutual information between two discrete integer sequences (in bits).
 pub trait DiscreteMutualInformationEstimator {
     fn estimate_mutual_information_bits(&self, x_values: &[i32], y_values: &[i32]) -> f64;
 }
 
-/// Baseline Shannon entropy estimator.
+// ── Continuous estimator traits ───────────────────────────────────────────────
+//
+// Continuous data is represented as &[Vec<f32>] (N points × D dimensions),
+// matching the kNN-based API in shannon.rs and mutual/calculation.rs.
+//
+// LP / multiplayer usage: plug in the right estimator per signal type without
+// changing call sites. Discrete handles genome symbols and material-type counts;
+// continuous handles position density fields, velocity distributions, and any
+// other physical observable sampled from the MPM grid.
+
+/// Compute differential entropy from a continuous multivariate sample (in bits).
+pub trait ContinuousEntropyEstimator {
+    fn estimate_entropy_bits(&self, data: &[Vec<f32>]) -> f64;
+}
+
+/// Compute mutual information between two continuous multivariate samples (in bits).
+pub trait ContinuousMutualInformationEstimator {
+    fn estimate_mutual_information_bits(&self, x: &[Vec<f32>], y: &[Vec<f32>]) -> f64;
+}
+
+// ── Discrete implementations ──────────────────────────────────────────────────
+
+/// Baseline maximum-likelihood Shannon entropy estimator for discrete data.
+/// Biased downward for small samples; prefer Miller-Madow correction for n < 100.
 #[derive(Debug, Default, Clone, Copy)]
 pub struct ShannonEstimator;
 
@@ -21,7 +46,7 @@ impl DiscreteEntropyEstimator for ShannonEstimator {
     }
 }
 
-/// Baseline empirical mutual-information estimator.
+/// Baseline empirical mutual-information estimator for discrete data.
 #[derive(Debug, Default, Clone, Copy)]
 pub struct EmpiricalMutualInformationEstimator;
 
@@ -31,10 +56,52 @@ impl DiscreteMutualInformationEstimator for EmpiricalMutualInformationEstimator
     }
 }
 
+// ── Continuous implementations (k-NN based) ───────────────────────────────────
+
+/// Kraskov et al. 2004 k-NN entropy estimator for continuous multivariate data.
+/// `k` controls the neighbourhood size; k=3..5 works well for most LP use cases.
+#[derive(Debug, Clone, Copy)]
+pub struct KnnEntropyEstimator {
+    pub k: usize,
+}
+
+impl Default for KnnEntropyEstimator {
+    fn default() -> Self {
+        Self { k: 3 }
+    }
+}
+
+impl ContinuousEntropyEstimator for KnnEntropyEstimator {
+    fn estimate_entropy_bits(&self, data: &[Vec<f32>]) -> f64 {
+        Shannon::continuous_entropy(data, self.k)
+    }
+}
+
+/// Kraskov et al. 2004 k-NN mutual information estimator for continuous data.
+/// `k` controls the neighbourhood size; k=3..5 works well for most LP use cases.
+#[derive(Debug, Clone, Copy)]
+pub struct KnnMutualInformationEstimator {
+    pub k: usize,
+}
+
+impl Default for KnnMutualInformationEstimator {
+    fn default() -> Self {
+        Self { k: 3 }
+    }
+}
+
+impl ContinuousMutualInformationEstimator for KnnMutualInformationEstimator {
+    fn estimate_mutual_information_bits(&self, x: &[Vec<f32>], y: &[Vec<f32>]) -> f64 {
+        MutualInfo::continuous_knn(x, y, self.k)
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
 
+    // ── discrete ─────────────────────────────────────────────────────────────
+
     #[test]
     fn shannon_entropy_for_fair_binary_is_one_bit() {
         let estimator = ShannonEstimator;
@@ -68,4 +135,33 @@ mod tests {
         let mi = estimator.estimate_mutual_information_bits(&x, &y);
         assert!(mi.abs() < 1e-9, "expected near 0, got {}", mi);
     }
+
+    // ── continuous ────────────────────────────────────────────────────────────
+
+    #[test]
+    fn knn_entropy_estimator_positive_for_spread_data() {
+        let estimator = KnnEntropyEstimator::default();
+        // 20 points spread across 1D — entropy should be clearly positive
+        let data: Vec<Vec<f32>> = (0..20).map(|i| vec![i as f32 * 0.5]).collect();
+        let h = estimator.estimate_entropy_bits(&data);
+        assert!(
+            h > 0.0,
+            "expected positive entropy for spread data, got {}",
+            h
+        );
+    }
+
+    #[test]
+    fn knn_mi_estimator_nonzero_for_correlated_data() {
+        let estimator = KnnMutualInformationEstimator::default();
+        // x and y perfectly correlated → MI should be positive
+        let x: Vec<Vec<f32>> = (0..20).map(|i| vec![i as f32]).collect();
+        let y: Vec<Vec<f32>> = (0..20).map(|i| vec![i as f32 * 2.0]).collect();
+        let mi = estimator.estimate_mutual_information_bits(&x, &y);
+        assert!(
+            mi > 0.0,
+            "expected positive MI for correlated data, got {}",
+            mi
+        );
+    }
 }