Can pyPTE reconstruct a known directed network of oscillators?

Two ground-truth topologies are simulated with delayed Kuramoto oscillators: a directed ring, where every node drives its successor, and two internally dense groups joined by a single one-way bridge.

Ground truth is sparse but dPTE gives every ordered pair a value, so the recovered matrix is dense no matter what. Scored with AUC rather than a hand-picked threshold: the chance a true edge outranks a non-edge.

Run with:

uv run python -m examples.kuramoto_network [--quick]

simulating: Directed ring (i -> i+1) ...
simulating: Two groups, one-way bridges ...

results
  Directed ring (i -> i+1)
    AUC, directed vs uncoupled   0.969   (0.5 = chance)
    mean dPTE on directed edges  0.755
    mean dPTE on uncoupled pairs 0.458
  Two groups, one-way bridges
    AUC, directed vs uncoupled   0.897   (0.5 = chance)
    mean dPTE on directed edges  0.592
    mean dPTE on uncoupled pairs 0.491
    mean dPTE on symmetric pairs 0.500   (expected ~0.5: coupled, but no net direction)

verified: directed edges outrank uncoupled pairs, while symmetrically
          coupled pairs correctly sit at 0.5. dPTE measures direction,
          not the presence of coupling, so a strong bidirectional link
          is invisible to it by design.

          The recovered matrix is dense while the truth is sparse, so
          turning these scores into a graph needs a significance test,
          not a threshold picked by eye.

import argparse

import matplotlib.pyplot as plt
import numpy as np
import numpy.typing as npt

from examples.models import kuramoto
from examples.utils import style
from examples.utils.plotting import FIGURES, plot_matrix, save, use_house_style
from pyPTE import PTE


def auc(scores: npt.NDArray, truth: npt.NDArray) -> float:
    """Area under the ROC curve, computed by rank comparison.

    Equivalent to the Mann-Whitney U statistic: the chance that a randomly
    chosen true edge scores above a randomly chosen non-edge. 0.5 is chance,
    1.0 is perfect separation.
    """
    positive = scores[truth]
    negative = scores[~truth]
    if positive.size == 0 or negative.size == 0:
        return float("nan")
    order = np.argsort(np.concatenate([positive, negative]), kind="mergesort")
    ranks = np.empty(order.size, dtype=float)
    ranks[order] = np.arange(1, order.size + 1)
    # average ranks over ties so exact duplicates do not bias the statistic
    values = np.concatenate([positive, negative])
    for value in np.unique(values):
        tied = values == value
        if tied.sum() > 1:
            ranks[tied] = ranks[tied].mean()
    rank_sum = ranks[: positive.size].sum()
    return float(
        (rank_sum - positive.size * (positive.size + 1) / 2)
        / (positive.size * negative.size)
    )


def off_diagonal(matrix: npt.NDArray) -> npt.NDArray:
    """Flatten a matrix, dropping the diagonal, which carries no information."""
    mask = ~np.eye(matrix.shape[0], dtype=bool)
    return matrix[mask]


def evaluate(coupling: npt.NDArray, *, n_samples: int, seed: int) -> dict:
    signal = kuramoto.simulate(coupling, n_samples=n_samples, seed=seed)
    dPTE, _ = PTE(signal)

    directed = kuramoto.net_flow(coupling)
    # symmetric pairs are coupled but carry no net direction, so they belong to
    # neither class: scoring them as edges would punish the correct answer
    symmetric = (coupling > 0) & (coupling == coupling.T)

    considered = off_diagonal(~symmetric)
    truth = off_diagonal(directed)[considered]
    scores = off_diagonal(dPTE)[considered]

    return {
        "dPTE": dPTE,
        "truth": directed,
        "symmetric": symmetric,
        "auc": auc(scores, truth),
        "edge_mean": scores[truth].mean(),
        "nonedge_mean": scores[~truth].mean(),
        "symmetric_mean": off_diagonal(dPTE)[off_diagonal(symmetric)].mean()
        if symmetric.any()
        else float("nan"),
    }


