[Vivian0108]

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.

[Vivian0108]

@kslote1 @chaewon-research

[kslote1]

Hi everyone!

This addition is great.

I get the sense that once this addition is merged, it will become our benchmark harness going forward. Is everyone happy with the design?

@jefish003 @chaewon-research @Parva-26 @VictorLiwentao

[chaewon-research]

The design looks good to me. The separation between the adapters, metrics, and result generation should make it easy to reuse for #55. I’d keep the k/T sweep as a separate benchmark entry point, but I don’t think #56 needs any changes for that.

[Parva-26]

@kslote1 the design looks good to me. Keeping the generator, adapters, metrics and runner separate makes it easy to reuse for #55. I ran it on ad6d35e and found three things worth fixing before it becomes the standard harness:

1. CE-only runs stop early without tigramite. In run_profile, the except ImportError branch does return rows, so the loop ends after the first CE case. With sizes=[(5, 200), (10, 200)] and seeds 0-2 it returns 1 CE row instead of 6, with only a warning. Setting skip_pcmci = True and continuing would keep the CE results.

2. The new test fails outside a source checkout. causationentropy/tests/test_benchmark_harness.py imports the top-level benchmarks package, but pyproject.toml only packages causationentropy*, and causationentropy/tests ships with the package. After pip install ., pytest --pyargs causationentropy.tests.test_benchmark_harness fails at collection with ModuleNotFoundError: No module named 'benchmarks'. A pytest.importorskip("benchmarks") at the top of the file avoids that.

3. Self-lags, for reuse with other generators. The standard method conditions on the target's own lags (Z_init), so it never reports a self-lag, while PCMCI does. This generator has no self-loops, so the current table is fine. On a VAR(1) with autoregression on every node, though, CE found 0/4 self-lags and PCMCI 4/4 in each of 3 seeds, with both getting 3/3 cross-links. Counting self-lags in the main score would penalise CE by construction there, so it may be worth reporting cross-links as the main metric and self-lags separately (question 3 in nil-vn's list on Performance, and the missing comparison to prior implementations #18).

Script for 1 and 3 (run from the repo root of this branch)

# Checks for #56. Run from the repo root of the PR branch: python pr56_review.py
import contextlib
import io

import numpy as np

import benchmarks.run as bench_run
from benchmarks.adapters import run_causation_entropy, run_pcmci_parcorr
from benchmarks.metrics import normalize_ce_edges, normalize_pcmci_edges
from causationentropy.graph.utils import pcmci_to_networkx


# [1] run_profile when tigramite can't be imported
def missing_tigramite(*args, **kwargs):
    raise ImportError("tigramite is required for the PCMCI baseline.")


bench_run.run_pcmci_parcorr = missing_tigramite
with contextlib.redirect_stdout(io.StringIO()), contextlib.redirect_stderr(io.StringIO()):
    rows = bench_run.run_profile(
        sizes=[(5, 200), (10, 200)],
        seeds=[0, 1, 2],
        max_lag=2,
        n_shuffles=10,
        pc_alpha=0.05,
        skip_pcmci=False,
        verbose=False,
    )
print(f"[1] without tigramite: {len(rows)} CE row(s) for 6 requested cases")

# [2] self-lags: VAR(1), autoregression 0.6 on every node, chain 0 -> 1 -> 2 -> 3
n, T = 4, 500
A = np.diag([0.6] * n)  # A[i, j] != 0 means j drives i
A[1, 0] = A[2, 1] = A[3, 2] = 0.4
truth = {(j, i, 1) for i in range(n) for j in range(n) if A[i, j] != 0}
for seed in range(3):
    rng = np.random.default_rng(seed)
    X = np.zeros((T, n))
    for t in range(1, T):
        X[t] = A @ X[t - 1] + 0.1 * rng.standard_normal(n)
    with contextlib.redirect_stdout(io.StringIO()):
        G, _ = run_causation_entropy(X, max_lag=1, seed=seed, n_shuffles=50)
    ce = normalize_ce_edges(G, [f"X{k}" for k in range(n)])
    res, _ = run_pcmci_parcorr(X, tau_max=1)
    pc = normalize_pcmci_edges(pcmci_to_networkx(res))
    for name, found in (("CE", ce), ("PCMCI", pc)):
        self_tp = sum(1 for e in found & truth if e[0] == e[1])
        cross_tp = sum(1 for e in found & truth if e[0] != e[1])
        print(f"[2] seed {seed} {name:5s}: self-lags {self_tp}/4, cross-links {cross_tp}/3")

I used an AI assistant's help to draft this comment. I have personally checked the details, but I thought it's better to just inform you beforehand.

[Vivian0108]

Thanks @Parva-26 for the review, all three fixed and verified:

1. Missing tigramite now warns once and continues CE-only (was returning 1/6 rows, now 6/6; regression test added).
2. Test file guarded with pytest.importorskip("benchmarks") for installed-package runs.
3. Cross links are now the main metric with self-lags in separate columns. Table regenerated from a fresh full run. discovery.py always conditions on the target's own lags, and the new table shows PCMCI's self-lag FPs separately.
   @chaewon-research — no changes needed for Benchmark kNN discovery sensitivity to k and time-series length #55 reuse. The score_edges breakdown carries over to the k/T sweep.