feat: transfer released physical drift profile to Rain EP

5 files changed, 306 insertions, 8 deletions

diff --git a/RAIN_EP_RELEASED_PROFILE.md b/RAIN_EP_RELEASED_PROFILE.md
new file mode 100644
index 0000000..d124cc3
--- /dev/null
+++ b/RAIN_EP_RELEASED_PROFILE.md
@@ -0,0 +1,57 @@
+# Rain EP Released-Drift Profile
+
+## Scope
+
+This protocol asks whether a local affine predictor adds value beyond a local
+intercept when hardware update bias changes with the adaptive parameter state.
+It is a controlled neural-network transfer of a shape measured in released
+physical traces. It is not a reproduction of the resistor hardware and is not
+labeled as a real-hardware neural-network result.
+
+The source is the committed analysis
+`results/physical_bias/p0_state_dependence.json`, derived from Zenodo record
+15692914, release v1.0.1. For each of four measured edges, the report supplies
+an affine offset, local gate-voltage slope, and the gate range visited in the
+retained drift traces. The transfer profile divides both the offset and the
+slope times observed gate range by the RMS of the four offsets. This freezes
+the following dimensionless values without a task-accuracy fit:
+
+- offsets: `0.93562, 1.60206, -0.10014, 0.74026`;
+- full-range state variations: `0.52595, 0.31317, 0.38320, -0.04782`.
+
+Neural parameters are deterministically assigned these four profiles. Their
+initial local update-offset RMS is set by `bias_ratio`; parameter displacement
+in units of that tensor's initial parameter RMS maps through `tanh` to the
+measured state-range coordinate. This mapping preserves the measured relative
+state dependence but does not claim that resistor gate volts equal neural
+weight units.
+
+## Local observation contract
+
+At a neutral probe, the teaching input to the local update circuit is disabled.
+The corrector receives the resulting measured circuit output at its current
+parameter state. It is not given an externally calibrated coefficient. An
+intercept-only predictor and affine SDIL receive identical probes. All
+parameters probe in parallel, so probe count is independent of parameter count;
+each parameter stores its own predictor coefficients.
+
+## R0 development screen
+
+R0 begins only after `dillavou_c0` finishes. It uses the same 10,000/2,000
+training-only split as development, batch size 128, one epoch, fixed bias ratio
+1.0, and the released affine profile. Eight cells run:
+
+1. clean positive EP;
+2. raw released-profile bias;
+3. intercept-only, one initial probe;
+4. affine SDIL, one initial probe;
+5. intercept-only, one initial probe then every 10 steps;
+6. affine SDIL with the same probes;
+7. intercept-only, one initial probe then every 50 steps;
+8. affine SDIL with the same probes.
+
+R0 is a mechanics/development screen. A useful result requires affine SDIL to
+have lower held-out residual bias and higher task accuracy than the matched
+intercept-only arm. Final evidence additionally requires multiple seeds, a
+frozen probe cadence, measurement noise/quantization, and the strong-clamp
+baseline.
diff --git a/experiments/rain_ep_bias_train.py b/experiments/rain_ep_bias_train.py
index fa62137..2ba3dd0 100644
--- a/experiments/rain_ep_bias_train.py
+++ b/experiments/rain_ep_bias_train.py
@@ -19,6 +19,7 @@ ROOT = Path(__file__).resolve().parents[1]
 sys.path.insert(0, str(ROOT))
 
 from sdil.rain_ep_adapter import (  # noqa: E402
+    DillavouBiasProfile,
     DillavouUpdateCorrector,
     RainGradientCorrector,
     RainLayerStateCorrector,
@@ -55,6 +56,9 @@ def parse_args() -> argparse.Namespace:
         "--dillavou-drift-ratio", type=float, default=0.0,
         help="zero is the exact fixed update-offset model from Dillavou et al.")
     parser.add_argument("--dillavou-calibration-steps", type=int, default=1)
