arco 0.3.0

//! Threshold calibration for emergence metrics.
//!
//! Per the Mathematical Constitution:
//!     Emergence thresholds must be calibrated per universe class
//!     using the 95th percentile of a null distribution. This ensures
//!     that "emergence" means "statistically distinguishable from
//!     noise."
//!
//! # Calibration procedure
//!
//! 1. Generate null trajectories using purely destructive rule sets.
//! 2. Compute the emergence metric for each null universe.
//! 3. Set the threshold to the specified percentile of the null
//!    distribution.
//! 4. Apply engineering floors to prevent degenerate thresholds.
//!
//! # Design
//!
//! Calibration is generic over any [`InformationUniverse`] type.
//! The null distribution is generated by the universe itself,
//! which must provide a way to create destructive rule sets.
//! This keeps calibration substrate-independent — each universe
//! defines what "destructive" means for its substrate.

use crate::metrics::{compute_memory, compute_persistence, compute_storage};
use crate::observation::Observation;
use crate::types::NullEnsembles;
use crate::universe::InformationUniverse;
use rand::RngExt;
use rand::SeedableRng;
use rand::rngs::StdRng;

/// Configuration for threshold calibration.
pub struct CalibrationConfig {
    pub percentile: f64,
    pub floor_persistence: f64,
    pub floor_storage: f64,
    pub floor_memory: f64,
    pub max_delta: usize,
    pub n_shuffles: usize,
    pub seed: u64,
}

impl Default for CalibrationConfig {
    fn default() -> Self {
        Self {
            percentile: 95.0,
            floor_persistence: 0.01,
            floor_storage: 0.01,
            floor_memory: 0.01,
            max_delta: 15,
            n_shuffles: 10,
            seed: 42,
        }
    }
}

// ===================================================================
// Null distribution generation
// ===================================================================

/// Generate trajectory ensembles for null-distribution calibration.
///
/// Each null universe uses a destructive rule set provided by the
/// universe's `null_rules()` method. Trajectories are generated
/// from random initial states sampled from the universe's state
/// space.
///
/// # Parameters
/// * `universe` — The universe to sample from. Provides state space,
///   observation operator, and schedule.
/// * `n_universes` — Number of null universes (≥ 30 recommended).
/// * `n_ensemble` — Ensemble size per universe.
/// * `steps` — Timesteps per trajectory.
/// * `seed` — Random seed for reproducibility.
///
/// # Returns
/// Vector of null trajectory ensembles, each a vector of trajectories,
/// each a vector of observation values.
pub fn generate_null_trajectories<U: InformationUniverse>(
    universe: &U,
    n_universes: usize,
    n_ensemble: usize,
    steps: usize,
    seed: u64,
) -> NullEnsembles<U>
where
    <U::Observation as Observation<U::State>>::Output: Clone,
{
    let mut rng = StdRng::seed_from_u64(seed);
    let state_space = universe.state_space();
    let schedule = universe.schedule();
    let observer = universe.observation();

    let mut null_ensembles = Vec::with_capacity(n_universes);

    for i in 0..n_universes {
        // Get a destructive rule set from the universe
        let rules = universe.null_rules(&mut rng);

        // Select random initial states
        let n_pool = state_space.len();
        let n_ensemble = n_ensemble.min(n_pool);
        let mut init_indices: Vec<usize> = (0..n_pool).collect();
        for j in 0..n_ensemble {
            let k = rng.random_range(j..n_pool);
            init_indices.swap(j, k);
        }
        let initial_states: Vec<U::State> = init_indices
            .iter()
            .take(n_ensemble)
            .map(|&idx| state_space[idx].clone())
            .collect();

        // Generate ensemble
        let ensemble = generate_trajectories(
            &initial_states,
            &rules,
            steps,
            schedule,
            observer,
            seed + i as u64 * 1000,
        );

        null_ensembles.push(ensemble);
    }

    null_ensembles
}

// ===================================================================
// Trajectory generation (generic)
// ===================================================================

/// Generate an ensemble of trajectories from distinct initial states.
///
/// This is the generic trajectory generator used by both calibration
/// and the scientific cycle. It is substrate-independent — it works
/// with any `State`, `Rule`, `Schedule`, and `Observation` types.
pub fn generate_trajectories<S, R, K, O>(
    initial_states: &[S],
    rules: &[R],
    steps: usize,
    schedule: &K,
    observer: &O,
    base_seed: u64,
) -> Vec<Vec<O::Output>>
where
    S: crate::state::State,
    R: crate::rules::Rule<S>,
    K: crate::schedule::Schedule<S, R>,
    O: crate::observation::Observation<S>,
    O::Output: Clone,
{
    let mut trajectories = Vec::with_capacity(initial_states.len());

    for (i, initial) in initial_states.iter().enumerate() {
        let seed = base_seed + i as u64 * 137;
        let mut rng = StdRng::seed_from_u64(seed);
        let mut observations = Vec::with_capacity(steps + 1);
        observations.push(observer.observe(initial));
        let mut current = initial.clone();

        for _ in 0..steps {
            current = schedule.step(&current, rules, &mut rng);
            observations.push(observer.observe(&current));
        }

        trajectories.push(observations);
    }

    trajectories
}

// ===================================================================
// Threshold calibration
// ===================================================================

