feat: add BP-free physical coupled-learning adapter

1 files changed, 246 insertions, 0 deletions

diff --git a/sdil/physical_coupled.py b/sdil/physical_coupled.py
new file mode 100644
index 0000000..0216048
--- /dev/null
+++ b/sdil/physical_coupled.py
@@ -0,0 +1,246 @@
+"""Backpropagation-free simulator for the Dillavou two-edge circuit.
+
+The circuit equations and standard/overclamping updates follow Appendix D/F
+of Dillavou et al. (arXiv:2505.22887v2).  Every adaptive operation is an
+explicit NumPy local rule; this module intentionally has no autodiff path.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import Callable, Iterable, Optional
+
+import numpy as np
+
+
+Array = np.ndarray
+
+
+@dataclass(frozen=True)
+class Circuit:
+    high: float = 0.4351
+    low: float = 0.0181
+    conductance_per_gate: float = 8.5e-4
+    threshold_voltage: float = 0.7
+    fixed_conductance: float = 1.0 / 500.0
+    measured_learning_rate: float = 2040.0
+    integration_step_seconds: float = 2.0e-4
+    gate_minimum: float = 1.0
+    gate_maximum: float = 5.2
+
+
+@dataclass(frozen=True)
+class Task:
+    name: str
+    input_voltage: float
+    label_voltage: float
+
+
+@dataclass(frozen=True)
+class LocalAffineBias:
+    reference_gate: Array
+    bias_at_reference: Array
+    local_slopes: Array
+
+    def __post_init__(self) -> None:
+        for value in (
+            self.reference_gate, self.bias_at_reference, self.local_slopes
+        ):
+            if np.asarray(value).shape != (2,):
+                raise ValueError("two-edge bias arrays must have shape (2,)")
+
+    def __call__(self, gates: Array, strength: float = 1.0) -> Array:
+        gates = np.asarray(gates, dtype=float)
+        if gates.shape != (2,):
+            raise ValueError("gates must have shape (2,)")
+        return strength * (
+            self.bias_at_reference
+            + self.local_slopes * (gates - self.reference_gate)
+        )
+
+
+@dataclass
+class LocalPredictor:
+    """Independent per-edge affine filters trained by normalized LMS."""
+
+    reference_gate: Array
+    feature_scale: Array
+    coefficients: Array
+    affine: bool
+
+    @classmethod
+    def zeros(
+        cls,
+        reference_gate: Array,
+        feature_scale: Array,
+        *,
+        affine: bool,
+    ) -> "LocalPredictor":
+        width = 2 if affine else 1
+        return cls(
+            reference_gate=np.asarray(reference_gate, dtype=float).copy(),
+            feature_scale=np.asarray(feature_scale, dtype=float).copy(),
+            coefficients=np.zeros((2, width), dtype=float),
+            affine=affine,
+        )
+
+    def features(self, gates: Array) -> Array:
+        gates = np.asarray(gates, dtype=float)
+        if gates.shape != (2,):
+            raise ValueError("gates must have shape (2,)")
+        if not self.affine:
+            return np.ones((2, 1), dtype=float)
+        normalized = (gates - self.reference_gate) / self.feature_scale
+        return np.column_stack((np.ones(2, dtype=float), normalized))
+
+    def predict(self, gates: Array) -> Array:
+        return np.sum(self.coefficients * self.features(gates), axis=1)
+
+    def update(self, gates: Array, neutral_measurement: Array, rate: float) -> Array:
+        """One local normalized-LMS update and its pre-update residual."""
+        neutral_measurement = np.asarray(neutral_measurement, dtype=float)
+        if neutral_measurement.shape != (2,):
+            raise ValueError("neutral measurement must have shape (2,)")
+        features = self.features(gates)
+        residual = neutral_measurement - np.sum(
+            self.coefficients * features, axis=1)
+        normalization = np.sum(features * features, axis=1, keepdims=True)
+        self.coefficients += (
+            rate * residual[:, None] * features / np.maximum(normalization, 1e-12)
+        )
+        return residual
+
+    def copy(self) -> "LocalPredictor":
+        return LocalPredictor(
+            reference_gate=self.reference_gate.copy(),
+            feature_scale=self.feature_scale.copy(),
+            coefficients=self.coefficients.copy(),
+            affine=self.affine,
+        )
+
+
+def free_output(circuit: Circuit, gates: Array, input_voltage: float) -> float:
+    """Appendix D, Eq. D13, with gates ordered (minus, plus)."""
+    gate_minus, gate_plus = np.asarray(gates, dtype=float)
+    scale = circuit.conductance_per_gate
+    numerator = (
+        input_voltage * circuit.fixed_conductance
+        + scale * (
+            gate_plus * circuit.high
+            + gate_minus * circuit.low
+            - (circuit.low + circuit.high) * circuit.threshold_voltage
+        )
+    )
+    denominator = (
+        circuit.fixed_conductance
+        + scale * (
+            gate_plus + gate_minus - 2.0 * circuit.threshold_voltage
+        )
+    )
+    if denominator <= 0.0:
+        raise ValueError("nonpositive effective conductance")
+    return float(numerator / denominator)
+
+
+def solution_line(circuit: Circuit, task: Task) -> tuple[float, float]:
+    """Return slope/intercept of gate_plus versus gate_minus at zero error."""
+    label = task.label_voltage
+    scale = circuit.conductance_per_gate
+    denominator = scale * (circuit.high - label)
+    if denominator == 0.0:
+        raise ValueError("label coincides with high boundary")
+    slope = -(circuit.low - label) / (circuit.high - label)
+    intercept = -(
+        circuit.fixed_conductance * (task.input_voltage - label)
+        + scale * circuit.threshold_voltage
+        * (2.0 * label - circuit.low - circuit.high)
+    ) / denominator
+    return float(slope), float(intercept)
+
+
+def joint_solution(circuit: Circuit, tasks: Iterable[Task]) -> Array:
+    tasks = tuple(tasks)
+    if len(tasks) != 2:
+        raise ValueError("joint_solution expects exactly two tasks")
+    slope_a, intercept_a = solution_line(circuit, tasks[0])
+    slope_b, intercept_b = solution_line(circuit, tasks[1])
+    if slope_a == slope_b:
+        raise ValueError("parallel solution lines have no unique joint solution")
+    gate_minus = (intercept_b - intercept_a) / (slope_a - slope_b)
+    return np.asarray(
+        (gate_minus, slope_a * gate_minus + intercept_a), dtype=float)
+
+
+def voltage_drop_squares(circuit: Circuit, output: float) -> Array:
+    return np.asarray(
+        ((output - circuit.low) ** 2, (circuit.high - output) ** 2),
+        dtype=float,
+    )
+
+
+def standard_clean_rate(
+    circuit: Circuit, gates: Array, task: Task, nudging: float = 1.0
+) -> tuple[Array, float, float]:
+    output_free = free_output(circuit, gates, task.input_voltage)
+    output_clamped = output_free + nudging * (
+        task.label_voltage - output_free)
+    rate = circuit.measured_learning_rate * (
+        voltage_drop_squares(circuit, output_free)
+        - voltage_drop_squares(circuit, output_clamped)
+    )
+    return rate, output_free, output_clamped
+
+
+def overclamped_clean_rate(
+    circuit: Circuit,
+    gates: Array,
+    task: Task,
+    *,
+    nudging: float = 0.25,
+    clamp_magnitude: Optional[float] = None,
+) -> tuple[Array, float, float]:
+    """Leading-order overclamping signal from Appendix F, Eq. F6--F8."""
+    output_free = free_output(circuit, gates, task.input_voltage)
+    error = task.label_voltage - output_free
+    magnitude = circuit.high if clamp_magnitude is None else clamp_magnitude
+    output_clamped = output_free + nudging * magnitude * np.sign(error)
+    rate = circuit.measured_learning_rate * (
+        voltage_drop_squares(circuit, output_free)
+        - voltage_drop_squares(circuit, output_clamped)
+    )
+    return rate, output_free, output_clamped
+
+
+def calibrate_predictor(
+    predictor: LocalPredictor,
+    states: Array,
+    measurement: Callable[[Array], Array],
+    *,
+    epochs: int,
+    learning_rate: float,
+) -> int:
+    """Sequential local LMS calibration; returns neutral observation count."""
+    states = np.asarray(states, dtype=float)
+    if states.ndim != 2 or states.shape[1] != 2:
+        raise ValueError("calibration states must have shape (observations, 2)")
+    if epochs < 1:
+        raise ValueError("epochs must be positive")
+    count = 0
+    for _ in range(epochs):
+        for gates in states:
+            predictor.update(gates, measurement(gates), learning_rate)
+            count += 1
+    return count
+
+
+def local_replay_update(
+    predictor: LocalPredictor,
+    gates: Array,
+    teaching_measurement: Array,
+    eligibility: Array,
+    learning_rate: float,
+) -> Array:
+    """The complete stored-tuple SDIL update, independent of any task/model."""
+    residual = np.asarray(teaching_measurement) - predictor.predict(gates)
+    return learning_rate * residual * np.asarray(eligibility)
+