From ad6d35ed2eb51c59a8d97840e06efdffeea55380 Mon Sep 17 00:00:00 2001 From: Vivek Pandit Date: Tue, 6 Oct 2026 20:40:42 -0700 Subject: [PATCH 1/2] benchmark: add reproducible CE vs PCMCI harness and issue #18 table Implements the benchmark harness proposed on #18 (nil-vn design, approved by kslote1): canonical linear-Gaussian generator (rho=0.7, p=0.2), exact (source, target, lag) matching, serial Gaussian oCSE vs PCMCI+ParCorr, median runtime + micro-averaged accuracy, profiles quick/small/full with CSV/JSON/Markdown outputs. Full table (n=5,10,20, T=500, 5 seeds) reproduces the ~32x gap at comparable accuracy. CI covers a tiny smoke path only. --- .gitignore | 3 + benchmarks/README.md | 47 +++ benchmarks/RESULTS_issue18.md | 14 + benchmarks/__init__.py | 4 + benchmarks/adapters.py | 83 ++++++ benchmarks/metrics.py | 66 +++++ benchmarks/run.py | 271 ++++++++++++++++++ .../tests/test_benchmark_harness.py | 89 ++++++ 8 files changed, 577 insertions(+) create mode 100644 benchmarks/README.md create mode 100644 benchmarks/RESULTS_issue18.md create mode 100644 benchmarks/__init__.py create mode 100644 benchmarks/adapters.py create mode 100644 benchmarks/metrics.py create mode 100644 benchmarks/run.py create mode 100644 causationentropy/tests/test_benchmark_harness.py diff --git a/.gitignore b/.gitignore index 3110646..78db695 100644 --- a/.gitignore +++ b/.gitignore @@ -83,3 +83,6 @@ venv.bak/ .Trashes ehthumbs.db Thumbs.db + +# Local benchmark artifacts (issue #18) +benchmark-results/ diff --git a/benchmarks/README.md b/benchmarks/README.md new file mode 100644 index 0000000..f3ef025 --- /dev/null +++ b/benchmarks/README.md @@ -0,0 +1,47 @@ +# Benchmark harness for issue #18 + +Reproducible runtime/accuracy comparison: CausationEntropy Gaussian oCSE +vs Tigramite PCMCI + ParCorr on identical synthetic data. + +## Quick start + +```bash +python -m benchmarks.run --profile quick --output benchmark-results +python -m benchmarks.run --profile full --output benchmark-results +``` + +## Conventions (answers to the open questions on #18) + +* **Generator:** `linear_stochastic_gaussian_process(rho=0.7, p=0.2)`, + explicit `erdos_renyi_graph(n, p, seed, directed=True)` ground truth. + With `p=0.2`, `n=20` yields ~77 true edges, matching the issue report. +* **Edge matching:** exact `(source, target, lag)`; all ground-truth + edges are lag 1, so lag-2 discoveries count as false positives. +* **Self-lags:** included in the main counts (the generator produces no + self-loops, so any self-lag discovery is a false positive). +* **CE config:** `method="standard"`, Gaussian, serial (`n_jobs=1`), + `alpha_forward=alpha_backward=0.05`, profile-dependent `n_shuffles`. +* **PCMCI config:** `ParCorr`, `tau_max=max_lag`, `pc_alpha=0.05`. + Only directed links at lag >= 1 count; undirected/conflicting and + contemporaneous links are ignored. +* **Significance comparability:** end-to-end defaults (permutation vs + analytical). An analytical Gaussian backend is out of scope here. +* **Repetitions:** fixed seeds; runtime = median, accuracy = + micro-averaged (TP/FP/FN summed across seeds). +* **Execution:** serial comparison is primary. Parallel (`n_jobs`) results + from #28 can be added as a separate table later. +* **Profiles:** `quick` (smoke), `small` (accuracy signal), + `full` (the issue table: `n=5,10,20`, `T=500`, 5 seeds). + +## Outputs + +`--output` receives `run_results.csv` (per-seed rows), `summary.json` +(config + aggregates), and `SUMMARY.md` (the comparison table). +Raw result directories are local artifacts (see `.gitignore`); copy the +final table to `benchmarks/RESULTS_issue18.md` when publishing. + +## CI + +The full benchmark does not run in CI. CI covers the generators, +adapters, and metrics on a tiny dataset via +`causationentropy/tests/test_benchmark_harness.py`. diff --git a/benchmarks/RESULTS_issue18.md b/benchmarks/RESULTS_issue18.md new file mode 100644 index 0000000..440dbe3 --- /dev/null +++ b/benchmarks/RESULTS_issue18.md @@ -0,0 +1,14 @@ +# Issue #18 benchmark: CausationEntropy vs Tigramite PCMCI + +Generator: `linear_stochastic_gaussian_process(rho=0.7, p=0.2)`, `max_lag`/`tau_max`=2. +CE: Gaussian oCSE (`method='standard'`), serial (`n_jobs=1`), `n_shuffles=200`. PCMCI: `ParCorr`, `pc_alpha=0.05`. +Matching: exact `(source, target, lag)`; runtime is the median across seeds, accuracy micro-averaged (TP/FP/FN summed). + +| method | n | T | seeds | runtime median s | TP | FP | FN | precision | recall | F1 | +|---|---|---|---|---|---|---|---|---|---|---| +| causationentropy-gaussian-oCSE | 5 | 500 | 5 | 1.6 | 14 | 7 | 3 | 0.667 | 0.824 | 0.737 | +| tigramite-pcmci-parcorr | 5 | 500 | 5 | 0.0 | 14 | 11 | 3 | 0.560 | 0.824 | 0.667 | +| causationentropy-gaussian-oCSE | 10 | 500 | 5 | 7.0 | 69 | 32 | 3 | 0.683 | 0.958 | 0.798 | +| tigramite-pcmci-parcorr | 10 | 500 | 5 | 0.2 | 68 | 47 | 4 | 0.591 | 0.944 | 0.727 | +| causationentropy-gaussian-oCSE | 20 | 500 | 5 | 31.1 | 329 | 135 | 55 | 0.709 | 0.857 | 0.776 | +| tigramite-pcmci-parcorr | 20 | 500 | 5 | 0.9 | 324 | 190 | 60 | 0.630 | 0.844 | 0.722 | diff --git a/benchmarks/__init__.py b/benchmarks/__init__.py new file mode 100644 index 0000000..7b94e1b --- /dev/null +++ b/benchmarks/__init__.py @@ -0,0 +1,4 @@ +"""Reproducible runtime/accuracy benchmarks (issue #18). + +See ``benchmarks/README.md`` for conventions and usage. +""" diff --git a/benchmarks/adapters.py b/benchmarks/adapters.py new file mode 100644 index 0000000..df29462 --- /dev/null +++ b/benchmarks/adapters.py @@ -0,0 +1,83 @@ +"""Dataset generation and method adapters for issue #18.""" + +from __future__ import annotations + +import time + +import networkx as nx +import numpy as np + +from causationentropy.datasets.synthetic import linear_stochastic_gaussian_process + +CANONICAL_RHO = 0.7 +CANONICAL_P = 0.2 +CANONICAL_EPSILON = 0.1 + + +def generate_case( + n: int, T: int, seed: int, rho: float = CANONICAL_RHO, p: float = CANONICAL_P +): + """Generate one benchmark case. + + Returns ``(data, truth, node_order)`` where ``truth`` is the exact-match + ``(source, target, lag)`` edge set (all ground-truth edges are lag 1) + and ``node_order`` maps column ``k`` to node label ``'X{k}'``. + """ + G_true = nx.erdos_renyi_graph(n, p, seed=seed, directed=True) + data, _ = linear_stochastic_gaussian_process( + rho, n=n, T=T, p=p, epsilon=CANONICAL_EPSILON, seed=seed, G=G_true + ) + truth = {(u, v, 1) for u, v in G_true.edges()} + node_order = [f"X{k}" for k in range(n)] + return data, truth, node_order + + +def run_causation_entropy( + data, + max_lag: int, + seed: int, + n_shuffles: int, + alpha_forward: float = 0.05, + alpha_backward: float = 0.05, +): + """Run Gaussian oCSE serially; return ``(graph, runtime_s)``.""" + from causationentropy.core.discovery import discover_network + + start = time.perf_counter() + graph = discover_network( + data, + method="standard", + information="gaussian", + max_lag=max_lag, + alpha_forward=alpha_forward, + alpha_backward=alpha_backward, + n_shuffles=n_shuffles, + n_jobs=1, + random_state=seed, + ) + return graph, time.perf_counter() - start + + +def run_pcmci_parcorr(data, tau_max: int, pc_alpha: float = 0.05): + """Run PCMCI + ParCorr; return ``(results_dict, runtime_s)``. + + Raises + ------ + ImportError + If ``tigramite`` is not installed. + """ + try: + from tigramite import data_processing as pp + from tigramite.independence_tests.parcorr import ParCorr + from tigramite.pcmci import PCMCI + except ImportError as exc: + raise ImportError( + "tigramite is required for the PCMCI baseline " + "(pip install 'causationentropy[dev]')." + ) from exc + + dataframe = pp.DataFrame(data, var_names=[f"X{i}" for i in range(data.shape[1])]) + start = time.perf_counter() + pcmci = PCMCI(dataframe=dataframe, cond_ind_test=ParCorr(), verbosity=0) + results = pcmci.run_pcmci(tau_max=tau_max, pc_alpha=pc_alpha) + return results, time.perf_counter() - start diff --git a/benchmarks/metrics.py b/benchmarks/metrics.py new file mode 100644 index 0000000..916c133 --- /dev/null +++ b/benchmarks/metrics.py @@ -0,0 +1,66 @@ +"""Edge-set normalization and graph-recovery metrics for issue #18. + +Canonical rule (see ``benchmarks/README.md``): a true positive requires an +exact ``(source, target, lag)`` match. Self-lag links are included in the +main counts. +""" + +from __future__ import annotations + + +def normalize_ce_edges(graph, node_order) -> set: + """Normalize a ``discover_network`` graph to ``(source, target, lag)`` ints. + + Parameters + ---------- + graph: + NetworkX MultiDiGraph from :func:`discover_network`. Nodes are + expected to be ``'X{k}'`` labels matching ``node_order``. + node_order: + Sequence mapping column index to node label. + """ + index = {node: k for k, node in enumerate(node_order)} + edges = set() + for u, v, data in graph.edges(data=True): + lag = int(data.get("lag", 0)) + if lag < 1: + continue + edges.add((int(index[u]), int(index[v]), lag)) + return edges + + +def normalize_pcmci_edges(graph) -> set: + """Normalize a :func:`pcmci_to_networkx` graph to ``(source, target, lag)``. + + Only directed links at lag >= 1 are counted as discoveries; + undirected / conflicting / contemporaneous links are ignored. + """ + edges = set() + for u, v, data in graph.edges(data=True): + if data.get("link_type", "directed") != "directed": + continue + lag = int(data.get("lag", 0)) + if lag < 1: + continue + edges.add((int(u), int(v), lag)) + return edges + + +def score_edges(truth: set, predicted: set) -> dict: + """Score ``predicted`` against ``truth`` with exact-match counting.""" + tp = len(truth & predicted) + fp = len(predicted - truth) + fn = len(truth - predicted) + precision = tp / (tp + fp) if (tp + fp) else 0.0 + recall = tp / (tp + fn) if (tp + fn) else 0.0 + f1 = 2 * precision * recall / (precision + recall) if (precision + recall) else 0.0 + return { + "TP": tp, + "FP": fp, + "FN": fn, + "precision": precision, + "recall": recall, + "F1": f1, + "n_truth": len(truth), + "n_predicted": len(predicted), + } diff --git a/benchmarks/run.py b/benchmarks/run.py new file mode 100644 index 0000000..31e4e0a --- /dev/null +++ b/benchmarks/run.py @@ -0,0 +1,271 @@ +"""Reproducible runtime/accuracy benchmark for issue #18. + +Usage:: + + python -m benchmarks.run --profile quick --output benchmark-results + python -m benchmarks.run --profile full --output benchmark-results + +The full benchmark runs locally (it takes several minutes); CI should only +exercise the tiny smoke path via ``causationentropy/tests/test_benchmark_harness.py``. +""" + +from __future__ import annotations + +import argparse +import csv +import json +import statistics +import sys +from pathlib import Path + +from benchmarks.adapters import generate_case, run_causation_entropy, run_pcmci_parcorr +from benchmarks.metrics import ( + normalize_ce_edges, + normalize_pcmci_edges, + score_edges, +) + +PROFILES = { + # Smoke: seconds, exercises the whole pipeline. + "quick": { + "sizes": [(5, 200)], + "seeds": [0, 1], + "n_shuffles": 30, + }, + # Reviewable accuracy signal without the n=20 cost. + "small": { + "sizes": [(5, 500), (10, 500)], + "seeds": [0, 1, 2], + "n_shuffles": 50, + }, + # The issue #18 table: n=5, 10, 20 at T=500. + "full": { + "sizes": [(5, 500), (10, 500), (20, 500)], + "seeds": [0, 1, 2, 3, 4], + "n_shuffles": 200, + }, +} + +CE_CONFIG = { + "method": "standard", + "information": "gaussian", + "alpha_forward": 0.05, + "alpha_backward": 0.05, +} +PCMCI_CONFIG = {"pc_alpha": 0.05} + + +def run_profile( + sizes, + seeds, + max_lag: int, + n_shuffles: int, + pc_alpha: float, + skip_pcmci: bool, + verbose: bool = True, +): + """Run all cases; return a list of per-seed result rows.""" + from causationentropy.graph.utils import pcmci_to_networkx + + rows = [] + for n, T in sizes: + for seed in seeds: + data, truth, node_order = generate_case(n, T, seed) + ce_graph, ce_time = run_causation_entropy( + data, + max_lag=max_lag, + seed=seed, + n_shuffles=n_shuffles, + alpha_forward=CE_CONFIG["alpha_forward"], + alpha_backward=CE_CONFIG["alpha_backward"], + ) + ce_scores = score_edges(truth, normalize_ce_edges(ce_graph, node_order)) + row = { + "n": n, + "T": T, + "max_lag": max_lag, + "seed": seed, + "method": "causationentropy-gaussian-oCSE", + "runtime_s": ce_time, + **ce_scores, + } + rows.append(row) + if verbose: + print( + f"CE n={n} T={T} seed={seed}: " + f"{ce_time:.1f}s TP={ce_scores['TP']} " + f"FP={ce_scores['FP']} FN={ce_scores['FN']}", + flush=True, + ) + + if skip_pcmci: + continue + try: + results, pcmci_time = run_pcmci_parcorr( + data, tau_max=max_lag, pc_alpha=pc_alpha + ) + except ImportError as exc: + print(f"WARNING: {exc} Skipping PCMCI baseline.", file=sys.stderr) + return rows + pcmci_graph = pcmci_to_networkx(results) + pcmci_scores = score_edges(truth, normalize_pcmci_edges(pcmci_graph)) + row = { + "n": n, + "T": T, + "max_lag": max_lag, + "seed": seed, + "method": "tigramite-pcmci-parcorr", + "runtime_s": pcmci_time, + **pcmci_scores, + } + rows.append(row) + if verbose: + print( + f"PCMCI n={n} T={T} seed={seed}: " + f"{pcmci_time:.1f}s TP={pcmci_scores['TP']} " + f"FP={pcmci_scores['FP']} FN={pcmci_scores['FN']}", + flush=True, + ) + return rows + + +def aggregate(rows): + """Aggregate per-seed rows to one row per (method, n, T). + + Runtime is the median across seeds; accuracy is micro-averaged + (TP/FP/FN summed, then precision/recall/F1 derived). + """ + groups = {} + for row in rows: + groups.setdefault((row["method"], row["n"], row["T"]), []).append(row) + summary = [] + for (method, n, T), group in sorted( + groups.items(), key=lambda kv: (kv[0][1], kv[0][0]) + ): + runtimes = sorted(r["runtime_s"] for r in group) + tp = sum(r["TP"] for r in group) + fp = sum(r["FP"] for r in group) + fn = sum(r["FN"] for r in group) + precision = tp / (tp + fp) if (tp + fp) else 0.0 + recall = tp / (tp + fn) if (tp + fn) else 0.0 + f1 = ( + 2 * precision * recall / (precision + recall) + if (precision + recall) + else 0.0 + ) + summary.append( + { + "method": method, + "n": n, + "T": T, + "n_seeds": len(group), + "runtime_median_s": statistics.median(runtimes), + "runtime_min_s": runtimes[0], + "runtime_max_s": runtimes[-1], + "TP": tp, + "FP": fp, + "FN": fn, + "n_truth": tp + fn, + "precision": precision, + "recall": recall, + "F1": f1, + } + ) + return summary + + +def markdown_table(summary, max_lag: int, n_shuffles: int, pc_alpha: float) -> str: + lines = [ + "# Issue #18 benchmark: CausationEntropy vs Tigramite PCMCI", + "", + f"Generator: `linear_stochastic_gaussian_process(rho=0.7, p=0.2)`, " + f"`max_lag`/`tau_max`={max_lag}.", + f"CE: Gaussian oCSE (`method='standard'`), serial (`n_jobs=1`), " + f"`n_shuffles={n_shuffles}`. PCMCI: `ParCorr`, `pc_alpha={pc_alpha}`.", + "Matching: exact `(source, target, lag)`; runtime is the median " + "across seeds, accuracy micro-averaged (TP/FP/FN summed).", + "", + "| method | n | T | seeds | runtime median s | TP | FP | FN | precision | recall | F1 |", + "|---|---|---|---|---|---|---|---|---|---|---|", + ] + for entry in summary: + lines.append( + f"| {entry['method']} | {entry['n']} | {entry['T']} | " + f"{entry['n_seeds']} | {entry['runtime_median_s']:.1f} | " + f"{entry['TP']} | {entry['FP']} | {entry['FN']} | " + f"{entry['precision']:.3f} | {entry['recall']:.3f} | {entry['F1']:.3f} |" + ) + return "\n".join(lines) + "\n" + + +def main(argv=None) -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--profile", choices=sorted(PROFILES), default="quick") + parser.add_argument( + "--output", + default="benchmark-results", + help="Directory for run_results.csv, summary.json, SUMMARY.md", + ) + parser.add_argument("--max-lag", type=int, default=2) + parser.add_argument( + "--n-shuffles", type=int, default=None, help="Override the profile default" + ) + parser.add_argument("--pc-alpha", type=float, default=PCMCI_CONFIG["pc_alpha"]) + parser.add_argument( + "--seeds", + type=int, + nargs="*", + default=None, + help="Override the profile seed list", + ) + parser.add_argument("--skip-pcmci", action="store_true") + args = parser.parse_args(argv) + + profile = PROFILES[args.profile] + seeds = args.seeds if args.seeds is not None else profile["seeds"] + n_shuffles = ( + args.n_shuffles if args.n_shuffles is not None else profile["n_shuffles"] + ) + + rows = run_profile( + profile["sizes"], + seeds, + args.max_lag, + n_shuffles, + args.pc_alpha, + args.skip_pcmci, + ) + summary = aggregate(rows) + + out = Path(args.output) + out.mkdir(parents=True, exist_ok=True) + with open(out / "run_results.csv", "w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=list(rows[0].keys())) + writer.writeheader() + writer.writerows(rows) + with open(out / "summary.json", "w") as fh: + json.dump( + { + "config": { + "profile": args.profile, + "max_lag": args.max_lag, + "n_shuffles": n_shuffles, + "pc_alpha": args.pc_alpha, + "seeds": seeds, + **CE_CONFIG, + }, + "summary": summary, + }, + fh, + indent=2, + ) + table = markdown_table(summary, args.max_lag, n_shuffles, args.pc_alpha) + with open(out / "SUMMARY.md", "w") as fh: + fh.write(table) + print() + print(table) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/causationentropy/tests/test_benchmark_harness.py b/causationentropy/tests/test_benchmark_harness.py new file mode 100644 index 0000000..5c0c1c7 --- /dev/null +++ b/causationentropy/tests/test_benchmark_harness.py @@ -0,0 +1,89 @@ +"""Smoke coverage for the issue #18 benchmark harness (CI-fast).""" + +import networkx as nx +import pytest + +from benchmarks.adapters import generate_case, run_causation_entropy +from benchmarks.metrics import normalize_ce_edges, score_edges +from benchmarks.run import aggregate, markdown_table + + +def test_generate_case_truth_matches_graph(): + data, truth, node_order = generate_case(4, 80, seed=0) + assert data.shape == (80, 4) + assert node_order == ["X0", "X1", "X2", "X3"] + G = nx.erdos_renyi_graph(4, 0.2, seed=0, directed=True) + assert truth == {(u, v, 1) for u, v in G.edges()} + + +def test_score_edges_exact_match(): + out = score_edges({(0, 1, 1), (1, 2, 1)}, {(0, 1, 1), (0, 2, 2)}) + assert (out["TP"], out["FP"], out["FN"]) == (1, 1, 1) + assert out["precision"] == pytest.approx(0.5) + assert out["recall"] == pytest.approx(0.5) + + +def test_tiny_end_to_end_ce_only(): + data, truth, node_order = generate_case(3, 60, seed=0) + graph, runtime = run_causation_entropy(data, max_lag=1, seed=0, n_shuffles=5) + assert runtime > 0 + predicted = normalize_ce_edges(graph, node_order) + assert all(isinstance(e[0], int) and e[2] >= 1 for e in predicted) + out = score_edges(truth, predicted) + assert out["TP"] + out["FN"] == len(truth) + + +def test_pcmci_adapter_smoke(): + pytest.importorskip("tigramite") + from benchmarks.adapters import run_pcmci_parcorr + from benchmarks.metrics import normalize_pcmci_edges + from causationentropy.graph.utils import pcmci_to_networkx + + data, truth, _ = generate_case(3, 60, seed=0) + results, runtime = run_pcmci_parcorr(data, tau_max=1, pc_alpha=0.05) + assert runtime > 0 + predicted = normalize_pcmci_edges(pcmci_to_networkx(results)) + out = score_edges(truth, predicted) + assert out["TP"] + out["FN"] == len(truth) + + +def test_aggregate_and_markdown(): + rows = [ + { + "method": "m", + "n": 5, + "T": 200, + "seed": 0, + "runtime_s": 2.0, + "TP": 4, + "FP": 1, + "FN": 1, + "precision": 0.8, + "recall": 0.8, + "F1": 0.8, + "n_truth": 5, + "n_predicted": 5, + }, + { + "method": "m", + "n": 5, + "T": 200, + "seed": 1, + "runtime_s": 4.0, + "TP": 3, + "FP": 3, + "FN": 2, + "precision": 0.5, + "recall": 0.6, + "F1": 0.545, + "n_truth": 5, + "n_predicted": 6, + }, + ] + summary = aggregate(rows) + assert len(summary) == 1 + entry = summary[0] + assert entry["runtime_median_s"] == pytest.approx(3.0) + assert (entry["TP"], entry["FP"], entry["FN"]) == (7, 4, 3) + table = markdown_table(summary, max_lag=2, n_shuffles=30, pc_alpha=0.05) + assert "| m | 5 | 200 | 2 | 3.0 | 7 | 4 | 3 |" in table From 7ae2ffd74b7305335679dff2b8bf56b1b84d6be4 Mon Sep 17 00:00:00 2001 From: Vivek Pandit Date: Mon, 12 Oct 2026 16:01:57 -0700 Subject: [PATCH 2/2] benchmark: address PR review (continue CE-only, packaged-install guard, self-lag split) 1. run_profile no longer returns early when tigramite is missing; it warns once and continues CE-only for all remaining cases. 2. test_benchmark_harness.py skips via importorskip when the top-level benchmarks package is not installed. 3. score_edges/aggregate/markdown report cross-links as the main metric with self-lags separately (standard oCSE conditions on the target's own lags and cannot discover them). Regenerated RESULTS_issue18.md from a fresh full-profile run. --- benchmarks/README.md | 6 +- benchmarks/RESULTS_issue18.md | 18 ++-- benchmarks/metrics.py | 41 ++++++- benchmarks/run.py | 60 ++++++++--- .../tests/test_benchmark_harness.py | 100 ++++++++++++------ 5 files changed, 161 insertions(+), 64 deletions(-) diff --git a/benchmarks/README.md b/benchmarks/README.md index f3ef025..6e1b92b 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -17,8 +17,10 @@ python -m benchmarks.run --profile full --output benchmark-results With `p=0.2`, `n=20` yields ~77 true edges, matching the issue report. * **Edge matching:** exact `(source, target, lag)`; all ground-truth edges are lag 1, so lag-2 discoveries count as false positives. -* **Self-lags:** included in the main counts (the generator produces no - self-loops, so any self-lag discovery is a false positive). +* **Self-lags:** reported separately from cross-links. Standard oCSE + conditions on the target's own lags, so it cannot discover self-lags + by construction; counting them in the main score would penalise it on + generators with autoregression. Cross-links are the main metric. * **CE config:** `method="standard"`, Gaussian, serial (`n_jobs=1`), `alpha_forward=alpha_backward=0.05`, profile-dependent `n_shuffles`. * **PCMCI config:** `ParCorr`, `tau_max=max_lag`, `pc_alpha=0.05`. diff --git a/benchmarks/RESULTS_issue18.md b/benchmarks/RESULTS_issue18.md index 440dbe3..9e565e2 100644 --- a/benchmarks/RESULTS_issue18.md +++ b/benchmarks/RESULTS_issue18.md @@ -2,13 +2,13 @@ Generator: `linear_stochastic_gaussian_process(rho=0.7, p=0.2)`, `max_lag`/`tau_max`=2. CE: Gaussian oCSE (`method='standard'`), serial (`n_jobs=1`), `n_shuffles=200`. PCMCI: `ParCorr`, `pc_alpha=0.05`. -Matching: exact `(source, target, lag)`; runtime is the median across seeds, accuracy micro-averaged (TP/FP/FN summed). +Matching: exact `(source, target, lag)`. Cross-links are the main metric; self-lags are reported separately because standard oCSE conditions on the target's own lags and cannot discover them. Runtime is the median across seeds, accuracy micro-averaged (TP/FP/FN summed). -| method | n | T | seeds | runtime median s | TP | FP | FN | precision | recall | F1 | -|---|---|---|---|---|---|---|---|---|---|---| -| causationentropy-gaussian-oCSE | 5 | 500 | 5 | 1.6 | 14 | 7 | 3 | 0.667 | 0.824 | 0.737 | -| tigramite-pcmci-parcorr | 5 | 500 | 5 | 0.0 | 14 | 11 | 3 | 0.560 | 0.824 | 0.667 | -| causationentropy-gaussian-oCSE | 10 | 500 | 5 | 7.0 | 69 | 32 | 3 | 0.683 | 0.958 | 0.798 | -| tigramite-pcmci-parcorr | 10 | 500 | 5 | 0.2 | 68 | 47 | 4 | 0.591 | 0.944 | 0.727 | -| causationentropy-gaussian-oCSE | 20 | 500 | 5 | 31.1 | 329 | 135 | 55 | 0.709 | 0.857 | 0.776 | -| tigramite-pcmci-parcorr | 20 | 500 | 5 | 0.9 | 324 | 190 | 60 | 0.630 | 0.844 | 0.722 | +| method | n | T | seeds | runtime median s | cross TP | cross FP | cross FN | cross precision | cross recall | cross F1 | self TP | self FP | self FN | +|---|---|---|---|---|---|---|---|---|---|---|---|---|---| +| causationentropy-gaussian-oCSE | 5 | 500 | 5 | 1.6 | 14 | 7 | 3 | 0.667 | 0.824 | 0.737 | 0 | 0 | 0 | +| tigramite-pcmci-parcorr | 5 | 500 | 5 | 0.0 | 14 | 8 | 3 | 0.636 | 0.824 | 0.718 | 0 | 3 | 0 | +| causationentropy-gaussian-oCSE | 10 | 500 | 5 | 6.9 | 69 | 32 | 3 | 0.683 | 0.958 | 0.798 | 0 | 0 | 0 | +| tigramite-pcmci-parcorr | 10 | 500 | 5 | 0.2 | 68 | 45 | 4 | 0.602 | 0.944 | 0.735 | 0 | 2 | 0 | +| causationentropy-gaussian-oCSE | 20 | 500 | 5 | 31.9 | 329 | 135 | 55 | 0.709 | 0.857 | 0.776 | 0 | 0 | 0 | +| tigramite-pcmci-parcorr | 20 | 500 | 5 | 0.9 | 324 | 179 | 60 | 0.644 | 0.844 | 0.731 | 0 | 11 | 0 | diff --git a/benchmarks/metrics.py b/benchmarks/metrics.py index 916c133..dfe2122 100644 --- a/benchmarks/metrics.py +++ b/benchmarks/metrics.py @@ -1,8 +1,9 @@ """Edge-set normalization and graph-recovery metrics for issue #18. Canonical rule (see ``benchmarks/README.md``): a true positive requires an -exact ``(source, target, lag)`` match. Self-lag links are included in the -main counts. +exact ``(source, target, lag)`` match. Cross-links are the main metric; +self-lags are reported separately because the standard oCSE method +conditions on the target's own lags and cannot discover them. """ from __future__ import annotations @@ -47,14 +48,18 @@ def normalize_pcmci_edges(graph) -> set: def score_edges(truth: set, predicted: set) -> dict: - """Score ``predicted`` against ``truth`` with exact-match counting.""" + """Score ``predicted`` against ``truth`` with exact-match counting. + + Returns overall counts plus a cross-link / self-lag breakdown + (``cross_*`` is the main metric; ``self_*`` is reported separately). + """ tp = len(truth & predicted) fp = len(predicted - truth) fn = len(truth - predicted) precision = tp / (tp + fp) if (tp + fp) else 0.0 recall = tp / (tp + fn) if (tp + fn) else 0.0 f1 = 2 * precision * recall / (precision + recall) if (precision + recall) else 0.0 - return { + out = { "TP": tp, "FP": fp, "FN": fn, @@ -64,3 +69,31 @@ def score_edges(truth: set, predicted: set) -> dict: "n_truth": len(truth), "n_predicted": len(predicted), } + truth_self = {e for e in truth if e[0] == e[1]} + predicted_self = {e for e in predicted if e[0] == e[1]} + cross_truth = truth - truth_self + cross_predicted = predicted - predicted_self + cross_tp = len(cross_truth & cross_predicted) + cross_fp = len(cross_predicted - cross_truth) + cross_fn = len(cross_truth - cross_predicted) + cross_precision = cross_tp / (cross_tp + cross_fp) if (cross_tp + cross_fp) else 0.0 + cross_recall = cross_tp / (cross_tp + cross_fn) if (cross_tp + cross_fn) else 0.0 + cross_f1 = ( + 2 * cross_precision * cross_recall / (cross_precision + cross_recall) + if (cross_precision + cross_recall) + else 0.0 + ) + out.update( + { + "cross_TP": cross_tp, + "cross_FP": cross_fp, + "cross_FN": cross_fn, + "cross_precision": cross_precision, + "cross_recall": cross_recall, + "cross_F1": cross_f1, + "self_TP": len(truth_self & predicted_self), + "self_FP": len(predicted_self - truth_self), + "self_FN": len(truth_self - predicted_self), + } + ) + return out diff --git a/benchmarks/run.py b/benchmarks/run.py index 31e4e0a..d7af84b 100644 --- a/benchmarks/run.py +++ b/benchmarks/run.py @@ -68,6 +68,7 @@ def run_profile( from causationentropy.graph.utils import pcmci_to_networkx rows = [] + pcmci_available = not skip_pcmci for n, T in sizes: for seed in seeds: data, truth, node_order = generate_case(n, T, seed) @@ -98,15 +99,16 @@ def run_profile( flush=True, ) - if skip_pcmci: + if not pcmci_available: continue try: results, pcmci_time = run_pcmci_parcorr( data, tau_max=max_lag, pc_alpha=pc_alpha ) except ImportError as exc: - print(f"WARNING: {exc} Skipping PCMCI baseline.", file=sys.stderr) - return rows + print(f"WARNING: {exc} Continuing with CE only.", file=sys.stderr) + pcmci_available = False + continue pcmci_graph = pcmci_to_networkx(results) pcmci_scores = score_edges(truth, normalize_pcmci_edges(pcmci_graph)) row = { @@ -133,8 +135,20 @@ def aggregate(rows): """Aggregate per-seed rows to one row per (method, n, T). Runtime is the median across seeds; accuracy is micro-averaged - (TP/FP/FN summed, then precision/recall/F1 derived). + (TP/FP/FN summed, then precision/recall/F1 derived). Cross-links are + the main metric; self-lag counts are reported separately. """ + + def _prf(tp, fp, fn): + precision = tp / (tp + fp) if (tp + fp) else 0.0 + recall = tp / (tp + fn) if (tp + fn) else 0.0 + f1 = ( + 2 * precision * recall / (precision + recall) + if (precision + recall) + else 0.0 + ) + return precision, recall, f1 + groups = {} for row in rows: groups.setdefault((row["method"], row["n"], row["T"]), []).append(row) @@ -146,13 +160,11 @@ def aggregate(rows): tp = sum(r["TP"] for r in group) fp = sum(r["FP"] for r in group) fn = sum(r["FN"] for r in group) - precision = tp / (tp + fp) if (tp + fp) else 0.0 - recall = tp / (tp + fn) if (tp + fn) else 0.0 - f1 = ( - 2 * precision * recall / (precision + recall) - if (precision + recall) - else 0.0 - ) + precision, recall, f1 = _prf(tp, fp, fn) + cross_tp = sum(r["cross_TP"] for r in group) + cross_fp = sum(r["cross_FP"] for r in group) + cross_fn = sum(r["cross_FN"] for r in group) + cross_precision, cross_recall, cross_f1 = _prf(cross_tp, cross_fp, cross_fn) summary.append( { "method": method, @@ -169,6 +181,15 @@ def aggregate(rows): "precision": precision, "recall": recall, "F1": f1, + "cross_TP": cross_tp, + "cross_FP": cross_fp, + "cross_FN": cross_fn, + "cross_precision": cross_precision, + "cross_recall": cross_recall, + "cross_F1": cross_f1, + "self_TP": sum(r["self_TP"] for r in group), + "self_FP": sum(r["self_FP"] for r in group), + "self_FN": sum(r["self_FN"] for r in group), } ) return summary @@ -182,18 +203,23 @@ def markdown_table(summary, max_lag: int, n_shuffles: int, pc_alpha: float) -> s f"`max_lag`/`tau_max`={max_lag}.", f"CE: Gaussian oCSE (`method='standard'`), serial (`n_jobs=1`), " f"`n_shuffles={n_shuffles}`. PCMCI: `ParCorr`, `pc_alpha={pc_alpha}`.", - "Matching: exact `(source, target, lag)`; runtime is the median " - "across seeds, accuracy micro-averaged (TP/FP/FN summed).", + "Matching: exact `(source, target, lag)`. Cross-links are the main " + "metric; self-lags are reported separately because standard oCSE " + "conditions on the target's own lags and cannot discover them. " + "Runtime is the median across seeds, accuracy micro-averaged " + "(TP/FP/FN summed).", "", - "| method | n | T | seeds | runtime median s | TP | FP | FN | precision | recall | F1 |", - "|---|---|---|---|---|---|---|---|---|---|---|", + "| method | n | T | seeds | runtime median s | cross TP | cross FP | cross FN | cross precision | cross recall | cross F1 | self TP | self FP | self FN |", + "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|", ] for entry in summary: lines.append( f"| {entry['method']} | {entry['n']} | {entry['T']} | " f"{entry['n_seeds']} | {entry['runtime_median_s']:.1f} | " - f"{entry['TP']} | {entry['FP']} | {entry['FN']} | " - f"{entry['precision']:.3f} | {entry['recall']:.3f} | {entry['F1']:.3f} |" + f"{entry['cross_TP']} | {entry['cross_FP']} | {entry['cross_FN']} | " + f"{entry['cross_precision']:.3f} | {entry['cross_recall']:.3f} | " + f"{entry['cross_F1']:.3f} | " + f"{entry['self_TP']} | {entry['self_FP']} | {entry['self_FN']} |" ) return "\n".join(lines) + "\n" diff --git a/causationentropy/tests/test_benchmark_harness.py b/causationentropy/tests/test_benchmark_harness.py index 5c0c1c7..a85d02c 100644 --- a/causationentropy/tests/test_benchmark_harness.py +++ b/causationentropy/tests/test_benchmark_harness.py @@ -3,9 +3,11 @@ import networkx as nx import pytest +pytest.importorskip("benchmarks") + from benchmarks.adapters import generate_case, run_causation_entropy from benchmarks.metrics import normalize_ce_edges, score_edges -from benchmarks.run import aggregate, markdown_table +from benchmarks.run import aggregate, markdown_table, run_profile def test_generate_case_truth_matches_graph(): @@ -23,6 +25,35 @@ def test_score_edges_exact_match(): assert out["recall"] == pytest.approx(0.5) +def test_score_edges_self_lag_breakdown(): + truth = {(0, 0, 1), (0, 1, 1), (1, 2, 1)} + predicted = {(0, 0, 1), (0, 1, 1), (2, 2, 1)} + out = score_edges(truth, predicted) + assert (out["cross_TP"], out["cross_FP"], out["cross_FN"]) == (1, 0, 1) + assert (out["self_TP"], out["self_FP"], out["self_FN"]) == (1, 1, 0) + assert out["cross_TP"] + out["self_TP"] == out["TP"] + + +def test_run_profile_continues_without_tigramite(monkeypatch): + def missing_tigramite(*args, **kwargs): + raise ImportError("tigramite is required for the PCMCI baseline.") + + import benchmarks.run as bench_run + + monkeypatch.setattr(bench_run, "run_pcmci_parcorr", missing_tigramite) + rows = run_profile( + sizes=[(3, 60), (4, 60)], + seeds=[0, 1], + max_lag=1, + n_shuffles=5, + pc_alpha=0.05, + skip_pcmci=False, + verbose=False, + ) + assert len(rows) == 4 + assert all(r["method"] == "causationentropy-gaussian-oCSE" for r in rows) + + def test_tiny_end_to_end_ce_only(): data, truth, node_order = generate_case(3, 60, seed=0) graph, runtime = run_causation_entropy(data, max_lag=1, seed=0, n_shuffles=5) @@ -47,43 +78,48 @@ def test_pcmci_adapter_smoke(): assert out["TP"] + out["FN"] == len(truth) +def _row( + seed, runtime, tp, fp, fn, cross_tp, cross_fp, cross_fn, self_tp, self_fp, self_fn +): + return { + "method": "m", + "n": 5, + "T": 200, + "seed": seed, + "runtime_s": runtime, + "TP": tp, + "FP": fp, + "FN": fn, + "precision": 0.0, + "recall": 0.0, + "F1": 0.0, + "n_truth": tp + fn, + "n_predicted": tp + fp, + "cross_TP": cross_tp, + "cross_FP": cross_fp, + "cross_FN": cross_fn, + "cross_precision": 0.0, + "cross_recall": 0.0, + "cross_F1": 0.0, + "self_TP": self_tp, + "self_FP": self_fp, + "self_FN": self_fn, + } + + def test_aggregate_and_markdown(): rows = [ - { - "method": "m", - "n": 5, - "T": 200, - "seed": 0, - "runtime_s": 2.0, - "TP": 4, - "FP": 1, - "FN": 1, - "precision": 0.8, - "recall": 0.8, - "F1": 0.8, - "n_truth": 5, - "n_predicted": 5, - }, - { - "method": "m", - "n": 5, - "T": 200, - "seed": 1, - "runtime_s": 4.0, - "TP": 3, - "FP": 3, - "FN": 2, - "precision": 0.5, - "recall": 0.6, - "F1": 0.545, - "n_truth": 5, - "n_predicted": 6, - }, + _row(0, 2.0, 4, 1, 1, 4, 1, 1, 0, 0, 0), + _row(1, 4.0, 3, 3, 2, 2, 2, 2, 1, 1, 0), ] summary = aggregate(rows) assert len(summary) == 1 entry = summary[0] assert entry["runtime_median_s"] == pytest.approx(3.0) assert (entry["TP"], entry["FP"], entry["FN"]) == (7, 4, 3) + assert (entry["cross_TP"], entry["cross_FP"], entry["cross_FN"]) == (6, 3, 3) + assert (entry["self_TP"], entry["self_FP"], entry["self_FN"]) == (1, 1, 0) + assert entry["cross_precision"] == pytest.approx(6 / 9) table = markdown_table(summary, max_lag=2, n_shuffles=30, pc_alpha=0.05) - assert "| m | 5 | 200 | 2 | 3.0 | 7 | 4 | 3 |" in table + assert "| m | 5 | 200 | 2 | 3.0 | 6 | 3 | 3 |" in table + assert "self TP" in table