/// Calibration result with thresholds and null distribution statistics.
#[derive(Debug, Clone)]
pub struct CalibrationResult {
    /// Calibrated persistence threshold.
    pub persistence_threshold: f64,
    /// Calibrated storage threshold.
    pub storage_threshold: f64,
    /// Calibrated memory threshold.
    pub memory_threshold: f64,
    /// Null distribution statistics for persistence.
    pub null_persistence: NullStats,
    /// Null distribution statistics for storage.
    pub null_storage: NullStats,
    /// Null distribution statistics for memory.
    pub null_memory: NullStats,
    /// The percentile used for thresholding.
    pub percentile: f64,
}

/// Statistics from a null distribution.
#[derive(Debug, Clone)]
pub struct NullStats {
    /// Mean of the null distribution.
    pub mean: f64,
    /// Standard deviation of the null distribution.
    pub std: f64,
    /// Raw scores from null universes.
    pub scores: Vec<f64>,
}

impl NullStats {
    /// Empirical p-value: fraction of null scores ≥ observed_value.
    pub fn empirical_p(&self, observed_value: f64) -> f64 {
        let count = self.scores.iter().filter(|&&s| s >= observed_value).count();
        if self.scores.is_empty() {
            1.0
        } else {
            count as f64 / self.scores.len() as f64
        }
    }
}

/// Calibrate emergence thresholds from null ensembles.
///
/// # Parameters
/// * `null_ensembles` — Null trajectory ensembles.
/// * `percentile` — Percentile for threshold (0–100).
/// * `floor_persistence` — Engineering floor for persistence.
/// * `floor_storage` — Engineering floor for storage.
/// * `floor_memory` — Engineering floor for memory.
/// * `max_delta` — Maximum timescale for storage/memory.
/// * `n_shuffles` — Number of shuffles for bias correction.
/// * `seed` — Seed for shuffle RNG.
pub fn calibrate_thresholds<T: Eq + std::hash::Hash + Clone>(
    null_ensembles: &[Vec<Vec<T>>],
    config: &CalibrationConfig,
) -> CalibrationResult {
    let mut persistence_scores = Vec::with_capacity(null_ensembles.len());
    let mut storage_scores = Vec::with_capacity(null_ensembles.len());
    let mut memory_scores = Vec::with_capacity(null_ensembles.len());

    for ensemble in null_ensembles {
        persistence_scores.push(compute_persistence(
            ensemble,
            1,
            config.n_shuffles,
            config.seed,
        ));
        storage_scores.push(compute_storage(
            ensemble,
            config.max_delta,
            config.n_shuffles,
            config.seed,
        ));
        memory_scores.push(compute_memory(
            ensemble,
            config.max_delta,
            config.n_shuffles,
            config.seed,
        ));
    }

    let persistence_threshold =
        percentile_value(&persistence_scores, config.percentile).max(config.floor_persistence);
    let storage_threshold =
        percentile_value(&storage_scores, config.percentile).max(config.floor_storage);
    let memory_threshold =
        percentile_value(&memory_scores, config.percentile).max(config.floor_memory);

    CalibrationResult {
        persistence_threshold,
        storage_threshold,
        memory_threshold,
        null_persistence: NullStats::from_scores(persistence_scores),
        null_storage: NullStats::from_scores(storage_scores),
        null_memory: NullStats::from_scores(memory_scores),
        percentile: config.percentile,
    }
}

/// Compute a percentile value with linear interpolation.
fn percentile_value(scores: &[f64], percentile: f64) -> f64 {
    if scores.is_empty() {
        return 0.0;
    }
    if scores.len() == 1 {
        return scores[0];
    }
    let mut sorted: Vec<f64> = scores.to_vec();
    sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));

    let k = (percentile / 100.0) * (sorted.len() - 1) as f64;
    let lo = k.floor() as usize;
    let hi = k.ceil() as usize;

    if lo == hi {
        sorted[lo]
    } else {
        let frac = k - lo as f64;
        sorted[lo] * (1.0 - frac) + sorted[hi] * frac
    }
}

impl NullStats {
    fn from_scores(scores: Vec<f64>) -> Self {
        let n = scores.len() as f64;
        let mean = if n > 0.0 {
            scores.iter().sum::<f64>() / n
        } else {
            0.0
        };
        let std = if n > 1.0 {
            let variance = scores.iter().map(|s| (s - mean).powi(2)).sum::<f64>() / (n - 1.0);
            variance.sqrt()
        } else {
            0.0
        };
        Self { mean, std, scores }
    }
}

// ===================================================================
// Convenience: run full calibration pipeline
// ===================================================================

/// Run the full calibration pipeline and return thresholds.
///
/// Convenience wrapper around `generate_null_trajectories` and
/// `calibrate_thresholds`.
pub fn calibrate<U: InformationUniverse>(
    universe: &U,
    n_null_universes: usize,
    n_ensemble: usize,
    steps: usize,
    config: &CalibrationConfig,
) -> CalibrationResult
where
    <U::Observation as Observation<U::State>>::Output: Eq + std::hash::Hash + Clone,
{
    let null_ensembles =
        generate_null_trajectories(universe, n_null_universes, n_ensemble, steps, config.seed);

    calibrate_thresholds(&null_ensembles, config)
}