dsfb_dscd/

paper.rs

//! DSCD paper-facing deterministic simulation and export pipeline.
//!
//! This module layers Deterministic Structural Causal Dynamics (DSCD) on top
//! of DSFB + ADD. The edge rule is explicitly local and deterministic:
//! temporal span bounds, structural radius bounds, and trust decay over
//! temporal/structural distance. This supports certification-friendly graph
//! complexity and finite-size threshold analysis in the DSCD paper.

use std::collections::{BTreeMap, BTreeSet, HashMap, HashSet, VecDeque};
use std::fs;
use std::path::{Path, PathBuf};

use anyhow::{anyhow, bail, Context, Result};
use csv::Writer;
use dsfb::sim::SimConfig;
use dsfb::DsfbParams;
use dsfb_add::aet::run_aet_sweep;
use dsfb_add::iwlt::run_iwlt_sweep;
use dsfb_add::SimulationConfig as AddSimulationConfig;
use rayon::prelude::*;
use serde::{Deserialize, Serialize};

use crate::integrations::{generate_dscd_events_from_dsfb, DscdObserverSample, ResidualState};

/// DSCD locality and threshold parameters.
///
/// These parameters are deterministic and directly control the local-causal
/// edge rule used for DSCD graph construction:
/// - Locality (temporal + structural) supports interval finiteness.
/// - Distance-based trust decay supports nontrivial threshold transition shape.
/// - Out-degree bound keeps complexity certification-friendly.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct DscdParams {
    pub num_events: usize,
    pub tau_min: f64,
    pub tau_max: f64,
    pub tau_steps: usize,
    pub max_temporal_lag: usize,
    pub max_structural_radius: f64,
    pub temporal_decay_gamma: f64,
    pub structural_decay_beta: f64,
    pub max_out_edges: usize,
}

impl DscdParams {
    /// Demo preset tuned for free Colab runtimes.
    pub fn demo() -> Self {
        Self {
            num_events: 10_000,
            tau_min: 0.0,
            tau_max: 1.0,
            tau_steps: 64,
            max_temporal_lag: 64,
            max_structural_radius: 0.1,
            temporal_decay_gamma: 0.0,
            structural_decay_beta: 2.0,
            max_out_edges: 8,
        }
    }

    /// Performance preset for workstation runs.
    pub fn performance() -> Self {
        Self {
            num_events: 100_000,
            tau_min: 0.0,
            tau_max: 1.0,
            tau_steps: 256,
            max_temporal_lag: 64,
            max_structural_radius: 0.1,
            temporal_decay_gamma: 0.0,
            structural_decay_beta: 2.0,
            max_out_edges: 8,
        }
    }

    pub fn tau_grid(&self) -> Vec<f64> {
        linspace(self.tau_min, self.tau_max, self.tau_steps)
    }
}

impl Default for DscdParams {
    fn default() -> Self {
        Self::demo()
    }
}

/// Runtime configuration for DSCD simulation + scaling outputs.
#[derive(Debug, Clone)]
pub struct DscdConfig {
    pub params: DscdParams,
    pub root_event_id: u64,
    pub output_dir: PathBuf,
    pub scaling_ns: Vec<usize>,
}

impl Default for DscdConfig {
    fn default() -> Self {
        Self {
            params: DscdParams::default(),
            root_event_id: 0,
            output_dir: PathBuf::from("output-dsfb-dscd"),
            scaling_ns: vec![5_000, 10_000, 20_000, 50_000, 100_000],
        }
    }
}

/// DSCD event node exported for causal graph and scaling analyses.
///
/// EVENT TYPE NOTE:
/// This is the canonical event struct written to `graph_events.csv`. It
/// carries event/time indices, trust, and ADD-derived structure.
#[derive(Debug, Clone, Serialize)]
pub struct DscdEvent {
    pub id: u64,
    pub time_index: u64,
    pub module_id: u32,
    pub residual_state: String,
    pub trust: f64,
    pub echo_slope: f64,
    pub entropy_density: f64,
    pub structural_level: f64,
}

/// DSCD directed edge candidate with locality + trust-decay provenance.
#[derive(Debug, Clone, Serialize)]
pub struct DscdEdge {
    pub edge_id: u64,
    pub src: u64,
    pub dst: u64,
    pub module_id: u32,
    pub rewrite_rule_id: Option<String>,
    pub base_trust: f64,
    pub effective_trust: f64,
    pub trust_at_creation: f64,
    pub temporal_lag: usize,
    pub structural_distance: f64,
    pub structural_level_source: f64,
    pub structural_level_target: f64,
}

/// Trust-gated DSCD graph at a fixed threshold.
#[derive(Debug, Clone, Default)]
pub struct DscdGraph {
    pub events: Vec<DscdEvent>,
    pub edges: Vec<DscdEdge>,
}

/// DSCD error type alias.
pub type DscdError = anyhow::Error;

#[derive(Debug, Clone)]
struct CurvePoint {
    tau: f64,
    expansion_ratio: f64,
}

#[derive(Debug, Clone)]
struct ScalingMetrics {
    num_events: usize,
    tau_star: f64,
    transition_width: f64,
    max_derivative: f64,
}

#[derive(Debug, Clone)]
struct EventAux {
    event: DscdEvent,
    rewrite_rule_id: String,
    observer_trust: HashMap<u32, f64>,
}

