mechanism: separate BCI role from performance velocity

2 files changed, 91 insertions, 0 deletions

diff --git a/experiments/bci_smoke.py b/experiments/bci_smoke.py
index 5cd5b65..4724d26 100644
--- a/experiments/bci_smoke.py
+++ b/experiments/bci_smoke.py
@@ -93,6 +93,55 @@ def check_phase_masks_and_velocity_reset():
     print("BCI online/plasticity phase masks and episode velocity reset: exact")
 
 
+def check_temporal_difference_role_vectorizer():
+    """Role perturbations times performance change give the Harnett sign."""
+    cfg = BCIConfig(
+        days=2, episodes_per_day=4096, steps_per_episode=2,
+        feedback="performance_velocity", vectorizer_eta=0.5)
+    model = BCISDIL(cfg, model_seed=12)
+    generator = torch.Generator().manual_seed(13)
+    soma = 0.2 * torch.randn(
+        cfg.episodes_per_day, cfg.n_neurons, generator=generator)
+    xi = torch.empty_like(soma).bernoulli_(
+        0.5, generator=generator).mul_(2).sub_(1)
+    target = model.causal_role_targets(soma, xi)
+    mean_target = target.mean(0)
+    cosine = torch.nn.functional.cosine_similarity(
+        mean_target, model.role, dim=0).item()
+    assert cosine > 0.99
+
+    active = torch.ones(cfg.episodes_per_day, dtype=torch.bool)
+    model.A.zero_()
+    model.update_vectorizer(
+        torch.ones(cfg.episodes_per_day, 1),
+        torch.zeros(cfg.episodes_per_day, cfg.n_neurons),
+        target, active)
+    learned_cosine = torch.nn.functional.cosine_similarity(
+        model.A[:, 0], model.role, dim=0).item()
+    assert learned_cosine > 0.99
+
+    # With the neutral predictor exact, r_i = role_i * delta_error.  Hence the
+    # P+ minus P- residual is positive in improving epochs and negative in
+    # worsening epochs, independent of instantaneous error magnitude.
+    model.A[:, 0].copy_(model.role)
+    model.P.copy_(model.coupling)
+    model.P_bias.zero_()
+    previous = torch.tensor([0.8, 0.4, 0.9, 0.3])
+    current = torch.tensor([0.6, 0.5, 0.5, 0.6])
+    delta = model.feedback_features(current, previous)
+    _, innovation, _, _ = model.apical_components(soma[:4], delta)
+    difference = (innovation[:, :cfg.n_plus].mean(1)
+                  - innovation[:, cfg.n_plus:cfg.n_plus + cfg.n_minus].mean(1))
+    improving = delta[:, 0] > 0
+    worsening = delta[:, 0] < 0
+    sign_index = 0.5 * (
+        difference[improving].mean() - difference[worsening].mean())
+    assert sign_index > 0
+    print(
+        "BCI temporal-difference role vectorizer: "
+        f"cos={cosine:.6f}, sign_index={sign_index.item():.6f}")
+
+
 def check_grouped_decoders_and_metrics():
     generator = torch.Generator().manual_seed(9)
     groups = torch.arange(200)
@@ -136,5 +185,6 @@ if __name__ == "__main__":
     check_causal_estimator()
     check_predictor_identification()
     check_phase_masks_and_velocity_reset()
+    check_temporal_difference_role_vectorizer()
     check_grouped_decoders_and_metrics()
     print("ALL BCI MECHANICS CHECKS PASSED")
diff --git a/experiments/verify_theory.py b/experiments/verify_theory.py
index 21a7c06..c00a091 100644
--- a/experiments/verify_theory.py
+++ b/experiments/verify_theory.py
@@ -348,6 +348,46 @@ def check_kolen_pollack_difference_dynamics():
     assert maximum_error < 3e-15
 
 
+def check_bci_error_derivative_structure():
+    """Separate causal role from temporal performance innovation."""
+    rng = np.random.default_rng(641)
+    role = np.concatenate((np.full(5, 0.2), np.full(5, -0.2),
+                           np.zeros(30)))
+    examples = 200_000
+    xi = 2.0 * rng.integers(0, 2, size=(examples, role.size)) - 1.0
+    cursor_direction = xi @ role
+    role_estimate = (cursor_direction[:, None] * xi).mean(axis=0)
+    role_error = float(np.abs(role_estimate - role).max())
+    role_cosine = float(role_estimate @ role / (
+        np.linalg.norm(role_estimate) * np.linalg.norm(role)))
+
+    # Positive target error throughout the BCI trial. Improving events have
+    # smaller current error; worsening events have larger current error.
+    previous = np.asarray((0.8, 0.7, 0.6, 0.5, 0.4, 0.3))
+    current = np.asarray((0.6, 0.5, 0.4, 0.7, 0.6, 0.5))
+    improvement = previous - current
+    improving = improvement > 0
+    worsening = improvement < 0
+
+    def sign_index(modulator):
+        residual = modulator[:, None] * role[None, :]
+        difference = residual[:, :5].mean(1) - residual[:, 5:10].mean(1)
+        return 0.5 * (difference[improving].mean()
+                      - difference[worsening].mean())
+
+    instantaneous_index = float(sign_index(current))
+    temporal_difference_index = float(sign_index(improvement))
+    print("\nBCI ERROR-DERIVATIVE STRUCTURE")
+    print(f"role_estimator_cosine={role_cosine:.9f} "
+          f"max_error={role_error:.3e} "
+          f"instantaneous_index={instantaneous_index:+.6f} "
+          f"td_index={temporal_difference_index:+.6f}")
+    assert role_cosine > 0.999
+    assert role_error < 0.005
+    assert instantaneous_index < 0.0
+    assert temporal_difference_index > 0.0
+
+
 def main():
     check_simultaneous_variance()
     check_sigma_bias()
@@ -359,6 +399,7 @@ def main():
     check_dynamic_neutral_projection()
     check_intermittent_feedback_tracking()
     check_kolen_pollack_difference_dynamics()
+    check_bci_error_derivative_structure()
     print("\nALL THEORY CHECKS PASSED")