#!/usr/bin/env python3
"""Phase 4M: isolated KC->MBON efficacy mechanisms on frozen MaleCNS inputs."""
from __future__ import annotations

import argparse
import gc
import json
import resource
import shutil
import subprocess
import sys
import time
from collections import defaultdict
from pathlib import Path
from typing import Any

import numpy as np
import pandas as pd
import torch
import yaml

from phase4h_audit import sha256, verify_checksum, write_json
from phase4i import Validator, phase4c_regression
from phase4j import apply_v2, effect, verify_frozen
from phase4k import dataset_checks, seal
from phase4l import OdorSimulator, pn_drive


MECHANISMS = ("additive", "multiplicative", "bounded_multiplicative")


def frozen_checks(cfg: dict[str, Any]) -> dict[str, Any]:
    earlier = verify_frozen(cfg); later = {}
    for key, control in (("phase4j_root", "manifests/checksums.sha256"), ("phase4j_model_root", "checksums.sha256"),
                         ("phase4k_root", "manifests/checksums.sha256"), ("phase4k_model_root", "checksums.sha256"),
                         ("phase4l_root", "manifests/checksums.sha256"), ("phase4l_model_root", "checksums.sha256")):
        root = Path(cfg["paths"][key]); later[key] = verify_checksum(root, root / control)
    report = {"phase4c_to_4i": earlier["all_passed"], "phase4j_to_4l": {k: v["passed"] for k, v in later.items()}}
    report["all_passed"] = bool(earlier["all_passed"] and all(v["passed"] for v in later.values()))
    if not report["all_passed"]: raise RuntimeError("frozen Phase 4C-4L checksum failure")
    return report


def frozen_hashes(cfg: dict[str, Any]) -> dict[str, str]:
    p = cfg["paths"]; root = Path(p["connectome_root"]); files = [root / "gpu/csr-indptr.int64.npy", root / "gpu/csr-indices.int32.npy", root / "gpu/csr-weights.int64.npy", root / "neuron-index/unified-metadata.parquet",
        Path(p["phase4h_root"]) / "plastic-mask/plastic-edge-candidates.parquet",
        Path(p["phase4j_root"]) / "audit/phase4j.py", Path(p["phase4j_model_root"]) / "phase4j.yaml",
        Path(p["phase4l_root"]) / "odor-encoder/frozen-odor-panel.parquet", Path(p["phase4l_root"]) / "odor-encoder/mapped-glomerular-vectors.npz",
        Path(p["phase4l_root"]) / "receptor-mapping/olfactory-mapping.parquet", Path(p["phase4l_root"]) / "manifests/checksums.sha256",
        Path(p["phase4l_model_root"]) / "checksums.sha256"]
    for f in files:
        if not f.is_file(): raise FileNotFoundError(f)
    return {str(f): sha256(f) for f in files}


def load_odor_inputs(cfg: dict[str, Any]) -> tuple[pd.DataFrame, dict[str, np.ndarray], dict[str, Any], pd.DataFrame]:
    exp = Path(cfg["paths"]["phase4l_root"]); panel = pd.read_parquet(exp / "odor-encoder/frozen-odor-panel.parquet"); vecfile = np.load(exp / "odor-encoder/mapped-glomerular-vectors.npz"); vectors = {str(oid): vecfile[str(oid)].copy() for oid in panel.odor_id}
    selection = json.loads((exp / "odor-encoder/selection.json").read_text()); mapping = pd.read_parquet(exp / "receptor-mapping/olfactory-mapping.parquet")
    if selection["ODOR_A"]["odor_id"] not in vectors or selection["ODOR_B"]["odor_id"] not in vectors: raise RuntimeError("frozen odor identity mismatch")
    return panel, vectors, selection, mapping


def anatomical_audit(v: Validator, cfg: dict[str, Any], exp: Path) -> tuple[np.ndarray, np.ndarray]:
    cand = v.cand; raw = np.asarray(v.rt.raw_np[v.edge_ids], dtype=np.int64); sign = v.rt.signs[cand.source_neuron_idx.to_numpy(np.int32)]
    w = (raw.astype(np.float32) * sign.astype(np.float32) * np.float32(cfg["model"]["w_syn_mv_per_anatomical_synapse"])).astype(np.float32)
    base = v.base_candidate.astype(np.float32)
    if not np.array_equal(raw, cand.raw_synapse_count.to_numpy(np.int64)): raise RuntimeError("plastic mask/raw CSR count mismatch")
    if not np.array_equal(w, base): raise RuntimeError("anatomical effective weight mismatch with frozen runtime")
    if not np.all(sign == 1): raise RuntimeError("presynaptic efficacy branch requires cholinergic positive KC candidates")
    if len(cand) != 14232 or len(np.unique(v.edge_ids)) != 14232: raise RuntimeError("plastic mask changed")
    report = {"plastic_mask_edges": len(cand), "raw_count_min": int(raw.min()), "raw_count_median": float(np.median(raw)), "raw_count_max": int(raw.max()),
              "raw_count_sum": int(raw.sum()), "anatomical_weight_min_mv": float(w.min()), "anatomical_weight_max_mv": float(w.max()),
              "all_cholinergic_positive": True, "anatomical_formula_verified_exactly": True, "base_CSR_read_only": True,
              "E_birth_state": 1.0, "other_brain_edges_E": 1.0}
    write_json(exp / "anatomical-weight/summary.json", report); return raw, w