#[derive(Debug, Clone)]
struct Dataset {
    events: Vec<DscdEvent>,
    candidate_edges: Vec<DscdEdge>,
}

#[derive(Debug, Clone)]
struct CandidateEdgeDraft {
    src: u64,
    dst: u64,
    module_id: u32,
    rewrite_rule_id: String,
    base_trust: f64,
    effective_trust: f64,
    temporal_lag: usize,
    structural_distance: f64,
    structural_level_source: f64,
    structural_level_target: f64,
}

/// Deterministic structural distance between two events.
pub fn structural_distance(e_i: &DscdEvent, e_j: &DscdEvent) -> f64 {
    (e_j.structural_level - e_i.structural_level).abs()
}

/// Deterministic trust attenuation over temporal + structural separation.
pub fn effective_trust(base_tau: f64, dt: usize, d_struct: f64, params: &DscdParams) -> f64 {
    let dt_f = dt as f64;
    let decay_temporal = (-params.temporal_decay_gamma * dt_f).exp();
    let decay_struct = (-params.structural_decay_beta * d_struct).exp();
    base_tau * decay_temporal * decay_struct
}

/// Build thresholded graph `G_tau` from precomputed local candidates.
///
/// EDGE SWEEP NOTE:
/// This is the threshold gate used by the tau-sweep writer. Monotonicity holds:
/// increasing `tau` can only remove edges.
pub fn build_graph_for_tau(
    events: &[DscdEvent],
    candidate_edges: &[DscdEdge],
    tau: f64,
) -> DscdGraph {
    let time_by_id: HashMap<u64, u64> = events
        .iter()
        .map(|event| (event.id, event.time_index))
        .collect();

    let mut edges = Vec::new();
    for edge in candidate_edges {
        if edge.effective_trust < tau {
            continue;
        }

        let Some(&src_time) = time_by_id.get(&edge.src) else {
            continue;
        };
        let Some(&dst_time) = time_by_id.get(&edge.dst) else {
            continue;
        };

        // Enforce acyclic time ordering.
        if src_time < dst_time {
            edges.push(edge.clone());
        }
    }

    DscdGraph {
        events: events.to_vec(),
        edges,
    }
}

/// Compute reachable component size from `root_id` using deterministic BFS.
pub fn compute_reachable_component_size(graph: &DscdGraph, root_id: u64) -> usize {
    if graph.events.is_empty() || !graph.events.iter().any(|event| event.id == root_id) {
        return 0;
    }

    let mut adjacency: HashMap<u64, Vec<u64>> = HashMap::new();
    for edge in &graph.edges {
        adjacency.entry(edge.src).or_default().push(edge.dst);
    }
    for neighbors in adjacency.values_mut() {
        neighbors.sort_unstable();
    }

    let mut queue = VecDeque::from([root_id]);
    let mut visited = HashSet::from([root_id]);
    while let Some(current) = queue.pop_front() {
        if let Some(neighbors) = adjacency.get(&current) {
            for next in neighbors {
                if visited.insert(*next) {
                    queue.push_back(*next);
                }
            }
        }
    }

    visited.len()
}

/// Reachable fraction for the root component in the trust-gated graph.
pub fn expansion_ratio(graph: &DscdGraph, root_id: u64) -> f64 {
    if graph.events.is_empty() {
        return 0.0;
    }
    compute_reachable_component_size(graph, root_id) as f64 / graph.events.len() as f64
}

/// Run deterministic DSCD simulation, threshold sweep, and scaling exports.
pub fn run_dscd_simulation(cfg: &DscdConfig) -> std::result::Result<(), DscdError> {
    validate_config(cfg)?;
    fs::create_dir_all(&cfg.output_dir).with_context(|| {
        format!(
            "failed to create DSCD output directory {}",
            cfg.output_dir.display()
        )
    })?;

    let mut run_ns = cfg.scaling_ns.clone();
    run_ns.push(cfg.params.num_events);
    run_ns.sort_unstable();
    run_ns.dedup();

    let max_n = run_ns
        .iter()
        .copied()
        .max()
        .ok_or_else(|| anyhow!("no event counts configured"))?;
    let dataset = generate_dataset(max_n, &cfg.params)?;
    let taus = cfg.params.tau_grid();

    let mut summary_rows = Vec::new();
    let mut curve_by_n: HashMap<usize, Vec<CurvePoint>> = HashMap::new();
    let mut main_events = Vec::new();
    let mut main_edges = Vec::new();

    for n in run_ns {
        let events = dataset.events[..n].to_vec();
        let candidate_edges: Vec<DscdEdge> = dataset
            .candidate_edges
            .iter()
            .filter(|edge| edge.src < n as u64 && edge.dst < n as u64)
            .cloned()
            .collect();

        let curve = compute_threshold_curve(&events, &candidate_edges, &taus, cfg.root_event_id);
        write_threshold_curve_csv(
            &cfg.output_dir.join(format!("threshold_curve_N_{n}.csv")),
            &curve,
        )?;

        let metrics = compute_scaling_metrics(n, &curve);
        summary_rows.push(metrics.clone());
        curve_by_n.insert(n, curve);

        if n == cfg.params.num_events {
            main_events = events;
            main_edges = candidate_edges;
        }
    }

    write_threshold_scaling_summary_csv(
        &cfg.output_dir.join("threshold_scaling_summary.csv"),
        &summary_rows,
    )?;
    write_finite_size_scaling_csv(
        &cfg.output_dir.join("finite_size_scaling.csv"),
        &summary_rows,
    )?;

    let Some(main_curve) = curve_by_n.get(&cfg.params.num_events) else {
        bail!(
            "missing threshold curve for main num_events={}",
            cfg.params.num_events
        );
    };
    let main_metrics = compute_scaling_metrics(cfg.params.num_events, main_curve);
    let main_tau = main_metrics.tau_star;

    write_graph_events_csv(&cfg.output_dir.join("graph_events.csv"), &main_events)?;
    let final_graph = build_graph_for_tau(&main_events, &main_edges, main_tau);
    write_graph_edges_csv(&cfg.output_dir.join("graph_edges.csv"), &final_graph.edges)?;
    write_degree_distribution_csv(
        &cfg.output_dir.join("degree_distribution.csv"),
        &final_graph,
    )?;
    write_interval_sizes_csv(
        &cfg.output_dir.join("interval_sizes.csv"),
        &final_graph,
        1_024,
    )?;
    write_path_lengths_csv(
        &cfg.output_dir.join("path_lengths.csv"),
        &final_graph,
        cfg.root_event_id,
    )?;

    let provenance_path = cfg.output_dir.join("edge_provenance_0.csv");
    write_edge_provenance_csv(&provenance_path, &final_graph, main_tau, 10)?;
    let compatibility_path = cfg.output_dir.join("edge_provenance.csv");
    write_edge_provenance_csv(&compatibility_path, &final_graph, main_tau, 10)?;

    Ok(())
}

