diff --git a/.github/workflows/reproduce-primary.yml b/.github/workflows/reproduce-primary.yml index 91754bb..22fe101 100644 --- a/.github/workflows/reproduce-primary.yml +++ b/.github/workflows/reproduce-primary.yml @@ -68,6 +68,8 @@ jobs: benchmarks/seven_method/pcp_faithful.py \ benchmarks/seven_method/wang_full.py \ benchmarks/seven_method/wang_sensitivity.py \ + benchmarks/seven_method/wang_fixed_candidate.py \ + benchmarks/seven_method/reproduce_external_simulations.py \ benchmarks/seven_method/swine_benchmark.py \ reproduce_simulations.py \ additional_noise_sensitivity.py \ @@ -75,7 +77,7 @@ jobs: synthetic_application_twin.py \ realdata_postselection_diagnostics.py \ candidate_contamination_simulations.py - python -c "import dagguard; import benchmarks.seven_method.pcp_faithful; import benchmarks.seven_method.wang_sensitivity" + python -c "import dagguard; import benchmarks.seven_method.pcp_faithful; import benchmarks.seven_method.wang_sensitivity; import benchmarks.seven_method.wang_fixed_candidate; import benchmarks.seven_method.reproduce_external_simulations" - name: Run deterministic unit tests run: python -m unittest discover -s tests -v - name: Run DAGGuard quickstart @@ -96,3 +98,18 @@ jobs: python candidate_contamination_simulations.py \ --replicates 1 --out results/contamination_smoke_test test -s results/contamination_smoke_test/exact_greedy_paired_summary.csv + - name: Verify Wang fixed-candidate manuscript table + run: | + python -m benchmarks.seven_method.wang_fixed_candidate \ + --replicates 100 --check-manuscript \ + --out results/wang_fixed_candidate_audit + test -s results/wang_fixed_candidate_audit/wang_fixed_candidate_replicates.csv + test -s results/wang_fixed_candidate_audit/wang_fixed_candidate_summary.csv + - name: Verify external seven-method comparator tables + run: | + python -m benchmarks.seven_method.reproduce_external_simulations \ + --replicates 20 --check-manuscript \ + --out results/external_comparator_audit + test -s results/external_comparator_audit/external_comparator_replicates.csv + test -s results/external_comparator_audit/external_comparator_summary.csv + test -s results/external_comparator_audit/wang_threshold_summary.csv diff --git a/benchmarks/seven_method/SOURCE_AUDIT.md b/benchmarks/seven_method/SOURCE_AUDIT.md index 1d72544..cab7f76 100644 --- a/benchmarks/seven_method/SOURCE_AUDIT.md +++ b/benchmarks/seven_method/SOURCE_AUDIT.md @@ -24,6 +24,12 @@ No official public implementation was located. The benchmark therefore labels th The source application is discrete and includes domain-specific root-node restrictions. For the generic continuous simulation benchmark, variables are discretized by empirical tertiles and no source-application root labels are transferred. The paper does not fully formalize the phrase describing exclusion of collider-related nodes in the Step-3 conditioning set, so the implementation uses a conservative documented interpretation. Step 4 is not benchmarked because it only orients the retained skeleton and the common endpoint is skeleton adjacency; its state-wise orientation rule is also not sufficiently specified for a generic continuous-variable adaptation. These choices are limitations of comparability, not claimed features of the authors' original implementation. +### Fixed-candidate Wang-style local-BIC comparison + +The manuscript's separate fixed-candidate comparison is not the end-to-end hybrid adaptation above. `wang_fixed_candidate.py` isolates the local-BIC backward-parent deletion logic in Wang et al. Algorithm 2 and applies it to the same oriented continuous-data candidates used in the 1,200 controlled DAGGuard experiments. Algorithm 2 evaluates each current parent against the same current-parent baseline, collects all individually improving deletions in a `ToRemove` set, removes that set, and repeats. The audit therefore implements batched deletion rather than DAGGuard-Greedy's single-best deletion. The published Wang score is higher-is-better; for the continuous simulations, the identical deletion logic is expressed using the Gaussian local BIC used by DAGGuard, where lower is better. + +The GitHub Actions reproducibility check regenerates all 12 settings x 100 replicates from the pinned deterministic seeds and verifies every rounded value reported in Supplementary Table S3, including the 68% versus 96% exact-agreement contrast and the 2.956 versus 0.145 mean BIC-gap contrast under combined contamination. This closes the reproduction path for the closest local-BIC precedent without representing the continuous-data score substitution as official Wang et al. software. + ## Ordinary PC The ordinary-PC baseline uses an original-style ordered-pair skeleton search with immediate graph updates and two-sided Gaussian Fisher-z tests at `alpha=0.05`. It is presented as a conventional PC skeleton baseline, not as a reproduction of a particular software package. diff --git a/benchmarks/seven_method/reproduce_external_simulations.py b/benchmarks/seven_method/reproduce_external_simulations.py new file mode 100644 index 0000000..bfcf760 --- /dev/null +++ b/benchmarks/seven_method/reproduce_external_simulations.py @@ -0,0 +1,201 @@ +#!/usr/bin/env python3 +"""Reproduce the external-method rows of the seven-method simulation benchmark. + +The 240 datasets are regenerated from the same deterministic seeds used by +``additional_noise_sensitivity.py``. This script reruns the four external +comparators (Wang adaptation, PC-FDR, PC-p, ordinary PC), writes replicate-level +metrics, and verifies the manuscript's primary and sensitivity summaries. + +NOTEARS/DAGGuard rows are produced by the separate pinned end-to-end workflow; +this audit deliberately does not refit NOTEARS, which is substantially more +expensive and already has archived replicate-level outputs. +""" +from __future__ import annotations + +import argparse +from pathlib import Path + +import numpy as np +import pandas as pd + +from .common import ( + BASE_SEED, + simulate_dag_seeded, + simulate_lsem_noise, + simulate_weights_seeded, + skeleton_metrics, +) +from .pc_fdr import pc_fdr_skeleton +from .pc_original import pc_original_skeleton +from .pcp_faithful import pc_p_skeleton_faithful +from .wang_sensitivity import THRESHOLDS, wang_full_skeleton_bins + + +REPO_ROOT = Path(__file__).resolve().parents[2] +PRIMARY_ARCHIVE = REPO_ROOT / "results" / "seven_method_benchmark" / "simulation_primary_seven_methods.csv" + +METHOD_CONFIGS = ( + ("Wang et al. adapted", "wang", (.008, .005, .009)), + ("PC-FDR (q=.05)", "pcfdr", (.05, False)), + ("PC-FDR (q=.10)", "pcfdr", (.10, False)), + ("PC-FDR BY (q=.05)", "pcfdr", (.05, True)), + ("PC-p (q=.05)", "pcp", .05), + ("PC-p (q=.10)", "pcp", .10), + ("Ordinary PC (alpha=.05)", "pc", .05), +) + +EXPECTED_SENSITIVITY = { + "PC-FDR (q=.05)": (0.191, 0.200, 0.467, 19.50), + "PC-FDR (q=.10)": (0.195, 0.205, 0.477, 19.39), + "PC-FDR BY (q=.05)": (0.190, 0.199, 0.448, 19.91), + "PC-p (q=.05)": (0.097, 0.104, 0.366, 20.28), + "PC-p (q=.10)": (0.106, 0.113, 0.381, 20.02), +} + +EXPECTED_WANG = { + "A": (0.531, 0.317, 31.27, 20.31), + "B": (0.531, 0.316, 31.24, 20.20), + "C": (0.533, 0.323, 31.36, 20.72), + "D": (0.532, 0.321, 31.32, 20.60), +} + + +def _dataset(d: int, s: float, noise: str, rep: int, n: int = 500): + noise_index = ["normal", "exponential", "gumbel"].index(noise) + seed = BASE_SEED + 100000 * d + 1000 * int(10 * s) + 10 * noise_index + rep + truth = simulate_dag_seeded(d, 2 * d, seed) + weights = simulate_weights_seeded(truth, s, seed + 1) + X = simulate_lsem_noise(weights, n, noise, seed + 2) + return seed, truth, X + + +def _row(method: str, d: int, s: float, noise: str, rep: int, + seed: int, truth: np.ndarray, estimate: np.ndarray, runtime: float): + return { + "method": method, + "d": d, + "s": s, + "noise": noise, + "rep": rep, + "seed": seed, + "runtime_seconds": runtime, + **skeleton_metrics(truth, estimate), + } + + +def regenerate(replicates: int = 20, n: int = 500): + rows = [] + wang_rows = [] + for d in (10, 20): + for s in (2, 5): + for noise in ("normal", "exponential", "gumbel"): + for rep in range(replicates): + seed, truth, X = _dataset(d, s, noise, rep, n) + + for label, kind, params in METHOD_CONFIGS: + if kind == "wang": + estimate, runtime, _ = wang_full_skeleton_bins(X, *params, bins=3) + elif kind == "pcfdr": + q, by = params + estimate, runtime, _ = pc_fdr_skeleton(X, q=q, by=by, heuristic=False) + elif kind == "pcp": + estimate, runtime, _ = pc_p_skeleton_faithful(X, q=params) + elif kind == "pc": + estimate, runtime, _ = pc_original_skeleton(X, alpha=params) + else: + raise AssertionError(kind) + rows.append(_row(label, d, s, noise, rep, seed, truth, estimate, runtime)) + + for idx, thresholds in enumerate(THRESHOLDS): + estimate, runtime, _ = wang_full_skeleton_bins(X, *thresholds, bins=3) + wang_rows.append(_row(chr(ord("A") + idx), d, s, noise, rep, + seed, truth, estimate, runtime)) + return pd.DataFrame(rows), pd.DataFrame(wang_rows) + + +def summarize(df: pd.DataFrame) -> pd.DataFrame: + out = [] + for method, g in df.groupby("method", sort=False): + selected = int(g["edges"].sum()) + fp = int(g["fp"].sum()) + tp = int(g["tp"].sum()) + fn = int(g["fn"].sum()) + out.append({ + "method": method, + "n": len(g), + "mean_fdp": float(g["fdr"].mean()), + "pooled_fdp": fp / selected if selected else np.nan, + "pooled_tpr": tp / (tp + fn) if tp + fn else np.nan, + "mean_shd": float(g["shd"].mean()), + "mean_adjacencies": float(g["edges"].mean()), + }) + return pd.DataFrame(out) + + +def _assert_close(name: str, observed, expected, digits): + failures = [] + for value, target, ndigits in zip(observed, expected, digits): + if round(float(value), ndigits) != round(float(target), ndigits): + failures.append((value, target, ndigits)) + if failures: + raise AssertionError(f"{name} mismatch: {failures}") + + +def check_manuscript(primary: pd.DataFrame, wang: pd.DataFrame): + # Primary external comparator rows are checked against the committed seven-method table. + archive = pd.read_csv(PRIMARY_ARCHIVE) + mapping = { + "Wang et al. adapted": "Wang et al. (2026) hybrid structural pipeline, tertile-adapted (published setting A)", + "PC-FDR (q=.05)": "Li & Wang (2009) PC-FDR, q=0.05", + "PC-p (q=.05)": "Strobl, Spirtes & Visweswaran (2019) PC-p, q=0.05", + "Ordinary PC (alpha=.05)": "Ordinary PC, alpha=0.05", + } + for ours, archived in mapping.items(): + got = primary.loc[primary.method == ours].iloc[0] + want = archive.loc[archive.method == archived].iloc[0] + _assert_close( + ours, + (got.pooled_fdp, got.pooled_tpr, got.mean_shd, got.mean_adjacencies), + (want.pooled_skeleton_false_discovery_proportion, + want.pooled_skeleton_tpr, want.mean_skeleton_shd, + want.mean_selected_adjacencies), + (12, 12, 12, 12), + ) + + for method, expected in EXPECTED_SENSITIVITY.items(): + got = primary.loc[primary.method == method].iloc[0] + _assert_close(method, + (got.mean_fdp, got.pooled_fdp, got.pooled_tpr, got.mean_shd), + expected, (3, 3, 3, 2)) + + for setting, expected in EXPECTED_WANG.items(): + got = wang.loc[wang.method == setting].iloc[0] + _assert_close("Wang " + setting, + (got.pooled_fdp, got.pooled_tpr, got.mean_shd, got.mean_adjacencies), + expected, (3, 3, 2, 2)) + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--replicates", type=int, default=20) + parser.add_argument("--out", type=Path, required=True) + parser.add_argument("--check-manuscript", action="store_true") + args = parser.parse_args() + args.out.mkdir(parents=True, exist_ok=True) + raw, wang_raw = regenerate(replicates=args.replicates) + summary = summarize(raw) + wang_summary = summarize(wang_raw) + raw.to_csv(args.out / "external_comparator_replicates.csv", index=False) + summary.to_csv(args.out / "external_comparator_summary.csv", index=False) + wang_raw.to_csv(args.out / "wang_threshold_replicates.csv", index=False) + wang_summary.to_csv(args.out / "wang_threshold_summary.csv", index=False) + print(summary.to_string(index=False)) + print(wang_summary.to_string(index=False)) + if args.check_manuscript: + if args.replicates != 20: + raise ValueError("Manuscript check requires 20 replicates per setting") + check_manuscript(summary, wang_summary) + + +if __name__ == "__main__": + main() diff --git a/benchmarks/seven_method/wang_fixed_candidate.py b/benchmarks/seven_method/wang_fixed_candidate.py new file mode 100644 index 0000000..53da3ba --- /dev/null +++ b/benchmarks/seven_method/wang_fixed_candidate.py @@ -0,0 +1,160 @@ +#!/usr/bin/env python3 +"""Audit the Wang et al. (2026) local-BIC pruning precedent on fixed candidates. + +This is distinct from the end-to-end Wang hybrid adaptation in ``wang_full.py``. +It applies Algorithm 2's batched backward-parent deletion logic to the same +oriented continuous-data candidate DAGs used in the controlled DAGGuard study. +The published Wang score is higher-is-better; here the identical deletion logic +uses the Gaussian local BIC used for these continuous simulations, where lower +is better. +""" +from __future__ import annotations + +import argparse +from pathlib import Path + +import numpy as np +import pandas as pd + +from candidate_contamination_simulations import ( + BASE_SEED, + SETTINGS, + assign_weights, + contaminate_candidate, + generate_dag, + population_standardize, + simulate_sem, +) +from local_bic_refinement import ( + edge_jaccard, + exact_refine_dag, + gaussian_local_bic, + greedy_refine_dag, + total_gaussian_bic, +) + + +REPO_ROOT = Path(__file__).resolve().parents[2] +ARCHIVED_SUMMARY = REPO_ROOT / "results" / "seven_method_benchmark" / "wang_fixed_candidate_summary.csv" + + +def wang_style_prune(X: np.ndarray, candidate: np.ndarray, + score_tolerance: float = 1e-10) -> np.ndarray: + """Algorithm-2-style batched backward deletion with Gaussian local BIC. + + Wang et al. Algorithm 2 evaluates every current parent relative to the same + current-parent score, collects all individually improving deletions in a + ``ToRemove`` set, removes that set, and repeats. This differs from + DAGGuard-Greedy, which removes only the single best improving parent per + iteration. + """ + selected = (np.asarray(candidate) != 0).astype(np.int8).copy() + for child in range(selected.shape[0]): + current = [int(p) for p in np.flatnonzero(selected[:, child])] + while current: + baseline = gaussian_local_bic(X, child, current) + to_remove = [] + for parent in current: + trial = [p for p in current if p != parent] + if gaussian_local_bic(X, child, trial) < baseline - score_tolerance: + to_remove.append(parent) + if not to_remove: + break + remove = set(to_remove) + current = [p for p in current if p not in remove] + selected[:, child] = 0 + if current: + selected[np.asarray(current, dtype=int), child] = 1 + return selected + + +def regenerate(replicates: int = 100, d: int = 20, n: int = 500) -> pd.DataFrame: + rows = [] + for setting_index, setting in enumerate(SETTINGS): + for rep in range(replicates): + seed = BASE_SEED + 100_000 * setting_index + rep + rng = np.random.default_rng(seed) + truth, _ = generate_dag(d, setting.true_edges, setting.true_max_indegree, rng) + W = assign_weights(truth, setting.weak_fraction, rng) + if setting.error_regime == "lowvar_heterogeneous": + W, error_scales = population_standardize(W) + else: + error_scales = np.ones(d) + X = simulate_sem(W, n, error_scales, rng) + candidate, _ = contaminate_candidate(truth, setting, rng) + + exact = exact_refine_dag(X, candidate) + greedy = greedy_refine_dag(X, candidate) + wang = wang_style_prune(X, candidate) + wang_gap = total_gaussian_bic(X, wang) - exact.total_bic + greedy_gap = greedy.total_bic - exact.total_bic + if wang_gap < -1e-7 or greedy_gap < -1e-7: + raise AssertionError("A comparator scored below the exact optimum") + rows.append({ + "setting": setting.name, + "rep": rep, + "seed": seed, + "wang_equals_exact": np.array_equal(wang, exact.adjacency), + "greedy_equals_exact": np.array_equal(greedy.adjacency, exact.adjacency), + "wang_exact_jaccard": edge_jaccard(wang, exact.adjacency), + "wang_gap": max(0.0, float(wang_gap)), + "greedy_gap": max(0.0, float(greedy_gap)), + }) + return pd.DataFrame(rows) + + +def summarize(raw: pd.DataFrame) -> pd.DataFrame: + return raw.groupby("setting", as_index=False).agg( + wang_exact=("wang_equals_exact", "mean"), + greedy_exact=("greedy_equals_exact", "mean"), + jaccard=("wang_exact_jaccard", "mean"), + wang_gap=("wang_gap", "mean"), + wang_max_gap=("wang_gap", "max"), + greedy_gap=("greedy_gap", "mean"), + ) + + +def check_archived_summary(summary: pd.DataFrame) -> None: + expected = pd.read_csv(ARCHIVED_SUMMARY).sort_values("setting").reset_index(drop=True) + observed = summary.sort_values("setting").reset_index(drop=True) + if list(observed["setting"]) != list(expected["setting"]): + raise AssertionError("Archived Wang fixed-candidate setting names do not match") + digits = { + "wang_exact": 2, + "greedy_exact": 2, + "jaccard": 4, + "wang_gap": 3, + "wang_max_gap": 3, + "greedy_gap": 3, + } + failures = [] + for column, ndigits in digits.items(): + got = observed[column].round(ndigits).to_numpy() + want = expected[column].round(ndigits).to_numpy() + if not np.array_equal(got, want): + failures.append((column, got.tolist(), want.tolist())) + if failures: + raise AssertionError("Supplement S7 mismatch: " + repr(failures)) + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("--replicates", type=int, default=100) + parser.add_argument("--out", type=Path) + parser.add_argument("--check-manuscript", action="store_true") + args = parser.parse_args() + raw = regenerate(args.replicates) + summary = summarize(raw) + print(summary.to_string(index=False)) + if args.out: + args.out.mkdir(parents=True, exist_ok=True) + raw.to_csv(args.out / "wang_fixed_candidate_replicates.csv", index=False) + summary.to_csv(args.out / "wang_fixed_candidate_summary.csv", index=False) + if args.check_manuscript: + if args.replicates != 100: + raise ValueError("Manuscript check requires 100 replicates per setting") + check_archived_summary(summary) + + +if __name__ == "__main__": + main() diff --git a/results/seven_method_benchmark/wang_fixed_candidate_summary.csv b/results/seven_method_benchmark/wang_fixed_candidate_summary.csv new file mode 100644 index 0000000..23e189f --- /dev/null +++ b/results/seven_method_benchmark/wang_fixed_candidate_summary.csv @@ -0,0 +1,13 @@ +setting,wang_exact,greedy_exact,jaccard,wang_gap,wang_max_gap,greedy_gap +clean_sparse,1.00,1.00,1.0000,0.000,0.000,0.000 +combined_contamination,0.68,0.96,0.9890,2.956,65.598,0.145 +dense_high_indegree,0.87,0.97,0.9959,0.316,8.380,0.065 +dense_moderate,0.96,1.00,0.9988,0.037,2.023,0.000 +fp_025,1.00,1.00,1.0000,0.000,0.000,0.000 +fp_050,0.99,1.00,0.9995,0.000,0.009,0.000 +fp_100,0.98,1.00,0.9991,0.010,0.626,0.000 +lowvar_heterogeneous,1.00,1.00,1.0000,0.000,0.000,0.000 +lowvar_weak_fp,0.93,0.99,0.9960,0.296,10.227,0.022 +missing_010,1.00,1.00,1.0000,0.000,0.000,0.000 +reversal_010,1.00,1.00,1.0000,0.000,0.000,0.000 +weak_fp,0.96,1.00,0.9980,0.406,30.943,0.000