class ReplayEngine:
    """Real KC spikes + existing KC->MBON edges + frozen MBON LIF, no invented edge."""

    def __init__(self, sim: OdorSimulator, targets: np.ndarray):
        self.sim = sim; self.targets = np.asarray(targets, np.int32); self.dt = sim.dt; self.rt = sim.rt; v = sim.v
        indices = v.rt.indices_np; sources = sim.sources; mask = np.isin(sources, v.kc) & np.isin(indices, self.targets); self.ids = np.flatnonzero(mask)
        self.src = sources[self.ids].astype(np.int32); self.target = indices[self.ids].astype(np.int32)
        pos = {int(x): i for i, x in enumerate(self.targets)}; self.target_local = np.fromiter((pos[int(x)] for x in self.target), np.int32, len(self.target))
        all_base = v.base.detach().cpu().numpy(); self.base = all_base[self.ids].astype(np.float32)
        can_index = {int(eid): i for i, eid in enumerate(v.edge_ids)}; self.candidate_local = np.fromiter((can_index.get(int(x), -1) for x in self.ids), np.int32, len(self.ids))
        self.by_source = {int(s): np.flatnonzero(self.src == s) for s in np.unique(self.src)}
        if not np.all(self.base > 0): raise RuntimeError("isolated KC replay contains nonpositive source weights")

    def run(self, spikes: pd.DataFrame, delta: np.ndarray, duration_ms: float) -> dict[str, Any]:
        weights = self.base.copy(); valid = self.candidate_local >= 0; weights[valid] += delta[self.candidate_local[valid]]
        grouped = {int(k): g.neuron_idx.to_numpy(np.int32) for k, g in spikes.groupby("step")}; nsteps = int(round(duration_ms / self.dt)); n = len(self.targets); q = np.zeros((self.rt.queue_len, n), np.float32)
        conductance = np.zeros(n, np.float32); voltage = np.full(n, self.rt.rest, np.float32); refractory = np.zeros(n, np.int32); counts = np.zeros(n, np.int32); first = np.full(n, -1, np.int32)
        integrated_g = np.zeros(n, np.float64); transmitted = 0.0; delivered = 0.0; wall0 = time.perf_counter()
        for step in range(nsteps):
            slot = step % self.rt.queue_len; conductance = conductance * np.float32(self.rt.g_decay) + q[slot]; delivered += float(q[slot].sum(dtype=np.float64)); q[slot].fill(0)
            integrated_g += conductance; refractory = np.maximum(refractory, 1) - 1; avail = refractory == 0
            voltage[avail] += np.float32(self.rt.alpha) * (conductance[avail] - (voltage[avail] - np.float32(self.rt.rest))); voltage[~avail] = np.float32(self.rt.reset_value)
            fired = avail & (voltage >= np.float32(self.rt.threshold)); counts += fired; first[(first < 0) & fired] = step; voltage[fired] = np.float32(self.rt.reset_value); refractory[fired] = self.rt.refractory_steps
            for source in grouped.get(step, np.empty(0, np.int32)):
                rel = self.by_source.get(int(source))
                if rel is not None:
                    np.add.at(q[(step + self.rt.delay_steps) % self.rt.queue_len], self.target_local[rel], weights[rel]); transmitted += float(weights[rel].sum(dtype=np.float64))
        return {"counts": counts, "first": first, "final_g": conductance, "final_v": voltage, "integrated_g": integrated_g,
                "transmitted_weighted_events_mv": transmitted, "delivered_input_mv": delivered, "wall_seconds": time.perf_counter() - wall0}


def replay_effect(pre: dict[str, Any], post: dict[str, Any], dt: float) -> dict[str, Any]:
    def first_ms(x: dict[str, Any]) -> float | None:
        a = x["first"]; return None if np.all(a < 0) else float(a[a >= 0].min() * dt)
    a, b = first_ms(pre), first_ms(post)
    return {"delta_g_max_abs_mv": float(np.abs(post["final_g"] - pre["final_g"]).max()), "delta_g_mean_abs_mv": float(np.abs(post["final_g"] - pre["final_g"]).mean()),
            "delta_v_max_abs_mv": float(np.abs(post["final_v"] - pre["final_v"]).max()), "delta_v_mean_abs_mv": float(np.abs(post["final_v"] - pre["final_v"]).mean()),
            "integrated_g_delta_sum_mv_steps": float((post["integrated_g"] - pre["integrated_g"]).sum()),
            "MBON_spikes_pre": int(pre["counts"].sum()), "MBON_spikes_post": int(post["counts"].sum()), "MBON_spike_count_delta": int(post["counts"].sum() - pre["counts"].sum()),
            "MBON_spike_L1": int(np.abs(post["counts"] - pre["counts"]).sum()), "first_latency_pre_ms": a, "first_latency_post_ms": b,
            "first_latency_delta_ms": None if a is None or b is None else float(b - a), "transmitted_delta_mv": float(post["transmitted_weighted_events_mv"] - pre["transmitted_weighted_events_mv"]),
            "transmitted_reduction_percent": float(100 * (pre["transmitted_weighted_events_mv"] - post["transmitted_weighted_events_mv"]) / pre["transmitted_weighted_events_mv"]) if pre["transmitted_weighted_events_mv"] else None,
            "MBON_total_KC_input_reduction_percent": float(100 * (pre["delivered_input_mv"] - post["delivered_input_mv"]) / pre["delivered_input_mv"]) if pre["delivered_input_mv"] else None,
            "isolated_input_scope": "all real KC->target-MBON edges, not all whole-brain MBON inputs"}


