benchmark: add reproducible CE vs PCMCI harness and issue #18 table

.gitignore

@@ -83,3 +83,6 @@ venv.bak/
 .Trashes
 ehthumbs.db
 Thumbs.db
+
+# Local benchmark artifacts (issue #18)
+benchmark-results/

benchmarks/README.md

@@ -0,0 +1,49 @@
+# 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:** 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`.
+  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`.

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)`. 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 | 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 |

benchmarks/__init__.py

@@ -0,0 +1,4 @@
+"""Reproducible runtime/accuracy benchmarks (issue #18).
+
+See ``benchmarks/README.md`` for conventions and usage.
+"""

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

benchmarks/metrics.py

@@ -0,0 +1,99 @@
+"""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. 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
+
+
+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.
+
+    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
+    out = {
+        "TP": tp,
+        "FP": fp,
+        "FN": fn,
+        "precision": precision,
+        "recall": recall,
+        "F1": f1,
+        "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