diff options
| -rw-r--r-- | experiments/physical_bias_p1_smoke.py | 87 | ||||
| -rw-r--r-- | sdil/physical_coupled.py | 246 |
2 files changed, 333 insertions, 0 deletions
diff --git a/experiments/physical_bias_p1_smoke.py b/experiments/physical_bias_p1_smoke.py new file mode 100644 index 0000000..dc9b6f1 --- /dev/null +++ b/experiments/physical_bias_p1_smoke.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python3 +"""Deterministic contract checks for the physical coupled-learning adapter.""" + +from __future__ import annotations + +from pathlib import Path +import sys + +import numpy as np + +ROOT = Path(__file__).resolve().parents[1] +sys.path.insert(0, str(ROOT)) + +from sdil.physical_coupled import ( + Circuit, + LocalAffineBias, + LocalPredictor, + Task, + calibrate_predictor, + free_output, + joint_solution, + local_replay_update, + standard_clean_rate, +) + + +def main() -> None: + circuit = Circuit() + tasks = ( + Task("alpha", circuit.high, 0.31), + Task("beta", circuit.low, 0.14), + ) + joint = joint_solution(circuit, tasks) + for task in tasks: + assert abs(free_output(circuit, joint, task.input_voltage) - task.label_voltage) < 1e-12 + rate, _, _ = standard_clean_rate(circuit, joint, task) + assert np.linalg.norm(rate) < 1e-10 + + field = LocalAffineBias( + reference_gate=np.asarray((3.0, 3.5)), + bias_at_reference=np.asarray((2.9, 4.7)), + local_slopes=np.asarray((0.6, -0.2)), + ) + states = np.column_stack(( + np.linspace(2.2, 4.8, 40), + np.linspace(2.8, 5.0, 40), + )) + scale = np.ptp(states, axis=0) + affine = LocalPredictor.zeros(field.reference_gate, scale, affine=True) + constant = LocalPredictor.zeros(field.reference_gate, scale, affine=False) + observations_affine = calibrate_predictor( + affine, states, field, epochs=30, learning_rate=0.2) + observations_constant = calibrate_predictor( + constant, states, field, epochs=30, learning_rate=0.2) + assert observations_affine == observations_constant + affine_error = np.mean([ + np.linalg.norm(field(state) - affine.predict(state)) for state in states + ]) + constant_error = np.mean([ + np.linalg.norm(field(state) - constant.predict(state)) for state in states + ]) + assert affine_error < 1e-4 + assert constant_error > 0.1 + + gates = states[7] + teaching = np.asarray((3.2, -1.7)) + eligibility = np.asarray((0.4, 0.8)) + update_a = local_replay_update(affine, gates, teaching, eligibility, 0.03) + update_b = 0.03 * (teaching - affine.predict(gates)) * eligibility + assert np.array_equal(update_a, update_b) + + downstream = np.random.default_rng(19).normal(size=(100, 100)) + downstream[:] = np.random.default_rng(29).normal(size=downstream.shape) + update_c = local_replay_update(affine, gates, teaching, eligibility, 0.03) + assert np.array_equal(update_a, update_c) + print({ + "joint_solution": joint.tolist(), + "affine_calibration_mae": float(affine_error), + "constant_calibration_mae": float(constant_error), + "neutral_observations_each": observations_affine, + "local_replay_exact": True, + "autodiff_used": False, + }) + + +if __name__ == "__main__": + main() 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) + |