def latency_direction(row: dict[str, Any]) -> bool:
    value = row.get("first_latency_delta_ms"); return bool(row["MBON_spike_count_delta"] < 0 or (value is not None and value > 0))


def compare_A_B(a: dict[str, Any], b: dict[str, Any]) -> bool:
    if not latency_direction(a): return False
    if a["MBON_spike_L1"] > b["MBON_spike_L1"]: return True
    al = 0.0 if a["first_latency_delta_ms"] is None else max(float(a["first_latency_delta_ms"]), 0.0)
    bl = 0.0 if b["first_latency_delta_ms"] is None else max(float(b["first_latency_delta_ms"]), 0.0)
    return bool(al > bl and a["MBON_spike_L1"] >= b["MBON_spike_L1"])


def main_experiment(v: Validator, sim: OdorSimulator, cfg: dict[str, Any], mcfg: dict[str, Any], exp: Path,
                    panel: pd.DataFrame, vectors: dict[str, np.ndarray], selection: dict[str, Any], mapping: pd.DataFrame,
                    anatomical: np.ndarray) -> tuple[dict[str, Any], dict[int, dict[int, dict[str, np.ndarray]]], ReplayEngine]:
    protocol = mcfg["protocol"]; points = list(map(int, protocol["pairing_counts"])); maxpair = max(points); seeds = list(map(int, protocol["train_seeds"])); eval_seeds = list(map(int, protocol["evaluation_seeds"]));
    if points != list(map(int, cfg["plasticity_accumulation"]["pairing_counts"])) or seeds != list(map(int, cfg["plasticity_accumulation"]["train_seeds"])) or eval_seeds != list(map(int, cfg["plasticity_accumulation"]["evaluation_seeds"])): raise RuntimeError("Phase 4L frozen pairing protocol mismatch")
    compartment = str(protocol["target_compartment_verbatim"]); dan = v.dan_targets(str(protocol["matched_dan_type"]), compartment); duration = float(protocol["duration_ms"])
    gloms = selection["glomeruli_order"]; matched_sum = float(np.median(panel.mapped_response_sum)); pn = {}
    for label in ("ODOR_A", "ODOR_B"):
        oid = selection[label]["odor_id"]; target, rate, _ = pn_drive(oid, vectors[oid], mapping, gloms, str(protocol["input_mode"]), matched_sum, cfg); pn[label] = (target, rate)
    src = v.cand.source_neuron_idx.to_numpy(np.int32); lookup = np.full(v.n, -1, np.int32); lookup[v.kc] = np.arange(len(v.kc), dtype=np.int32); loc = lookup[src]
    if np.any(loc < 0): raise RuntimeError("non-KC source in plastic mask")
    comp = v.cand.target_compartment_verbatim.astype(str).to_numpy() == compartment; k = float(cfg["eligibility_v2"]["eta_effective_mv"]) / float(cfg["model"]["w_syn_mv_per_anatomical_synapse"])
    if k <= 0 or k >= 1: raise RuntimeError("invalid frozen eta/Wsyn dimensionless conversion")
    bound = 1.0 - float(cfg["eligibility_v2"]["max_fractional_depression"]); checkpoints: dict[int, dict[int, dict[str, np.ndarray]]] = {}; histories = []; update_times = []; training_walls = []; reproduction = []
    zero = np.zeros(len(v.edge_ids), np.float32); target_mbon = np.unique(v.cand.loc[comp, "target_neuron_idx"].to_numpy(np.int32)); replay = ReplayEngine(sim, target_mbon)
    # Baseline replay uses the same deterministic evaluation KC spike trains for all mechanisms.
    evaluation = {}; pre_replay = {}
    for es in eval_seeds:
        for label, (targets, rates) in pn.items():
            run = sim.run_rates(targets, rates, es, duration, collect_kc_spikes=True); evaluation[(es, label)] = run["kc_spikes"]
            run["kc_spikes"].to_parquet(exp / "isolated-replay" / f"KC-spikes-{label}-seed-{es}.parquet", index=False)
            pre_replay[(es, label)] = replay.run(run["kc_spikes"], zero, duration)
    for ts in seeds:
        delta_add = zero.copy(); em = np.ones(len(v.edge_ids), np.float32); eb = np.ones(len(v.edge_ids), np.float32); checkpoints[ts] = {}; p4l = np.load(Path(cfg["paths"]["phase4l_root"]) / "plasticity-magnitude" / f"seed-{ts}-overlays.npz")
        if not np.array_equal(p4l["edge_id"], v.edge_ids): raise RuntimeError("Phase 4L additive checkpoint edge IDs differ")
        for pairing in range(1, maxpair + 1):
            run = sim.run_rates(*pn["ODOR_A"], ts + pairing, duration, delta=delta_add, dan_targets=dan, dan_pulses=[float(protocol["dan_time_ms"])])
            training_walls.append(run["wall_seconds"]); trace = run["kc_trace_at_dan"]; elig = trace[loc].astype(np.float32); active = comp & (elig > 0)
            old_m, old_b = em[active].copy(), eb[active].copy(); t0 = time.perf_counter()
            delta_add, _ = apply_v2(v, trace, compartment, 1.0, delta_add, cfg)
            factor = np.float32(1.0) - np.float32(k) * elig[active]; em[active] = np.maximum(np.float32(0.0), em[active] * factor)
            eb[active] = np.maximum(np.float32(bound), eb[active] * factor); update_times.append(time.perf_counter() - t0)
            for local, om, ob in zip(np.flatnonzero(active), old_m, old_b):
                histories.append({"train_seed": ts, "pairing": pairing, "edge_id": int(v.edge_ids[local]), "eligibility": float(elig[local]),
                                  "E_multiplicative_before": float(om), "E_multiplicative_after": float(em[local]), "E_bounded_before": float(ob), "E_bounded_after": float(eb[local])})
            if pairing in points:
                saved = p4l[f"pairing_{pairing}"].astype(np.float32); error = float(np.max(np.abs(saved - delta_add)))
                reproduction.append({"train_seed": ts, "pairing": pairing, "max_abs_error_mv": error, "exact": bool(np.array_equal(saved, delta_add))})
                if error > 1e-6: raise RuntimeError(f"Phase 4L additive training reproduction divergence seed={ts} pairing={pairing}: {error}")
                checkpoints[ts][pairing] = {"additive": delta_add.copy(), "multiplicative": em.copy(), "bounded_multiplicative": eb.copy(), "eligibility_at_last_DAN": elig.copy()}
                np.savez_compressed(exp / "efficacy-overlay" / f"E-seed-{ts}-pairings-{pairing}.npz", E_multiplicative=em, E_bounded_multiplicative=eb)
        p4l.close()
    pd.DataFrame(histories).to_parquet(exp / "efficacy-overlay/update-history.parquet", index=False)
    write_json(exp / "regression/phase4l-additive-reproduction.json", {"all_passed": all(x["max_abs_error_mv"] <= 1e-6 for x in reproduction), "checks": reproduction})
    if not np.all(np.isfinite([h["E_multiplicative_after"] for h in histories])) or not np.all(np.isfinite([h["E_bounded_after"] for h in histories])):
        raise RuntimeError("nonfinite efficacy")
    rows = []; replay_times = []
    for ts, es in zip(seeds, eval_seeds):
        for pairing in points:
            state = checkpoints[ts][pairing]
            for mechanism in MECHANISMS:
                if mechanism == "additive": delta = state["additive"]; efficacy = 1.0 + delta / anatomical
                else:
                    efficacy = state[mechanism]; delta = anatomical * (efficacy - 1.0)
                if not np.all(np.isfinite(efficacy)) or np.any(efficacy < -1e-6) or np.any(efficacy > 1.0 + 1e-6): raise RuntimeError("invalid E state")
                changed = np.abs(delta) > 1e-8; e_changed = efficacy[changed]; weighted_E = float(np.sum(anatomical[changed] * e_changed) / np.sum(anatomical[changed])) if np.any(changed) else 1.0
                pairrow = {"train_seed": ts, "eval_seed": es, "pairings": pairing, "mechanism": mechanism, "updated_edges": int(changed.sum()),
                    "E_mean_all_mask": float(np.mean(efficacy)), "E_median_all_mask": float(np.median(efficacy)), "E_min_all_mask": float(np.min(efficacy)),
                    "E_mean_changed": float(np.mean(e_changed)) if len(e_changed) else 1.0, "E_median_changed": float(np.median(e_changed)) if len(e_changed) else 1.0,
                    "E_anatomical_weighted_changed": weighted_E, "effective_weight_delta_sum_mv": float(delta.sum()),
                    "effective_weight_reduction_percent_changed": float(100 * (1 - weighted_E)), "saturation_fraction_mask": float(np.count_nonzero(efficacy <= (bound + 1e-6)) / len(efficacy))}
                result = {}
                for label in ("ODOR_A", "ODOR_B"):
                    t0 = time.perf_counter(); post = replay.run(evaluation[(es, label)], delta, duration); replay_times.append(time.perf_counter() - t0)
                    result[label] = replay_effect(pre_replay[(es, label)], post, sim.dt)
                pairrow.update({f"A_{k}": value for k, value in result["ODOR_A"].items() if k != "isolated_input_scope"})
                pairrow.update({f"B_{k}": value for k, value in result["ODOR_B"].items() if k != "isolated_input_scope"})
                pairrow["A_specific_functional"] = compare_A_B(result["ODOR_A"], result["ODOR_B"])
                rows.append(pairrow)
    table = pd.DataFrame(rows); table.to_parquet(exp / "isolated-replay/mechanism-comparison.parquet", index=False)
    curve = []
    for (mechanism, pairing), frame in table.groupby(["mechanism", "pairings"]):
        curve.append({"mechanism": mechanism, "pairings": int(pairing), "seeds": int(len(frame)), "A_specific_functional_seeds": int(frame.A_specific_functional.sum()),
            "E_mean_changed": float(frame.E_mean_changed.mean()), "E_median_changed": float(frame.E_median_changed.mean()), "E_anatomical_weighted_changed": float(frame.E_anatomical_weighted_changed.mean()),
            "E_min_mask": float(frame.E_min_all_mask.min()), "effective_weight_reduction_percent_changed": float(frame.effective_weight_reduction_percent_changed.mean()),
            "A_transmission_reduction_percent": float(frame.A_transmitted_reduction_percent.mean()), "B_transmission_reduction_percent": float(frame.B_transmitted_reduction_percent.mean()),
            "A_total_KC_input_reduction_percent": float(frame.A_MBON_total_KC_input_reduction_percent.mean()), "B_total_KC_input_reduction_percent": float(frame.B_MBON_total_KC_input_reduction_percent.mean()),
            "A_MBON_delta_g_max_abs_mv": float(frame.A_delta_g_max_abs_mv.mean()), "A_MBON_delta_v_max_abs_mv": float(frame.A_delta_v_max_abs_mv.mean()),
            "A_MBON_spike_count_delta_mean": float(frame.A_MBON_spike_count_delta.mean()), "A_MBON_spike_L1_mean": float(frame.A_MBON_spike_L1.mean()),
            "A_MBON_first_latency_delta_ms_mean": float(frame.A_first_latency_delta_ms.mean()) if frame.A_first_latency_delta_ms.notna().any() else None,
            "B_MBON_spike_L1_mean": float(frame.B_MBON_spike_L1.mean()), "saturation_fraction_mask": float(frame.saturation_fraction_mask.mean())})
    gate = {}; majority = int(mcfg["gates"]["seed_majority"])
    for mech in MECHANISMS:
        successes = [x for x in curve if x["mechanism"] == mech and x["A_specific_functional_seeds"] >= majority]
        gate[mech] = {"pass": bool(successes), "first_passing_pairing": min((x["pairings"] for x in successes), default=None), "seed_majority": majority,
                      "endpoint_100_specific_functional_seeds": next(x["A_specific_functional_seeds"] for x in curve if x["mechanism"] == mech and x["pairings"] == 100)}
    primary = gate["multiplicative"]["pass"] or gate["bounded_multiplicative"]["pass"]
    report = {"status": "PASS-KC-MBON-FUNCTIONAL-TRANSFER-v2" if primary else "FAIL-KC-MBON-FUNCTIONAL-TRANSFER-v2", "presynaptic_functional_transfer_pass": primary,
        "mechanism_gates": gate, "curve": curve, "same_KC_spikes_DAN_timing_mask_and_seeds": True, "phase4l_additive_reproduced": True,
        "postsynaptic_branch_run": False, "no_dense_time_by_edge_state": True,
        "performance": {"mean_training_brain_wall_seconds_per_pairing": float(np.mean(training_walls)), "mean_update_seconds_per_pairing": float(np.mean(update_times)),
                        "mean_isolated_replay_seconds": float(np.mean(replay_times)), "E_overlay_memory_bytes_two_float32_arrays": 2 * len(v.edge_ids) * 4}}
    write_json(exp / "isolated-replay/summary.json", report)
    return report, checkpoints, replay


