feat: add sparse digital CLLN scaling ladder

1 files changed, 336 insertions, 0 deletions

diff --git a/sdil/coupled_ladder.py b/sdil/coupled_ladder.py
new file mode 100644
index 0000000..a9b8014
--- /dev/null
+++ b/sdil/coupled_ladder.py
@@ -0,0 +1,336 @@
+"""Sparse digital coupled-learning grids for controlled scaling experiments.
+
+The network is a linear resistor lattice.  Its local learning signal is the
+same free-minus-clamped voltage-square difference used by coupled learning.
+Component imperfections reuse the per-edge measurement model used by the
+hardware-realistic simulator in :mod:`sdil.physical_grid`.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+from scipy.sparse import coo_matrix
+from scipy.sparse.linalg import spsolve
+
+from sdil.physical_grid import (
+    GridCircuit,
+    GridSquareLawImperfection,
+    RingClassificationDataset,
+    edge_voltage_drops,
+    output_difference,
+)
+
+
+Array = np.ndarray
+
+
+def make_scaled_grid(side: int) -> GridCircuit:
+    """Return a square grid whose side-4 boundary layout matches Figure 5."""
+    if side < 4 or side % 4:
+        raise ValueError("grid side must be a positive multiple of four")
+
+    def node(row: int, column: int) -> int:
+        return row * side + column
+
+    return GridCircuit(
+        rows=side,
+        columns=side,
+        source_nodes=(
+            node(3 * side // 4, 3 * side // 4),
+            node(3 * side // 4, side // 4),
+            node(side // 4, 3 * side // 4),
+            node(side // 4, side // 4),
+        ),
+        target_nodes=(
+            node(side // 2, side // 2),
+            node(side // 2, 0),
+        ),
+    )
+
+
+def tile_figure5_gates(base_gates: Array, side: int) -> Array:
+    """Tile a released 4-by-4 horizontal/vertical gate pattern."""
+    gates = np.asarray(base_gates, dtype=float)
+    if gates.shape != (32,):
+        raise ValueError("the released base gate vector must have 32 entries")
+    if side < 4 or side % 4:
+        raise ValueError("grid side must be a positive multiple of four")
+    horizontal = gates[:16].reshape(4, 4)
+    vertical = gates[16:].reshape(4, 4)
+    tiled_horizontal = np.asarray([
+        horizontal[row % 4, column % 4]
+        for row in range(side)
+        for column in range(side)
+    ])
+    tiled_vertical = np.asarray([
+        vertical[row % 4, column % 4]
+        for row in range(side)
+        for column in range(side)
+    ])
+    return np.concatenate((tiled_horizontal, tiled_vertical))
+
+
+def solve_linear_grid_state(
+    circuit: GridCircuit,
+    gates: Array,
+    source_values: Array,
+    *,
+    target_values: Array | None = None,
+) -> Array:
+    """Solve the sparse linear Kirchhoff system with fixed boundary nodes."""
+    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")
+    conductances = circuit.conductance_scale * (
+        gates - circuit.threshold_voltage)
+    if np.any(conductances <= 0.0):
+        raise ValueError("all digital conductances must be positive")
+
+    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_nodes = [
+        node for node in range(circuit.node_count) if node not in fixed
+    ]
+    unknown_index = {node: index for index, node in enumerate(unknown_nodes)}
+    matrix_rows: list[int] = []
+    matrix_columns: list[int] = []
+    matrix_values: list[float] = []
+    rhs = np.zeros(len(unknown_nodes), dtype=float)
+
+    for conductance, (first, second) in zip(
+        conductances, circuit.edge_pairs
+    ):
+        for node, neighbour in ((first, second), (second, first)):
+            if node in fixed:
+                continue
+            row = unknown_index[node]
+            matrix_rows.append(row)
+            matrix_columns.append(row)
+            matrix_values.append(float(conductance))
+            if neighbour in fixed:
+                rhs[row] += conductance * fixed[neighbour]
+            else:
+                matrix_rows.append(row)
+                matrix_columns.append(unknown_index[neighbour])
+                matrix_values.append(float(-conductance))
+
+    matrix = coo_matrix(
+        (matrix_values, (matrix_rows, matrix_columns)),
+        shape=(len(unknown_nodes), len(unknown_nodes)),
+    ).tocsr()
+    unknown_voltages = np.asarray(spsolve(matrix, rhs), dtype=float)
+    if not np.all(np.isfinite(unknown_voltages)):
+        raise RuntimeError("linear grid solve returned a nonfinite state")
+
+    voltages = np.empty(circuit.node_count, dtype=float)
+    for node, voltage in fixed.items():
+        voltages[node] = voltage
+    voltages[unknown_nodes] = unknown_voltages
+    return voltages
+
+
+def evaluate_digital_grid(
+    circuit: GridCircuit,
+    gates: Array,
+    dataset: RingClassificationDataset,
+) -> dict:
+    outputs = []
+    for inputs in dataset.inputs_v:
+        state = solve_linear_grid_state(
+            circuit, gates, circuit.source_values(*inputs))
+        outputs.append(output_difference(circuit, state))
+    outputs_array = np.asarray(outputs)
+    labels = np.asarray(dataset.labels_v)
+    errors = labels - outputs_array
+    active = labels * errors > 0.0
+    return {
+        "classification_error": float(np.mean(
+            np.sign(outputs_array) != np.sign(labels))),
+        "hinge_loss_v2": float(np.mean(np.where(
+            active, np.square(errors), 0.0))),
+        "margin_success_fraction": float(np.mean(~active)),
+        "outputs_v": outputs_array.tolist(),
+    }
+
+
+@dataclass(frozen=True)
+class DigitalTrainingConfig:
+    epochs: int = 600
+    record_every: int = 10
+    standard_nudging: float = 128.0 / 129.0
+    learning_time_seconds: float = 1.0e-3
+    overclamp_nudging: float = 32.0 / 129.0
+    overclamp_time_seconds_per_v: float = 2.5e-3
+    overclamp_target_magnitude_v: float | None = None
+
+    def __post_init__(self) -> None:
+        if self.epochs < 1 or self.record_every < 1:
+            raise ValueError("epochs and record interval must be positive")
+        if self.learning_time_seconds <= 0.0:
+            raise ValueError("learning time must be positive")
+
+
+def train_digital_grid(
+    circuit: GridCircuit,
+    initial_gates: Array,
+    dataset: RingClassificationDataset,
+    imperfection: GridSquareLawImperfection,
+    *,
+    method: str,
+    config: DigitalTrainingConfig,
+    constant_bias_v_per_s: Array | None = None,
+    noise_seed: int = 0,
+) -> dict:
+    """Train with a hand-written local coupled-learning update."""
+    allowed = {
+        "clean",
+        "matched_noise",
+        "raw",
+        "constant",
+        "sdil",
+        "overclamp",
+        "overclamp_sdil",
+    }
+    if method not in allowed:
+        raise ValueError(f"unrecognized method {method}")
+    if method == "constant" and constant_bias_v_per_s is None:
+        raise ValueError("constant calibration requires a per-edge baseline")
+
+    gates = np.asarray(initial_gates, dtype=float).copy()
+    if gates.shape != (circuit.edge_count,):
+        raise ValueError("initial gate vector has the wrong shape")
+    if constant_bias_v_per_s is not None:
+        constant_bias = np.asarray(constant_bias_v_per_s, dtype=float)
+        if constant_bias.shape != gates.shape:
+            raise ValueError("constant baseline has the wrong shape")
+    else:
+        constant_bias = None
+
+    rng = np.random.default_rng(noise_seed)
+    local_updates = 0
+    neutral_observations = 0
+    clipped_updates = 0
+    cumulative_learning_time = 0.0
+    trace = [{
+        "epoch": 0,
+        "local_updates": 0,
+        **evaluate_digital_grid(circuit, gates, dataset),
+    }]
+
+    for epoch in range(1, config.epochs + 1):
+        for inputs, label in zip(dataset.inputs_v, dataset.labels_v):
+            sources = circuit.source_values(*inputs)
+            free_state = solve_linear_grid_state(circuit, gates, sources)
+            output_free = output_difference(circuit, free_state)
+            error = label - output_free
+            if label * error <= 0.0:
+                continue
+
+            free_drops = edge_voltage_drops(circuit, free_state)
+            is_overclamp = method.startswith("overclamp")
+            if is_overclamp:
+                target_magnitude = (
+                    circuit.high_voltage
+                    if config.overclamp_target_magnitude_v is None
+                    else config.overclamp_target_magnitude_v
+                )
+                output_clamped = output_free + config.overclamp_nudging * (
+                    target_magnitude * np.sign(error) - output_free)
+                duration = (
+                    config.overclamp_time_seconds_per_v * abs(error))
+            else:
+                output_clamped = output_free + (
+                    config.standard_nudging * error)
+                duration = config.learning_time_seconds
+
+            target_mean = float(np.mean(
+                free_state[list(circuit.target_nodes)]))
+            target_values = np.asarray((
+                target_mean + 0.5 * output_clamped,
+                target_mean - 0.5 * output_clamped,
+            ))
+            clamped_state = solve_linear_grid_state(
+                circuit, gates, sources, target_values=target_values)
+            clamped_drops = edge_voltage_drops(circuit, clamped_state)
+            ideal_rate = imperfection.ideal_rate(
+                circuit.measured_learning_rate, free_drops, clamped_drops)
+            observed_rate = imperfection.observed_rate(
+                circuit.measured_learning_rate, free_drops, clamped_drops)
+
+            if method == "clean":
+                applied_rate = ideal_rate
+            elif method == "matched_noise":
+                measurement_error = observed_rate - ideal_rate
+                random_sign = rng.choice((-1.0, 1.0), size=len(gates))
+                applied_rate = ideal_rate + random_sign * np.abs(
+                    measurement_error)
+            elif method == "constant":
+                applied_rate = observed_rate - constant_bias
+            elif method in {"sdil", "overclamp_sdil"}:
+                neutral_rate = imperfection.neutral_bias(
+                    circuit.measured_learning_rate, free_drops)
+                applied_rate = observed_rate - neutral_rate
+                neutral_observations += 1
+            else:
+                applied_rate = observed_rate
+
+            proposed = gates + duration * applied_rate
+            clipped = np.clip(
+                proposed, circuit.gate_minimum, circuit.gate_maximum)
+            clipped_updates += int(np.any(clipped != proposed))
+            gates = clipped
+            local_updates += 1
+            cumulative_learning_time += duration
+
+        if epoch % config.record_every == 0 or epoch == config.epochs:
+            trace.append({
+                "epoch": epoch,
+                "local_updates": local_updates,
+                **evaluate_digital_grid(circuit, gates, dataset),
+            })
+
+    zero_records = [
+        record for record in trace
+        if record["classification_error"] == 0.0
+    ]
+    if zero_records:
+        epochs_to_zero = int(zero_records[0]["epoch"])
+        updates_to_zero = int(zero_records[0]["local_updates"])
+        reached_zero = True
+    else:
+        epochs_to_zero = config.epochs
+        updates_to_zero = local_updates
+        reached_zero = False
+    epoch_axis = np.asarray([record["epoch"] for record in trace])
+    error_axis = np.asarray([
+        record["classification_error"] for record in trace])
+    error_auc = float(np.trapezoid(error_axis, epoch_axis) / config.epochs)
+    final = trace[-1]
+    return {
+        "method": method,
+        "classification_error": final["classification_error"],
+        "hinge_loss_v2": final["hinge_loss_v2"],
+        "margin_success_fraction": final["margin_success_fraction"],
+        "outputs_v": final["outputs_v"],
+        "reached_zero_error": reached_zero,
+        "restricted_epochs_to_zero_error": epochs_to_zero,
+        "restricted_updates_to_zero_error": updates_to_zero,
+        "classification_error_auc": error_auc,
+        "local_updates": local_updates,
+        "neutral_observations": neutral_observations,
+        "cumulative_learning_time_seconds": float(cumulative_learning_time),
+        "clipped_updates": clipped_updates,
+        "final_gates_v": gates.tolist(),
+        "trace": trace,
+    }
+