/// Build evenly spaced deterministic threshold grid.
pub fn linspace(start: f64, end: f64, count: usize) -> Vec<f64> {
    if count <= 1 {
        return vec![start];
    }
    let step = (end - start) / (count.saturating_sub(1) as f64);
    (0..count).map(|idx| start + step * idx as f64).collect()
}

fn validate_config(cfg: &DscdConfig) -> Result<()> {
    let params = &cfg.params;
    if params.num_events == 0 {
        bail!("num_events must be greater than zero");
    }
    if params.tau_steps < 2 {
        bail!("tau_steps must be at least 2");
    }
    if !params.tau_min.is_finite() || !params.tau_max.is_finite() {
        bail!("tau bounds must be finite");
    }
    if params.tau_max < params.tau_min {
        bail!("tau_max must be >= tau_min");
    }
    if params.max_temporal_lag == 0 {
        bail!("max_temporal_lag must be > 0");
    }
    if !params.max_structural_radius.is_finite() || params.max_structural_radius <= 0.0 {
        bail!("max_structural_radius must be > 0");
    }
    if !params.temporal_decay_gamma.is_finite() || params.temporal_decay_gamma < 0.0 {
        bail!("temporal_decay_gamma must be >= 0");
    }
    if !params.structural_decay_beta.is_finite() || params.structural_decay_beta < 0.0 {
        bail!("structural_decay_beta must be >= 0");
    }
    if params.max_out_edges == 0 {
        bail!("max_out_edges must be > 0");
    }
    if cfg.scaling_ns.iter().any(|n| *n == 0) {
        bail!("scaling_ns values must be greater than zero");
    }
    Ok(())
}