def anatomy_vs_efficacy(v: Validator, sim: OdorSimulator, cfg: dict[str, Any], exp: Path, anatomical: np.ndarray,
                        checkpoints: dict[int, dict[int, dict[str, np.ndarray]]], mcfg: dict[str, Any]) -> dict[str, Any]:
    targets = v.cand.target_neuron_idx.to_numpy(np.int32); base = v.base.detach().cpu().numpy(); positive = np.maximum(base, 0.0)
    target_total = np.bincount(v.rt.indices_np, weights=positive.astype(np.float64), minlength=v.n)
    frames = []
    for ts in map(int, mcfg["protocol"]["train_seeds"]):
        state = checkpoints[ts][100]
        for mechanism in ("multiplicative", "bounded_multiplicative"):
            E = state[mechanism]; frame = v.cand[["edge_id", "source_neuron_idx", "source_body_id", "target_neuron_idx", "target_body_id", "raw_synapse_count", "target_compartment_verbatim"]].copy()
            frame["train_seed"] = ts; frame["mechanism"] = mechanism; frame["anatomical_weight_mv"] = anatomical; frame["E_baseline"] = 1.0; frame["E_current"] = E
            frame["effective_weight_mv"] = anatomical * E; frame["target_total_positive_anatomical_input_mv"] = target_total[targets]
            frame["target_MBON_total_input_fraction_static"] = anatomical * E / np.maximum(target_total[targets], 1e-12)
            frame["eligibility_at_last_DAN"] = state["eligibility_at_last_DAN"]
            frame["DAN_channel"] = mcfg["protocol"]["matched_dan_type"]; frames.append(frame)
    out = pd.concat(frames, ignore_index=True); out.to_parquet(exp / "anatomy-vs-efficacy/anatomy-vs-efficacy.parquet", index=False)
    from scipy.stats import spearmanr
    by = out[(out.train_seed == int(mcfg["protocol"]["train_seeds"][0])) & (out.mechanism == "bounded_multiplicative")]
    rho = spearmanr(by.raw_synapse_count, by.effective_weight_mv).statistic
    report = {"rows": int(len(out)), "unique_edges": int(out.edge_id.nunique()), "raw_count_vs_effective_weight_spearman_seed1_bounded": float(rho),
        "high_anatomical_count_implies_high_current_effective_in_this_assumed_model_only": bool(rho > 0.9),
        "physiological_influence_inferred_from_EM_count": False, "static_fraction_is_not_dynamic_MBON_input": True}
    write_json(exp / "anatomy-vs-efficacy/summary.json", report); return report


