#!/usr/bin/env python3 """Phase 4F: official-data retinotopic visual encoding for frozen MaleCNS-LIF.""" from __future__ import annotations import argparse import gc import hashlib import json import math import os import platform import shutil import subprocess import sys import time from collections import Counter from pathlib import Path from typing import Any import numpy as np import pandas as pd import pyarrow as pa import pyarrow.parquet as pq import torch import yaml from scipy.spatial import cKDTree from brain_runtime import BrainRuntime, InputBatch, InputEvent def np_json(value: Any) -> Any: if isinstance(value, (np.generic,)): return value.item() if isinstance(value, np.ndarray): return value.tolist() if isinstance(value, Path): return str(value) raise TypeError(type(value).__name__) def write_json(path: Path, value: Any) -> None: path.parent.mkdir(parents=True, exist_ok=True) path.write_text(json.dumps(value, indent=2, sort_keys=True, default=np_json) + "\n", encoding="utf-8") def sha256(path: Path, chunk: int = 8 << 20) -> str: h = hashlib.sha256() with path.open("rb") as f: while data := f.read(chunk): h.update(data) return h.hexdigest() def verify_checksum_file(root: Path, checksum_file: Path) -> dict[str, Any]: checked, failures = 0, [] for line in checksum_file.read_text(encoding="utf-8").splitlines(): if not line.strip(): continue expected, rel = line.split(None, 1) path = root / rel.strip().lstrip("*") checked += 1 if not path.is_file() or sha256(path) != expected: failures.append(str(path)) return {"checksum_file": str(checksum_file), "checked": checked, "failures": failures, "passed": not failures} def immutable_audit(cfg: dict[str, Any]) -> tuple[list[Path], dict[str, str], dict[str, Any]]: p = cfg["paths"] controls = [ Path(p["reference_root"]) / "checksums.sha256", Path(p["runtime_root"]) / "checksums.sha256", Path(p["io_root"]) / "checksums.sha256", Path(p["phase4c_root"]) / "manifests/checksums.sha256", Path(p["phase4d_root"]) / "manifests/checksums.sha256", Path(p["phase4e_root"]) / "manifests/checksums.sha256", Path(p["model_root"]) / "visual-interface-v2/checksums.sha256", ] missing = [str(x) for x in controls if not x.is_file()] if missing: raise FileNotFoundError(f"immutable controls missing: {missing}") snapshots = {str(x): sha256(x) for x in controls} checks = [] for x in controls: root = x.parent.parent if x.parent.name == "manifests" else x.parent checks.append(verify_checksum_file(root, x)) if not all(x["passed"] for x in checks): raise RuntimeError(f"frozen input checksum failure: {checks}") return controls, snapshots, {"all_passed": True, "controls": checks} def git_commit(path: Path) -> str: return subprocess.check_output(["/usr/bin/git", "-C", str(path), "rev-parse", "HEAD"], text=True).strip() def external_resource_audit(cfg: dict[str, Any], exp: Path) -> dict[str, Any]: root = Path(cfg["paths"]["external_root"]) specs = { "visualpathways": (cfg["official_resources"]["visualpathways_commit"], "GPL-3.0"), "celltype-explorer": (cfg["official_resources"]["celltype_explorer_commit"], "CC-BY-4.0"), "flyvis": (cfg["official_resources"]["flyvis_commit"], "MIT"), } repos = {} for name, (expected, license_id) in specs.items(): path = root / name; actual = git_commit(path) if actual != expected: raise RuntimeError(f"{name} commit {actual} != {expected}") license_files = [x for x in path.iterdir() if x.is_file() and x.name.lower().startswith(("license", "copying"))] repos[name] = {"url": {"visualpathways": "https://github.com/reiserlab/visualpathways", "celltype-explorer": "https://github.com/reiserlab/celltype-explorer-drosophila-male-cns", "flyvis": "https://github.com/TuragaLab/flyvis"}[name], "commit": actual, "expected_commit": expected, "license": license_id, "license_files": [str(x) for x in license_files]} source = root / "retinotopy-tbars-body-summary.parquet" source_summary = json.loads((root / "retinotopy-tbars-body-summary.json").read_text()) result = {"repositories": repos, "official_tbar_snapshot": {**source_summary, "sha256": sha256(source)}, "paper": {"doi": cfg["official_resources"]["paper_doi"], "role": "official visual-pathway definition reference"}} write_json(exp / "audit/external-resource-provenance.json", result) return result def _hierarchy_labels(cfg: dict[str, Any], n: int) -> np.ndarray: labels = np.full(n, None, dtype=object) root = Path(cfg["paths"]["phase4c_root"]) / "visual-audit" for name in ("visual_sensory", "first_order", "second_order", "higher_visual", "central_brain"): path = root / f"hierarchy-{name}.parquet" if path.is_file(): idx = pq.read_table(path, columns=["neuron_idx"])["neuron_idx"].to_numpy().astype(np.int64) labels[idx] = name return labels def _read_eye_map(path: Path, eye: str) -> pd.DataFrame: workbook = pd.ExcelFile(path) candidates = [name for name in workbook.sheet_names if "column" in name.lower()] if not candidates and len(workbook.sheet_names) == 1: candidates = workbook.sheet_names if len(candidates) != 1: raise RuntimeError(f"cannot uniquely identify curated-column sheet in {path}: {workbook.sheet_names}") data = pd.read_excel(path, sheet_name=candidates[0]) id_cols = [x for x in data.columns if x not in {"hex1_id", "hex2_id"}] rows = [] for cell_type in id_cols: for body, h1, h2 in zip(data[cell_type], data["hex1_id"], data["hex2_id"]): if pd.notna(body): rows.append((int(body), float(h1), float(h2), eye, str(cell_type))) out = pd.DataFrame(rows, columns=["body_id", "hex1", "hex2", "eye", "eye_map_cell_type"]) return out.drop_duplicates(["body_id", "hex1", "hex2", "eye", "eye_map_cell_type"]) def build_retinotopy(cfg: dict[str, Any], runtime: BrainRuntime, exp: Path) -> tuple[pd.DataFrame, dict[str, Any]]: ext = Path(cfg["paths"]["external_root"]); vp_root = ext / "visualpathways" md = runtime.table.to_pandas(); n = len(md) body = md["body_id"].to_numpy(np.uint64); body_to_idx = {int(b): i for i, b in enumerate(body)} h1 = md["annotation_assignedOlHex1"].to_numpy(float); h2 = md["annotation_assignedOlHex2"].to_numpy(float) direct = np.isfinite(h1) & np.isfinite(h2) coord1, coord2 = h1.copy(), h2.copy() kind = np.full(n, None, dtype=object); kind[direct] = "official_annotation_assignedOlHex" provenance = np.full(n, None, dtype=object); provenance[direct] = "MaleCNS v1.0 annotation assignedOlHex1/assignedOlHex2" reliability = np.full(n, None, dtype=object); reliability[direct] = "official-direct" contribution = np.full(n, np.nan); rf_r = np.full(n, np.nan); rf_theta = np.full(n, np.nan); rf_vic = np.full(n, np.nan) eye = pd.concat([ _read_eye_map(vp_root / "params/ME_columnar-cells_location.xlsx", "R"), _read_eye_map(vp_root / "params/ME_L_columnar-cells_location.xlsx", "L")], ignore_index=True) eye = eye[eye.body_id.isin(body_to_idx)].copy() eye_consistency = [] for row in eye.itertuples(index=False): i = body_to_idx[int(row.body_id)] if direct[i]: eye_consistency.append(abs(coord1[i] - row.hex1) < 1e-6 and abs(coord2[i] - row.hex2) < 1e-6) else: coord1[i], coord2[i], kind[i], reliability[i] = row.hex1, row.hex2, "official_visualpathways_curated_column", "official-direct" provenance[i] = f"reiserlab/visualpathways curated column {row.eye} {row.eye_map_cell_type}" tbar = pq.read_table(ext / "retinotopy-tbars-body-summary.parquet").to_pandas() tbar_canonical = 0 for row in tbar.itertuples(index=False): i = body_to_idx.get(int(row.body_id)) if i is None: continue tbar_canonical += 1 contribution[i] = row.visual_input_contribution_sum; rf_r[i] = row.tbar_r_weighted_mean rf_theta[i] = row.tbar_theta_weighted_mean; rf_vic[i] = row.tbar_vic_weighted_mean if not np.isfinite(coord1[i]) and np.isfinite(row.tbar_hex1_weighted_mean) and np.isfinite(row.tbar_hex2_weighted_mean): coord1[i], coord2[i] = row.tbar_hex1_weighted_mean, row.tbar_hex2_weighted_mean kind[i], reliability[i] = "official_retinotopy_tbar_weighted_mean", "official-derived-field" provenance[i] = cfg["official_resources"]["retinotopy_tbar_url"] # Anatomical receptive-field centroids for projection neurons are derived # only where no official body coordinate exists. They are explicitly # labelled derived and never counted as direct official coordinate coverage. official = np.isfinite(coord1) & np.isfinite(coord2) src = np.flatnonzero(official) deg = runtime.indptr_np[src + 1] - runtime.indptr_np[src] edge_ids = np.concatenate([np.arange(runtime.indptr_np[s], runtime.indptr_np[s + 1]) for s in src]) edge_src = np.repeat(src, deg); dst = runtime.indices_np[edge_ids]; w = runtime.raw_np[edge_ids].astype(np.float64) sw = np.bincount(dst, weights=w, minlength=n) sx = np.bincount(dst, weights=w * coord1[edge_src], minlength=n) sy = np.bincount(dst, weights=w * coord2[edge_src], minlength=n) superclass = np.asarray(md["annotation_superclass"].tolist(), dtype=object) vp = np.isin(superclass, ["visual_projection", "visual_projection_tbc"]) derived = vp & ~official & (sw > 0) coord1[derived] = sx[derived] / sw[derived]; coord2[derived] = sy[derived] / sw[derived] contribution[derived] = sw[derived]; kind[derived] = "connectome_weighted_centroid_from_official_eye_map" reliability[derived] = "connectome-derived"; provenance[derived] = "official eye-map coordinates + unchanged MaleCNS v1.0 raw synapse-count CSR" mapped = np.isfinite(coord1) & np.isfinite(coord2) xmin, xmax = np.nanmin(coord1[mapped]), np.nanmax(coord1[mapped]); ymin, ymax = np.nanmin(coord2[mapped]), np.nanmax(coord2[mapped]) x = 2 * (coord1 - xmin) / (xmax - xmin) - 1; y = 2 * (coord2 - ymin) / (ymax - ymin) - 1 hierarchy = _hierarchy_labels(cfg, n); hierarchy[vp] = "visual_projection" instance = np.asarray(md["annotation_instance"].tolist(), object) path_class, path_source = np.full(n, None, object), np.full(n, None, object) for filename, field, source_name in (("TableS1_VPN_in_out_summary.csv", "oc_class", "official VPN table"), ("TableS3_VCBN_in_out_summary.csv", "cc_class", "official VCBN table")): table = pd.read_csv(vp_root / "results/tables" / filename) lookup = dict(zip(table["instance"].astype(str), table[field])) for i, val in enumerate(instance): if val is not None and str(val) in lookup and pd.notna(lookup[str(val)]): path_class[i], path_source[i] = str(lookup[str(val)]), source_name eye_by_body = eye.groupby("body_id").agg(eye=("eye", lambda z: ";".join(sorted(set(z)))), eye_map_cell_type=("eye_map_cell_type", lambda z: ";".join(sorted(set(z))))).to_dict("index") rows = pd.DataFrame({ "body_id": body, "neuron_idx": np.arange(n, dtype=np.int32), "visual_layer": hierarchy, "cell_type": md["annotation_type"], "instance": md["annotation_instance"], "side": md["side"], "neurotransmitter": md["nt_consensus_nt"], "official_hex1": coord1, "official_hex2": coord2, "retinotopic_x": x, "retinotopic_y": y, "coordinate_kind": kind, "coordinate_reliability": reliability, "visual_input_contribution": contribution, "receptive_field_r": rf_r, "receptive_field_theta": rf_theta, "receptive_field_vic": rf_vic, "pathway_class": path_class, "pathway_class_source": path_source, "provenance": provenance, "eye": [eye_by_body.get(int(b), {}).get("eye") for b in body], "eye_map_cell_type": [eye_by_body.get(int(b), {}).get("eye_map_cell_type") for b in body], }) rows = rows[mapped].copy().reset_index(drop=True) out = exp / "retinotopy/visual-retinotopy.parquet"; pq.write_table(pa.Table.from_pandas(rows, preserve_index=False), out, compression="zstd") official_mask = rows.coordinate_reliability.isin(["official-direct", "official-derived-field"]) summary = { "canonical_neurons": n, "mapped_total_including_connectome_derived": len(rows), "mapped_total_coverage_percent": 100 * len(rows) / n, "reliable_official_coordinate_neurons": int(official_mask.sum()), "reliable_official_coordinate_coverage_percent": 100 * official_mask.sum() / n, "coordinate_kind_counts": rows.coordinate_kind.value_counts(dropna=False).to_dict(), "official_annotation_both_coordinates": int(direct.sum()), "official_tbar_unique_bodies": int(len(tbar)), "official_tbar_canonical_bodies": tbar_canonical, "visual_projection_total": int(vp.sum()), "visual_projection_mapped": int(np.count_nonzero(vp & mapped)), "connectome_derived_projection": int(derived.sum()), "visualpathways_eye_map_canonical_rows": int(len(eye)), "eye_map_vs_annotation_checked": len(eye_consistency), "eye_map_vs_annotation_exact": int(sum(eye_consistency)), "normalization": {"method": "independent min-max of preserved official-equivalent hex coordinates to [-1,1]", "hex1_min_max": [xmin, xmax], "hex2_min_max": [ymin, ymax], "no_cartesian_claim": True}, "table": str(out), "schema": str(pq.read_schema(out)), } write_json(exp / "retinotopy/summary.json", summary) return rows, summary def hex_grid_for_count(count: int) -> tuple[np.ndarray, np.ndarray, int]: extent = 0 while 1 + 3 * extent * (extent + 1) < count: extent += 1 u, v = [], [] for q in range(-extent, extent + 1): for r in range(max(-extent, -extent - q), min(extent, extent - q) + 1): u.append(q); v.append(r) u, v = np.asarray(u), np.asarray(v) # Exact FlyVis public default convention from flyvis.utils.hex_utils. x, y = 1.5 * v, -np.sqrt(3) * (u + v / 2) x /= max(abs(x).max(), 1); y /= max(abs(y).max(), 1) return x, y, extent def luminance_frames(pattern: str, x: np.ndarray, y: np.ndarray, segments: int) -> np.ndarray: frames = [] for t in range(segments): z = t / max(segments - 1, 1) if pattern == "full-field": f = np.ones_like(x) elif pattern == "left-spot": cx, cy, s = -0.65, 0.0, 0.22; f = np.exp(-((x-cx)**2+(y-cy)**2)/(2*s*s)) elif pattern == "right-spot": cx, cy, s = 0.65, 0.0, 0.22; f = np.exp(-((x-cx)**2+(y-cy)**2)/(2*s*s)) elif pattern == "upper-spot": cx, cy, s = 0.0, 0.65, 0.22; f = np.exp(-((x-cx)**2+(y-cy)**2)/(2*s*s)) elif pattern == "lower-spot": cx, cy, s = 0.0, -0.65, 0.22; f = np.exp(-((x-cx)**2+(y-cy)**2)/(2*s*s)) elif pattern in {"left-to-right", "right-to-left", "up-to-down", "down-to-up"}: p = -0.8 + 1.6*z; p = -p if pattern in {"right-to-left", "down-to-up"} else p cx, cy = (p, 0.0) if pattern in {"left-to-right", "right-to-left"} else (0.0, -p) f = np.exp(-((x-cx)**2+(y-cy)**2)/(2*0.18**2)) elif pattern == "looming": s = 0.08 + 0.48*z; f = np.exp(-(x*x+y*y)/(2*s*s)) elif pattern == "sparse-moving-object": cx, cy = -0.75+1.5*z, 0.55-1.1*z; f = np.exp(-((x-cx)**2+(y-cy)**2)/(2*0.10**2)) elif pattern == "static-object": f = np.exp(-((x+0.25)**2+(y-0.20)**2)/(2*0.10**2)) else: raise KeyError(pattern) frames.append(f.astype(np.float32)) return np.asarray(frames) def feature_channels(frames: np.ndarray, coords: np.ndarray) -> dict[str, np.ndarray]: previous = np.vstack([np.zeros((1, frames.shape[1]), np.float32), frames[:-1]]) delta = frames - previous tree = cKDTree(coords); neighbours = tree.query(coords, k=min(7, len(coords)))[1][:, 1:] dx = np.abs(coords[:, None, 0] - coords[neighbours, 0]); dy = np.abs(coords[:, None, 1] - coords[neighbours, 1]) h, v = [], [] for frame in frames: diff = np.abs(frame[:, None] - frame[neighbours]) h.append((diff * dx).sum(1) / np.maximum(dx.sum(1), 1e-8)); v.append((diff * dy).sum(1) / np.maximum(dy.sum(1), 1e-8)) return {"luminance": frames, "on_change": np.maximum(delta, 0), "off_change": np.maximum(-delta, 0), "temporal_motion": np.abs(delta), "horizontal_motion_tendency": np.asarray(h), "vertical_motion_tendency": np.asarray(v)} def encode_pattern(cfg: dict[str, Any], route: str, pattern: str, pool: pd.DataFrame) -> tuple[InputBatch, np.ndarray, dict[str, Any]]: t0 = time.perf_counter(); enc = cfg["encoder"]; segments = int(enc["temporal_segments"]); k = int(enc["targets_per_segment"]) coords = pool[["retinotopic_x", "retinotopic_y"]].to_numpy(float) if len(pool) < k: raise RuntimeError(f"retinotopic visual-projection pool {len(pool)} < target count {k}") if route == "route-a-malecns-native": source_coords = coords; frames = luminance_frames(pattern, source_coords[:, 0], source_coords[:, 1], segments) channels = feature_channels(frames, source_coords); mapped = channels frontend = "MaleCNS official/derived retinotopic coordinates" elif route == "route-b-flyvis-assisted": official_columns = int(pool.attrs["official_unique_eye_columns"]) fx, fy, extent = hex_grid_for_count(official_columns); fly_coords = np.c_[fx, fy] channels = feature_channels(luminance_frames(pattern, fx, fy, segments), fly_coords) nearest = cKDTree(fly_coords).query(coords, k=1)[1] mapped = {name: value[:, nearest] for name, value in channels.items()} frontend = f"FlyVis axial hex geometry extent={extent}, {len(fx)} hexals; nearest-coordinate functional bridge" else: raise KeyError(route) weights = enc["flyvis_bridge"]["channel_weights"] score = sum(float(weights[name]) * mapped[name] for name in weights) dt = float(enc["trial_duration_ms"]) / segments; events, unions, selections = [], [], [] for i in range(segments): order = np.lexsort((pool.neuron_idx.to_numpy(), -score[i])); chosen = np.sort(pool.neuron_idx.to_numpy(np.int64)[order[:k]]) unions.append(chosen); selections.append(chosen) events.append(InputEvent(target_neuron_idx=chosen, stimulation_mode="external_synaptic_drive", value_kind="poisson_rate_hz", value=float(enc["poisson_rate_hz"]), start_ms=i*dt, duration_ms=dt, source_id=route, channel_id=f"{pattern}:{i}", amplitude_mv=float(enc["external_drive_mv"]), metadata={"route": route, "pattern": pattern, "segment": i, "engineering_feature_channels": list(weights), "connectome_modified": False, "brain_gain_modified": False})) union = np.unique(np.concatenate(unions)) audit = {"route": route, "pattern": pattern, "frontend": frontend, "encoder_wall_seconds": time.perf_counter()-t0, "segments": segments, "targets_per_segment": k, "target_slots": k*segments, "rate_hz": enc["poisson_rate_hz"], "amplitude_mv": enc["external_drive_mv"], "union_target_count": len(union), "feature_channels": {name: {"min": float(val.min()), "max": float(val.max()), "mean": float(val.mean())} for name, val in mapped.items()}} return InputBatch(events=tuple(events), batch_id=f"{route}:{pattern}", metadata={"phase": "4F"}), union, audit def group_members(cfg: dict[str, Any], runtime: BrainRuntime) -> dict[str, np.ndarray]: root = Path(cfg["paths"]["phase4c_root"]) / "visual-audit" optic_parts = [pq.read_table(root/f"hierarchy-{x}.parquet", columns=["neuron_idx"])["neuron_idx"].to_numpy() for x in ("visual_sensory", "first_order", "second_order", "higher_visual")] superc = np.asarray(runtime.superclass, object) return {"optic_lobe": np.unique(np.concatenate(optic_parts)).astype(np.int32), "visual_projection": np.flatnonzero(np.isin(superc, ["visual_projection", "visual_projection_tbc"])).astype(np.int32), "central_brain": pq.read_table(root/"hierarchy-central_brain.parquet", columns=["neuron_idx"])["neuron_idx"].to_numpy().astype(np.int32), "descending": runtime.output_populations["descending"], "motor_related": runtime.output_populations["motor_related"]} def group_counts(nodes: np.ndarray, members: np.ndarray, n: int) -> np.ndarray: lookup = np.full(n, -1, np.int32); lookup[members] = np.arange(len(members), dtype=np.int32) pos = lookup[nodes]; pos = pos[pos >= 0] return np.bincount(pos, minlength=len(members)).astype(np.int32) def peak_fraction(steps: np.ndarray, n: int) -> float: return float(np.bincount(steps.astype(np.int64)).max(initial=0) / n) if len(steps) else 0.0 def run_trial(cfg: dict[str, Any], runtime: BrainRuntime, batch: InputBatch, union: np.ndarray, seed: int, groups: dict[str, np.ndarray], exp: Path, route: str, pattern: str) -> tuple[dict[str, Any], dict[str, np.ndarray]]: runtime.reset(seed); t0 = time.perf_counter(); runtime.step(batch, float(cfg["encoder"]["trial_duration_ms"])); brain_wall = time.perf_counter()-t0 spikes = runtime.read_spikes(); nodes, steps = spikes["neuron_idx"], spikes["step"]; unique = np.unique(nodes) distance = runtime.runtime_module.p4b.directed_distances(runtime.indptr_np, runtime.indices_np, union.astype(np.int64), runtime.n) dvals = distance[unique] if len(unique) else np.empty(0, np.int16); dvals = dvals[dvals >= 0] counts = {name: group_counts(nodes, members, runtime.n) for name, members in groups.items()} masks = {name: np.isin(nodes, members) for name, members in groups.items()} active_fraction = len(unique)/runtime.n; bio_s = runtime.simulation_time_ms/1000.0 mean_rate = len(nodes)/(runtime.n*bio_s); peak = peak_fraction(steps, runtime.n) runaway_cfg = cfg["encoder"]["runaway"] runaway = active_fraction >= float(runaway_cfg["active_fraction_gte"]) or mean_rate >= float(runaway_cfg["mean_firing_rate_hz_gte"]) or peak >= float(runaway_cfg["peak_step_fraction_gte"]) target_mask = np.zeros(runtime.n, bool); target_mask[union] = True record = {"route": route, "pattern": pattern, "seed": seed, "biological_ms": runtime.simulation_time_ms, "brain_wall_seconds": brain_wall, "total_spikes": len(nodes), "active_neurons": len(unique), "active_fraction": active_fraction, "mean_firing_rate_hz": mean_rate, "peak_step_fraction": peak, "runaway": runaway, "external_event_count": int(sum(len(x) for x in runtime.drive_nodes)), "external_total_amplitude_mv": float(sum(float(x.sum()) for x in runtime.drive_values)), "target_union_count": len(union), "target_active": int(np.count_nonzero(target_mask[unique])) if len(unique) else 0, "non_target_active": int(np.count_nonzero(~target_mask[unique])) if len(unique) else 0, "propagation_depth": int(dvals.max(initial=0)), "numerically_finite": bool(torch.isfinite(runtime.v).all() and torch.isfinite(runtime.g).all())} for name in groups: record[f"{name}_active"] = int(np.count_nonzero(counts[name])); record[f"{name}_spikes"] = int(counts[name].sum()) record[f"first_{name}_latency_ms"] = float(spikes["timestamp_ms"][masks[name]].min()) if np.any(masks[name]) else None record["real_anatomical_multi_hop"] = bool(record["non_target_active"] > 0 and record["propagation_depth"] >= int(cfg["success_gate"]["minimum_graph_hops"])) trial_dir = exp / "propagation" / route; trial_dir.mkdir(parents=True, exist_ok=True) np.savez_compressed(trial_dir/f"{pattern}-seed{seed}.npz", step=steps, neuron_idx=nodes, **{f"counts_{k}": v for k,v in counts.items()}) write_json(trial_dir/f"{pattern}-seed{seed}.json", record) return record, counts def metric(a: np.ndarray, b: np.ndarray) -> dict[str, float]: a, b = a.astype(np.float64), b.astype(np.float64); an, bn = np.linalg.norm(a), np.linalg.norm(b) cosine = 1-(float(np.dot(a,b)/(an*bn)) if an and bn else (1.0 if an == bn else 0.0)) if np.array_equal(a,b): corr = 0.0 elif np.std(a) and np.std(b): corr = 1-float(np.corrcoef(a,b)[0,1]) else: corr = 1.0 aa, bb = a>0, b>0; union = np.count_nonzero(aa|bb) return {"cosine_distance": cosine, "correlation_distance": corr, "jaccard_distance": 1-float(np.count_nonzero(aa&bb)/union) if union else 0.0, "euclidean_distance": float(np.linalg.norm(a-b))} def separability(cfg: dict[str, Any], features: dict[str, dict[str, dict[int, dict[str, np.ndarray]]]], exp: Path) -> dict[str, Any]: result = {} for route, patterns in features.items(): route_result = {"matched_pairs": {}, "readouts": {}} for readout in ("visual_projection", "central_brain", "descending", "motor_related"): within = {} for pattern, seeds in patterns.items(): vals = [seeds[s][readout] for s in sorted(seeds)] d = [metric(vals[i], vals[j]) for i in range(len(vals)) for j in range(i+1,len(vals))] within[pattern] = {k: float(np.mean([x[k] for x in d])) for k in d[0]} if d else {k:0.0 for k in metric(vals[0],vals[0])} route_result["readouts"][readout] = {"within_condition_mean": within} for a,b in cfg["encoder"]["matched_pairs"]: key=f"{a}__{b}"; sa, sb=patterns[a],patterns[b] same=[metric(sa[s][readout],sb[s][readout]) for s in sorted(set(sa)&set(sb))] all_between=[metric(x[readout],y[readout]) for x in sa.values() for y in sb.values()] mean_between={k:float(np.mean([z[k] for z in all_between])) for k in all_between[0]} max_within=max(within[a]["euclidean_distance"],within[b]["euclidean_distance"]) repeat_nonzero=all(x["euclidean_distance"]>0 for x in same) repeatable=bool(repeat_nonzero and mean_between["euclidean_distance"]>max_within) route_result["matched_pairs"].setdefault(key,{})[readout]={"same_seed":same,"between_condition_mean":mean_between, "within_condition_euclidean_max":max_within,"repeatably_separable":repeatable} route_result["separable_pairs_in_descending_or_motor"] = [k for k,v in route_result["matched_pairs"].items() if v["descending"]["repeatably_separable"] or v["motor_related"]["repeatably_separable"]] route_result["separable_pair_count_in_descending_or_motor"] = len(route_result["separable_pairs_in_descending_or_motor"]) result[route]=route_result; write_json(exp/f"separability/{route}.json",route_result) write_json(exp/"separability/summary.json",result); return result def regression_suite(cfg: dict[str, Any], runtime: BrainRuntime, exp: Path) -> dict[str, Any]: sys.path.insert(0, cfg["paths"]["io_root"]); import phase4d as p4d reg=p4d.regression(exp,cfg,runtime); api=p4d.api_tests(exp,runtime); stream=p4d.stream_contract(exp,runtime) result={"phase4c_five_trials_passed":all(x["passed"] for x in reg.values()),"phase4d_brainio_passed":api["all_passed"], "activity_stream_passed":stream["replay_exact"],"phase4e_input_contract_checksum_verified":True, "trials":reg,"api":api,"activity_stream":stream} result["all_passed"]=all(result[k] for k in ("phase4c_five_trials_passed","phase4d_brainio_passed","activity_stream_passed","phase4e_input_contract_checksum_verified")) write_json(exp/"regression/summary.json",result); return result def manifest_tree(root: Path, output: Path, exclude: set[Path] | None=None) -> dict[str, Any]: exclude=exclude or set(); files=[] for path in sorted(root.rglob("*")): if path.is_file() and path not in exclude and ".git" not in path.parts: files.append({"path":str(path.relative_to(root)),"bytes":path.stat().st_size,"sha256":sha256(path)}) value={"root":str(root),"file_count":len(files),"total_bytes":sum(x["bytes"] for x in files),"files":files} write_json(output,value); return value def main() -> None: ap=argparse.ArgumentParser(); ap.add_argument("--config",type=Path,required=True); args=ap.parse_args() cfg=yaml.safe_load(args.config.read_text()); exp=Path(cfg["paths"]["experiment_root"]); model=Path(cfg["paths"]["encoder_root"]) for sub in ("external-resources","retinotopy","encoders","stimuli","propagation","separability","latency","regression","audit","manifests","compiler-reports/extension-build"): (exp/sub).mkdir(parents=True,exist_ok=True) shutil.copy2(args.config,exp/"encoders/phase4f.yaml") controls,before,immutable=immutable_audit(cfg); write_json(exp/"audit/immutable-inputs-before.json",immutable) external=external_resource_audit(cfg,exp) sys.path.insert(0,cfg["paths"]["model_root"]); runtime=BrainRuntime(args.config,exp/"compiler-reports/extension-build") retino,retino_summary=build_retinotopy(cfg,runtime,exp) superc=np.asarray(runtime.superclass,object); pool=retino[retino.neuron_idx.isin(np.flatnonzero(np.isin(superc,["visual_projection","visual_projection_tbc"])))].copy() # Prefer direct official receptive-field records, then retain audited connectome-derived centroids. rank={"official-direct":0,"official-derived-field":1,"connectome-derived":2}; pool["_rank"]=pool.coordinate_reliability.map(rank) pool=pool.sort_values(["_rank","neuron_idx"]).drop(columns="_rank").reset_index(drop=True) official_unique=int(retino.loc[retino.coordinate_reliability=="official-direct",["official_hex1","official_hex2"]].drop_duplicates().shape[0]) pool.attrs["official_unique_eye_columns"]=official_unique pq.write_table(pa.Table.from_pandas(pool,preserve_index=False),exp/"retinotopy/visual-projection-entry-pool.parquet",compression="zstd") write_json(exp/"encoders/functional-channel-bridge.json",{"route_a":"features evaluated directly at MaleCNS coordinates", "route_b":"features evaluated on FlyVis public axial-hex geometry then nearest-coordinate bridged; no neuron IDs shared", "flyvis_and_malecns_are_distinct_connectomes":True,"target_expression":cfg["encoder"]["flyvis_bridge"]["target_expression"], "channel_weights":cfg["encoder"]["flyvis_bridge"]["channel_weights"],"biological_claim":"engineering preprocessing only; not phototransduction"}) groups=group_members(cfg,runtime); features={}; trial_records=[]; encoder_audits=[] for route in ("route-a-malecns-native","route-b-flyvis-assisted"): features[route]={} for pattern in cfg["encoder"]["patterns"]: batch,union,enc_audit=encode_pattern(cfg,route,pattern,pool); encoder_audits.append(enc_audit) write_json(exp/f"stimuli/{route}-{pattern}.json",enc_audit); features[route][pattern]={} for seed in map(int,cfg["encoder"]["seeds"]): rec,counts=run_trial(cfg,runtime,batch,union,seed,groups,exp,route,pattern) rec["encoder_wall_seconds"]=enc_audit["encoder_wall_seconds"]; rec["total_control_wall_seconds"]=rec["brain_wall_seconds"]+rec["encoder_wall_seconds"] trial_records.append(rec); features[route][pattern][seed]=counts # Exact repeatability with the same seed, input, and initial state. deterministic={} for route in features: batch,union,_=encode_pattern(cfg,route,"left-spot",pool); seed=int(cfg["encoder"]["seeds"][0]) runtime.reset(seed); runtime.step(batch,float(cfg["encoder"]["trial_duration_ms"])); a=runtime.read_spikes() runtime.reset(seed); runtime.step(batch,float(cfg["encoder"]["trial_duration_ms"])); b=runtime.read_spikes() deterministic[route]=bool(np.array_equal(a["step"],b["step"]) and np.array_equal(a["neuron_idx"],b["neuron_idx"])) sep=separability(cfg,features,exp) matched={} for route in features: matched[route]={} for seed in map(int,cfg["encoder"]["seeds"]): vals=[x["external_event_count"] for x in trial_records if x["route"]==route and x["seed"]==seed] amps=[x["external_total_amplitude_mv"] for x in trial_records if x["route"]==route and x["seed"]==seed] matched[route][str(seed)]={"event_counts":sorted(set(vals)),"amplitudes_mv":sorted(set(amps)),"matched":len(set(vals))==1 and len(set(amps))==1} write_json(exp/"audit/matched-input-energy.json",matched); write_json(exp/"audit/determinism.json",deterministic) latency={"trials":trial_records,"encoder_wall_seconds":{"min":min(x["encoder_wall_seconds"] for x in trial_records),"max":max(x["encoder_wall_seconds"] for x in trial_records)}, "brain_wall_seconds":{"min":min(x["brain_wall_seconds"] for x in trial_records),"max":max(x["brain_wall_seconds"] for x in trial_records)}, "total_control_wall_seconds":{"min":min(x["total_control_wall_seconds"] for x in trial_records),"max":max(x["total_control_wall_seconds"] for x in trial_records)}} write_json(exp/"latency/summary.json",latency); write_json(exp/"propagation/summary.json",trial_records) del runtime; gc.collect(); torch.cuda.empty_cache() reg_runtime=BrainRuntime(args.config,exp/"compiler-reports/extension-build"); regression=regression_suite(cfg,reg_runtime,exp) del reg_runtime; gc.collect(); torch.cuda.empty_cache() after={str(x):sha256(x) for x in controls}; unchanged=before==after all_matched=all(v["matched"] for r in matched.values() for v in r.values()); all_deterministic=all(deterministic.values()) non_runaway=all(not x["runaway"] and x["numerically_finite"] for x in trial_records) multi_hop=any(x["real_anatomical_multi_hop"] for x in trial_records) output_latencies=[z for x in trial_records for z in (x["first_descending_latency_ms"],x["first_motor_related_latency_ms"]) if z is not None] within_100=bool(output_latencies and min(output_latencies)<=float(cfg["success_gate"]["maximum_output_latency_ms"])) separable_pairs={r:v["separable_pair_count_in_descending_or_motor"] for r,v in sep.items()} separability_gate=any(v>=int(cfg["success_gate"]["minimum_separable_matched_pairs_in_descending_or_motor"]) for v in separable_pairs.values()) passed=all([unchanged,regression["all_passed"],all_matched,all_deterministic,non_runaway,multi_hop,within_100,separability_gate]) status="PASS-CANDIDATE" if passed else "FAIL-NO-V1-CANDIDATE" summary={"status":status,"candidate_generated":passed,"retinotopy":retino_summary,"entry_pool_count":len(pool), "routes":{"route-a-malecns-native":{"trial_count":sum(x["route"]=="route-a-malecns-native" for x in trial_records),"separability":sep["route-a-malecns-native"]}, "route-b-flyvis-assisted":{"trial_count":sum(x["route"]=="route-b-flyvis-assisted" for x in trial_records),"separability":sep["route-b-flyvis-assisted"]}}, "matched_input_energy":all_matched,"same_seed_exact_reproducibility":deterministic,"real_anatomical_multi_hop":multi_hop, "non_runaway_and_finite":non_runaway,"minimum_descending_or_motor_latency_ms":min(output_latencies) if output_latencies else None, "separable_pair_counts":separable_pairs,"separability_gate":separability_gate,"regression":regression, "frozen_inputs_unchanged":unchanged,"external_resources":external,"game_rl_decoder_web_3d_started":False} write_json(exp/"audit/summary.json",summary) model.mkdir(parents=True,exist_ok=True) for src in (Path(__file__).resolve(),args.config.resolve(),Path(__file__).with_name("phase4f_extract_retinotopy.py")): if src.is_file() and src.resolve() != (model/src.name).resolve(): shutil.copy2(src,model/src.name) shutil.copy2(Path(cfg["paths"]["io_root"])/"brain_runtime.py",model/"brain_runtime.py") write_json(model/"VERSION.json",{"name":cfg["names"]["encoder"],"status":status,"candidate_name":cfg["names"]["candidate"] if passed else None, "brain_frozen":True,"connectome_changed":False,"lif_semantics_changed":False,"transmitter_sign_changed":False}) external_manifest=manifest_tree(Path(cfg["paths"]["external_root"]),exp/"manifests/external-resources-manifest.json") model_manifest=manifest_tree(model,model/"MANIFEST.json",{model/"MANIFEST.json",model/"checksums.sha256"}) model_files=sorted(p for p in model.rglob("*") if p.is_file() and p.name not in {"MANIFEST.json","checksums.sha256"})+[model/"MANIFEST.json"] (model/"checksums.sha256").write_text("".join(f"{sha256(p)} {p.relative_to(model)}\n" for p in model_files),encoding="utf-8") exp_manifest=manifest_tree(exp,exp/"manifests/manifest.json",{exp/"manifests/manifest.json",exp/"manifests/checksums.sha256"}) outputs=sorted(p for p in exp.rglob("*") if p.is_file() and p != exp/"manifests/checksums.sha256") (exp/"manifests/checksums.sha256").write_text("".join(f"{sha256(p)} {p.relative_to(exp)}\n" for p in outputs),encoding="utf-8") print(json.dumps({"status":status,"candidate_generated":passed,"mapped":retino_summary["mapped_total_including_connectome_derived"], "official_mapped":retino_summary["reliable_official_coordinate_neurons"],"separable_pairs":separable_pairs, "min_output_latency_ms":summary["minimum_descending_or_motor_latency_ms"]},indent=2)) if __name__=="__main__": main()