fn generate_dataset(max_events: usize, params: &DscdParams) -> Result<Dataset> {
    let dsfb_cfg = SimConfig {
        steps: max_events,
        ..SimConfig::default()
    };
    let event_batch = generate_dscd_events_from_dsfb(&dsfb_cfg, DsfbParams::default(), max_events)?;

    let mut add_cfg = AddSimulationConfig::default();
    let add_steps = if max_events <= 20_000 {
        max_events.min(2_048)
    } else {
        max_events.min(4_096)
    }
    .max(128);
    add_cfg.steps_per_run = add_steps;
    add_cfg.multi_steps_per_run = vec![add_steps];
    let lambda_grid = add_cfg.lambda_grid();
    let aet = run_aet_sweep(&add_cfg, &lambda_grid)
        .map_err(|error| anyhow!("failed to compute AET sweep for DSCD: {error}"))?;
    let iwlt = run_iwlt_sweep(&add_cfg, &lambda_grid)
        .map_err(|error| anyhow!("failed to compute IWLT sweep for DSCD: {error}"))?;
    let echo_len = aet.echo_slope.len().max(1);
    let entropy_len = iwlt.entropy_density.len().max(1);

    let mut samples_by_event: BTreeMap<u64, Vec<&DscdObserverSample>> = BTreeMap::new();
    for sample in &event_batch.observer_samples {
        samples_by_event
            .entry(sample.event_id)
            .or_default()
            .push(sample);
    }

    let mut aux_events = Vec::with_capacity(event_batch.events.len());
    let mut raw_levels = Vec::with_capacity(event_batch.events.len());

    for (idx, event) in event_batch.events.iter().enumerate() {
        let event_id = event.id.0;
        let samples = samples_by_event.get(&event_id).cloned().unwrap_or_default();

        let mut trust_sum = 0.0;
        let mut residual_sum = 0.0;
        let mut observer_trust = HashMap::new();
        let mut best_observer = 0_u32;
        let mut best_trust = f64::NEG_INFINITY;
        let mut rewrite_rule_id = String::from("stable_envelope");

        for sample in &samples {
            trust_sum += sample.trust;
            residual_sum += sample.residual_summary;
            observer_trust.insert(sample.observer_id, sample.trust);
            if sample.trust > best_trust {
                best_trust = sample.trust;
                best_observer = sample.observer_id;
                rewrite_rule_id = sample.rewrite_rule_label.to_string();
            }
        }

        let sample_count = samples.len().max(1) as f64;
        let avg_residual = residual_sum / sample_count;
        let avg_trust = (trust_sum / sample_count).clamp(0.0, 1.0);
        let residual_state = ResidualState::from_residual(avg_residual)
            .as_str()
            .to_string();

        let echo_slope = aet.echo_slope[idx % echo_len];
        let entropy_density = iwlt.entropy_density[idx % entropy_len];
        let raw_structural = 0.5 * echo_slope + 0.5 * entropy_density;

        let event = DscdEvent {
            id: event_id,
            time_index: idx as u64,
            module_id: best_observer,
            residual_state,
            trust: avg_trust,
            echo_slope,
            entropy_density,
            structural_level: 0.0,
        };

        aux_events.push(EventAux {
            event,
            rewrite_rule_id,
            observer_trust,
        });
        raw_levels.push(raw_structural);
    }

    let (min_level, max_level) = raw_levels
        .iter()
        .fold((f64::INFINITY, f64::NEG_INFINITY), |(mn, mx), value| {
            (mn.min(*value), mx.max(*value))
        });
    let level_span = (max_level - min_level).abs();

    for (idx, aux) in aux_events.iter_mut().enumerate() {
        aux.event.structural_level = if level_span <= f64::EPSILON {
            0.5
        } else {
            ((raw_levels[idx] - min_level) / level_span).clamp(0.0, 1.0)
        };
    }

    // EDGE CONSTRUCTION NOTE:
    // Local candidate generation + bounded out-degree happen here.
    let candidate_edges = build_candidate_edges(&aux_events, params);
    let events = aux_events.into_iter().map(|entry| entry.event).collect();
    Ok(Dataset {
        events,
        candidate_edges,
    })
}

fn build_candidate_edges(events: &[EventAux], params: &DscdParams) -> Vec<DscdEdge> {
    if events.len() < 2 {
        return Vec::new();
    }

    let mut edges = Vec::new();

    for src_idx in 0..(events.len() - 1) {
        let src = &events[src_idx];
        let mut local_candidates = Vec::<CandidateEdgeDraft>::new();

        let max_dst_idx = src_idx
            .saturating_add(params.max_temporal_lag)
            .min(events.len() - 1);

        for dst_idx in (src_idx + 1)..=max_dst_idx {
            let dst = &events[dst_idx];
            let dt = dst_idx - src_idx;
            if dt == 0 {
                continue;
            }

            // Locality constraints: temporal lag + structural neighborhood.
            let d_struct = structural_distance(&src.event, &dst.event);
            if d_struct > params.max_structural_radius {
                continue;
            }

            let base_tau = dst
                .observer_trust
                .get(&dst.event.module_id)
                .copied()
                .unwrap_or(dst.event.trust);
            let tau_eff = effective_trust(base_tau, dt, d_struct, params);
            if !tau_eff.is_finite() || tau_eff <= 0.0 {
                continue;
            }

            local_candidates.push(CandidateEdgeDraft {
                src: src.event.id,
                dst: dst.event.id,
                module_id: dst.event.module_id,
                rewrite_rule_id: src.rewrite_rule_id.clone(),
                base_trust: base_tau,
                effective_trust: tau_eff,
                temporal_lag: dt,
                structural_distance: d_struct,
                structural_level_source: src.event.structural_level,
                structural_level_target: dst.event.structural_level,
            });
        }

        // Bounded out-degree selection: best local edges first.
        local_candidates.sort_by(|a, b| {
            b.effective_trust
                .total_cmp(&a.effective_trust)
                .then(a.temporal_lag.cmp(&b.temporal_lag))
                .then(a.structural_distance.total_cmp(&b.structural_distance))
                .then(a.dst.cmp(&b.dst))
        });

        for draft in local_candidates.into_iter().take(params.max_out_edges) {
            let edge_id = edges.len() as u64;
            edges.push(DscdEdge {
                edge_id,
                src: draft.src,
                dst: draft.dst,
                module_id: draft.module_id,
                rewrite_rule_id: Some(draft.rewrite_rule_id),
                base_trust: draft.base_trust,
                effective_trust: draft.effective_trust,
                trust_at_creation: draft.effective_trust,
                temporal_lag: draft.temporal_lag,
                structural_distance: draft.structural_distance,
                structural_level_source: draft.structural_level_source,
                structural_level_target: draft.structural_level_target,
            });
        }
    }

    edges
}

fn compute_threshold_curve(
    events: &[DscdEvent],
    candidate_edges: &[DscdEdge],
    taus: &[f64],
    root_event_id: u64,
) -> Vec<CurvePoint> {
    // THRESHOLD SWEEP NOTE:
    // This loop writes rho(tau) through repeated threshold gating of the same
    // local candidate set.
    taus.par_iter()
        .map(|tau| {
            let graph = build_graph_for_tau(events, candidate_edges, *tau);
            CurvePoint {
                tau: *tau,
                expansion_ratio: expansion_ratio(&graph, root_event_id),
            }
        })
        .collect()
}

