results: test SDIL on measured physical surrogate

1 files changed, 155 insertions, 0 deletions

diff --git a/sdil/physical_coupled.py b/sdil/physical_coupled.py
index 0216048..f451947 100644
--- a/sdil/physical_coupled.py
+++ b/sdil/physical_coupled.py
@@ -244,3 +244,158 @@ def local_replay_update(
     residual = np.asarray(teaching_measurement) - predictor.predict(gates)
     return learning_rate * residual * np.asarray(eligibility)
 
+
+def task_errors(circuit: Circuit, gates: Array, tasks: Iterable[Task]) -> Array:
+    return np.asarray([
+        (task.label_voltage - free_output(
+            circuit, gates, task.input_voltage)) ** 2
+        for task in tasks
+    ], dtype=float)
+
+
+def simulate_alternating_tasks(
+    circuit: Circuit,
+    tasks: Iterable[Task],
+    bias_field: LocalAffineBias,
+    *,
+    method: str,
+    period_seconds: float,
+    cycles: int,
+    initial_gates: Array,
+    bias_strength: float = 1.0,
+    predictor: Optional[LocalPredictor] = None,
+    online_predictor_rate: float = 0.05,
+    noise_standard_deviation: Optional[Array] = None,
+    seed: int = 0,
+    summary_cycles: int = 20,
+    record_history: bool = False,
+) -> dict:
+    """Alternate two tasks using explicit local circuit updates.
+
+    `online_constant` and `online_sdil` take one neutral observation at the
+    beginning of each half-cycle.  Frozen predictors take none during task
+    learning.  The overclamping implementation uses the leading-order
+    constant-displacement signal of Eq. F6 and the error-proportional update
+    duration of Eq. F8; it is an analogue for these regression tasks, not a
+    reproduction of the paper's classification experiment.
+    """
+    allowed = {
+        "raw", "frozen_constant", "frozen_sdil", "online_constant",
+        "online_sdil", "oracle", "same_rms_noise", "overclamp",
+    }
+    if method not in allowed:
+        raise ValueError(f"unrecognized method {method}")
+    tasks = tuple(tasks)
+    if len(tasks) != 2:
+        raise ValueError("exactly two alternating tasks are required")
+    if period_seconds <= 0.0 or cycles < 1:
+        raise ValueError("period and cycles must be positive")
+    if method in {
+        "frozen_constant", "frozen_sdil", "online_constant", "online_sdil"
+    } and predictor is None:
+        raise ValueError(f"{method} requires a predictor")
+    if method == "same_rms_noise" and noise_standard_deviation is None:
+        raise ValueError("same_rms_noise requires a standard deviation")
+
+    active_predictor = predictor.copy() if predictor is not None else None
+    gates = np.asarray(initial_gates, dtype=float).copy()
+    if gates.shape != (2,):
+        raise ValueError("initial gates must have shape (2,)")
+    nominal_step = circuit.integration_step_seconds
+    half_steps = max(1, int(round(period_seconds / (2.0 * nominal_step))))
+    rng = np.random.default_rng(seed)
+    initial_error_scale = float(np.mean([
+        abs(task.label_voltage - free_output(
+            circuit, gates, task.input_voltage))
+        for task in tasks
+    ]))
+    initial_error_scale = max(initial_error_scale, 1e-6)
+
+    combined_error_history = []
+    cycle_span_history = []
+    gate_history = []
+    task_error_history = []
+    learning_on_time = 0.0
+    neutral_observations = 0
+    clipped_updates = 0
+
+    for _ in range(cycles):
+        half_endpoints = []
+        half_task_errors = []
+        for task in tasks:
+            if method in {"online_constant", "online_sdil"}:
+                neutral = bias_field(gates, bias_strength)
+                active_predictor.update(
+                    gates, neutral, online_predictor_rate)
+                neutral_observations += 1
+            for _ in range(half_steps):
+                physical_bias = bias_field(gates, bias_strength)
+                if method == "overclamp":
+                    clean_rate, output_free, _ = overclamped_clean_rate(
+                        circuit, gates, task)
+                    duration = nominal_step * abs(
+                        task.label_voltage - output_free) / initial_error_scale
+                    residual_bias = physical_bias
+                else:
+                    clean_rate, _, _ = standard_clean_rate(circuit, gates, task)
+                    duration = nominal_step
+                    if method == "raw":
+                        residual_bias = physical_bias
+                    elif method == "oracle":
+                        residual_bias = np.zeros(2, dtype=float)
+                    elif method == "same_rms_noise":
+                        residual_bias = rng.normal(
+                            loc=0.0,
+                            scale=np.asarray(noise_standard_deviation, dtype=float),
+                            size=2,
+                        )
+                    else:
+                        residual_bias = (
+                            physical_bias - active_predictor.predict(gates)
+                        )
+                proposed = gates + duration * (clean_rate + residual_bias)
+                clipped = np.clip(
+                    proposed, circuit.gate_minimum, circuit.gate_maximum)
+                clipped_updates += int(np.any(clipped != proposed))
+                gates = clipped
+                learning_on_time += duration
+            half_endpoints.append(gates.copy())
+            half_task_errors.append(task_errors(circuit, gates, tasks))
+        half_task_errors_array = np.asarray(half_task_errors)
+        combined_error_history.append(float(np.mean(half_task_errors_array)))
+        cycle_span_history.append(float(np.linalg.norm(
+            half_endpoints[1] - half_endpoints[0])))
+        gate_history.append(np.asarray(half_endpoints).tolist())
+        task_error_history.append(half_task_errors_array.tolist())
+
+    summary_count = min(summary_cycles, cycles)
+    combined = np.asarray(combined_error_history[-summary_count:])
+    spans = np.asarray(cycle_span_history[-summary_count:])
+    result = {
+        "method": method,
+        "period_seconds": period_seconds,
+        "cycles": cycles,
+        "half_steps": half_steps,
+        "initial_gates": np.asarray(initial_gates, dtype=float).tolist(),
+        "final_gates": gates.tolist(),
+        "bias_strength": bias_strength,
+        "mean_combined_error": float(np.mean(combined)),
+        "std_combined_error": float(np.std(combined)),
+        "mean_cycle_span": float(np.mean(spans)),
+        "std_cycle_span": float(np.std(spans)),
+        "neutral_observations_during_learning": neutral_observations,
+        "learning_on_time_seconds": float(learning_on_time),
+        "clipped_updates": clipped_updates,
+        "final_predictor_coefficients": (
+            None if active_predictor is None
+            else active_predictor.coefficients.tolist()
+        ),
+    }
+    if record_history:
+        result.update({
+            "combined_error_history": combined_error_history,
+            "cycle_span_history": cycle_span_history,
+            "half_cycle_gate_history": gate_history,
+            "half_cycle_task_error_history": task_error_history,
+        })
+    return result