Add TE-directed adaptive coupling + scaling benchmark (measured)

te_adapt_coupling(): K_ij += lr·TE(i→j) — strengthens coupling along
causal information channels (Lizier 2012).

Scaling benchmark (measured on GTX 1060 host, Python 3.12, numpy):
  N=16:   0.03 ms/step (29K steps/s)
  N=64:   0.09 ms/step (10K steps/s)
  N=256:  7.0  ms/step (143 steps/s)
  N=1000: 47   ms/step (21 steps/s, 16 MB)
Bottleneck: O(N²) coupling matrix multiply.

Co-Authored-By: Arcane Sapience <protoscience@anulum.li>

Author: anulum

benchmarks/scaling_benchmark.py

@@ -0,0 +1,87 @@
+# SPDX-License-Identifier: AGPL-3.0-or-later | Commercial license available
+# © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
+# © Code 2020–2026 Miroslav Šotek. All rights reserved.
+# ORCID: 0009-0009-3560-0851
+# Contact: www.anulum.li | protoscience@anulum.li
+# SCPN Phase Orchestrator — Scaling benchmark
+
+"""Measure performance scaling from N=16 to N=4000.
+
+Usage:
+    python benchmarks/scaling_benchmark.py
+    python benchmarks/scaling_benchmark.py --output benchmarks/scaling_results.json
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import time
+
+import numpy as np
+
+from scpn_phase_orchestrator.upde.engine import UPDEEngine
+from scpn_phase_orchestrator.upde.order_params import compute_order_parameter
+
+
+def bench_one(n: int, n_steps: int = 100, dt: float = 0.01) -> dict:
+    rng = np.random.default_rng(42)
+    phases = rng.uniform(0, 2 * np.pi, n)
+    omegas = np.ones(n)
+    knm = np.full((n, n), 0.5 / n)
+    np.fill_diagonal(knm, 0.0)
+    alpha = np.zeros((n, n))
+
+    engine = UPDEEngine(n, dt=dt)
+
+    t0 = time.perf_counter()
+    for _ in range(n_steps):
+        phases = engine.step(phases, omegas, knm, 0.0, 0.0, alpha)
+    elapsed = time.perf_counter() - t0
+
+    R, _ = compute_order_parameter(phases)
+    mem_mb = (knm.nbytes + alpha.nbytes + phases.nbytes) / 1e6
+
+    return {
+        "N": n,
+        "n_steps": n_steps,
+        "wall_time_s": round(elapsed, 4),
+        "ms_per_step": round(elapsed / n_steps * 1000, 3),
+        "steps_per_sec": round(n_steps / elapsed, 1),
+        "final_R": round(float(R), 4),
+        "memory_MB": round(mem_mb, 2),
+    }
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--output", type=str, default=None)
+    parser.add_argument(
+        "--sizes",
+        type=int,
+        nargs="+",
+        default=[16, 64, 256, 1000],
+    )
+    args = parser.parse_args()
+
+    print(f"{'N':>6} {'ms/step':>10} {'steps/s':>10} {'R':>8} {'mem MB':>8}")
+    print("-" * 50)
+
+    results = []
+    for n in args.sizes:
+        r = bench_one(n)
+        results.append(r)
+        print(
+            f"{r['N']:>6} {r['ms_per_step']:>10.3f} "
+            f"{r['steps_per_sec']:>10.1f} {r['final_R']:>8.4f} "
+            f"{r['memory_MB']:>8.2f}"
+        )
+
+    if args.output:
+        with open(args.output, "w") as f:
+            json.dump(results, f, indent=2)
+        print(f"\nSaved to {args.output}")
+
+
+if __name__ == "__main__":
+    main()

benchmarks/scaling_results.json

@@ -0,0 +1,38 @@
+[
+  {
+    "N": 16,
+    "n_steps": 100,
+    "wall_time_s": 0.0034,
+    "ms_per_step": 0.034,
+    "steps_per_sec": 29155.4,
+    "final_R": 0.2323,
+    "memory_MB": 0.0
+  },
+  {
+    "N": 64,
+    "n_steps": 100,
+    "wall_time_s": 0.0093,
+    "ms_per_step": 0.093,
+    "steps_per_sec": 10798.7,
+    "final_R": 0.1441,
+    "memory_MB": 0.07
+  },
+  {
+    "N": 256,
+    "n_steps": 100,
+    "wall_time_s": 0.7003,
+    "ms_per_step": 7.003,
+    "steps_per_sec": 142.8,
+    "final_R": 0.0237,
+    "memory_MB": 1.05
+  },
+  {
+    "N": 1000,
+    "n_steps": 100,
+    "wall_time_s": 4.6746,
+    "ms_per_step": 46.746,
+    "steps_per_sec": 21.4,
+    "final_R": 0.0254,
+    "memory_MB": 16.01
+  }
+]

src/scpn_phase_orchestrator/coupling/te_adaptive.py

@@ -0,0 +1,48 @@
+# SPDX-License-Identifier: AGPL-3.0-or-later | Commercial license available
+# © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
+# © Code 2020–2026 Miroslav Šotek. All rights reserved.
+# ORCID: 0009-0009-3560-0851
+# Contact: www.anulum.li | protoscience@anulum.li
+# SCPN Phase Orchestrator — Transfer entropy directed adaptive coupling
+
+from __future__ import annotations
+
+import numpy as np
+from numpy.typing import NDArray
+
+from scpn_phase_orchestrator.monitor.transfer_entropy import (
+    transfer_entropy_matrix,
+)
+
+__all__ = ["te_adapt_coupling"]
+
+
+def te_adapt_coupling(
+    knm: NDArray,
+    phase_history: NDArray,
+    lr: float = 0.01,
+    decay: float = 0.0,
+    n_bins: int = 8,
+) -> NDArray:
+    """Adapt coupling matrix using transfer entropy as learning signal.
+
+    K_ij(t+1) = (1-decay) * K_ij(t) + lr * TE(i→j)
+
+    Strengthens coupling along causal information flow channels.
+    Weakens where there is no causal influence.
+
+    Lizier 2012, "Local Information Transfer as a Spatiotemporal Filter
+    for Complex Systems," Physical Review E 77(2):026110.
+
+    Args:
+        knm: current (n, n) coupling matrix.
+        phase_history: (n, T) recent phase trajectories.
+        lr: learning rate for TE-based update.
+        decay: coupling decay rate per update (0 = no decay).
+        n_bins: histogram bins for TE estimation.
+    """
+    te = transfer_entropy_matrix(phase_history, n_bins=n_bins)
+    knm_new = (1.0 - decay) * knm + lr * te
+    np.fill_diagonal(knm_new, 0.0)
+    result: NDArray = np.maximum(knm_new, 0.0)
+    return result

src/scpn_phase_orchestrator/supervisor/predictive.py

@@ -14,8 +14,6 @@
 
 from scpn_phase_orchestrator.actuation.mapper import ControlAction
 from scpn_phase_orchestrator.monitor.boundaries import BoundaryState
-from scpn_phase_orchestrator.supervisor.regimes import Regime
-from scpn_phase_orchestrator.upde.engine import UPDEEngine
 from scpn_phase_orchestrator.upde.metrics import UPDEState
 from scpn_phase_orchestrator.upde.order_params import compute_order_parameter
 from scpn_phase_orchestrator.upde.reduction import OttAntonsenReduction

tests/test_coupling_est.py

@@ -98,19 +98,16 @@ def test_single_harmonic_sign_consistent(self):
         phases = rng.uniform(0, 2 * np.pi, (3, 500))
         omegas = np.ones(3)
         knm_std = estimate_coupling(phases, omegas, dt=0.01)
-        knm_harm = estimate_coupling_harmonics(
-            phases, omegas, dt=0.01, n_harmonics=1
-        )
+        knm_harm = estimate_coupling_harmonics(phases, omegas, dt=0.01, n_harmonics=1)
         # Sign of dominant coupling should match
         mask = np.abs(knm_std) > 0.1
         if np.any(mask):
-            assert np.sum(np.sign(knm_harm["sin_1"][mask]) == np.sign(knm_std[mask])) > 0
+            signs_match = np.sign(knm_harm["sin_1"][mask]) == np.sign(knm_std[mask])
+            assert np.sum(signs_match) > 0
 
     def test_custom_n_harmonics(self):
         rng = np.random.default_rng(42)
         phases = rng.uniform(0, 2 * np.pi, (3, 100))
-        result = estimate_coupling_harmonics(
-            phases, np.ones(3), dt=0.01, n_harmonics=3
-        )
+        result = estimate_coupling_harmonics(phases, np.ones(3), dt=0.01, n_harmonics=3)
         assert "sin_3" in result
         assert "cos_3" in result

tests/test_predictive_supervisor.py

@@ -80,9 +80,7 @@ def test_hard_violation_overrides(self):
         knm = np.eye(4)
         alpha = np.zeros((4, 4))
         bs = BoundaryState(hard_violations=["test"])
-        actions = ps.decide(
-            phases, omegas, knm, alpha, _make_state(0.9), bs
-        )
+        actions = ps.decide(phases, omegas, knm, alpha, _make_state(0.9), bs)
         assert len(actions) == 1
         assert actions[0].knob == "zeta"
 

tests/test_te_adaptive.py

@@ -0,0 +1,55 @@
+# SPDX-License-Identifier: AGPL-3.0-or-later | Commercial license available
+# © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
+# © Code 2020–2026 Miroslav Šotek. All rights reserved.
+# ORCID: 0009-0009-3560-0851
+# Contact: www.anulum.li | protoscience@anulum.li
+# SCPN Phase Orchestrator — TE adaptive coupling tests
+
+from __future__ import annotations
+
+import numpy as np
+
+from scpn_phase_orchestrator.coupling.te_adaptive import te_adapt_coupling
+
+
+class TestTEAdaptive:
+    def test_output_shape(self):
+        knm = np.full((4, 4), 0.5)
+        np.fill_diagonal(knm, 0.0)
+        rng = np.random.default_rng(42)
+        history = rng.uniform(0, 2 * np.pi, (4, 100))
+        result = te_adapt_coupling(knm, history)
+        assert result.shape == (4, 4)
+
+    def test_zero_diagonal(self):
+        knm = np.full((3, 3), 0.5)
+        np.fill_diagonal(knm, 0.0)
+        rng = np.random.default_rng(42)
+        history = rng.uniform(0, 2 * np.pi, (3, 100))
+        result = te_adapt_coupling(knm, history)
+        np.testing.assert_array_equal(np.diag(result), 0.0)
+
+    def test_non_negative(self):
+        knm = np.full((3, 3), 0.1)
+        np.fill_diagonal(knm, 0.0)
+        rng = np.random.default_rng(42)
+        history = rng.uniform(0, 2 * np.pi, (3, 200))
+        result = te_adapt_coupling(knm, history, lr=0.01)
+        assert np.all(result >= 0)
+
+    def test_coupling_increases(self):
+        knm = np.full((4, 4), 0.1)
+        np.fill_diagonal(knm, 0.0)
+        rng = np.random.default_rng(42)
+        history = rng.uniform(0, 2 * np.pi, (4, 200))
+        result = te_adapt_coupling(knm, history, lr=1.0)
+        # With lr=1.0, TE contribution should increase off-diagonal
+        assert float(result.sum()) >= float(knm.sum())
+
+    def test_decay_reduces(self):
+        knm = np.full((3, 3), 1.0)
+        np.fill_diagonal(knm, 0.0)
+        rng = np.random.default_rng(42)
+        history = rng.uniform(0, 2 * np.pi, (3, 100))
+        result = te_adapt_coupling(knm, history, lr=0.0, decay=0.5)
+        assert float(result.sum()) < float(knm.sum())