From 3844564b2761ae29ac90ef4144015cfa4a445822 Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:33:47 -0500 Subject: [PATCH 1/9] Add Wang fixed-candidate reproducibility audit --- .../seven_method/wang_fixed_candidate.py | 143 ++++++++++++++++++ 1 file changed, 143 insertions(+) create mode 100644 benchmarks/seven_method/wang_fixed_candidate.py diff --git a/benchmarks/seven_method/wang_fixed_candidate.py b/benchmarks/seven_method/wang_fixed_candidate.py new file mode 100644 index 0000000..18a8dee --- /dev/null +++ b/benchmarks/seven_method/wang_fixed_candidate.py @@ -0,0 +1,143 @@ +#!/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 SETTINGS, run_one +from local_bic_refinement import edge_jaccard, gaussian_local_bic, total_gaussian_bic + + +EXPECTED = { + "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), +} + + +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.""" + 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): + base = run_one(setting, setting_index, rep, d, n) + # run_one is deterministic; rebuild its candidate/data through the archived helper + # by importing the simulation primitives only when needed below. + from candidate_contamination_simulations import ( + BASE_SEED, assign_weights, contaminate_candidate, generate_dag, + population_standardize, simulate_sem, + ) + 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) + from local_bic_refinement import exact_refine_dag, greedy_refine_dag + 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 + rows.append({ + "setting": setting.name, + "rep": rep, + "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.total_bic - exact.total_bic)), + }) + 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_expected(summary: pd.DataFrame) -> None: + failures = [] + for row in summary.itertuples(index=False): + observed = (round(row.wang_exact, 2), round(row.greedy_exact, 2), + round(row.jaccard, 4), round(row.wang_gap, 3), + round(row.wang_max_gap, 3), round(row.greedy_gap, 3)) + if observed != EXPECTED[row.setting]: + failures.append((row.setting, observed, EXPECTED[row.setting])) + 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_expected(summary) + + +if __name__ == "__main__": + main() From 65a6f1e97b2430a854506cc861bad8dc81d61ae6 Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:34:13 -0500 Subject: [PATCH 2/9] Tighten Wang fixed-candidate reproduction script --- .../seven_method/wang_fixed_candidate.py | 53 +++++++++++++------ 1 file changed, 38 insertions(+), 15 deletions(-) diff --git a/benchmarks/seven_method/wang_fixed_candidate.py b/benchmarks/seven_method/wang_fixed_candidate.py index 18a8dee..c75cf6e 100644 --- a/benchmarks/seven_method/wang_fixed_candidate.py +++ b/benchmarks/seven_method/wang_fixed_candidate.py @@ -16,8 +16,22 @@ import numpy as np import pandas as pd -from candidate_contamination_simulations import SETTINGS, run_one -from local_bic_refinement import edge_jaccard, gaussian_local_bic, total_gaussian_bic +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, +) EXPECTED = { @@ -38,7 +52,14 @@ 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.""" + """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])] @@ -63,13 +84,6 @@ 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): - base = run_one(setting, setting_index, rep, d, n) - # run_one is deterministic; rebuild its candidate/data through the archived helper - # by importing the simulation primitives only when needed below. - from candidate_contamination_simulations import ( - BASE_SEED, assign_weights, contaminate_candidate, generate_dag, - population_standardize, simulate_sem, - ) 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) @@ -80,19 +94,23 @@ def regenerate(replicates: int = 100, d: int = 20, n: int = 500) -> pd.DataFrame error_scales = np.ones(d) X = simulate_sem(W, n, error_scales, rng) candidate, _ = contaminate_candidate(truth, setting, rng) - from local_bic_refinement import exact_refine_dag, greedy_refine_dag + 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.total_bic - exact.total_bic)), + "greedy_gap": max(0.0, float(greedy_gap)), }) return pd.DataFrame(rows) @@ -111,9 +129,14 @@ def summarize(raw: pd.DataFrame) -> pd.DataFrame: def check_expected(summary: pd.DataFrame) -> None: failures = [] for row in summary.itertuples(index=False): - observed = (round(row.wang_exact, 2), round(row.greedy_exact, 2), - round(row.jaccard, 4), round(row.wang_gap, 3), - round(row.wang_max_gap, 3), round(row.greedy_gap, 3)) + observed = ( + round(float(row.wang_exact), 2), + round(float(row.greedy_exact), 2), + round(float(row.jaccard), 4), + round(float(row.wang_gap), 3), + round(float(row.wang_max_gap), 3), + round(float(row.greedy_gap), 3), + ) if observed != EXPECTED[row.setting]: failures.append((row.setting, observed, EXPECTED[row.setting])) if failures: From 0280f1f1d2034d4fbd6f28aecd43751fce74f00a Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:34:30 -0500 Subject: [PATCH 3/9] CI-check Wang fixed-candidate manuscript table --- .github/workflows/reproduce-primary.yml | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/.github/workflows/reproduce-primary.yml b/.github/workflows/reproduce-primary.yml index 91754bb..2bbb768 100644 --- a/.github/workflows/reproduce-primary.yml +++ b/.github/workflows/reproduce-primary.yml @@ -68,6 +68,7 @@ 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/swine_benchmark.py \ reproduce_simulations.py \ additional_noise_sensitivity.py \ @@ -75,7 +76,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" - name: Run deterministic unit tests run: python -m unittest discover -s tests -v - name: Run DAGGuard quickstart @@ -96,3 +97,10 @@ 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 benchmarks/seven_method/wang_fixed_candidate.py \ + --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 From 03696129a16075f69b7cbdc836233c6c01087d95 Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:37:12 -0500 Subject: [PATCH 4/9] Run Wang manuscript audit as a module --- .github/workflows/reproduce-primary.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/reproduce-primary.yml b/.github/workflows/reproduce-primary.yml index 2bbb768..c7065ce 100644 --- a/.github/workflows/reproduce-primary.yml +++ b/.github/workflows/reproduce-primary.yml @@ -99,7 +99,7 @@ jobs: test -s results/contamination_smoke_test/exact_greedy_paired_summary.csv - name: Verify Wang fixed-candidate manuscript table run: | - python benchmarks/seven_method/wang_fixed_candidate.py \ + 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 From 25d13061328c920172840e5abdf1f013b3287f00 Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:41:41 -0500 Subject: [PATCH 5/9] Document reproducible Wang fixed-candidate audit --- benchmarks/seven_method/SOURCE_AUDIT.md | 6 ++++++ 1 file changed, 6 insertions(+) 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. From 0e635ab8bde12a1ab68c8a710e31f2c4410d689e Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:42:51 -0500 Subject: [PATCH 6/9] Archive Wang fixed-candidate summary --- .../wang_fixed_candidate_summary.csv | 13 +++++++++++++ 1 file changed, 13 insertions(+) create mode 100644 results/seven_method_benchmark/wang_fixed_candidate_summary.csv 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 From 032425452e9c494ae7e4fd91d5539ba54678a772 Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:43:15 -0500 Subject: [PATCH 7/9] Use archived Wang summary as CI source of truth --- .../seven_method/wang_fixed_candidate.py | 50 ++++++++----------- 1 file changed, 22 insertions(+), 28 deletions(-) diff --git a/benchmarks/seven_method/wang_fixed_candidate.py b/benchmarks/seven_method/wang_fixed_candidate.py index c75cf6e..53da3ba 100644 --- a/benchmarks/seven_method/wang_fixed_candidate.py +++ b/benchmarks/seven_method/wang_fixed_candidate.py @@ -34,20 +34,8 @@ ) -EXPECTED = { - "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), -} +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, @@ -56,7 +44,7 @@ def wang_style_prune(X: np.ndarray, candidate: np.ndarray, 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 + ``ToRemove`` set, removes that set, and repeats. This differs from DAGGuard-Greedy, which removes only the single best improving parent per iteration. """ @@ -126,19 +114,25 @@ def summarize(raw: pd.DataFrame) -> pd.DataFrame: ) -def check_expected(summary: pd.DataFrame) -> None: +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 row in summary.itertuples(index=False): - observed = ( - round(float(row.wang_exact), 2), - round(float(row.greedy_exact), 2), - round(float(row.jaccard), 4), - round(float(row.wang_gap), 3), - round(float(row.wang_max_gap), 3), - round(float(row.greedy_gap), 3), - ) - if observed != EXPECTED[row.setting]: - failures.append((row.setting, observed, EXPECTED[row.setting])) + 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)) @@ -159,7 +153,7 @@ def main() -> None: if args.check_manuscript: if args.replicates != 100: raise ValueError("Manuscript check requires 100 replicates per setting") - check_expected(summary) + check_archived_summary(summary) if __name__ == "__main__": From 01c68ae8c6920a401de2c2d4470fb9a5d0de3392 Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:46:56 -0500 Subject: [PATCH 8/9] Add seven-method external simulation reproduction --- .../reproduce_external_simulations.py | 201 ++++++++++++++++++ 1 file changed, 201 insertions(+) create mode 100644 benchmarks/seven_method/reproduce_external_simulations.py 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() From 9348c5a36f485395470c5827d3d99f8f292decd0 Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:47:12 -0500 Subject: [PATCH 9/9] CI-check all external comparator simulation tables --- .github/workflows/reproduce-primary.yml | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/.github/workflows/reproduce-primary.yml b/.github/workflows/reproduce-primary.yml index c7065ce..22fe101 100644 --- a/.github/workflows/reproduce-primary.yml +++ b/.github/workflows/reproduce-primary.yml @@ -69,6 +69,7 @@ jobs: 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 \ @@ -76,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; import benchmarks.seven_method.wang_fixed_candidate" + 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 @@ -104,3 +105,11 @@ jobs: --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