fn compute_scaling_metrics(num_events: usize, curve: &[CurvePoint]) -> ScalingMetrics {
    let max_rho = curve
        .iter()
        .map(|point| point.expansion_ratio)
        .fold(0.0_f64, f64::max);
    let high_target = 0.9 * max_rho;
    let low_target = 0.1 * max_rho;
    let critical_target = 0.5 * max_rho;

    let tau_star = find_tau_at_or_below(curve, critical_target)
        .unwrap_or_else(|| curve.last().map(|point| point.tau).unwrap_or_default());
    let tau_high = find_tau_at_or_below(curve, high_target)
        .unwrap_or_else(|| curve.first().map(|point| point.tau).unwrap_or_default());
    let tau_low = find_tau_at_or_below(curve, low_target)
        .unwrap_or_else(|| curve.last().map(|point| point.tau).unwrap_or(tau_high));

    let max_derivative = curve
        .windows(2)
        .filter_map(|pair| {
            let dt = pair[1].tau - pair[0].tau;
            if dt.abs() <= f64::EPSILON {
                return None;
            }
            Some(((pair[1].expansion_ratio - pair[0].expansion_ratio) / dt).abs())
        })
        .fold(0.0_f64, f64::max);

    ScalingMetrics {
        num_events,
        tau_star,
        transition_width: (tau_low - tau_high).max(0.0),
        max_derivative,
    }
}

fn find_tau_at_or_below(curve: &[CurvePoint], ratio_threshold: f64) -> Option<f64> {
    if curve.is_empty() {
        return None;
    }
    if curve[0].expansion_ratio <= ratio_threshold {
        return Some(curve[0].tau);
    }

    for pair in curve.windows(2) {
        let prev = &pair[0];
        let curr = &pair[1];
        if curr.expansion_ratio > ratio_threshold {
            continue;
        }
        let drho = curr.expansion_ratio - prev.expansion_ratio;
        if drho.abs() <= f64::EPSILON {
            return Some(curr.tau);
        }
        let alpha = (ratio_threshold - prev.expansion_ratio) / drho;
        return Some(prev.tau + alpha * (curr.tau - prev.tau));
    }

    None
}

fn write_threshold_curve_csv(path: &Path, curve: &[CurvePoint]) -> Result<()> {
    let mut writer = Writer::from_path(path)?;
    writer.write_record(["tau", "expansion_ratio"])?;
    for point in curve {
        writer.write_record([point.tau.to_string(), point.expansion_ratio.to_string()])?;
    }
    writer.flush()?;
    Ok(())
}

fn write_threshold_scaling_summary_csv(path: &Path, rows: &[ScalingMetrics]) -> Result<()> {
    let mut writer = Writer::from_path(path)?;
    writer.write_record([
        "num_events",
        "tau_star",
        "transition_width",
        "max_derivative",
    ])?;
    for row in rows {
        writer.write_record([
            row.num_events.to_string(),
            row.tau_star.to_string(),
            row.transition_width.to_string(),
            row.max_derivative.to_string(),
        ])?;
    }
    writer.flush()?;
    Ok(())
}

fn write_finite_size_scaling_csv(path: &Path, rows: &[ScalingMetrics]) -> Result<()> {
    let mut writer = Writer::from_path(path)?;
    writer.write_record([
        "num_events",
        "tau_star",
        "transition_width",
        "max_derivative",
        "width_0_1_to_0_9",
    ])?;
    for row in rows {
        writer.write_record([
            row.num_events.to_string(),
            row.tau_star.to_string(),
            row.transition_width.to_string(),
            row.max_derivative.to_string(),
            row.transition_width.to_string(),
        ])?;
    }
    writer.flush()?;
    Ok(())
}

fn write_graph_events_csv(path: &Path, events: &[DscdEvent]) -> Result<()> {
    let mut writer = Writer::from_path(path)?;
    writer.write_record([
        "id",
        "time_index",
        "module_id",
        "residual_state",
        "trust",
        "echo_slope",
        "entropy_density",
        "structural_level",
    ])?;
    for event in events {
        writer.write_record([
            event.id.to_string(),
            event.time_index.to_string(),
            event.module_id.to_string(),
            event.residual_state.clone(),
            event.trust.to_string(),
            event.echo_slope.to_string(),
            event.entropy_density.to_string(),
            event.structural_level.to_string(),
        ])?;
    }
    writer.flush()?;
    Ok(())
}

fn write_graph_edges_csv(path: &Path, edges: &[DscdEdge]) -> Result<()> {
    let mut writer = Writer::from_path(path)?;
    writer.write_record([
        "edge_id",
        "src",
        "dst",
        "module_id",
        "rewrite_rule_id",
        "base_trust",
        "effective_trust",
        "trust_at_creation",
        "temporal_lag",
        "structural_distance",
        "structural_level_source",
        "structural_level_target",
    ])?;
    for edge in edges {
        writer.write_record([
            edge.edge_id.to_string(),
            edge.src.to_string(),
            edge.dst.to_string(),
            edge.module_id.to_string(),
            edge.rewrite_rule_id.clone().unwrap_or_default(),
            edge.base_trust.to_string(),
            edge.effective_trust.to_string(),
            edge.trust_at_creation.to_string(),
            edge.temporal_lag.to_string(),
            edge.structural_distance.to_string(),
            edge.structural_level_source.to_string(),
            edge.structural_level_target.to_string(),
        ])?;
    }
    writer.flush()?;
    Ok(())
}

