antecedent_validate/stability/
orientation.rs

1//! Orientation stability via PCMCI+ block bootstrap.
2//!
3//! SPDX-License-Identifier: MIT OR Apache-2.0
4
5#![allow(clippy::cast_possible_truncation, clippy::cast_precision_loss)]
6
7use std::collections::BTreeMap;
8use std::sync::Arc;
9
10use antecedent_core::{ExecutionContext, Lag, VariableId};
11use antecedent_data::{ResamplingPlan, TableView, TimeSeriesData, resample_timeseries};
12use antecedent_discovery::{DiscoveryWorkspace, LaggedLink, PcmciPlus};
13use antecedent_graph::{DenseNodeId, Endpoint, NodeRef};
14
15use crate::error::ValidationError;
16
17use super::pcmci_grid::{LinkStability, report_from_counts};
18
19/// Undirected contemporaneous edge retention frequency.
20#[derive(Clone, Copy, Debug, PartialEq)]
21pub struct UndirectedLinkStability {
22    /// Endpoints with `a.raw() <= b.raw()`.
23    pub a: VariableId,
24    /// Other endpoint.
25    pub b: VariableId,
26    /// Fraction of replicates retaining an undirected contemp edge.
27    pub frequency: f64,
28}
29
30/// Report from [`OrientationStability`].
31#[derive(Clone, Debug)]
32pub struct OrientationStabilityReport {
33    /// Directed contemporaneous edge frequencies.
34    pub directed: Arc<[LinkStability]>,
35    /// Undirected contemporaneous edge frequencies.
36    pub undirected: Arc<[UndirectedLinkStability]>,
37    /// Fraction of replicates that produced ≥1 conflict edge among contemp pairs.
38    pub conflict_rate: f64,
39    /// Bootstrap replicates.
40    pub replicates: u32,
41    /// Moving-block length.
42    pub block_size: usize,
43}
44
45/// Block-bootstrap orientation stability around a [`PcmciPlus`] configuration.
46#[derive(Clone, Debug)]
47pub struct OrientationStability {
48    /// PCMCI+ configuration to re-run.
49    pub pcmci_plus: PcmciPlus,
50    /// Bootstrap replicates.
51    pub replicates: u32,
52    /// Block length.
53    pub block_size: usize,
54}
55
56impl Default for OrientationStability {
57    fn default() -> Self {
58        Self::new()
59    }
60}
61
62impl OrientationStability {
63    /// Defaults: 20 replicates, block size 20, FDR off.
64    #[must_use]
65    pub fn new() -> Self {
66        Self { pcmci_plus: PcmciPlus::new().with_fdr(false), replicates: 20, block_size: 20 }
67    }
68
69    /// Run orientation stability assessment.
70    ///
71    /// # Errors
72    ///
73    /// Data or discovery failures.
74    pub fn run(
75        &self,
76        data: &TimeSeriesData,
77        variables: &[VariableId],
78        workspace: &mut DiscoveryWorkspace,
79        ctx: &ExecutionContext,
80    ) -> Result<OrientationStabilityReport, ValidationError> {
81        if self.replicates == 0 || self.block_size == 0 {
82            return Err(ValidationError::NotApplicable {
83                message: "orientation stability requires positive replicates and block_size",
84            });
85        }
86        if self.block_size > data.row_count() {
87            return Err(ValidationError::NotApplicable {
88                message: "block_size exceeds series length",
89            });
90        }
91        let mut directed_counts: BTreeMap<LaggedLink, u32> = BTreeMap::new();
92        let mut undirected_counts: BTreeMap<(VariableId, VariableId), u32> = BTreeMap::new();
93        let mut conflict_reps = 0u32;
94        let mut rng = ctx.rng.stream(0x0E1E_u64);
95        let mut index_scratch = Vec::new();
96        for _ in 0..self.replicates {
97            let boot = resample_timeseries(
98                data,
99                ResamplingPlan::MovingBlock { length: self.block_size },
100                &mut rng,
101                &mut index_scratch,
102            )
103            .map_err(ValidationError::from)?;
104            let result = self
105                .pcmci_plus
106                .run(&boot, variables, workspace, ctx)
107                .map_err(ValidationError::from)?;
108            let mut had_conflict = false;
109            for edge in result.evidence.graph.edges() {
110                let (Some(va), Some(vb)) = (
111                    lagged_var_lag0(result.evidence.graph.nodes(), edge.a),
112                    lagged_var_lag0(result.evidence.graph.nodes(), edge.b),
113                ) else {
114                    continue;
115                };
116                if edge.is_conflict() {
117                    had_conflict = true;
118                    continue;
119                }
120                if edge.is_undirected() {
121                    let (a, b) = if va.raw() <= vb.raw() { (va, vb) } else { (vb, va) };
122                    *undirected_counts.entry((a, b)).or_insert(0) += 1;
123                    continue;
124                }
125                if edge.is_dag_directed() {
126                    let (src, tgt) = match (edge.at_a, edge.at_b) {
127                        (Endpoint::Tail, Endpoint::Arrow) => (va, vb),
128                        (Endpoint::Arrow, Endpoint::Tail) => (vb, va),
129                        _ => continue,
130                    };
131                    let link = LaggedLink {
132                        source: src,
133                        source_lag: Lag::CONTEMPORANEOUS,
134                        target: tgt,
135                        target_lag: Lag::CONTEMPORANEOUS,
136                    };
137                    *directed_counts.entry(link).or_insert(0) += 1;
138                }
139            }
140            if had_conflict {
141                conflict_reps += 1;
142            }
143        }
144        let directed = report_from_counts(directed_counts, self.replicates, self.block_size);
145        let mut undirected = Vec::with_capacity(undirected_counts.len());
146        for ((a, b), c) in undirected_counts {
147            undirected.push(UndirectedLinkStability {
148                a,
149                b,
150                frequency: f64::from(c) / f64::from(self.replicates),
151            });
152        }
153        undirected.sort_by(|x, y| {
154            y.frequency.partial_cmp(&x.frequency).unwrap_or(std::cmp::Ordering::Equal)
155        });
156        Ok(OrientationStabilityReport {
157            directed: directed.frequencies,
158            undirected: Arc::from(undirected),
159            conflict_rate: f64::from(conflict_reps) / f64::from(self.replicates),
160            replicates: self.replicates,
161            block_size: self.block_size,
162        })
163    }
164}
165
166fn lagged_var_lag0(nodes: &[NodeRef], id: DenseNodeId) -> Option<VariableId> {
167    let node = nodes.get(id.as_usize())?;
168    match node {
169        NodeRef::Lagged { variable, lag } if lag.raw() == 0 => Some(*variable),
170        _ => None,
171    }
172}
173
174#[cfg(test)]
175#[allow(clippy::cast_precision_loss)]
176mod tests {
177    use std::sync::Arc;
178
179    use antecedent_core::{
180        CausalSchemaBuilder, ExecutionContext, Lag, MeasurementSpec, RoleHint, SmallRoleSet,
181        ValueType, VariableId,
182    };
183    use antecedent_data::{
184        Float64Column, OwnedColumn, OwnedColumnarStorage, SamplingRegularity, TimeIndex,
185        TimeSeriesData, ValidityBitmap,
186    };
187    use antecedent_discovery::{DiscoveryConstraints, DiscoveryWorkspace, TemporalConstraints};
188
189    use super::*;
190
191    fn contemp_chain() -> (TimeSeriesData, Vec<VariableId>) {
192        let n = 250usize;
193        let mut b = CausalSchemaBuilder::new();
194        for name in ["x", "y"] {
195            b.add_variable(
196                name,
197                ValueType::Continuous,
198                SmallRoleSet::from_hint(RoleHint::Context),
199                None,
200                None,
201                MeasurementSpec::default(),
202            )
203            .unwrap();
204        }
205        let schema = b.build().unwrap();
206        let mut x = vec![0.0; n];
207        let mut y = vec![0.0; n];
208        for t in 0..n {
209            x[t] = (t as f64 * 0.03).sin();
210            y[t] = 0.85 * x[t] + 0.05 * (t as f64 * 0.07).cos();
211        }
212        let cols = vec![
213            OwnedColumn::Float64(
214                Float64Column::new(
215                    VariableId::from_raw(0),
216                    Arc::from(x),
217                    ValidityBitmap::all_valid(n),
218                )
219                .unwrap(),
220            ),
221            OwnedColumn::Float64(
222                Float64Column::new(
223                    VariableId::from_raw(1),
224                    Arc::from(y),
225                    ValidityBitmap::all_valid(n),
226                )
227                .unwrap(),
228            ),
229        ];
230        let storage = OwnedColumnarStorage::try_new(schema, cols, None, None).unwrap();
231        let data = TimeSeriesData::try_new(
232            storage,
233            TimeIndex { regularity: SamplingRegularity::Regular { interval_ns: 1 }, length: n },
234        )
235        .unwrap();
236        (data, vec![VariableId::from_raw(0), VariableId::from_raw(1)])
237    }
238
239    #[test]
240    fn contemp_dependence_appears_in_orientation_report() {
241        let (data, vars) = contemp_chain();
242        let constraints = DiscoveryConstraints {
243            temporal: TemporalConstraints {
244                max_lag: Lag::from_raw(1),
245                min_lag: Lag::CONTEMPORANEOUS,
246            },
247            max_cond_size: 1,
248            alpha: 0.1,
249            ..Default::default()
250        };
251        let stab = OrientationStability {
252            pcmci_plus: PcmciPlus::new().with_fdr(false).with_constraints(constraints),
253            replicates: 6,
254            block_size: 30,
255        };
256        let mut ws = DiscoveryWorkspace::default();
257        let ctx = ExecutionContext::for_tests(9);
258        let report = stab.run(&data, &vars, &mut ws, &ctx).unwrap();
259        assert_eq!(report.replicates, 6);
260        let any = report.directed.iter().any(|l| l.frequency > 0.0)
261            || report.undirected.iter().any(|l| l.frequency > 0.0);
262        assert!(
263            any,
264            "expected contemp edge retention; directed={:?} undirected={:?}",
265            report.directed, report.undirected
266        );
267    }
268}
antecedent_validate/stability/orientation.rs

antecedent_validate/stability/
orientation.rs
