diff options
| -rw-r--r-- | experiments/physical_grid_smoke.py | 121 | ||||
| -rw-r--r-- | sdil/physical_grid.py | 314 |
2 files changed, 435 insertions, 0 deletions
diff --git a/experiments/physical_grid_smoke.py b/experiments/physical_grid_smoke.py new file mode 100644 index 0000000..a84ed84 --- /dev/null +++ b/experiments/physical_grid_smoke.py @@ -0,0 +1,121 @@ +#!/usr/bin/env python3 +"""Mechanics checks for the reconstructed Figure-5 physical grid.""" + +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_grid import ( # noqa: E402 + EdgePolynomialPredictor, + GridCircuit, + GridSquareLawImperfection, + _residual_and_jacobian, + edge_voltage_drops, + fit_edge_predictor, + output_difference, + solve_grid_state, +) + + +def sample_local_states( + circuit: GridCircuit, count: int, seed: int +) -> np.ndarray: + rng = np.random.default_rng(seed) + states = [] + for _ in range(count): + gates = rng.normal(2.33, 0.08, circuit.edge_count) + sources = circuit.source_values( + rng.uniform(circuit.low_voltage, circuit.high_voltage), + rng.uniform(circuit.low_voltage, circuit.high_voltage), + ) + voltages = solve_grid_state(circuit, gates, sources) + states.append(edge_voltage_drops(circuit, voltages)) + return np.asarray(states) + + +def main() -> None: + circuit = GridCircuit() + assert circuit.edge_count == 32 + assert circuit.source_nodes == (15, 13, 7, 5) + assert circuit.target_nodes == (10, 8) + gates = np.full(circuit.edge_count, 2.33) + sources = circuit.source_values(0.2266, 0.2266) + free = solve_grid_state(circuit, gates, sources) + assert np.allclose(free[list(circuit.source_nodes)], sources) + assert abs(output_difference(circuit, free)) < 1e-12 + residual, _ = _residual_and_jacobian(circuit, gates, free) + unknown = [ + node for node in range(circuit.node_count) + if node not in circuit.source_nodes + ] + assert np.linalg.norm(residual[unknown], ord=np.inf) < 1e-10 + + target_difference = 0.018 + target_mean = float(np.mean(free[list(circuit.target_nodes)])) + target_values = np.asarray(( + target_mean + 0.5 * target_difference, + target_mean - 0.5 * target_difference, + )) + clamped = solve_grid_state( + circuit, gates, sources, target_values=target_values, + initial_state=free) + assert np.allclose(clamped[list(circuit.target_nodes)], target_values) + assert abs(output_difference(circuit, clamped) - target_difference) < 1e-12 + + ideal = GridSquareLawImperfection.ideal(circuit.edge_count) + free_drops = edge_voltage_drops(circuit, free) + clamped_drops = edge_voltage_drops(circuit, clamped) + assert np.array_equal( + ideal.observed_rate( + circuit.measured_learning_rate, free_drops, clamped_drops), + ideal.ideal_rate( + circuit.measured_learning_rate, free_drops, clamped_drops), + ) + + hardware = GridSquareLawImperfection.sample_appendix_c( + circuit.edge_count, 20260829) + calibration_states = sample_local_states(circuit, 48, 20260830) + neutral = np.asarray([ + hardware.neutral_bias(circuit.measured_learning_rate, state) + for state in calibration_states + ]) + center = np.mean(calibration_states, axis=0) + scale = np.maximum(np.ptp(calibration_states, axis=0), 1e-3) + constant = EdgePolynomialPredictor.zeros(center, scale, degree=0) + quadratic = EdgePolynomialPredictor.zeros(center, scale, degree=2) + assert fit_edge_predictor(constant, calibration_states, neutral) == 48 + assert fit_edge_predictor(quadratic, calibration_states, neutral) == 48 + heldout_states = sample_local_states(circuit, 16, 20260831) + heldout_neutral = np.asarray([ + hardware.neutral_bias(circuit.measured_learning_rate, state) + for state in heldout_states + ]) + constant_rmse = float(np.sqrt(np.mean([ + np.square(measurement - constant.predict(state)) + for state, measurement in zip(heldout_states, heldout_neutral) + ]))) + quadratic_rmse = float(np.sqrt(np.mean([ + np.square(measurement - quadratic.predict(state)) + for state, measurement in zip(heldout_states, heldout_neutral) + ]))) + assert quadratic_rmse < 1e-7 + assert constant_rmse > 1e-3 + print({ + "nodes": circuit.node_count, + "edges": circuit.edge_count, + "source_nodes": circuit.source_nodes, + "target_nodes": circuit.target_nodes, + "constant_neutral_rmse_v_per_s": constant_rmse, + "quadratic_neutral_rmse_v_per_s": quadratic_rmse, + "autodiff_used": False, + }) + + +if __name__ == "__main__": + main() diff --git a/sdil/physical_grid.py b/sdil/physical_grid.py new file mode 100644 index 0000000..10d667b --- /dev/null +++ b/sdil/physical_grid.py @@ -0,0 +1,314 @@ +"""Circuit-faithful 4x4 physical learning network used by Dillavou et al. + +The nonlinear conductance, periodic topology, and local voltage-square update +follow Eqs. (2)--(3) of arXiv:2505.22887v2. Source and target node locations +match the released Figure-5 experiment objects. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np + + +Array = np.ndarray + + +@dataclass(frozen=True) +class GridCircuit: + rows: int = 4 + columns: int = 4 + conductance_scale: float = 8.0e-4 + threshold_voltage: float = 0.7 + measured_learning_rate: float = 2.5e3 + low_voltage: float = 0.0181 + high_voltage: float = 0.4351 + gate_minimum: float = 1.0 + gate_maximum: float = 5.2 + source_nodes: tuple[int, ...] = (15, 13, 7, 5) + target_nodes: tuple[int, int] = (10, 8) + + @property + def node_count(self) -> int: + return self.rows * self.columns + + @property + def edge_pairs(self) -> tuple[tuple[int, int], ...]: + horizontal = [] + vertical = [] + for row in range(self.rows): + for column in range(self.columns): + node = row * self.columns + column + horizontal.append(( + node, + row * self.columns + (column + 1) % self.columns, + )) + vertical.append(( + node, + ((row + 1) % self.rows) * self.columns + column, + )) + return tuple(horizontal + vertical) + + @property + def edge_count(self) -> int: + return len(self.edge_pairs) + + def source_values(self, input_one: float, input_two: float) -> Array: + return np.asarray(( + input_one, + input_two, + self.low_voltage, + self.high_voltage, + ), dtype=float) + + +def edge_voltage_drops(circuit: GridCircuit, node_voltages: Array) -> Array: + voltages = np.asarray(node_voltages, dtype=float) + if voltages.shape != (circuit.node_count,): + raise ValueError("node voltage vector has the wrong shape") + return np.asarray([ + voltages[first] - voltages[second] + for first, second in circuit.edge_pairs + ]) + + +def output_difference(circuit: GridCircuit, node_voltages: Array) -> float: + positive, negative = circuit.target_nodes + return float(node_voltages[positive] - node_voltages[negative]) + + +def _residual_and_jacobian( + circuit: GridCircuit, gates: Array, voltages: Array +) -> tuple[Array, Array]: + residual = np.zeros(circuit.node_count, dtype=float) + jacobian = np.zeros( + (circuit.node_count, circuit.node_count), dtype=float) + scale = circuit.conductance_scale + threshold = circuit.threshold_voltage + for gate, (first, second) in zip(gates, circuit.edge_pairs): + voltage_first = voltages[first] + voltage_second = voltages[second] + conductance = scale * ( + gate - threshold - 0.5 * (voltage_first + voltage_second)) + current = conductance * (voltage_first - voltage_second) + residual[first] += current + residual[second] -= current + derivative_first = scale * (gate - threshold - voltage_first) + derivative_second = scale * (-gate + threshold + voltage_second) + jacobian[first, first] += derivative_first + jacobian[first, second] += derivative_second + jacobian[second, first] -= derivative_first + jacobian[second, second] -= derivative_second + return residual, jacobian + + +def solve_grid_state( + circuit: GridCircuit, + gates: Array, + source_values: Array, + *, + target_values: Array | None = None, + initial_state: Array | None = None, + tolerance: float = 1e-11, + maximum_iterations: int = 20, +) -> Array: + """Solve Kirchhoff's laws by Newton iteration with an analytic Jacobian.""" + gates = np.asarray(gates, dtype=float) + sources = np.asarray(source_values, dtype=float) + if gates.shape != (circuit.edge_count,): + raise ValueError("gate vector has the wrong shape") + if sources.shape != (len(circuit.source_nodes),): + raise ValueError("source voltage vector has the wrong shape") + fixed = dict(zip(circuit.source_nodes, sources)) + if target_values is not None: + targets = np.asarray(target_values, dtype=float) + if targets.shape != (2,): + raise ValueError("target voltage vector must have shape (2,)") + fixed.update(zip(circuit.target_nodes, targets)) + unknown = np.asarray([ + node for node in range(circuit.node_count) if node not in fixed + ]) + voltages = np.full( + circuit.node_count, float(np.mean(sources)), dtype=float) + if initial_state is not None: + initial = np.asarray(initial_state, dtype=float) + if initial.shape != (circuit.node_count,): + raise ValueError("initial state vector has the wrong shape") + voltages[:] = initial + for node, value in fixed.items(): + voltages[node] = value + + for _ in range(maximum_iterations): + residual, jacobian = _residual_and_jacobian( + circuit, gates, voltages) + unknown_residual = residual[unknown] + if np.linalg.norm(unknown_residual, ord=np.inf) <= tolerance: + return voltages + unknown_jacobian = jacobian[np.ix_(unknown, unknown)] + step = np.linalg.solve(unknown_jacobian, unknown_residual) + voltages[unknown] -= step + residual, _ = _residual_and_jacobian(circuit, gates, voltages) + raise RuntimeError( + "grid state did not converge; residual=" + f"{np.linalg.norm(residual[unknown], ord=np.inf):.3e}") + + +@dataclass(frozen=True) +class GridSquareLawImperfection: + free_gain: Array + clamped_gain: Array + free_input_offset_v: Array + clamped_input_offset_v: Array + multiplier_output_offset_v_per_s: Array + + def __post_init__(self) -> None: + shapes = { + np.asarray(value).shape + for value in ( + self.free_gain, + self.clamped_gain, + self.free_input_offset_v, + self.clamped_input_offset_v, + self.multiplier_output_offset_v_per_s, + ) + } + if len(shapes) != 1: + raise ValueError("grid imperfection arrays disagree") + shape = next(iter(shapes)) + if len(shape) != 1 or shape[0] < 1: + raise ValueError("grid imperfection arrays must be nonempty vectors") + + @classmethod + def ideal(cls, edge_count: int) -> "GridSquareLawImperfection": + return cls( + free_gain=np.ones(edge_count), + clamped_gain=np.ones(edge_count), + free_input_offset_v=np.zeros(edge_count), + clamped_input_offset_v=np.zeros(edge_count), + multiplier_output_offset_v_per_s=np.zeros(edge_count), + ) + + @classmethod + def sample_appendix_c( + cls, + edge_count: int, + seed: int, + *, + gain_standard_deviation: float = 0.01, + twin_mismatch_standard_deviation_v: float = 0.001, + multiplier_offset_standard_deviation_v_per_s: float = 2.3, + ) -> "GridSquareLawImperfection": + rng = np.random.default_rng(seed) + common_gain = 1.0 + rng.normal( + 0.0, gain_standard_deviation, edge_count) + differential_gain = rng.normal( + 0.0, gain_standard_deviation, edge_count) + common_offset = rng.normal( + 0.0, twin_mismatch_standard_deviation_v, edge_count) + differential_offset = rng.normal( + 0.0, twin_mismatch_standard_deviation_v, edge_count) + return cls( + free_gain=common_gain + 0.5 * differential_gain, + clamped_gain=common_gain - 0.5 * differential_gain, + free_input_offset_v=common_offset + 0.5 * differential_offset, + clamped_input_offset_v=common_offset - 0.5 * differential_offset, + multiplier_output_offset_v_per_s=rng.normal( + 0.0, multiplier_offset_standard_deviation_v_per_s, + edge_count), + ) + + def observed_rate( + self, + learning_rate: float, + free_drops: Array, + clamped_drops: Array, + ) -> Array: + measured_free = ( + self.free_gain * free_drops + self.free_input_offset_v) + measured_clamped = ( + self.clamped_gain * clamped_drops + self.clamped_input_offset_v) + return ( + learning_rate + * (np.square(measured_free) - np.square(measured_clamped)) + + self.multiplier_output_offset_v_per_s + ) + + @staticmethod + def ideal_rate( + learning_rate: float, free_drops: Array, clamped_drops: Array + ) -> Array: + return learning_rate * ( + np.square(free_drops) - np.square(clamped_drops)) + + def neutral_bias(self, learning_rate: float, free_drops: Array) -> Array: + return self.observed_rate(learning_rate, free_drops, free_drops) + + +@dataclass +class EdgePolynomialPredictor: + feature_center: Array + feature_scale: Array + coefficients: Array + + @classmethod + def zeros( + cls, feature_center: Array, feature_scale: Array, *, degree: int + ) -> "EdgePolynomialPredictor": + center = np.asarray(feature_center, dtype=float) + scale = np.asarray(feature_scale, dtype=float) + if center.ndim != 1 or scale.shape != center.shape: + raise ValueError("edge feature metadata disagree") + if degree < 0 or np.any(scale <= 0.0): + raise ValueError("invalid polynomial degree or feature scale") + return cls( + feature_center=center.copy(), + feature_scale=scale.copy(), + coefficients=np.zeros((len(center), degree + 1), dtype=float), + ) + + @property + def degree(self) -> int: + return int(self.coefficients.shape[1] - 1) + + def features(self, local_state: Array) -> Array: + state = np.asarray(local_state, dtype=float) + if state.shape != self.feature_center.shape: + raise ValueError("edge local state has the wrong shape") + normalized = (state - self.feature_center) / self.feature_scale + return np.stack([ + normalized ** power for power in range(self.degree + 1) + ], axis=1) + + def predict(self, local_state: Array) -> Array: + return np.sum(self.coefficients * self.features(local_state), axis=1) + + def copy(self) -> "EdgePolynomialPredictor": + return EdgePolynomialPredictor( + feature_center=self.feature_center.copy(), + feature_scale=self.feature_scale.copy(), + coefficients=self.coefficients.copy(), + ) + + +def fit_edge_predictor( + predictor: EdgePolynomialPredictor, + local_states: Array, + neutral_measurements: Array, + *, + ridge: float = 1e-12, +) -> int: + states = np.asarray(local_states, dtype=float) + measurements = np.asarray(neutral_measurements, dtype=float) + if states.ndim != 2 or measurements.shape != states.shape: + raise ValueError("edge calibration matrices disagree") + if states.shape[1] != len(predictor.feature_center): + raise ValueError("edge calibration width changed") + features = np.asarray([predictor.features(state) for state in states]) + for edge in range(states.shape[1]): + design = features[:, edge, :] + gram = design.T @ design + rhs = design.T @ measurements[:, edge] + predictor.coefficients[edge] = np.linalg.solve( + gram + ridge * np.eye(gram.shape[0]), rhs) + return int(len(states)) |