+    parser.add_argument(
+        "--dillavou-profile-json", type=Path,
+        help="released physical state-dependence report; omitted means constant B")
     parser.add_argument("--predictor-rate", type=float, default=0.1)
     parser.add_argument(
         "--neutral-cadence", type=int, default=1,
@@ -230,12 +234,18 @@ def main() -> None:
             raise ValueError(
                 "Dillavou calibration probes the local update circuit and "
                 "does not require equilibrium batches")
+        empirical_profile = None
+        if args.dillavou_profile_json is not None:
+            profile_path = args.dillavou_profile_json.resolve()
+            empirical_profile = DillavouBiasProfile.from_state_dependence_report(
+                json.loads(profile_path.read_text()), source=str(profile_path))
         corrector = DillavouUpdateCorrector(
             mode=args.mode,
             bias_ratio=args.bias_ratio,
             predictor_rate=args.predictor_rate,
             neutral_cadence=args.neutral_cadence,
             drift_ratio=args.dillavou_drift_ratio,
+            empirical_profile=empirical_profile,
             seed=args.seed + 1729,
         )
         attach_dillavou_to_rain_estimator(estimator, corrector)
@@ -354,6 +364,10 @@ def main() -> None:
             "bias_ratio": args.bias_ratio,
             "dillavou_drift_ratio": args.dillavou_drift_ratio,
             "dillavou_calibration_steps": args.dillavou_calibration_steps,
+            "dillavou_profile": (
+                None if args.adapter != "dillavou"
+                or corrector.empirical_profile is None
+                else corrector.empirical_profile.as_dict()),
             "predictor_rate": args.predictor_rate,
             "neutral_cadence": args.neutral_cadence,
             "layer_calibration_steps": args.layer_calibration_steps,
diff --git a/experiments/rain_ep_dillavou_smoke.py b/experiments/rain_ep_dillavou_smoke.py
index 9519210..39bf5fd 100644
--- a/experiments/rain_ep_dillavou_smoke.py
+++ b/experiments/rain_ep_dillavou_smoke.py
@@ -13,6 +13,7 @@ ROOT = Path(__file__).resolve().parents[1]
 sys.path.insert(0, str(ROOT))
 
 from sdil.rain_ep_adapter import (  # noqa: E402
+    DillavouBiasProfile,
     DillavouUpdateCorrector,
     attach_dillavou_to_rain_estimator,
 )
@@ -21,6 +22,9 @@ from sdil.rain_ep_adapter import (  # noqa: E402
 def parse_args() -> argparse.Namespace:
     parser = argparse.ArgumentParser()
     parser.add_argument("--author-root", type=Path, required=True)
+    parser.add_argument(
+        "--profile-json", type=Path,
+        default=ROOT / "results/physical_bias/p0_state_dependence.json")
     return parser.parse_args()
 
 
@@ -92,6 +96,58 @@ def main() -> None:
     assert constant.debiaser.neutral_observations == 1
     assert innovation.debiaser.neutral_observations == 1
 
+    # Build the state-dependent shape only from the committed analysis of the
+    # released physical traces. With matched neutral observations, an affine
+    # predictor must generalize across local parameter states better than an
+    # intercept-only predictor.
+    import json
+    profile = DillavouBiasProfile.from_state_dependence_report(
+        json.loads(args.profile_json.read_text()),
+        source=str(args.profile_json.resolve()),
+    )
+    profile_constant = DillavouUpdateCorrector(
+        mode="constant", bias_ratio=0.2, predictor_rate=0.2,
+        calibration_steps=1, neutral_cadence=1,
+        empirical_profile=profile, seed=67)
+    profile_innovation = DillavouUpdateCorrector(
+        mode="innovation", bias_ratio=0.2, predictor_rate=0.2,
+        calibration_steps=1, neutral_cadence=1,
+        empirical_profile=profile, seed=67)
+    parameter_scale = [
+        value.square().mean().sqrt().clamp_min(1e-6)
+        for value in parameters_a
+    ]
+    for _ in range(12):
+        for displacement in torch.linspace(-1.0, 1.0, 21):
+            state = [
+                value + displacement * scale
+                for value, scale in zip(parameters_a, parameter_scale)
+            ]
+            profile_constant.apply(clean_a, state)
+            profile_innovation.apply(clean_a, state)
+    assert (
+        profile_constant.debiaser.neutral_observations
+        == profile_innovation.debiaser.neutral_observations
+    )
+    profile_constant.neutral_cadence = 0
+    profile_innovation.neutral_cadence = 0
+    held_state = [
+        value - 0.55 * scale
+        for value, scale in zip(parameters_a, parameter_scale)
+    ]
+    held_constant = profile_constant.apply(clean_b, held_state)
+    held_innovation = profile_innovation.apply(clean_b, held_state)
+    held_constant_error = sum(
+        float((actual - target).square().sum())
+        for actual, target in zip(held_constant, clean_b)
+    )
+    held_innovation_error = sum(
+        float((actual - target).square().sum())
+        for actual, target in zip(held_innovation, clean_b)
+    )
+    assert held_innovation_error < 0.05 * held_constant_error, (
+        held_innovation_error, held_constant_error)
+
     # Integration check: the corruption is attached after Rain's hand-written
     # local EP estimator and introduces no autograd graph.
     energy, network, cost, augmented, minimizer, estimator = build_estimator(
@@ -122,6 +178,11 @@ def main() -> None:
         "innovation_relative_error": innovation_error,
         "integrated_bias_to_clean_update_rms": observed_ratio,
         "neutral_observations": constant.debiaser.neutral_observations,
+        "released_profile_normalized_offsets": profile.normalized_offsets,
+        "released_profile_normalized_state_variations": (
+            profile.normalized_state_variations),
+        "released_profile_heldout_mse_ratio_affine_over_constant": (
+            held_innovation_error / held_constant_error),
         "autodiff_used_for_learning": False,
     })
 
diff --git a/experiments/rain_ep_released_profile_r0.sh b/experiments/rain_ep_released_profile_r0.sh
new file mode 100755
index 0000000..7529fc5
--- /dev/null
+++ b/experiments/rain_ep_released_profile_r0.sh
@@ -0,0 +1,63 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+ROOT=/home/yurenh2/sdil
+AUTHOR=/scratch/yurenh2/energy-based-learning
+PYTHON=/scratch/yurenh2/venvs/burstccn/bin/python
+PROFILE="$ROOT/results/physical_bias/p0_state_dependence.json"
+OUT="$ROOT/results/ep_bias/released_profile_r0"
+mkdir -p "$OUT"
+cd "$AUTHOR"
+
+run_cell() {
+  local gpu=$1
+  local tag=$2
+  local mode=$3
+  local cadence=$4
+  CUBLAS_WORKSPACE_CONFIG=:4096:8 CUDA_VISIBLE_DEVICES="$gpu" "$PYTHON" \
+    "$ROOT/experiments/rain_ep_bias_train.py" \
+    --author-root "$AUTHOR" --device cuda \
+    --adapter dillavou --network-protocol comparative32 \
+    --beta-policy fixed_positive --beta-value 0.25 \
+    --mode "$mode" --bias-ratio 1 \
+    --dillavou-drift-ratio 0 --dillavou-profile-json "$PROFILE" \
+    --predictor-rate 1 --dillavou-calibration-steps 1 \
+    --neutral-cadence "$cadence" --epochs 1 --schedule-epochs 100 \
+    --train-limit 10000 --test-limit 2000 --batch-size 128 \
+    --training-iterations 15 --inference-iterations 60 \
+    --evaluation-split train_holdout --data-seed 6200 \
+    --seed 1988 --beta-seed 7100 --deterministic \
+    --output "$OUT/$tag.json" > "$OUT/$tag.log" 2>&1
+}
+
+run_cell 0 pep_clean clean 0 &
+run_cell 1 profile_raw raw 0 &
+run_cell 2 intercept_initial constant 0 &
+run_cell 3 sdil_initial innovation 0 &
+run_cell 4 intercept_c10 constant 10 &
+run_cell 5 sdil_c10 innovation 10 &
+run_cell 6 intercept_c50 constant 50 &
+run_cell 7 sdil_c50 innovation 50 &
+wait
+
+"$PYTHON" - "$OUT" <<'PY'
+import json
+from pathlib import Path
+import sys
+
+root = Path(sys.argv[1])
+summary = {}
+for path in sorted(root.glob("*.json")):
+    if path.name == "summary.json":
+        continue
+    report = json.loads(path.read_text())
+    summary[path.stem] = {
+        "test_accuracy": report["final"]["test_accuracy"],
+        "test_cost": report["final"]["test_cost"],
+        "finite": report["final"]["finite"],
+        "wall_seconds": report["final"]["wall_seconds"],
+        "corrector": report["final"]["corrector"],
+    }
+(root / "summary.json").write_text(json.dumps(summary, indent=2) + "\n")
+print(json.dumps(summary, indent=2))
+PY
diff --git a/sdil/rain_ep_adapter.py b/sdil/rain_ep_adapter.py
index 3711a6a..b16e96d 100644
--- a/sdil/rain_ep_adapter.py
+++ b/sdil/rain_ep_adapter.py
@@ -8,8 +8,10 @@ mode is introduced here.
 
 from __future__ import annotations
 
+from dataclasses import dataclass
+import math
 from types import MethodType
-from typing import Iterable
+from typing import Iterable, Mapping
 
 import torch
 
@@ -22,6 +24,70 @@ from sdil.two_state_debias import (
 Tensor = torch.Tensor
 
 
+@dataclass(frozen=True)
+class DillavouBiasProfile:
+    """Dimensionless per-edge affine shapes fitted from released drift data."""
+
+    normalized_offsets: tuple[float, ...]
+    normalized_state_variations: tuple[float, ...]
+    source: str
+
+    def __post_init__(self) -> None:
+        if not self.normalized_offsets:
+            raise ValueError("Dillavou profile must contain at least one edge")
+        if len(self.normalized_offsets) != len(
+            self.normalized_state_variations
+        ):
+            raise ValueError("Dillavou profile arrays disagree")
+        values = self.normalized_offsets + self.normalized_state_variations
+        if not all(math.isfinite(value) for value in values):
+            raise ValueError("Dillavou profile contains a nonfinite value")
+
+    @classmethod
+    def from_state_dependence_report(
+        cls, report: Mapping, *, source: str
+    ) -> "DillavouBiasProfile":
+        offsets = []
+        state_variations = []
+        pairs = report.get("pairs", {})
+        for pair_name in sorted(pairs):
+            pair = pairs[pair_name]
+            reference = pair["reference_gate"]
+            affine = pair["local_affine_model"]
+            pair_offsets = affine["bias_at_reference_v_per_s"]
+            pair_slopes = affine["local_slopes_per_s"]
+            gates_by_edge = [[], []]
+            for trace in pair["traces"]:
+                gates_by_edge[0].extend(trace["retained_gate_minus"])
+                gates_by_edge[1].extend(trace["retained_gate_plus"])
+            for edge in range(2):
+                observed_range = max(
+                    abs(float(value) - float(reference[edge]))
+                    for value in gates_by_edge[edge]
+                )
+                offsets.append(float(pair_offsets[edge]))
+                state_variations.append(
+                    float(pair_slopes[edge]) * observed_range)
+        offset_rms = math.sqrt(
+            sum(value * value for value in offsets) / len(offsets))
+        if offset_rms <= 0.0:
+            raise ValueError("released Dillavou offsets have zero RMS")
+        return cls(
+            normalized_offsets=tuple(value / offset_rms for value in offsets),
+            normalized_state_variations=tuple(
+                value / offset_rms for value in state_variations),
+            source=source,
+        )
+
+    def as_dict(self) -> dict:
+        return {
+            "normalized_offsets": list(self.normalized_offsets),
+            "normalized_state_variations": list(
+                self.normalized_state_variations),
+            "source": self.source,
+        }
+
+
 class LocalStructuredBias:
     """Fixed per-element bias affine in a local first-state measurement."""
 
@@ -239,12 +305,16 @@ class DillavouUpdateCorrector:
         calibration_steps: int = 1,
         neutral_cadence: int = 1,
         drift_ratio: float = 0.0,
+        empirical_profile: DillavouBiasProfile | None = None,
         seed: int = 1729,
     ) -> None:
         if mode not in self.MODES:
             raise ValueError(f"unrecognized correction mode {mode}")
         if bias_ratio < 0.0 or drift_ratio < 0.0:
             raise ValueError("Dillavou bias ratios must be nonnegative")
+        if empirical_profile is not None and drift_ratio != 0.0:
+            raise ValueError(
+                "free drift ratio cannot be combined with an empirical profile")
         if not 0.0 < predictor_rate <= 1.0:
             raise ValueError("predictor rate must lie in (0, 1]")
         if neutral_cadence < 0:
@@ -257,6 +327,7 @@ class DillavouUpdateCorrector:
         self.calibration_steps = calibration_steps
         self.neutral_cadence = neutral_cadence
         self.drift_ratio = drift_ratio
+        self.empirical_profile = empirical_profile
         self.seed = seed
         self.steps = 0
         self.debiaser: LocalAffineDebiaser | None = None
@@ -293,14 +364,41 @@ class DillavouUpdateCorrector:
         ):
             gradient_scale = gradient.square().mean().sqrt().clamp_min(1e-12)
             parameter_scale = parameter.square().mean().sqrt().clamp_min(1e-6)