fn write_degree_distribution_csv(path: &Path, graph: &DscdGraph) -> Result<()> {
    let mut in_degree: HashMap<u64, usize> =
        graph.events.iter().map(|event| (event.id, 0)).collect();
    let mut out_degree: HashMap<u64, usize> =
        graph.events.iter().map(|event| (event.id, 0)).collect();

    for edge in &graph.edges {
        *out_degree.entry(edge.src).or_insert(0) += 1;
        *in_degree.entry(edge.dst).or_insert(0) += 1;
    }

    let mut writer = Writer::from_path(path)?;
    writer.write_record(["event_id", "in_degree", "out_degree"])?;
    for event in &graph.events {
        writer.write_record([
            event.id.to_string(),
            in_degree.get(&event.id).copied().unwrap_or(0).to_string(),
            out_degree.get(&event.id).copied().unwrap_or(0).to_string(),
        ])?;
    }
    writer.flush()?;
    Ok(())
}

fn write_interval_sizes_csv(path: &Path, graph: &DscdGraph, sample_count: usize) -> Result<()> {
    let mut ids_by_time: Vec<(u64, u64)> = graph
        .events
        .iter()
        .map(|event| (event.time_index, event.id))
        .collect();
    ids_by_time.sort_unstable();
    let ordered_ids: Vec<u64> = ids_by_time.into_iter().map(|(_, id)| id).collect();

    let mut adjacency: HashMap<u64, Vec<u64>> = HashMap::new();
    let mut reverse_adjacency: HashMap<u64, Vec<u64>> = HashMap::new();
    for edge in &graph.edges {
        adjacency.entry(edge.src).or_default().push(edge.dst);
        reverse_adjacency
            .entry(edge.dst)
            .or_default()
            .push(edge.src);
    }

    let pairs = deterministic_pairs(&ordered_ids, sample_count);
    let mut forward_cache: HashMap<u64, HashSet<u64>> = HashMap::new();
    let mut reverse_cache: HashMap<u64, HashSet<u64>> = HashMap::new();

    let mut writer = Writer::from_path(path)?;
    writer.write_record(["src", "dst", "interval_size"])?;
    for (src, dst) in pairs {
        let forward = reachability_set(src, &adjacency, &mut forward_cache);
        let interval_size = if !forward.contains(&dst) {
            0
        } else {
            let backward = reachability_set(dst, &reverse_adjacency, &mut reverse_cache);
            forward
                .intersection(&backward)
                .filter(|node| **node != src && **node != dst)
                .count()
        };

        writer.write_record([src.to_string(), dst.to_string(), interval_size.to_string()])?;
    }
    writer.flush()?;
    Ok(())
}

fn reachability_set(
    start: u64,
    adjacency: &HashMap<u64, Vec<u64>>,
    cache: &mut HashMap<u64, HashSet<u64>>,
) -> HashSet<u64> {
    if let Some(cached) = cache.get(&start) {
        return cached.clone();
    }

    let mut queue = VecDeque::from([start]);
    let mut visited = HashSet::from([start]);
    while let Some(current) = queue.pop_front() {
        if let Some(neighbors) = adjacency.get(&current) {
            for next in neighbors {
                if visited.insert(*next) {
                    queue.push_back(*next);
                }
            }
        }
    }

    cache.insert(start, visited.clone());
    visited
}

fn deterministic_pairs(ordered_ids: &[u64], sample_count: usize) -> Vec<(u64, u64)> {
    if ordered_ids.len() < 2 || sample_count == 0 {
        return Vec::new();
    }

    let n = ordered_ids.len();
    let max_pairs = n.saturating_mul(n.saturating_sub(1)) / 2;
    let target = sample_count.min(max_pairs);

    let mut state = 0x9E37_79B9_7F4A_7C15_u64;
    let mut seen = BTreeSet::new();
    let mut pairs = Vec::with_capacity(target);
    let mut attempts = 0_usize;
    let max_attempts = target.saturating_mul(20).max(128);

    while pairs.len() < target && attempts < max_attempts {
        attempts += 1;
        state = state
            .wrapping_mul(6_364_136_223_846_793_005)
            .wrapping_add(1);
        let i = (state as usize) % (n - 1);
        state = state
            .wrapping_mul(6_364_136_223_846_793_005)
            .wrapping_add(1);
        let gap = ((state as usize) % (n - i - 1)) + 1;
        let j = i + gap;
        if seen.insert((i, j)) {
            pairs.push((ordered_ids[i], ordered_ids[j]));
        }
    }

    if pairs.len() < target {
        for i in 0..(n - 1) {
            for j in (i + 1)..n {
                if seen.insert((i, j)) {
                    pairs.push((ordered_ids[i], ordered_ids[j]));
                    if pairs.len() == target {
                        return pairs;
                    }
                }
            }
        }
    }

    pairs
}

