|
| 1 | +# SPDX-License-Identifier: AGPL-3.0-or-later | Commercial license available |
| 2 | +# © Concepts 1996–2026 Miroslav Šotek. All rights reserved. |
| 3 | +# © Code 2020–2026 Miroslav Šotek. All rights reserved. |
| 4 | +# ORCID: 0009-0009-3560-0851 |
| 5 | +# Contact: www.anulum.li | protoscience@anulum.li |
| 6 | +# SCPN Phase Orchestrator — SPO supervision on neurolib output |
| 7 | +# |
| 8 | +# Demonstrates SPO's unique value: regime detection, TCBO consciousness |
| 9 | +# gate, and predictive supervision on top of neurolib's neural dynamics. |
| 10 | + |
| 11 | +"""Feed neurolib ALN output into SPO's supervision layer. |
| 12 | +
|
| 13 | +Shows what neurolib CAN'T do: regime classification, boundary monitoring, |
| 14 | +TCBO consciousness boundary, predictive control. |
| 15 | +
|
| 16 | +Usage: |
| 17 | + python experiments/regime_on_neurolib.py |
| 18 | +""" |
| 19 | + |
| 20 | +from __future__ import annotations |
| 21 | + |
| 22 | +import json |
| 23 | +import time |
| 24 | + |
| 25 | +import numpy as np |
| 26 | +from scipy.signal import hilbert |
| 27 | + |
| 28 | +from neurolib.models.aln import ALNModel |
| 29 | +from neurolib.utils.loadData import Dataset |
| 30 | + |
| 31 | +from scpn_phase_orchestrator.monitor.boundaries import BoundaryObserver, BoundaryState |
| 32 | +from scpn_phase_orchestrator.monitor.npe import compute_npe |
| 33 | +from scpn_phase_orchestrator.ssgf.tcbo import TCBOObserver |
| 34 | +from scpn_phase_orchestrator.supervisor.events import EventBus |
| 35 | +from scpn_phase_orchestrator.supervisor.regimes import Regime, RegimeManager |
| 36 | +from scpn_phase_orchestrator.upde.metrics import LayerState, UPDEState |
| 37 | +from scpn_phase_orchestrator.upde.order_params import compute_order_parameter |
| 38 | + |
| 39 | + |
| 40 | +def extract_phases_from_rates(rates: np.ndarray, window: int = 100) -> np.ndarray: |
| 41 | + """Extract instantaneous phases from neural firing rates via Hilbert.""" |
| 42 | + n_regions, n_t = rates.shape |
| 43 | + phases = np.zeros((n_regions, n_t)) |
| 44 | + for i in range(n_regions): |
| 45 | + analytic = hilbert(rates[i]) |
| 46 | + phases[i] = np.angle(analytic) % (2 * np.pi) |
| 47 | + return phases |
| 48 | + |
| 49 | + |
| 50 | +def main() -> None: |
| 51 | + print("Loading HCP + running neurolib ALN...") |
| 52 | + ds = Dataset("hcp") |
| 53 | + model = ALNModel(Cmat=ds.Cmat, Dmat=ds.Dmat) |
| 54 | + model.params["duration"] = 30000 # 30s |
| 55 | + model.params["Ke_gl"] = 2.0 # optimal from baseline |
| 56 | + |
| 57 | + t0 = time.perf_counter() |
| 58 | + model.run() |
| 59 | + sim_time = time.perf_counter() - t0 |
| 60 | + print(f" neurolib simulation: {sim_time:.1f}s") |
| 61 | + |
| 62 | + rates = model.rates_exc # (80, T) |
| 63 | + print(f" rates shape: {rates.shape}") |
| 64 | + |
| 65 | + # Extract phases |
| 66 | + phases = extract_phases_from_rates(rates) |
| 67 | + n_regions, n_t = phases.shape |
| 68 | + |
| 69 | + # Downsample for SPO processing (every 100 steps = 10ms windows) |
| 70 | + ds_factor = 100 |
| 71 | + n_windows = n_t // ds_factor |
| 72 | + |
| 73 | + # SPO supervision components |
| 74 | + event_bus = EventBus() |
| 75 | + regime_manager = RegimeManager(event_bus=event_bus) |
| 76 | + tcbo = TCBOObserver(tau_h1=0.72, window_size=30, embed_dim=2, embed_delay=1) |
| 77 | + |
| 78 | + regime_history = [] |
| 79 | + R_history = [] |
| 80 | + npe_history = [] |
| 81 | + tcbo_history = [] |
| 82 | + |
| 83 | + print(f"\nRunning SPO supervision on {n_windows} windows...") |
| 84 | + for w in range(n_windows): |
| 85 | + start = w * ds_factor |
| 86 | + end = start + ds_factor |
| 87 | + window_phases = phases[:, end - 1] # phase snapshot at window end |
| 88 | + |
| 89 | + # Kuramoto R across all regions |
| 90 | + R, psi = compute_order_parameter(window_phases) |
| 91 | + R_history.append(float(R)) |
| 92 | + |
| 93 | + # NPE |
| 94 | + npe = compute_npe(window_phases) |
| 95 | + npe_history.append(float(npe)) |
| 96 | + |
| 97 | + # TCBO |
| 98 | + tcbo_state = tcbo.observe(window_phases) |
| 99 | + tcbo_history.append(tcbo_state.p_h1) |
| 100 | + |
| 101 | + # Regime classification |
| 102 | + layer_states = [LayerState(R=R, psi=psi)] |
| 103 | + upde_state = UPDEState( |
| 104 | + layers=layer_states, |
| 105 | + cross_layer_alignment=np.eye(1), |
| 106 | + stability_proxy=R, |
| 107 | + regime_id=regime_manager.current_regime.value, |
| 108 | + ) |
| 109 | + proposed = regime_manager.evaluate(upde_state, BoundaryState()) |
| 110 | + regime_manager.transition(proposed) |
| 111 | + regime_history.append(regime_manager.current_regime.value) |
| 112 | + |
| 113 | + # Report |
| 114 | + regime_counts = {} |
| 115 | + for r in regime_history: |
| 116 | + regime_counts[r] = regime_counts.get(r, 0) + 1 |
| 117 | + |
| 118 | + transitions = len(regime_manager.transition_history) |
| 119 | + |
| 120 | + print(f"\n{'='*60}") |
| 121 | + print(f"SPO Supervision Results on neurolib ALN output") |
| 122 | + print(f"{'='*60}") |
| 123 | + print(f"Windows analyzed: {n_windows}") |
| 124 | + print(f"R range: [{min(R_history):.4f}, {max(R_history):.4f}]") |
| 125 | + print(f"R mean: {np.mean(R_history):.4f}") |
| 126 | + print(f"NPE range: [{min(npe_history):.4f}, {max(npe_history):.4f}]") |
| 127 | + print(f"TCBO p_h1 mean: {np.mean(tcbo_history):.4f}") |
| 128 | + print(f"Regime distribution: {regime_counts}") |
| 129 | + print(f"Regime transitions: {transitions}") |
| 130 | + print(f"Events logged: {event_bus.count}") |
| 131 | + |
| 132 | + results = { |
| 133 | + "n_windows": n_windows, |
| 134 | + "R_mean": round(float(np.mean(R_history)), 4), |
| 135 | + "R_std": round(float(np.std(R_history)), 4), |
| 136 | + "R_min": round(float(min(R_history)), 4), |
| 137 | + "R_max": round(float(max(R_history)), 4), |
| 138 | + "NPE_mean": round(float(np.mean(npe_history)), 4), |
| 139 | + "TCBO_p_h1_mean": round(float(np.mean(tcbo_history)), 4), |
| 140 | + "regime_counts": regime_counts, |
| 141 | + "n_transitions": transitions, |
| 142 | + "n_events": event_bus.count, |
| 143 | + "sim_duration_ms": 30000, |
| 144 | + "K": 2.0, |
| 145 | + } |
| 146 | + |
| 147 | + with open("experiments/regime_on_neurolib_results.json", "w") as f: |
| 148 | + json.dump(results, f, indent=2) |
| 149 | + print(f"\nSaved to experiments/regime_on_neurolib_results.json") |
| 150 | + |
| 151 | + |
| 152 | +if __name__ == "__main__": |
| 153 | + main() |
0 commit comments