-            offset_pattern = self._unit_pattern(
-                gradient, float(self.seed + 97 * index))
-            slope_pattern = self._unit_pattern(
-                gradient, float(self.seed + 193 * index + 41))
+            if self.empirical_profile is None:
+                offset_pattern = self._unit_pattern(
+                    gradient, float(self.seed + 97 * index))
+                slope_pattern = self._unit_pattern(
+                    gradient, float(self.seed + 193 * index + 41))
+                offset_gain = self.bias_ratio
+                slope_gain = self.drift_ratio
+            else:
+                flat_index = torch.arange(
+                    gradient.numel(), dtype=torch.int64,
+                    device=gradient.device).reshape(gradient.shape)
+                profile_index = torch.remainder(
+                    flat_index + self.seed + 97 * index,
+                    len(self.empirical_profile.normalized_offsets))
+                offset_values = torch.as_tensor(
+                    self.empirical_profile.normalized_offsets,
+                    dtype=gradient.dtype, device=gradient.device)
+                slope_values = torch.as_tensor(
+                    self.empirical_profile.normalized_state_variations,
+                    dtype=gradient.dtype, device=gradient.device)
+                offset_pattern = offset_values[profile_index]
+                slope_pattern = slope_values[profile_index]
+                # Small tensors need not contain the four profiles equally.
+                # Renormalize only their common amplitude; the measured
+                # per-profile slope/offset ratios remain unchanged.
+                pattern_rms = offset_pattern.square().mean().sqrt().clamp_min(
+                    1e-30)
+                offset_pattern = offset_pattern / pattern_rms
+                slope_pattern = slope_pattern / pattern_rms
+                offset_gain = self.bias_ratio
+                slope_gain = self.bias_ratio
             self._offsets.append(
-                self.bias_ratio * gradient_scale * offset_pattern)
+                offset_gain * gradient_scale * offset_pattern)
             self._slopes.append(
-                self.drift_ratio * gradient_scale * slope_pattern)
+                slope_gain * gradient_scale * slope_pattern)
             self._centers.append(parameter.clone())
             self._scales.append(parameter_scale.clone())
             feature_scales.append(torch.ones_like(parameter_scale))
@@ -395,10 +493,15 @@ class DillavouUpdateCorrector:
         self.last_diagnostics = {
             "step": self.steps,
             "bias_model": (
-                "dillavou_constant_update"
+                "dillavou_released_affine_update"
+                if self.empirical_profile is not None
+                else "dillavou_constant_update"
                 if self.drift_ratio == 0.0
                 else "dillavou_plus_local_state_drift"
             ),
+            "bias_profile_source": (
+                None if self.empirical_profile is None
+                else self.empirical_profile.source),
             "clean_update_rms": clean_rms,
             "bias_update_rms": bias_rms,
             "residual_update_rms": residual_rms,