def run_wholebrain(v: Validator, sim: OdorSimulator, cfg: dict[str, Any], mcfg: dict[str, Any], exp: Path, panel: pd.DataFrame,
                   vectors: dict[str, np.ndarray], selection: dict[str, Any], mapping: pd.DataFrame,
                   anatomical: np.ndarray, checkpoints: dict[int, dict[int, dict[str, np.ndarray]]], isolated: dict[str, Any]) -> dict[str, Any]:
    passed = [x for x in ("multiplicative", "bounded_multiplicative") if isolated["mechanism_gates"][x]["pass"]]
    if not passed:
        out = {"status": "NOT-RUN-ISOLATED-ENDPOINT-GATE-BLOCKED", "tested_mechanisms": [], "feedback_masks_memory": None}; write_json(exp / "whole-brain-replay/summary.json", out); return out
    duration = float(mcfg["protocol"]["duration_ms"]); gloms = selection["glomeruli_order"]; matched = float(np.median(panel.mapped_response_sum)); pn = {}
    for label in ("ODOR_A", "ODOR_B"):
        oid = selection[label]["odor_id"]; pn[label] = pn_drive(oid, vectors[oid], mapping, gloms, str(mcfg["protocol"]["input_mode"]), matched, cfg)[:2]
    rows = []; walls = []; zero = np.zeros(len(v.edge_ids), np.float32)
    for ts, es in zip(map(int, mcfg["protocol"]["train_seeds"]), map(int, mcfg["protocol"]["evaluation_seeds"])):
        baseline = {label: sim.run_rates(*pn[label], es, duration, delta=zero) for label in pn}
        for mechanism in passed:
            E = checkpoints[ts][100][mechanism]; delta = anatomical * (E - 1)
            effects = {}
            for label in pn:
                post = sim.run_rates(*pn[label], es, duration, delta=delta); walls.append(post["wall_seconds"]); effects[label] = effect(baseline[label], post)
                rows.append({"train_seed": ts, "eval_seed": es, "mechanism": mechanism, "odor": label, "MBON_spike_L1": effects[label]["MBON_spike_l1"],
                    "MBON_g_max_abs_mv": effects[label]["MBON_g_max_abs"], "MBON_v_max_abs_mv": effects[label]["MBON_v_max_abs"],
                    "MBON_latency_delta_ms": effects[label]["MBON_first_latency_delta_ms"], "downstream_spike_L1": effects[label]["downstream_spike_l1"],
                    "KC_spike_L1": int(np.abs(post["kc_counts"] - baseline[label]["kc_counts"]).sum()),
                    "finite": bool(post["finite"]), "active_neuron_fraction": float(np.count_nonzero(post["counts"]) / v.n), "wall_seconds": post["wall_seconds"]})
    frame = pd.DataFrame(rows); frame.to_parquet(exp / "whole-brain-replay/results.parquet", index=False); gates = {}
    for mechanism in passed:
        sub = frame[frame.mechanism == mechanism]; a = sub[sub.odor == "ODOR_A"].sort_values("train_seed"); b = sub[sub.odor == "ODOR_B"].sort_values("train_seed")
        wins = int(((a.MBON_spike_L1.to_numpy() > b.MBON_spike_L1.to_numpy()) & (a.MBON_spike_L1.to_numpy() > 0)).sum())
        gates[mechanism] = {"A_specific_functional_seeds": wins, "pass": bool(wins >= int(mcfg["gates"]["seed_majority"]) and sub.finite.all() and sub.active_neuron_fraction.max() < .5),
            "A_MBON_spike_L1_mean": float(a.MBON_spike_L1.mean()), "B_MBON_spike_L1_mean": float(b.MBON_spike_L1.mean()),
            "A_downstream_spike_L1_mean": float(a.downstream_spike_L1.mean()), "B_downstream_spike_L1_mean": float(b.downstream_spike_L1.mean()),
            "A_KC_spike_L1_mean": float(a.KC_spike_L1.mean()), "B_KC_spike_L1_mean": float(b.KC_spike_L1.mean())}
    status = "PASS-WHOLE-BRAIN-KCMBON-TRANSFER" if any(x["pass"] for x in gates.values()) else "WHOLE-BRAIN-FEEDBACK-MASKS-MEMORY"
    out = {"status": status, "tested_mechanisms": passed, "gates": gates, "feedback_masks_memory": status != "PASS-WHOLE-BRAIN-KCMBON-TRANSFER", "mean_post_brain_wall_seconds": float(np.mean(walls))}
    write_json(exp / "whole-brain-replay/summary.json", out); return out