fn write_path_lengths_csv(path: &Path, graph: &DscdGraph, root_event_id: u64) -> Result<()> {
    let mut events = graph.events.clone();
    events.sort_by_key(|event| event.time_index);

    let mut adjacency: HashMap<u64, Vec<u64>> = HashMap::new();
    for edge in &graph.edges {
        adjacency.entry(edge.src).or_default().push(edge.dst);
    }

    let mut distances: HashMap<u64, Option<usize>> =
        events.iter().map(|event| (event.id, None)).collect();
    if distances.contains_key(&root_event_id) {
        distances.insert(root_event_id, Some(0));
    }

    for event in &events {
        let Some(Some(dist)) = distances.get(&event.id).copied() else {
            continue;
        };

        if let Some(neighbors) = adjacency.get(&event.id) {
            for next in neighbors {
                let candidate = dist + 1;
                let current = distances.get(next).copied().flatten();
                if current.map_or(true, |value| candidate > value) {
                    distances.insert(*next, Some(candidate));
                }
            }
        }
    }

    let mut writer = Writer::from_path(path)?;
    writer.write_record(["event_id", "longest_path_length"])?;
    for event in &events {
        writer.write_record([
            event.id.to_string(),
            distances
                .get(&event.id)
                .copied()
                .flatten()
                .unwrap_or(0)
                .to_string(),
        ])?;
    }
    writer.flush()?;
    Ok(())
}

fn write_edge_provenance_csv(
    path: &Path,
    graph: &DscdGraph,
    tau_star: f64,
    limit: usize,
) -> Result<()> {
    let event_by_id: HashMap<u64, &DscdEvent> =
        graph.events.iter().map(|event| (event.id, event)).collect();
    let mut writer = Writer::from_path(path)?;
    writer.write_record([
        "edge_id",
        "source_event",
        "target_event",
        "tau_star",
        "base_trust",
        "effective_trust",
        "dt",
        "d_struct",
        "structural_level_source",
        "structural_level_target",
        "module_id",
        "rewrite_rule_id",
        "residual_state_source",
        "residual_state_target",
    ])?;

    for edge in graph.edges.iter().take(limit) {
        let Some(src) = event_by_id.get(&edge.src) else {
            continue;
        };
        let Some(dst) = event_by_id.get(&edge.dst) else {
            continue;
        };
        writer.write_record([
            edge.edge_id.to_string(),
            edge.src.to_string(),
            edge.dst.to_string(),
            tau_star.to_string(),
            edge.base_trust.to_string(),
            edge.effective_trust.to_string(),
            edge.temporal_lag.to_string(),
            edge.structural_distance.to_string(),
            edge.structural_level_source.to_string(),
            edge.structural_level_target.to_string(),
            edge.module_id.to_string(),
            edge.rewrite_rule_id.clone().unwrap_or_default(),
            src.residual_state.clone(),
            dst.residual_state.clone(),
        ])?;
    }
    writer.flush()?;
    Ok(())
}

