Digital Drosophila Archive
529 lines
37 KiB
Python
529 lines
37 KiB
Python
#!/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()
|