def conditioning_gate(cfg: dict[str, Any], mcfg: dict[str, Any], exp: Path, isolated: dict[str, Any], whole: dict[str, Any], identity: dict[str, Any]) -> dict[str, Any]:
    encoder = json.loads((Path(cfg["paths"]["phase4l_root"]) / "kc-representation/summary.json").read_text()); temporal = json.loads((Path(cfg["paths"]["phase4j_root"]) / "eligibility-v2/unit-tests.json").read_text())
    gates = {"OlfactoryEncoder_v1": bool(encoder["candidate_generated"]), "TemporalEligibility_v2": bool(temporal["all_passed"]),
             "KC_sparse_code": bool(encoder["modes"]["physiological"]["gate"] and encoder["modes"]["matched_energy"]["gate"]),
             "A_B_specific_edges": bool(identity["A_specific_eligible_edges"] > 0 and identity["B_specific_eligible_edges"] > 0),
             "isolated_functional_transfer_v2": bool(isolated["presynaptic_functional_transfer_pass"]),
             "whole_brain_transfer": whole["status"] == "PASS-WHOLE-BRAIN-KCMBON-TRANSFER"}
    report = {"status": "OPEN" if all(gates.values()) else "NOT-RUN-GATE-BLOCKED", "gates": gates, "formal_conditioning_run": False,
              "acquisition": "NOT-RUN" if not all(gates.values()) else "PENDING", "controls": "NOT-RUN" if not all(gates.values()) else "PENDING"}
    write_json(exp / "conditioning/gate.json", report); return report