def main(quick: bool = False) -> None:
    use_house_style()
    n_samples = 4000 if quick else 12000

    topologies = {
        "Directed ring (i -> i+1)": kuramoto.ring_with_shortcut(8, strength=3.0),
        "Two groups, one-way bridges": kuramoto.two_groups(4, within=2.0, across=3.0),
    }

    results = {}
    for name, coupling in topologies.items():
        print(f"simulating: {name} ...")
        results[name] = evaluate(coupling, n_samples=n_samples, seed=1)

    fig, axes = plt.subplots(
        2, 3, figsize=(12, 7.2), constrained_layout=True, width_ratios=[1, 1, 1.15]
    )

    for row, (name, res) in enumerate(results.items()):
        plot_matrix(
            axes[row, 0],
            res["truth"].astype(float),
            title=f"Ground truth\n{name}",
            cmap="Greys",
            vmin=0,
            vmax=1,
            colorbar=False,
        )
        plot_matrix(
            axes[row, 1],
            res["dPTE"],
            title=f"Recovered dPTE\nAUC = {res['auc']:.3f}",
            center_on_half=True,
            mask_diagonal=True,
        )

        ax = axes[row, 2]
        flat_dPTE = off_diagonal(res["dPTE"])
        flat_directed = off_diagonal(res["truth"])
        flat_symmetric = off_diagonal(res["symmetric"])
        bins = np.linspace(0.2, 0.8, 25)

        ax.hist(
            flat_dPTE[~flat_directed & ~flat_symmetric],
            bins=bins,
            alpha=0.65,
            label="no coupling",
            color=style.NULL,
        )
        if flat_symmetric.any():
            ax.hist(
                flat_dPTE[flat_symmetric],
                bins=bins,
                alpha=0.7,
                label="symmetric (no net direction)",
                color=style.AMBER,
            )
        ax.hist(
            flat_dPTE[flat_directed],
            bins=bins,
            alpha=0.8,
            label="directed edge",
            color=style.OBS,
        )
        ax.axvline(0.5, ls="--", lw=1, color="0.4")
        ax.set(xlabel="dPTE", ylabel="count", title="Score separation")
        ax.legend(fontsize=7)

    save(fig, FIGURES / "kuramoto_network.png")

    print("\nresults")
    for name, res in results.items():
        print(f"  {name}")
        print(f"    AUC, directed vs uncoupled   {res['auc']:.3f}   (0.5 = chance)")
        print(f"    mean dPTE on directed edges  {res['edge_mean']:.3f}")
        print(f"    mean dPTE on uncoupled pairs {res['nonedge_mean']:.3f}")
        if not np.isnan(res["symmetric_mean"]):
            print(
                f"    mean dPTE on symmetric pairs {res['symmetric_mean']:.3f}"
                "   (expected ~0.5: coupled, but no net direction)"
            )

    for name, res in results.items():
        assert res["auc"] > 0.75, f"{name}: directed edges should outrank non-edges"
        assert res["edge_mean"] > res["nonedge_mean"], (
            f"{name}: directed edges must score above uncoupled pairs"
        )
        if not np.isnan(res["symmetric_mean"]):
            assert abs(res["symmetric_mean"] - 0.5) < 0.05, (
                f"{name}: symmetrically coupled pairs carry no net direction and "
                "should sit at 0.5"
            )

    print(
        "\nverified: directed edges outrank uncoupled pairs, while symmetrically"
        "\n          coupled pairs correctly sit at 0.5. dPTE measures direction,"
        "\n          not the presence of coupling, so a strong bidirectional link"
        "\n          is invisible to it by design."
        "\n\n          The recovered matrix is dense while the truth is sparse, so"
        "\n          turning these scores into a graph needs a significance test,"
        "\n          not a threshold picked by eye."
    )


if __name__ == "__main__":
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--quick", action="store_true", help="fewer samples")
    main(**vars(parser.parse_args()))