/// Create a timestamped DSCD run directory under an output root.
pub fn create_dscd_run_dir(output_root: &Path) -> Result<PathBuf> {
    fs::create_dir_all(output_root).with_context(|| {
        format!(
            "failed to create DSCD output root {}",
            output_root.display()
        )
    })?;

    let stamp = chrono::Local::now().format("%Y%m%d_%H%M%S").to_string();

    let mut run_dir = output_root.join(&stamp);
    let mut suffix = 1_usize;
    while run_dir.exists() {
        run_dir = output_root.join(format!("{stamp}_{suffix:02}"));
        suffix += 1;
    }

    fs::create_dir_all(&run_dir)
        .with_context(|| format!("failed to create {}", run_dir.display()))?;
    Ok(run_dir)
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::time::{SystemTime, UNIX_EPOCH};

    fn toy_events() -> Vec<DscdEvent> {
        (0..8_u64)
            .map(|idx| DscdEvent {
                id: idx,
                time_index: idx,
                module_id: (idx % 2) as u32,
                residual_state: "M".to_string(),
                trust: 0.9 - 0.05 * idx as f64,
                echo_slope: idx as f64 * 0.1,
                entropy_density: idx as f64 * 0.1,
                structural_level: idx as f64 * 0.05,
            })
            .collect()
    }

    fn toy_candidates(events: &[DscdEvent], params: &DscdParams) -> Vec<DscdEdge> {
        let aux: Vec<EventAux> = events
            .iter()
            .map(|event| {
                let mut observer_trust = HashMap::new();
                observer_trust.insert(event.module_id, event.trust);
                EventAux {
                    event: event.clone(),
                    rewrite_rule_id: "stable_envelope".to_string(),
                    observer_trust,
                }
            })
            .collect();
        build_candidate_edges(&aux, params)
    }

    #[test]
    fn trust_gate_respects_tau_and_time_order() {
        let events = vec![
            DscdEvent {
                id: 0,
                time_index: 0,
                module_id: 0,
                residual_state: "L".to_string(),
                trust: 0.9,
                echo_slope: 0.1,
                entropy_density: 0.2,
                structural_level: 0.1,
            },
            DscdEvent {
                id: 1,
                time_index: 1,
                module_id: 0,
                residual_state: "M".to_string(),
                trust: 0.7,
                echo_slope: 0.2,
                entropy_density: 0.3,
                structural_level: 0.2,
            },
        ];
        let edges = vec![
            DscdEdge {
                edge_id: 0,
                src: 0,
                dst: 1,
                module_id: 0,
                rewrite_rule_id: Some("rule_a".to_string()),
                base_trust: 0.8,
                effective_trust: 0.8,
                trust_at_creation: 0.8,
                temporal_lag: 1,
                structural_distance: 0.1,
                structural_level_source: 0.1,
                structural_level_target: 0.2,
            },
            DscdEdge {
                edge_id: 1,
                src: 1,
                dst: 0,
                module_id: 0,
                rewrite_rule_id: Some("rule_b".to_string()),
                base_trust: 0.9,
                effective_trust: 0.9,
                trust_at_creation: 0.9,
                temporal_lag: 1,
                structural_distance: 0.1,
                structural_level_source: 0.2,
                structural_level_target: 0.1,
            },
        ];

        let graph = build_graph_for_tau(&events, &edges, 0.75);
        assert_eq!(graph.edges.len(), 1);
        assert_eq!(graph.edges[0].src, 0);
        assert_eq!(graph.edges[0].dst, 1);
    }

    #[test]
    fn local_rule_enforces_lag_radius_and_outdegree() {
        let params = DscdParams {
            num_events: 128,
            tau_min: 0.0,
            tau_max: 1.0,
            tau_steps: 8,
            max_temporal_lag: 2,
            max_structural_radius: 0.08,
            temporal_decay_gamma: 0.0,
            structural_decay_beta: 2.0,
            max_out_edges: 2,
        };
        let events = toy_events();
        let candidates = toy_candidates(&events, &params);

        let mut out_degree = HashMap::<u64, usize>::new();
        for edge in &candidates {
            let src_event = events
                .iter()
                .find(|event| event.id == edge.src)
                .expect("src");
            let dst_event = events
                .iter()
                .find(|event| event.id == edge.dst)
                .expect("dst");

            let lag = dst_event.time_index - src_event.time_index;
            assert!(lag as usize <= params.max_temporal_lag);

            let d_struct = structural_distance(src_event, dst_event);
            assert!(d_struct <= params.max_structural_radius + 1.0e-12);

            *out_degree.entry(edge.src).or_insert(0) += 1;
        }

        assert!(out_degree
            .values()
            .all(|degree| *degree <= params.max_out_edges));
    }

    #[test]
    fn threshold_monotonicity_holds() {
        let params = DscdParams {
            num_events: 128,
            tau_min: 0.0,
            tau_max: 1.0,
            tau_steps: 8,
            max_temporal_lag: 4,
            max_structural_radius: 1.0,
            temporal_decay_gamma: 0.01,
            structural_decay_beta: 2.0,
            max_out_edges: 3,
        };
        let events = toy_events();
        let candidates = toy_candidates(&events, &params);

        let tau_low = 0.2;
        let tau_high = 0.6;
        let graph_low = build_graph_for_tau(&events, &candidates, tau_low);
        let graph_high = build_graph_for_tau(&events, &candidates, tau_high);

        let low_set: HashSet<(u64, u64)> = graph_low
            .edges
            .iter()
            .map(|edge| (edge.src, edge.dst))
            .collect();
        for edge in &graph_high.edges {
            assert!(low_set.contains(&(edge.src, edge.dst)));
        }
    }

    #[test]
    fn tiny_scaling_run_writes_required_outputs() {
        let out_dir = std::env::temp_dir().join(format!(
            "dsfb-dscd-sim-test-{}",
            SystemTime::now()
                .duration_since(UNIX_EPOCH)
                .expect("clock")
                .as_nanos()
        ));

        let cfg = DscdConfig {
            params: DscdParams {
                num_events: 128,
                tau_min: 0.0,
                tau_max: 1.0,
                tau_steps: 3,
                max_temporal_lag: 8,
                max_structural_radius: 1.0,
                temporal_decay_gamma: 0.0,
                structural_decay_beta: 1.0,
                max_out_edges: 4,
            },
            root_event_id: 0,
            output_dir: out_dir.clone(),
            scaling_ns: vec![64, 128],
        };
        run_dscd_simulation(&cfg).expect("simulation should run");

        assert!(out_dir.join("threshold_scaling_summary.csv").exists());
        assert!(out_dir.join("finite_size_scaling.csv").exists());
        assert!(out_dir.join("threshold_curve_N_64.csv").exists());
        assert!(out_dir.join("threshold_curve_N_128.csv").exists());
        assert!(out_dir.join("graph_events.csv").exists());
        assert!(out_dir.join("graph_edges.csv").exists());
        assert!(out_dir.join("edge_provenance_0.csv").exists());

        let mut reader = csv::Reader::from_path(out_dir.join("threshold_scaling_summary.csv"))
            .expect("summary csv");
        for row in reader.records() {
            let row = row.expect("csv row");
            let tau_star: f64 = row.get(1).expect("tau_star").parse().expect("parse tau");
            let width: f64 = row.get(2).expect("width").parse().expect("parse width");
            assert!((0.0..=1.0).contains(&tau_star));
            assert!(width >= 0.0);
        }

        let _ = fs::remove_dir_all(out_dir);
    }
}