def integrity(cfg: dict[str, Any], mcfg: dict[str, Any], base_config: Path, exp: Path, before_hash: dict[str, str]) -> dict[str, Any]:
    del mcfg
    gc.collect(); torch.cuda.empty_cache(); checks = frozen_checks(cfg); data = dataset_checks(cfg); regression = phase4c_regression(cfg, base_config, exp); after = frozen_hashes(cfg)
    report = {"all_passed": bool(checks["all_passed"] and data["all_passed"] and regression["all_passed"] and before_hash == after),
              "frozen_checks": checks, "phase1_to_3b_dataset": data, "phase4c_reference_regression": regression,
              "frozen_hashes_identical": before_hash == after, "before": before_hash, "after": after,
              "base_CSR_modified": False, "APL_modified": False, "eligibility_v2_modified": False, "plastic_mask_expanded": False,
              "OlfactoryEncoder_modified": False, "PlasticityAccumulation_modified": False}
    write_json(exp / "audit/integrity.json", report); return report


def main() -> None:
    ap = argparse.ArgumentParser(); ap.add_argument("--config", type=Path, required=True); args = ap.parse_args(); mcfg = yaml.safe_load(args.config.read_text()); base_config = Path(mcfg["frozen_phase4l_config"])
    cfg = yaml.safe_load(base_config.read_text()); cfg["paths"]["phase4l_root"] = cfg["paths"]["experiment_root"]; cfg["paths"]["phase4l_model_root"] = mcfg["frozen_phase4l_config"].rsplit("/", 1)[0]
    exp, model = Path(mcfg["experiment_root"]), Path(mcfg["model_root"])
    for sub in ("literature-audit", "anatomical-weight", "efficacy-overlay", "presynaptic", "postsynaptic", "anatomy-vs-efficacy", "isolated-replay", "whole-brain-replay", "conditioning", "controls", "performance", "regression", "audit", "manifests", "compiler-reports/extension-build"):
        (exp / sub).mkdir(parents=True, exist_ok=True)
    model.mkdir(parents=True, exist_ok=True)
    if not (exp / "literature-audit/phase4m_literature_audit.json").is_file(): raise RuntimeError("literature audit must predate model implementation")
    prior = frozen_checks(cfg); before_hash = frozen_hashes(cfg); write_json(exp / "audit/frozen-before.json", {"checks": prior, "hashes": before_hash})
    shutil.copy2(args.config, exp / "audit/phase4m.yaml"); shutil.copy2(Path(__file__).resolve(), exp / "audit/phase4m.py")
    panel, vectors, selection, mapping = load_odor_inputs(cfg); v = Validator(cfg, base_config, exp); sim = OdorSimulator(v, cfg, exp); raw, anatomical = anatomical_audit(v, cfg, exp)
    isolated, checkpoints, replay = main_experiment(v, sim, cfg, mcfg, exp, panel, vectors, selection, mapping, anatomical)
    anatomy = anatomy_vs_efficacy(v, sim, cfg, exp, anatomical, checkpoints, mcfg)
    post = {"status": mcfg["postsynaptic"]["current_status"], "branch_defined": True, "run": False,
        "reason": "postsynaptic nAChR evidence exists in adult Drosophila, but a quantitative rule for exact MaleCNS y5B'2a target and its receptor map is unavailable",
        "combined_presynaptic_postsynaptic_run": False, "global_MBON_threshold_modified": False}
    write_json(exp / "postsynaptic/summary.json", post)
    whole = run_wholebrain(v, sim, cfg, mcfg, exp, panel, vectors, selection, mapping, anatomical, checkpoints, isolated)
    edge_identity = json.loads((Path(cfg["paths"]["phase4l_root"]) / "eligible-edges/summary.json").read_text())
    formal = conditioning_gate(cfg, mcfg, exp, isolated, whole, edge_identity)
    # Formal conditioning implementation is intentionally gated separately: this
    # diagnostic run never creates a memory claim from the efficacy comparison.
    if formal["status"] == "OPEN":
        formal["status"] = "NOT-RUN-IMPLEMENTATION-REQUIRES-CONTROLLED-FORMAL-PROTOCOL"
        formal["formal_conditioning_run"] = False
        formal["acquisition"] = "NOT-RUN"
        formal["controls"] = "NOT-RUN"
        write_json(exp / "conditioning/gate.json", formal)
    controls = {"isolated_neutral_B": "RUN", "matched_seed_same_KC_spikes": "RUN", "D_zero": "analytic identity E=1, no dopamine update", "random_reward": "NOT-RUN-FORMAL-CONDITIONING-BLOCKED", "unrelated_DAN": "NOT-RUN-FORMAL-CONDITIONING-BLOCKED", "frozen_no_plasticity": "isolated E=1 reference replay"}
    write_json(exp / "controls/summary.json", controls)
    checkpoints_bytes = sum(p.stat().st_size for p in (exp / "efficacy-overlay").glob("E-seed-*.npz")); performance = isolated["performance"] | {"overlay_checkpoint_bytes_total": checkpoints_bytes,
        "peak_RAM_kib": int(resource.getrusage(resource.RUSAGE_SELF).ru_maxrss), "peak_VRAM_bytes": int(torch.cuda.max_memory_allocated()),
        "whole_brain_mean_post_wall_seconds": whole.get("mean_post_brain_wall_seconds"), "dense_time_by_edge_state": False}
    write_json(exp / "performance/summary.json", performance)
    candidates = {"KCMBON-SynapticEfficacy-v1": True, "PresynapticEfficacy-v1": True,
        "KCMBON-FunctionalTransfer-v2": bool(isolated["presynaptic_functional_transfer_pass"]), "MaleCNS-OlfactoryAssociativeMemory-v1-candidate": False}
    del sim, v, replay, checkpoints; integrity_report = integrity(cfg, mcfg, base_config, exp, before_hash)
    summary = {"phase": "4M", "literature_audit": "literature-audit/phase4m_literature_audit.json", "E_math": mcfg["mechanisms"], "anatomical": json.loads((exp / "anatomical-weight/summary.json").read_text()),
               "isolated": isolated, "postsynaptic": post, "anatomy_vs_efficacy": anatomy, "whole_brain": whole, "conditioning": formal, "controls": controls,
               "performance": performance, "candidates": candidates, "integrity": integrity_report, "phase4n_started": False}
    write_json(exp / "audit/final-summary.json", summary); shutil.copy2(args.config, model / "phase4m.yaml"); shutil.copy2(Path(__file__).resolve(), model / "phase4m.py"); shutil.copy2(exp / "literature-audit/phase4m_literature_audit.json", model / "literature-audit.json")
    write_json(model / "VERSION.json", {"name": mcfg["name"], "candidates": candidates, "summary_sha256": sha256(exp / "audit/final-summary.json"), "base_CSR_read_only": True})
    seal(model, model / "MANIFEST.json", model / "checksums.sha256"); seal(exp, exp / "manifests/manifest.json", exp / "manifests/checksums.sha256")
    print(json.dumps({"phase": "4M", "isolated": isolated["status"], "whole_brain": whole["status"], "conditioning": formal["status"], "integrity": integrity_report["all_passed"]}, indent=2)); sys.stdout.flush()


if __name__ == "__main__": main()