Seal Phase4AE same-precision GPU numerical gate and audit limits

This commit is contained in:
Digital Drosophila Archive
2026-09-30 07:50:57 +08:00
parent 92241aea02
commit 1eb4c52cf7
648 changed files with 51363 additions and 22 deletions
@@ -0,0 +1,68 @@
# Phase 4AE — GPU author-semantics numerical gate
User authorized an independent model branch. This is not conditioning and does
not replace any Phase 1–4AD candidate/artifact. Learning remains unachieved.
The fixed author equations and all physical parameters come from the validated
4AD oracle, with dt=.1ms, delay=18 steps, refractory=22 steps and W=.275mV.
Canonical CSR, raw contact counts and frozen transmitter signs are read-only.
## Predeclared tests
Five synthetic numerical fixtures: ordered cancellation (including deliberately
nonphysiological large signed integers), 256-source fan-in bursts, refractory
release/self-loop collisions, representable threshold neighbours, and a fixed
seed signed random graph. Synthetic negative node labels are NOT MaleCNS IDs;
duplicate/zero edges occur only in labelled numerical fixtures.
Two 128-neuron real induced graphs are selected before simulation: the highest
in-degree target and its 127 largest-contact presynaptic partners; separately the
highest-in-degree target with >=64 positive and >=63 negative presynaptic inputs,
64 positive and 63 negative largest-contact partners. Ties break by canonical
neuron_idx. Include every existing edge within the selected node set, without
changing raw counts/signs or adding edges. This is numerical stress, not an input
interface or a physiological validation. Preserve identities and selection audit.
Also run the same five 4AD motifs for 100ms with one common seed request (or no
request for no_input). Report every outcome; do not select by activity/separability.
These tests do not exhaust all full-brain numerical corner cases.
## Exact arithmetic and schedule
Execute unchanged pinned author create_model in Brian2 2.5.1/NumPy1.22.3 using
the isolated 4AD environment. Export raw SI variables, actual coefficients and
weights. GPU float32 coefficients round to float32 before each multiply/add;
float64 uses double. No FMA or fast-math. Match generated expression order.
Update refractory and integrate; apply only currently permitted voltage seed
requests; strict threshold closes write access immediately; then accept delayed
anatomical arrivals in author order within each target, external synthetic input,
and reset. A fused reset is equivalent only if all arrivals reject current spikes.
Never replace sequential g+=w with g+sum(w). A queue proxy observes real Brian
arrival order without changing values; compare proxy and unwrapped runs exactly.
Use an exact spike ring buffer, no approximated delay or dense edge scan.
Before-threshold requests use the existing common requested-event convention,
not the author's distinct poi interface that changes target refractory to zero.
Report requested/accepted/rejected input separately. No graded APL in these isolated
tests; its future integration and frozen coefficients require explicit audit.
For each case/precision compare GPU to the same-precision real Brian oracle:
exact events, counts, refractory, input delivery and per-target arrival order/guard;
report full-history max/mean V/g error and bitwise state comparison. Fixed tolerance
is 1e-11V float64 / 1e-6V float32. Any event divergence fails, regardless of tolerance.
Each GPU case repeats three times with all arrays exactly equal. Cross-precision
differences are separately descriptive, not grounds to change thresholds.
## Integrity and progression
Preregister sources, old input hashes and deterministic fixture hashes before any
model execution. Preserve failures; use a separately recorded new attempt for
implementation repairs. Do not change frozen artifacts. Save manifests/SHA-256.
Old five deterministic-v2 full-brain reference regressions must pass again.
Only then consider a short full-graph author comparison, followed by DoOR CAL;
never claim full-brain equivalence or learning from isolated test PASS.
Author source/license/provenance are inherited read-only from 4AD. Supporting
primary documentation: https://numpy.org/doc/1.22/reference/generated/numpy.ufunc.at.html
and https://brian2.readthedocs.io/en/2.5.1/user/refractoriness.html . The fixed
generated source and actual boundary fixtures, not simplified prose, define timing.
@@ -0,0 +1,18 @@
#include <torch/extension.h>
void author_state_cuda(torch::Tensor v, torch::Tensor g, torch::Tensor ref,
torch::Tensor requests, torch::Tensor spikes, torch::Tensor counts,
double gd, double c0, double cg, double cv, double threshold,
double reset, double seed_volt, int64_t refractory_steps);
void author_fill_edges_cuda(torch::Tensor indptr, torch::Tensor sources,
torch::Tensor offsets, torch::Tensor output);
void author_arrivals_cuda(torch::Tensor edges, torch::Tensor starts, torch::Tensor targets,
torch::Tensor weights, torch::Tensor g, torch::Tensor ref);
void author_external_cuda(torch::Tensor external, torch::Tensor g, torch::Tensor ref);
PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
m.def("state_step", &author_state_cuda, "Author NumPy-ordered integration/threshold/reset");
m.def("fill_edges", &author_fill_edges_cuda, "Expand only active-source canonical CSR edges");
m.def("apply_arrivals", &author_arrivals_cuda, "Sequential per-target arrival additions");
m.def("apply_external", &author_external_cuda, "Post-anatomical external arrival additions");
}
@@ -0,0 +1,213 @@
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAException.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <cmath>
#include <cstdint>
#include <limits>
namespace {
constexpr int BLOCK = 128;
void check_tensor(const torch::Tensor& x, const torch::Tensor& like,
at::ScalarType dtype, const char* name) {
TORCH_CHECK(x.is_cuda() && x.is_contiguous() && x.dim() == 1,
name, " must be a contiguous CUDA vector");
TORCH_CHECK(x.device() == like.device(), name, " is on a different CUDA device");
TORCH_CHECK(x.scalar_type() == dtype, name, " has the wrong dtype");
}
void check_state(const torch::Tensor& g, const torch::Tensor& ref) {
TORCH_CHECK(g.scalar_type() == torch::kFloat32 || g.scalar_type() == torch::kFloat64,
"state must be float32 or float64");
check_tensor(g, g, g.scalar_type(), "g");
check_tensor(ref, g, torch::kInt32, "ref");
TORCH_CHECK(g.numel() == ref.numel(), "g/ref length mismatch");
TORCH_CHECK(g.numel() <= std::numeric_limits<int>::max(), "too many neurons");
}
void index_range(const torch::Tensor& x, int64_t upper, const char* name) {
if (x.numel()) {
TORCH_CHECK(x.min().item<int64_t>() >= 0 && x.max().item<int64_t>() < upper,
name, " index out of bounds");
}
}
template <typename T> __device__ __forceinline__ T add_rn(T a, T b);
template <> __device__ __forceinline__ float add_rn(float a, float b) { return __fadd_rn(a, b); }
template <> __device__ __forceinline__ double add_rn(double a, double b) { return __dadd_rn(a, b); }
template <typename T> __device__ __forceinline__ T mul_rn(T a, T b);
template <> __device__ __forceinline__ float mul_rn(float a, float b) { return __fmul_rn(a, b); }
template <> __device__ __forceinline__ double mul_rn(double a, double b) { return __dmul_rn(a, b); }
template <typename T>
__global__ void state_kernel(T* v, T* g, int32_t* ref, const uint8_t* requests,
uint8_t* spikes, int32_t* counts, int n, T gd, T c0,
T cg, T cv, T threshold, T reset, T seed_volt, int ref_steps) {
const int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i >= n) return;
int r = ref[i];
r = r > 0 ? r - 1 : 0;
const bool available = r == 0;
T vi = v[i];
if (available) {
const T oldg = g[i];
const T oldv = vi;
// Match NumPy float-array expression evaluation: each operation rounds in T.
g[i] = mul_rn(gd, oldg);
const T driven = mul_rn(cg, oldg);
const T leaked = mul_rn(cv, oldv);
vi = add_rn(c0, add_rn(driven, leaked));
if (requests[i]) vi = seed_volt;
}
const bool fired = available && vi > threshold;
if (fired) {
vi = reset;
g[i] = T(0);
r = ref_steps;
counts[i] += 1;
}
// Nonavailable neurons retain V/g. Resetting fired cells before arrivals is
// equivalent to reset-after-arrivals because all their arrivals are rejected.
v[i] = vi;
ref[i] = r;
spikes[i] = fired ? 1 : 0;
}
__global__ void fill_kernel(const int64_t* indptr, const int64_t* sources,
const int64_t* offsets, int64_t* output, int n) {
const int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i >= n) return;
const int64_t source = sources[i];
const int64_t begin = indptr[source];
const int64_t end = indptr[source + 1];
for (int64_t j = begin; j < end; ++j) output[offsets[i] + j - begin] = j;
}
template <typename T>
__global__ void arrival_kernel(const int64_t* edges, const int64_t* starts,
const int32_t* targets, const T* weights, T* g,
const int32_t* ref, int segments, int64_t edge_count) {
const int segment = blockIdx.x * blockDim.x + threadIdx.x;
if (segment >= segments) return;
const int target = targets[segment];
if (ref[target] != 0) return;
const int64_t end = segment + 1 < segments ? starts[segment + 1] : edge_count;
// Start at existing g, never sum weights separately: cancellation can differ.
T accumulator = g[target];
for (int64_t j = starts[segment]; j < end; ++j)
accumulator = add_rn(accumulator, weights[edges[j]]);
g[target] = accumulator;
}
template <typename T>
__global__ void external_kernel(const T* external, T* g, const int32_t* ref, int n) {
const int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n && ref[i] == 0 && external[i] != T(0))
g[i] = add_rn(g[i], external[i]);
}
} // namespace
void author_state_cuda(torch::Tensor v, torch::Tensor g, torch::Tensor ref,
torch::Tensor requests, torch::Tensor spikes, torch::Tensor counts,
double gd, double c0, double cg, double cv, double threshold,
double reset, double seed_volt, int64_t refractory_steps) {
check_state(g, ref);
const c10::cuda::CUDAGuard guard(g.device());
check_tensor(v, g, g.scalar_type(), "v");
check_tensor(requests, g, torch::kUInt8, "requests");
check_tensor(spikes, g, torch::kUInt8, "spikes");
check_tensor(counts, g, torch::kInt32, "counts");
const int64_t n = g.numel();
TORCH_CHECK(v.numel() == n && requests.numel() == n && spikes.numel() == n && counts.numel() == n,
"state vector length mismatch");
TORCH_CHECK(refractory_steps > 0 && refractory_steps <= std::numeric_limits<int>::max(),
"invalid refractory length");
for (double x : {gd, c0, cg, cv, threshold, reset, seed_volt})
TORCH_CHECK(std::isfinite(x), "nonfinite coefficient");
if (!n) return;
const auto stream = at::cuda::getCurrentCUDAStream(g.get_device());
AT_DISPATCH_FLOATING_TYPES(g.scalar_type(), "author_state", [&] {
state_kernel<scalar_t><<<(n + BLOCK - 1) / BLOCK, BLOCK, 0, stream>>>(
v.data_ptr<scalar_t>(), g.data_ptr<scalar_t>(), ref.data_ptr<int32_t>(),
requests.data_ptr<uint8_t>(), spikes.data_ptr<uint8_t>(), counts.data_ptr<int32_t>(),
static_cast<int>(n), static_cast<scalar_t>(gd), static_cast<scalar_t>(c0),
static_cast<scalar_t>(cg), static_cast<scalar_t>(cv), static_cast<scalar_t>(threshold),
static_cast<scalar_t>(reset), static_cast<scalar_t>(seed_volt), static_cast<int>(refractory_steps));
});
C10_CUDA_KERNEL_LAUNCH_CHECK();
}
void author_fill_edges_cuda(torch::Tensor indptr, torch::Tensor sources,
torch::Tensor offsets, torch::Tensor output) {
check_tensor(indptr, indptr, torch::kInt64, "indptr");
const c10::cuda::CUDAGuard guard(indptr.device());
check_tensor(sources, indptr, torch::kInt64, "sources");
check_tensor(offsets, indptr, torch::kInt64, "offsets");
check_tensor(output, indptr, torch::kInt64, "output");
TORCH_CHECK(indptr.numel() >= 1 && sources.numel() == offsets.numel(), "CSR shape mismatch");
TORCH_CHECK(sources.numel() <= std::numeric_limits<int>::max(), "too many sources");
TORCH_CHECK(indptr[0].item<int64_t>() == 0, "CSR indptr must begin at zero");
if (indptr.numel() > 1)
TORCH_CHECK((indptr.slice(0, 1) >= indptr.slice(0, 0, -1)).all().item<bool>(), "CSR indptr not monotonic");
index_range(sources, indptr.numel() - 1, "source");
if (!sources.numel()) {
TORCH_CHECK(output.numel() == 0, "nonempty expansion for no active sources");
return;
}
const auto degrees = indptr.index_select(0, sources + 1) - indptr.index_select(0, sources);
const auto expected_offsets = degrees.cumsum(0) - degrees;
TORCH_CHECK(torch::equal(offsets, expected_offsets), "incorrect output offsets");
TORCH_CHECK(degrees.sum().item<int64_t>() == output.numel(), "incorrect expansion length");
const auto stream = at::cuda::getCurrentCUDAStream(indptr.get_device());
fill_kernel<<<(sources.numel() + BLOCK - 1) / BLOCK, BLOCK, 0, stream>>>(
indptr.data_ptr<int64_t>(), sources.data_ptr<int64_t>(), offsets.data_ptr<int64_t>(),
output.data_ptr<int64_t>(), static_cast<int>(sources.numel()));
C10_CUDA_KERNEL_LAUNCH_CHECK();
}
void author_arrivals_cuda(torch::Tensor edges, torch::Tensor starts, torch::Tensor targets,
torch::Tensor weights, torch::Tensor g, torch::Tensor ref) {
check_state(g, ref);
const c10::cuda::CUDAGuard guard(g.device());
check_tensor(edges, g, torch::kInt64, "edges");
check_tensor(starts, g, torch::kInt64, "starts");
check_tensor(targets, g, torch::kInt32, "targets");
check_tensor(weights, g, g.scalar_type(), "weights");
TORCH_CHECK(starts.numel() == targets.numel(), "segment count mismatch");
TORCH_CHECK(starts.numel() <= std::numeric_limits<int>::max(), "too many target segments");
if (!starts.numel()) {
TORCH_CHECK(!edges.numel(), "edges without target segments");
return;
}
index_range(edges, weights.numel(), "edge");
index_range(starts, edges.numel(), "segment start");
index_range(targets, g.numel(), "target");
TORCH_CHECK(starts[0].item<int64_t>() == 0, "first segment must start at zero");
if (starts.numel() > 1) {
TORCH_CHECK((starts.slice(0, 1) > starts.slice(0, 0, -1)).all().item<bool>(), "segment starts not increasing");
TORCH_CHECK((targets.slice(0, 1) > targets.slice(0, 0, -1)).all().item<bool>(), "targets must be unique and increasing");
}
const auto stream = at::cuda::getCurrentCUDAStream(g.get_device());
AT_DISPATCH_FLOATING_TYPES(g.scalar_type(), "author_arrivals", [&] {
arrival_kernel<scalar_t><<<(starts.numel() + BLOCK - 1) / BLOCK, BLOCK, 0, stream>>>(
edges.data_ptr<int64_t>(), starts.data_ptr<int64_t>(), targets.data_ptr<int32_t>(),
weights.data_ptr<scalar_t>(), g.data_ptr<scalar_t>(), ref.data_ptr<int32_t>(),
static_cast<int>(starts.numel()), edges.numel());
});
C10_CUDA_KERNEL_LAUNCH_CHECK();
}
void author_external_cuda(torch::Tensor external, torch::Tensor g, torch::Tensor ref) {
check_state(g, ref);
const c10::cuda::CUDAGuard guard(g.device());
check_tensor(external, g, g.scalar_type(), "external");
TORCH_CHECK(external.numel() == g.numel(), "external/state length mismatch");
if (!g.numel()) return;
const auto stream = at::cuda::getCurrentCUDAStream(g.get_device());
AT_DISPATCH_FLOATING_TYPES(g.scalar_type(), "author_external", [&] {
external_kernel<scalar_t><<<(g.numel() + BLOCK - 1) / BLOCK, BLOCK, 0, stream>>>(
external.data_ptr<scalar_t>(), g.data_ptr<scalar_t>(), ref.data_ptr<int32_t>(), static_cast<int>(g.numel()));
});
C10_CUDA_KERNEL_LAUNCH_CHECK();
}
@@ -0,0 +1,12 @@
72601f46dde3202d62e71161958f620a3ce5b8f9c24cf1f5b7158289b5f9a3f0 PROTOCOL.md
208a937d17bbb041b9bf6b25b89264ee566830b60609e58491b44583d8541f15 author_ext.cpp
669f166e18ce9a3f63fc723db8b0a8b74f7fa4757194b1bd65b4cf6008a7787e author_ext.cu
356b6161576f6f6d79745af1927bb43b12a7c3d50922bf9e1bd29e473fcf2106 config.json
9eea9a592241cab7f4303850e969def7887eee3f4ec18421103605cafd258cf8 context.py
cb75d512ddf9a5ff344f3ac82eefb112406593bf7b414dfac79d7c5d730371ee finalize.py
03758d07572f4b44d272059d9308fac534268ad201b57bdd77d69fee95e67278 gpu_backend.py
32eb1c42383068ce52701d2139a376ce155ad7072cde019949f8ca2d55492409 oracle_stress.py
2519eaa2c50b2d91fd0c4f623aeabd6f5d212c5b537951d0a66bc81c569db78c prepare.py
92733a48921f95c3c0935bf16a704634b89b6540352f31cc1169ec6d3211b57a run_suite.py
720d04b379b6f1f64cae28a2203444a6d4c4e72ecb121cc69e93c2740fae69ab verify_frozen.py
7e8b0ac4adcd48eceabca7c9cc4f2c42a8746ebc673aee4490bae5a4746bbe12 manifest.json
@@ -0,0 +1,19 @@
{
"phase": "4AE",
"name": "author-runtime-lab-v1",
"human_authorization": "允许独立新模型分支测试动力学假设",
"seed": 410301,
"gpu_uuid": "GPU-79f3a85d-adce-7833-7003-6488140ca596",
"precisions": ["float64", "float32"],
"gpu_repeats": 3,
"max_abs_volt": {"float64": 1e-11, "float32": 1e-6},
"exact_events_required": true,
"exact_refractory_required": true,
"exact_arrival_order_per_tick_target_required": true,
"training": false,
"whole_brain_author_backend_claimed": false,
"candidate_generation": false,
"synthetic_cases": ["ordered_cancellation", "fanin_bursts", "release_selfloop", "threshold_rounding", "signed_random"],
"real_cases": ["highest_indegree", "mixed_high_indegree"],
"prior_motifs": ["pair", "trio", "five", "feedforward", "no_input"]
}
@@ -0,0 +1,40 @@
import hashlib
import json
from pathlib import Path
MODEL=Path(__file__).resolve().parent
ROOT=Path('/data/digital-drosophila')
EXP=ROOT/'experiments/phase4ae-author-runtime'
AD=ROOT/'experiments/phase4ad-author-semantics'
CFG=json.loads((MODEL/'config.json').read_text())
def digest(path):
h=hashlib.sha256()
with Path(path).open('rb') as f:
for b in iter(lambda:f.read(8<<20),b''):h.update(b)
return h.hexdigest()
def save(path,value):
path=Path(path);path.parent.mkdir(parents=True,exist_ok=True)
with path.open('x',encoding='utf-8') as f:json.dump(value,f,indent=2,ensure_ascii=False,allow_nan=False);f.write('\n')
def sources():
return {p.name:digest(p) for p in sorted(MODEL.iterdir()) if p.is_file() and p.name not in ('manifest.json','checksums.sha256')}
def check():
p=json.loads((EXP/'protocol/preregistered.json').read_text())
assert sources()==p['sources'],'Source changed after preregistration'
for row in p['inputs']:assert digest(ROOT/row['path'])==row['sha256'],row['path']
for row in json.loads((EXP/'protocol/case-hashes.json').read_text()):assert digest(EXP/row['path'])==row['sha256'],row['path']
return p
def integrity():
import sys
sys.path.insert(0,str(ROOT/'models/male-cns-lif/lif-semantics-audit-v1'))
import audit_semantics as old
r=old.integrity()
for rel in ('models/male-cns-lif/lif-semantics-audit-v1','experiments/phase4ac-lif-semantics',
'registry/feedback-context-review-20260930','models/male-cns-lif/author-semantics-lab-v1','experiments/phase4ad-author-semantics'):
r[rel]=old.manifest_check(ROOT/rel,digest)
assert all(x['pass'] for x in r.values())
return r
@@ -0,0 +1,52 @@
"""Gate the diagnostic prototype, never promote it as learned/full-brain runtime."""
import json
from pathlib import Path
import sys
from context import ROOT,EXP,MODEL,CFG,save,check,integrity,digest
from run_suite import load,compare
def seal(root):
rows=[]
for p in sorted(root.rglob('*')):
if '__pycache__' in p.parts or '.venv' in p.parts:continue
if p.is_file():rows.append(dict(path=str(p.relative_to(root)),bytes=p.stat().st_size,sha256=digest(p)))
assert not (root/'manifest.json').exists() and not (root/'checksums.sha256').exists()
save(root/'manifest.json',dict(root=str(root),files=rows))
with (root/'checksums.sha256').open('x') as f:
for r in rows:f.write(r['sha256']+' '+r['path']+'\n')
f.write(digest(root/'manifest.json')+' manifest.json\n')
def main():
check();after=integrity();save(EXP/'regression/old-integrity-after.json',after)
summary=json.loads((EXP/'regression/gpu-complete.json').read_text());ref=json.loads((EXP/'regression/frozen-five.json').read_text())
expected=(len(CFG['synthetic_cases'])+len(CFG['real_cases'])+len(CFG['prior_motifs']))*len(CFG['precisions'])*CFG['gpu_repeats']
assert len(summary['trials'])==expected
rows=[]
for record in summary['trials']:
name=record['case'];precision=record['precision'];repeat=record['repeat'];case=EXP/'cases'/name
oracle=load(EXP/'oracle'/(name+'-'+precision)/'state.npz');gpu=load(EXP/'gpu'/(name+'-'+precision+f'-r{repeat}')/'state.npz')
actual=compare(oracle,gpu,case,precision);assert actual==record['comparison'],'Stored report differs from raw states'
first=load(EXP/'gpu'/(name+'-'+precision+'-r0')/'state.npz')
import numpy as np
equal={k:bool(np.array_equal(v,gpu[k])) for k,v in first.items()};assert equal==record['repeat_arrays']
rows.append(record)
passed=all(r['pass_all'] for r in rows) and ref['pass_all'] and all(x['pass'] for x in after.values())
status='PASS-ISOLATED-GPU-NUMERICAL-GATE' if passed else 'FAIL-ISOLATED-GPU-NUMERICAL-GATE'
save(EXP/'audit/final-status.json',dict(status=status,pass_all=passed,cases=expected//6,gpu_runs=expected,
full_brain_author_equivalence=False,conditioning=False,learning_achieved=False,candidate_generated=False,
next_gate='Short full-graph same-precision author comparison, then fixed DoOR CAL only if numerical gates pass'))
lines=['# Phase 4AE GPU作者语义数值门禁','',f'状态:`{status}`。','',
'仅验证固定小网络,尚未证明全脑等价、生物学正确或学习可用。未训练、未修改冻结输入或旧模型。','',
'| case | precision | events | max V/g error V | exact whole history | pass |','|---|---|---:|---:|---|---|']
for r in rows:
if r['repeat']:continue
c=r['comparison'];lines.append(f"| {r['case']} | {r['precision']} | {c['oracle_events']} | {c['max_abs_volt']} | {c['bitwise_history']} | {r['pass_all']} |")
lines+=['','每项GPU执行3次,检查原始事件、输入请求实际送达、逐目标到达顺序及guard、计数、不应期和完整状态。',
'大整数、随机图、重复/零边属于纯数值夹具,不是新增MaleCNS连接或生理假设。真实子图使用全部既有诱导边、原contact count与冻结transmitter sign。',
'原始冻结完整性与旧五组全脑deterministic-v2回归再次核验;该旧回归PASS不等于新动力学全脑PASS。',
'没有根据活动强度选择参数;尚须短时全图对照及嗅觉信息保持测试。','']
with (EXP/'REPORT.md').open('x',encoding='utf-8') as f:f.write('\n'.join(lines))
check();seal(MODEL);seal(EXP);print(status,flush=True)
if not passed:raise SystemExit(1)
if __name__=='__main__':main()
@@ -0,0 +1,275 @@
"""Small-case diagnostic CUDA backend; SI states, explicit dtype rounding.
run_gpu is serial, records CPU diagnostics, and is not a performance runtime.
Coefficients and optional weights must come from the same-precision Brian2 oracle.
"""
from functools import lru_cache
import hashlib
import json
import math
import os
from pathlib import Path
import time
import traceback
import numpy as np
import torch
from torch.utils.cpp_extension import load
SOURCE_ROOT = Path(__file__).resolve().parent
COEFFICIENT_NAMES = ("gd", "c0", "cg", "cv", "threshold", "reset", "seed_volt")
def require(condition, message):
if not condition:
raise ValueError(message)
def write_json(path, data):
with Path(path).open("x", encoding="utf-8") as output:
json.dump(data, output, indent=2, allow_nan=False)
output.write("\n")
def digest(path):
checksum = hashlib.sha256()
with Path(path).open("rb") as source:
for block in iter(lambda: source.read(1 << 20), b""):
checksum.update(block)
return checksum.hexdigest()
@lru_cache(maxsize=None)
def _load_extension(build_path, signature):
build = Path(build_path) / ("author_cuda_" + signature[:16])
build.mkdir(parents=True, exist_ok=True)
module = load(name="author_cuda_" + signature[:16],
sources=[str(SOURCE_ROOT / "author_ext.cpp"), str(SOURCE_ROOT / "author_ext.cu")],
build_directory=str(build), extra_cflags=["-O2"],
extra_cuda_cflags=["-O2", "--fmad=false", "--ftz=false", "--prec-div=true", "--prec-sqrt=true"],
verbose=False)
return module
def build_extension(build_root):
provenance = {"sources": {name: digest(SOURCE_ROOT / name)
for name in ("author_ext.cpp", "author_ext.cu")},
"torch": torch.__version__, "torch_cuda": torch.version.cuda,
"arch_list": os.environ.get("TORCH_CUDA_ARCH_LIST"),
"cuda_flags": ["-O2", "--fmad=false", "--ftz=false", "--prec-div=true", "--prec-sqrt=true"]}
signature = hashlib.sha256(json.dumps(provenance, sort_keys=True).encode()).hexdigest()
return _load_extension(str(Path(build_root).resolve()), signature), provenance
def read_case(case_dir, precision, coefficients, weights):
case_dir = Path(case_dir)
require(precision in ("float32", "float64"), "precision must be float32 or float64")
dtype = np.dtype(precision)
with (case_dir / "spec.json").open(encoding="utf-8") as source:
spec = json.load(source)
for key in ("dt_ms", "steps", "duration_ms", "model"):
require(key in spec, "Missing spec field: " + key)
model = spec["model"]
for key in ("threshold_mv", "reset_potential_mv", "resting_potential_mv",
"synaptic_delay_ms", "refractory_ms", "w_syn_mv_per_anatomical_synapse"):
require(key in model and math.isfinite(float(model[key])), "Missing/nonfinite model field: " + key)
steps = spec["steps"]
require(type(steps) is int and 0 < steps < np.iinfo(np.int32).max, "Invalid step count")
dt = float(spec["dt_ms"])
require(math.isfinite(dt) and dt > 0, "Invalid dt")
require(math.isfinite(float(spec["duration_ms"])) and
round(float(spec["duration_ms"]) / dt) == steps, "Duration/steps mismatch")
delay = round(float(model["synaptic_delay_ms"]) / dt)
refractory = round(float(model["refractory_ms"]) / dt)
require(delay == 18 and refractory == 22, "This declared backend requires delay=18 and refractory=22")
require(float(model["w_syn_mv_per_anatomical_synapse"]) == 0.275, "W_syn differs from frozen value")
coeff = {}
for key in COEFFICIENT_NAMES:
require(key in coefficients and math.isfinite(float(coefficients[key])), "Missing/nonfinite coefficient: " + key)
coeff[key] = float(coefficients[key])
require(bool(np.isfinite(dtype.type(coeff[key]))), "Coefficient overflows requested dtype: " + key)
require(dtype.type(coeff["threshold"]) == dtype.type(float(model["threshold_mv"]) * .001),
"Oracle threshold differs from case spec")
require(dtype.type(coeff["reset"]) == dtype.type(float(model["reset_potential_mv"]) * .001),
"Oracle reset differs from case spec")
with np.load(case_dir / "arrays.npz", allow_pickle=False) as inputs:
arrays = {key: inputs[key].copy() for key in inputs.files}
required_types = {"nodes": "int32", "edge_id": "int64", "src": "int32", "dst": "int32",
"signed_raw": "int64", "requests": "uint8", "external_mv": "float64",
"initial_v_mv": "float64", "initial_g_mv": "float64"}
for key, expected in required_types.items():
require(key in arrays and arrays[key].dtype == np.dtype(expected), "Missing/wrong array dtype: " + key)
nodes = arrays["nodes"]
require(nodes.ndim == 1 and len(np.unique(nodes)) == len(nodes), "Invalid neuron identities")
n = len(nodes)
require(n < 512, "This diagnostic backend is restricted to fewer than 512 cells")
edge_id = arrays["edge_id"]
require(edge_id.ndim == 1 and len(np.unique(edge_id)) == len(edge_id), "Invalid edge identities")
e = len(edge_id)
for key in ("src", "dst", "signed_raw"):
require(arrays[key].shape == (e,), "Edge vector shape mismatch: " + key)
for key in ("src", "dst"):
require(bool(np.all((arrays[key] >= 0) & (arrays[key] < n))), "Local edge endpoint out of bounds")
require(np.array_equal(np.lexsort((edge_id, arrays["src"])), np.arange(e)),
"Edges must already be canonical source/edge_id sorted; indices will not be silently reordered")
for key in ("requests", "external_mv"):
require(arrays[key].shape == (steps, n), "Trial array shape mismatch: " + key)
require(bool(np.all(arrays["requests"] <= 1)), "Request array is not Boolean uint8")
for key in ("initial_v_mv", "initial_g_mv"):
require(arrays[key].shape == (n,), "Initial state shape mismatch: " + key)
for key in ("external_mv", "initial_v_mv", "initial_g_mv"):
require(bool(np.isfinite(arrays[key]).all()), "Nonfinite case array: " + key)
if weights is None:
weight_array = np.asarray(arrays["signed_raw"] * (.275 * .001), dtype=dtype)
weight_origin = "signed_raw * (0.275 * 0.001), then cast to requested dtype"
else:
weight_array = np.asarray(weights, dtype=dtype)
weight_origin = "caller supplied oracle SI weight array"
require(weight_array.shape == (e,) and bool(np.isfinite(weight_array).all()), "Invalid SI weights")
for key in ("initial_v_mv", "initial_g_mv", "external_mv"):
converted = np.asarray(arrays[key] * .001, dtype=dtype)
require(bool(np.isfinite(converted).all()), "Conversion overflow: " + key)
return spec, arrays, coeff, np.ascontiguousarray(weight_array), weight_origin
def run_gpu(case_dir, out_dir, precision, coefficients, build_root, weights=None):
"""Run one case; return state dict and save SI state plus GPU-selected arrivals.
``weights`` is optional SI volts in the original arrays.npz edge-row order.
Repeated calls share a cached extension; each needs a new exclusive out_dir.
"""
case_dir, out_dir = Path(case_dir), Path(out_dir)
require(not out_dir.exists(), "Refusing existing output directory: " + str(out_dir))
spec, arrays, coeff, weight_array, weight_origin = read_case(case_dir, precision, coefficients, weights)
require(torch.cuda.is_available(), "CUDA unavailable")
require(torch.cuda.device_count() == 1 and "4090" in torch.cuda.get_device_name(), "Only the single declared RTX4090 is permitted")
out_dir.mkdir(parents=True, exist_ok=False)
start = time.perf_counter()
try:
ext, build_provenance = build_extension(build_root)
build_done = time.perf_counter()
device = torch.device("cuda", torch.cuda.current_device())
torch.cuda.reset_peak_memory_stats(device)
dtype = torch.float32 if precision == "float32" else torch.float64
numpy_dtype = np.dtype(precision)
n, steps = len(arrays["nodes"]), spec["steps"]
def upload(value, target_dtype=None):
return torch.as_tensor(np.ascontiguousarray(value), device=device, dtype=target_dtype).contiguous()
v = upload(np.asarray(arrays["initial_v_mv"] * .001, dtype=numpy_dtype), dtype)
g = upload(np.asarray(arrays["initial_g_mv"] * .001, dtype=numpy_dtype), dtype)
ref = torch.zeros(n, dtype=torch.int32, device=device)
counts = torch.zeros(n, dtype=torch.int32, device=device)
requests = upload(arrays["requests"], torch.uint8)
accepted_requests = torch.empty_like(requests)
rejected_requests = torch.empty_like(requests)
external = upload(np.asarray(arrays["external_mv"] * .001, dtype=numpy_dtype), dtype)
synaptic_weights = upload(weight_array, dtype)
dst = upload(arrays["dst"], torch.int32)
degrees_np = np.bincount(arrays["src"], minlength=n).astype(np.int64)
indptr_np = np.concatenate((np.zeros(1, np.int64), np.cumsum(degrees_np, dtype=np.int64)))
indptr = upload(indptr_np, torch.int64)
ring = torch.zeros((19, n), dtype=torch.uint8, device=device)
history = torch.empty((steps, 3, n), dtype=dtype, device=device)
event_steps, event_nodes = [], []
arrival_steps, arrival_indices, arrival_accepted = [], [], []
torch.cuda.synchronize(device)
simulation_start = time.perf_counter()
for step in range(steps):
spikes = ring[step % 19]
# The state kernel decrements refractory before considering a request.
# These diagnostics use that exact pre-kernel eligibility on the GPU.
available_after_decrement = ref <= 1
accepted_requests[step].copy_(requests[step] & available_after_decrement.to(torch.uint8))
rejected_requests[step].copy_(requests[step] & (~available_after_decrement).to(torch.uint8))
ext.state_step(v, g, ref, requests[step], spikes, counts,
*[coeff[key] for key in COEFFICIENT_NAMES], 22)
if step >= 18:
due = ring[(step - 18) % 19]
sources = torch.nonzero(due, as_tuple=False).flatten().contiguous()
if sources.numel():
degrees = indptr[sources + 1] - indptr[sources]
offsets = (torch.cumsum(degrees, dim=0) - degrees).contiguous()
edge_count = int(degrees.sum().item())
expanded = torch.empty(edge_count, device=device, dtype=torch.int64)
ext.fill_edges(indptr, sources, offsets, expanded)
if edge_count:
destination = dst[expanded]
order = torch.argsort(destination, stable=True)
edge_indices = expanded[order].contiguous()
ordered_targets = destination[order]
boundary = torch.ones(edge_count, dtype=torch.bool, device=device)
boundary[1:] = ordered_targets[1:] != ordered_targets[:-1]
starts = torch.nonzero(boundary, as_tuple=False).flatten().contiguous()
targets = ordered_targets[starts].contiguous()
accepted = (ref[ordered_targets.long()] == 0).to(torch.uint8)
# Capture the actual CUDA-expanded and GPU-stable-sorted edge stream.
arrival_steps.append(np.full(edge_count, step, dtype=np.int32))
arrival_indices.append(edge_indices.cpu().numpy().copy())
arrival_accepted.append(accepted.cpu().numpy().copy())
ext.apply_arrivals(edge_indices, starts, targets, synaptic_weights, g, ref)
ext.apply_external(external[step], g, ref)
history[step, 0].copy_(v)
history[step, 1].copy_(g)
history[step, 2].copy_(ref)
fired = torch.nonzero(spikes, as_tuple=False).flatten().cpu().numpy()
if len(fired):
event_steps.append(np.full(len(fired), step, dtype=np.int32))
event_nodes.append(arrays["nodes"][fired])
torch.cuda.synchronize(device)
simulation_done = time.perf_counter()
def joined(chunks, dt):
return np.concatenate(chunks).astype(dt, copy=False) if chunks else np.empty(0, dtype=dt)
result = {"nodes": arrays["nodes"], "history": history.cpu().numpy(),
"counts": counts.cpu().numpy(), "v": v.cpu().numpy(), "g": g.cpu().numpy(),
"ref": ref.cpu().numpy(), "weights": synaptic_weights.cpu().numpy(),
"step": joined(event_steps, np.int32), "neuron_idx": joined(event_nodes, np.int32),
"request_attempted": arrays["requests"].copy(),
"request_accepted": accepted_requests.cpu().numpy(),
"request_rejected": rejected_requests.cpu().numpy()}
arrival_index = joined(arrival_indices, np.int64)
arrivals = {"step": joined(arrival_steps, np.int32), "edge_index": arrival_index,
"edge_id": arrays["edge_id"][arrival_index],
"target_local": arrays["dst"][arrival_index],
"accepted": joined(arrival_accepted, np.uint8)}
result.update(arrival_step=arrivals["step"], arrival_edge_index=arrivals["edge_index"],
arrival_accepted=arrivals["accepted"])
with (out_dir / "state.npz").open("xb") as target:
np.savez_compressed(target, **result)
with (out_dir / "arrival-order.npz").open("xb") as target:
np.savez_compressed(target, **arrivals)
finite = bool(all(np.isfinite(value).all() for value in result.values()))
metrics = {"precision": precision, "state_units": "SI volts; refractory is steps",
"coefficients_double_input": coeff,
"coefficients_effective_dtype": {key: float(numpy_dtype.type(value)) for key, value in coeff.items()},
"weights_origin": weight_origin, "nodes": n, "edges": len(weight_array), "steps": steps,
"total_spikes": len(result["step"]), "arriving_edges": len(arrival_index),
"accepted_arriving_edges": int(np.count_nonzero(arrivals["accepted"])),
"requested": int(np.count_nonzero(result["request_attempted"])),
"accepted_requests": int(np.count_nonzero(result["request_accepted"])),
"rejected_requests": int(np.count_nonzero(result["request_rejected"])),
"all_finite": finite, "delay_steps": 18, "refractory_steps": 22,
"build_or_cache_seconds": build_done - start,
"diagnostic_simulation_seconds": simulation_done - simulation_start,
"wall_seconds": time.perf_counter() - start,
"gpu": torch.cuda.get_device_name(device), "device": str(device),
"peak_allocated_vram_bytes": int(torch.cuda.max_memory_allocated(device)),
"peak_reserved_vram_bytes": int(torch.cuda.max_memory_reserved(device)),
"build": build_provenance,
"input_sha256": {name: digest(case_dir / name) for name in ("spec.json", "arrays.npz")},
"arrival_order_scope": "GPU expanded edges, stable target grouping; compare relative order per tick and target",
"performance_runtime_claim": False, "oracle_validation_claim": False,
"whole_brain_run": False, "learning": False,
"cpu_transfer_each_step_for_diagnostics": True}
write_json(out_dir / "metrics.json", metrics)
require(finite, "Nonfinite GPU state; raw trace preserved")
return result
except Exception as error:
write_json(out_dir / "error.json", {"type": type(error).__name__, "message": str(error),
"traceback": traceback.format_exc(), "precision": precision,
"validation_pass": False})
raise
@@ -0,0 +1,60 @@
{
"root": "/data/digital-drosophila/models/male-cns-lif/author-runtime-lab-v1",
"files": [
{
"path": "PROTOCOL.md",
"bytes": 4188,
"sha256": "72601f46dde3202d62e71161958f620a3ce5b8f9c24cf1f5b7158289b5f9a3f0"
},
{
"path": "author_ext.cpp",
"bytes": 1203,
"sha256": "208a937d17bbb041b9bf6b25b89264ee566830b60609e58491b44583d8541f15"
},
{
"path": "author_ext.cu",
"bytes": 10422,
"sha256": "669f166e18ce9a3f63fc723db8b0a8b74f7fa4757194b1bd65b4cf6008a7787e"
},
{
"path": "config.json",
"bytes": 805,
"sha256": "356b6161576f6f6d79745af1927bb43b12a7c3d50922bf9e1bd29e473fcf2106"
},
{
"path": "context.py",
"bytes": 1707,
"sha256": "9eea9a592241cab7f4303850e969def7887eee3f4ec18421103605cafd258cf8"
},
{
"path": "finalize.py",
"bytes": 3814,
"sha256": "cb75d512ddf9a5ff344f3ac82eefb112406593bf7b414dfac79d7c5d730371ee"
},
{
"path": "gpu_backend.py",
"bytes": 16277,
"sha256": "03758d07572f4b44d272059d9308fac534268ad201b57bdd77d69fee95e67278"
},
{
"path": "oracle_stress.py",
"bytes": 29235,
"sha256": "32eb1c42383068ce52701d2139a376ce155ad7072cde019949f8ca2d55492409"
},
{
"path": "prepare.py",
"bytes": 8828,
"sha256": "2519eaa2c50b2d91fd0c4f623aeabd6f5d212c5b537951d0a66bc81c569db78c"
},
{
"path": "run_suite.py",
"bytes": 5318,
"sha256": "92733a48921f95c3c0935bf16a704634b89b6540352f31cc1169ec6d3211b57a"
},
{
"path": "verify_frozen.py",
"bytes": 2157,
"sha256": "720d04b379b6f1f64cae28a2203444a6d4c4e72ecb121cc69e93c2740fae69ab"
}
]
}
@@ -0,0 +1,502 @@
"""Small, fixed-input Brian2 2.5.1 stress oracle for Phase 4AE.
The archived author source, equations, reset, parameters, and synapses are left
unchanged. Diagnostic initial states and the old common request interface are
explicit inputs. This file does not install packages or run whole-brain cases.
"""
from __future__ import annotations
import hashlib
import importlib.metadata
import json
import platform
import sys
import types
from pathlib import Path
import numpy as np
import pandas as pd
AUTHOR_PATH = Path("/data/digital-drosophila/experiments/phase4ad-author-semantics/external-resources/model.py")
AUTHOR_SHA256 = "fc45837d7122c6ce2a7f3f2f23c515992e4b232aadb919efabb72337fac88e4e"
AUTHOR_COMMIT = "91bdd1e7dcf193f3e7ca5a8933497fcef63b7960"
LIO_NAMES = {"gd": "_lio_2", "c0": "_lio_8", "cg": "_lio_9", "cv": "_lio_10"}
def _save(path: Path, value) -> None:
with path.open("x", encoding="utf-8") as handle:
json.dump(value, handle, indent=2, sort_keys=True, allow_nan=False)
handle.write("\n")
def _save_arrays(path: Path, arrays: dict) -> None:
with path.open("xb") as handle:
np.savez_compressed(handle, **arrays)
def _digest(path: Path) -> str:
h = hashlib.sha256()
with path.open("rb") as handle:
for block in iter(lambda: handle.read(8 << 20), b""):
h.update(block)
return h.hexdigest()
def _check(condition, message: str) -> None:
if not condition:
raise ValueError(message)
def _load_case(case_dir: Path) -> tuple[dict, dict]:
spec = json.loads((case_dir / "spec.json").read_text(encoding="utf-8"))
with np.load(case_dir / "arrays.npz", allow_pickle=False) as archive:
arrays = {name: archive[name].copy() for name in archive.files}
expected_types = {"nodes": "int32", "edge_id": "int64", "src": "int32",
"dst": "int32", "signed_raw": "int64", "requests": "uint8",
"external_mv": "float64", "initial_v_mv": "float64",
"initial_g_mv": "float64"}
for name, dtype in expected_types.items():
_check(name in arrays, f"Missing input array: {name}")
_check(arrays[name].dtype == np.dtype(dtype), f"{name} must have dtype {dtype}")
nodes = arrays["nodes"]
_check(nodes.ndim == 1 and 1 <= len(nodes) <= 512, "Stress cases require 1..512 nodes")
_check(np.all(np.diff(nodes.astype(np.int64)) > 0), "nodes must be sorted and unique")
steps = int(spec["steps"])
_check(steps == spec["steps"] and 1 <= steps <= 1000, "steps must be an integer in 1..1000")
_check(float(spec["dt_ms"]) == .1, "The fixed stress protocol requires dt_ms=0.1")
if "duration_ms" in spec:
_check(np.isclose(float(spec["duration_ms"]), steps * .1, rtol=0, atol=1e-12),
"duration_ms disagrees with steps and dt_ms")
if "dt_ms" in spec["model"]:
_check(float(spec["model"]["dt_ms"]) == .1, "model dt_ms disagrees with case dt_ms")
edge_count = len(arrays["edge_id"])
# The unchanged author calls Synapses.connect(i=[], j=[]), which Brian2
# 2.5.1 does not support. Empty-network fixtures must explicitly provide
# a documented zero-weight synthetic edge in the input case.
_check(edge_count > 0, "Use an explicit zero-weight fixture edge instead of empty connectivity")
for name in ("edge_id", "src", "dst", "signed_raw"):
_check(arrays[name].shape == (edge_count,), f"{name} has an inconsistent shape")
_check(len(np.unique(arrays["edge_id"])) == edge_count, "edge_id must be unique")
for name in ("src", "dst"):
_check(np.all((arrays[name] >= 0) & (arrays[name] < len(nodes))), f"{name} is out of local range")
_check(np.array_equal(np.lexsort((arrays["edge_id"], arrays["src"])), np.arange(edge_count)),
"Edges must be in canonical source-then-edge order")
for name in ("requests", "external_mv"):
_check(arrays[name].shape == (steps, len(nodes)), f"{name} must have shape (steps, N)")
_check(np.all(arrays["requests"] <= 1), "requests must contain only 0 or 1")
for name in ("initial_v_mv", "initial_g_mv"):
_check(arrays[name].shape == (len(nodes),), f"{name} must have shape (N,)")
for name in ("external_mv", "initial_v_mv", "initial_g_mv"):
_check(np.isfinite(arrays[name]).all(), f"{name} contains a nonfinite value")
return spec, arrays
def _code_records(network) -> list:
records = []
for obj in network.sorted_objects:
row = {"name": obj.name, "type": type(obj).__name__,
"when": obj.when, "order": int(obj.order), "active": bool(obj.active)}
if hasattr(obj, "abstract_code"):
abstract = obj.abstract_code
row["abstract_code"] = abstract if isinstance(abstract, (str, dict)) else str(abstract)
row["code_objects"] = [
{"name": code.name, "generated_code": str(code.code)}
for code in getattr(obj, "code_objects", ())]
records.append(row)
return records
class _QueueReadProxy:
"""Forward the queue unchanged; copy only observations from peek()."""
def __init__(self, queue, get_tick, get_available, targets, observations, counters):
self._original_queue = queue
self._get_tick = get_tick
self._get_available = get_available
self._targets = targets
self._observations = observations
self._counters = counters
def peek(self, *args, **kwargs):
value = self._original_queue.peek(*args, **kwargs)
edge_indices = np.asarray(value, dtype=np.int64).copy()
if np.any(edge_indices < 0) or np.any(edge_indices >= len(self._targets)):
raise AssertionError("The anatomical queue returned an invalid local edge index")
tick = self._get_tick()
accepted = np.asarray(self._get_available(), dtype=bool)[self._targets[edge_indices]].astype(np.uint8)
self._observations.append((tick, edge_indices, accepted))
self._counters["peek_calls"] += 1
return value
def advance(self, *args, **kwargs):
self._counters["advance_calls"] += 1
return self._original_queue.advance(*args, **kwargs)
def __getattr__(self, name):
return getattr(self._original_queue, name)
def _scalar_namespace(code_runner, name: str):
for code in getattr(code_runner, "code_objects", ()):
if name in code.namespace:
value = np.asarray(code.namespace[name])
if value.ndim == 0 and np.isfinite(value):
return float(value), {"source": "CodeObject.namespace", "code_object": code.name,
"name": name, "namespace_dtype": str(value.dtype)}
return None, None
def _coefficients(neu, parameters, b2, dt_seconds: float) -> dict:
generated = "\n".join(str(code.code) for code in neu.state_updater.code_objects)
_check("_lio_2 * g" in generated and
"_lio_8 + ((_lio_9 * g) + (_lio_10 * v))" in generated,
"Generated state expression changed; the recorded _lio-to-GPU mapping is not valid")
dt = np.float64(dt_seconds)
tau = np.float64(parameters["tau"] / b2.second)
membrane = np.float64(parameters["t_mbr"] / b2.second)
rest = np.float64(parameters["v_0"] / b2.volt)
# This is the actual scalar expression order emitted by the pinned Brian
# state updater, not an algebraically rearranged analytic solution.
lio_1 = (-dt) / tau
lio_3 = tau / (membrane - tau)
lio_4 = dt / membrane
lio_5 = dt / tau
lio_6 = (-dt) / membrane
fallback = {
"gd": np.exp(lio_1), "c0": rest - (rest * np.exp(lio_6)),
"cg": ((lio_3 * ((-np.exp(lio_4)) + np.exp(lio_5))) * np.exp(lio_6)) * np.exp(lio_1),
"cv": np.exp(lio_6),
}
result, provenance = {}, {}
for key, name in LIO_NAMES.items():
value, record = _scalar_namespace(neu.state_updater, name)
if value is None:
value = float(fallback[key])
record = {"source": "same generated scalar expression evaluated with numpy.float64", "name": name}
result[key], provenance[key] = value, record
for key, runner, name, author_key in (
("threshold", neu.thresholder["spike"], "v_th", "v_th"),
("reset", neu.resetter["spike"], "v_rst", "v_rst"),
):
value, record = _scalar_namespace(runner, name)
if value is None:
value = float(parameters[author_key] / b2.volt)
record = {"source": "unchanged author parameter", "name": name}
result[key], provenance[key] = value, record
result["seed_volt"] = float((parameters["v_th"] + 1 * b2.mV) / b2.volt)
provenance["seed_volt"] = {"source": "common request interface: unchanged v_th + 1 mV"}
result["metadata"] = {
"units": "SI volts and dimensionless coefficients", "names": LIO_NAMES,
"provenance": provenance, "hex": {key: value.hex() for key, value in result.items()},
"voltage_expression": "c0 + ((cg * old_g) + (cv * old_v))",
"g_expression": "gd * old_g",
"constant_values_are_double": True,
"float32_execution_note": "Brian NumPy scalar/array promotion and operation order must be matched separately; double coefficients alone do not assert float64 intermediate array arithmetic.",
}
return result
def _parameters_match(parameters, model, b2) -> dict:
rows = {}
for config_key, author_key, unit in (
("membrane_tau_ms", "t_mbr", b2.ms), ("synaptic_tau_ms", "tau", b2.ms),
("refractory_ms", "t_rfc", b2.ms), ("synaptic_delay_ms", "t_dly", b2.ms),
("resting_potential_mv", "v_0", b2.mV), ("reset_potential_mv", "v_rst", b2.mV),
("threshold_mv", "v_th", b2.mV), ("w_syn_mv_per_anatomical_synapse", "w_syn", b2.mV),
):
requested, actual = float(model[config_key]), float(parameters[author_key] / unit)
_check(np.isfinite(requested) and np.isclose(requested, actual, rtol=0, atol=1e-12),
f"Model parameter {config_key} differs from the archived author")
rows[config_key] = {"case": requested, "author": actual}
return rows
def _run_once(spec, arrays, precision, output, source_bytes, trace_queue):
import brian2 as b2
from brian2.core.base import BrianObject
output.mkdir(parents=True, exist_ok=False)
saved = {key: b2.prefs[key] for key in (
"codegen.target", "core.default_float_dtype", "legacy.refractory_timing")}
old_dt = b2.defaultclock.dt
dtype = np.dtype(precision)
steps, nodes = int(spec["steps"]), arrays["nodes"]
dt_ms, dt_seconds = float(spec["dt_ms"]), float(spec["dt_ms"]) * .001
try:
b2.start_scope()
b2.prefs["codegen.target"] = "numpy"
b2.prefs["core.default_float_dtype"] = dtype.type
b2.prefs["legacy.refractory_timing"] = False
b2.defaultclock.dt = dt_ms * b2.ms
author = types.ModuleType("phase4ae_author_" + precision)
author.__file__ = str(AUTHOR_PATH)
exec(compile(source_bytes, str(AUTHOR_PATH), "exec"), author.__dict__)
parameters = dict(author.default_params)
parameter_checks = _parameters_match(parameters, spec["model"], b2)
completeness = pd.DataFrame({"local_index": np.arange(len(nodes), dtype=np.int32)},
index=pd.Index(nodes, name="global_neuron_idx"))
connectivity = pd.DataFrame({
"Presynaptic_Index": arrays["src"], "Postsynaptic_Index": arrays["dst"],
"Excitatory x Connectivity": arrays["signed_raw"], "edge_id": arrays["edge_id"],
"source_neuron_idx": nodes[arrays["src"]], "target_neuron_idx": nodes[arrays["dst"]],
})
comp_path, con_path = output / "completeness.csv", output / "connectivity.parquet"
completeness.to_csv(comp_path)
connectivity.to_parquet(con_path, engine="pyarrow", index=False)
neu, syn, spikes = author.create_model(comp_path, con_path, parameters)
_check(np.dtype(neu.variables["v"].dtype) == dtype and np.dtype(neu.variables["g"].dtype) == dtype,
"Brian state dtype disagrees with the precision arm")
_check(np.dtype(syn.variables["w"].dtype) == dtype, "Brian synaptic weight dtype disagrees")
np.testing.assert_array_equal(np.asarray(syn.i[:]), arrays["src"])
np.testing.assert_array_equal(np.asarray(syn.j[:]), arrays["dst"])
neu.v = arrays["initial_v_mv"] * b2.mV
neu.g = arrays["initial_g_mv"] * b2.mV
initial_v = np.asarray(neu.v_[:], dtype=dtype).copy()
initial_g = np.asarray(neu.g_[:], dtype=dtype).copy()
attempted = arrays["requests"].copy()
accepted = np.zeros_like(attempted)
rejected = np.zeros_like(attempted)
request_slots_seen = np.zeros(steps, dtype=np.uint8)
def current_tick():
return int(round(float(neu.clock.t / b2.second) / dt_seconds))
def requests_before_thresholds():
tick = current_tick()
if tick < 0 or tick >= steps or request_slots_seen[tick]:
raise AssertionError("Common request operation executed outside its unique tick")
request_slots_seen[tick] = 1
request = attempted[tick].astype(bool)
available = np.asarray(neu.not_refractory[:], dtype=bool)
permitted = request & available
accepted[tick] = permitted
rejected[tick] = request & ~available
if permitted.any():
neu.v[np.flatnonzero(permitted)] = parameters["v_th"] + 1 * b2.mV
request_op = b2.NetworkOperation(requests_before_thresholds, clock=neu.clock,
when="before_thresholds", order=0, name="oracle_common_requests")
objects = [neu, syn, spikes, request_op]
ext_ticks, ext_targets = np.nonzero(arrays["external_mv"])
external_weights = np.empty(0, dtype=dtype)
if len(ext_ticks):
ext_ids = np.arange(len(ext_ticks), dtype=np.int32)
ext_source = b2.SpikeGeneratorGroup(len(ext_ticks), ext_ids,
ext_ticks.astype(np.float64) * dt_ms * b2.ms, sorted=True,
clock=neu.clock, name="zz_external_sources")
ext_syn = b2.Synapses(ext_source, neu, "w : volt", on_pre="g_post += w",
delay=0 * b2.ms, name="zz_external_arrivals")
ext_syn.connect(i=ext_ids, j=ext_targets)
ext_syn.w = arrays["external_mv"][ext_ticks, ext_targets] * b2.mV
ext_syn.pre.order = 0
_check(syn.pre.when == ext_syn.pre.when == "synapses" and syn.pre.order < 0,
"External delivery must follow the original anatomical pathway")
external_weights = np.asarray(ext_syn.w_[:], dtype=dtype).copy()
objects.extend([ext_source, ext_syn])
observations = []
queue_counters = {"peek_calls": 0, "advance_calls": 0, "installations": []}
if trace_queue:
class QueueInstaller(BrianObject):
def __init__(self):
super().__init__(clock=neu.clock, when="end", order=1_000_000,
name="zz_oracle_queue_installer")
def before_run(self, run_namespace):
# Network.before_run traverses sorted objects. All original
# code runners have compiled before this last end-slot object.
installed = []
for code in syn.pre.code_objects:
if "_queue" not in code.namespace or "_queue.peek()" not in str(code.code):
continue
original = code.namespace["_queue"]
if isinstance(original, _QueueReadProxy):
original = original._original_queue
code.namespace["_queue"] = _QueueReadProxy(
original, current_tick, lambda: neu.not_refractory[:],
arrays["dst"], observations, queue_counters)
installed.append(code.name)
if len(installed) != 1:
raise AssertionError(f"Expected exactly one anatomical queue consumer, got {installed}")
queue_counters["installations"].append(installed[0])
objects.append(QueueInstaller())
trace = b2.StateMonitor(neu, ["v", "g", "lastspike", "not_refractory"],
record=True, when="end", order=0, name="oracle_end_state")
objects.append(trace)
network = b2.Network(*objects)
_save(output / "schedule.json", {"schedule": list(network.schedule),
"summary": str(network.scheduling_summary())})
_save(output / "input-audit.json", {
"source_sha256": AUTHOR_SHA256, "parameter_checks": parameter_checks,
"source_parameters_modified": False, "precision": precision,
"dt_ms": dt_ms, "steps": steps, "request_interface": "before_thresholds; v=v_th+1mV only if not_refractory",
"poisson_helper_used": False, "all_neurons_keep_original_rfc": True,
"queue_trace_enabled": trace_queue, "state_units": "SI volts, original dtype",
"completeness_sha256": _digest(comp_path), "connectivity_sha256": _digest(con_path),
})
network.run(0 * b2.ms, namespace={})
_save(output / "compiled-code.json", _code_records(network))
if request_slots_seen.any() or observations:
raise AssertionError("Zero-duration compilation unexpectedly advanced the case")
network.run(steps * dt_ms * b2.ms, namespace={})
_save(output / "generated-code.json", _code_records(network))
_check(np.all(request_slots_seen == 1), "Not every request tick was observed exactly once")
event_float = np.asarray(spikes.t[:] / b2.second, dtype=np.float64) / dt_seconds
event_ticks = np.rint(event_float).astype(np.int32)
_check(np.allclose(event_float, event_ticks, rtol=0, atol=1e-6), "Spike time is off the integer clock")
local_indices = np.asarray(spikes.i[:], dtype=np.int32).copy()
order = np.lexsort((nodes[local_indices], event_ticks))
event_ticks, local_indices = event_ticks[order], local_indices[order]
fired = np.zeros_like(attempted)
fired[event_ticks, local_indices] = 1
counts = np.asarray(spikes.count[:], dtype=np.int32).copy()
_check(np.array_equal(fired.sum(axis=0), counts), "Events and per-neuron counts disagree")
_check(not np.any(accepted & ~fired), "An applicable common request did not produce a real spike")
_check(np.array_equal(attempted, accepted + rejected), "Request accounting is incomplete")
state_ticks = np.asarray(trace.t[:] / b2.second, dtype=np.float64) / dt_seconds
_check(state_ticks.shape == (steps,) and np.allclose(state_ticks, np.arange(steps), rtol=0, atol=1e-6),
"StateMonitor does not cover the exact requested ticks")
voltage = np.asarray(trace.v_, dtype=dtype).T.copy()
conductance = np.asarray(trace.g_, dtype=dtype).T.copy()
last = np.asarray(trace.lastspike_, dtype=np.float64).T
seen = np.isfinite(last) & (last >= 0)
last_ticks = np.zeros(last.shape, dtype=np.int64)
last_ticks[seen] = np.rint(last[seen] / dt_seconds).astype(np.int64)
refractory_steps = int(round(float(parameters["t_rfc"] / b2.ms) / dt_ms))
remaining = np.where(seen, np.maximum(refractory_steps - (np.arange(steps)[:, None] - last_ticks), 0), 0).astype(np.int32)
available_end = np.asarray(trace.not_refractory, dtype=bool).T
_check(np.array_equal(remaining == 0, available_end), "Refractory flag and integer countdown disagree")
history = np.stack([voltage, conductance, remaining.astype(dtype)], axis=1)
_check(np.isfinite(history).all(), "Author state became nonfinite")
arrival_step = np.empty(0, dtype=np.int32)
arrival_edge = np.empty(0, dtype=np.int64)
arrival_accepted = np.empty(0, dtype=np.uint8)
if trace_queue:
_check(queue_counters["peek_calls"] == steps and queue_counters["advance_calls"] == steps,
"Queue proxy did not observe exactly one peek and advance per tick")
_check([entry[0] for entry in observations] == list(range(steps)), "Queue peek ticks are not complete and ordered")
arrival_step = np.concatenate([np.full(len(ids), tick, dtype=np.int32) for tick, ids, _ in observations])
arrival_edge = np.concatenate([ids for _, ids, _ in observations]).astype(np.int64, copy=False)
arrival_accepted = np.concatenate([mask for _, _, mask in observations]).astype(np.uint8, copy=False)
delay_steps = int(round(float(parameters["t_dly"] / b2.ms) / dt_ms))
if len(arrival_step):
_check(np.all(arrival_step >= delay_steps), "An anatomical event arrived before the fixed delay")
_check(np.all(fired[arrival_step - delay_steps, arrays["src"][arrival_edge]]),
"A queued arrival has no matching presynaptic spike at the author delay")
_check(np.array_equal(arrival_accepted,
available_end[arrival_step, arrays["dst"][arrival_edge]].astype(np.uint8)),
"Observed arrival gating disagrees with the post-threshold refractory state")
source_counts_due = fired[:max(0, steps - delay_steps)].sum(axis=0, dtype=np.int64)
expected_arrivals = int(source_counts_due[arrays["src"]].sum(dtype=np.int64))
_check(len(arrival_edge) == expected_arrivals, "Queue trace omitted or duplicated an anatomical arrival")
result = {
"nodes": nodes.copy(), "step": event_ticks, "local_idx": local_indices,
"neuron_idx": nodes[local_indices], "counts": counts,
"v": voltage[-1].copy(), "g": conductance[-1].copy(), "ref": remaining[-1].copy(),
"history": history, "ref_history": remaining,
"not_refractory": available_end.astype(np.uint8),
"weights": np.asarray(syn.w_[:], dtype=dtype).copy(),
"initial_v": initial_v, "initial_g": initial_g,
"request_attempted": attempted, "request_accepted": accepted, "request_rejected": rejected,
"request_count": attempted.sum(axis=1, dtype=np.int32),
"applicable_count": accepted.sum(axis=1, dtype=np.int32),
"rejected_count": rejected.sum(axis=1, dtype=np.int32), "fired": fired,
"arrival_step": arrival_step, "arrival_edge_index": arrival_edge,
"arrival_accepted": arrival_accepted,
"external_step": ext_ticks.astype(np.int32), "external_target": ext_targets.astype(np.int32),
"external_weights": external_weights,
"external_accepted": available_end[ext_ticks, ext_targets].astype(np.uint8),
}
coefficients = _coefficients(neu, parameters, b2, dt_seconds)
_save(output / "coefficients.json", coefficients)
_save(output / "queue-audit.json", {
**queue_counters, "trace_enabled": trace_queue, "due_events": len(arrival_edge),
"accepted_due_events": int(arrival_accepted.sum()),
"all_due_events_preserved": True if trace_queue else None,
"ordering": "raw anatomical queue.peek local edge-index sequence; no sorting or filtering",
"proxy_contract": "peek returns original object; advance and all other methods forward unchanged",
})
_save_arrays(output / "state.npz", result)
return result, coefficients
except Exception as error:
_save(output / "error.json", {"type": type(error).__name__, "message": str(error)})
raise
finally:
for key, value in saved.items():
b2.prefs[key] = value
b2.defaultclock.dt = old_dt
def run_case(case_dir: Path, out_dir: Path, precision: str) -> dict:
"""Run one <=512-node, <=1000-step fixed case with exclusive outputs.
Main outputs are ``state.npz`` and ``coefficients.json``. Queue tracing is
enabled by default. ``spec.trace_queue`` can disable it; when enabled,
``spec.validate_queue_proxy`` defaults to True and repeats the same case
without instrumentation. That repeat must match every non-arrival array
exactly. An opt-out is recorded as unverified, never as a successful check.
This function changes global Brian preferences temporarily; run it serially.
"""
_check(precision in ("float64", "float32"), "precision must be float64 or float32")
case_dir, out_dir = Path(case_dir).resolve(strict=True), Path(out_dir).resolve()
spec, arrays = _load_case(case_dir)
source_bytes = AUTHOR_PATH.read_bytes()
source_hash = hashlib.sha256(source_bytes).hexdigest()
_check(source_hash == AUTHOR_SHA256, "Archived author source hash mismatch")
trace_queue = spec.get("trace_queue", True)
validate_proxy = spec.get("validate_queue_proxy", trace_queue)
_check(isinstance(trace_queue, bool) and isinstance(validate_proxy, bool), "Queue switches must be JSON booleans")
_check(not validate_proxy or trace_queue, "Proxy validation requires queue tracing")
out_dir.mkdir(parents=True, exist_ok=False)
try:
versions = {}
for name in ("brian2", "numpy", "pandas", "scipy", "pyarrow", "sympy", "cython"):
versions[name] = importlib.metadata.version(name)
_save(out_dir / "environment.json", {"python": sys.version, "platform": platform.platform(),
"packages": versions, "backend": "Brian2 numpy runtime"})
_check(versions["brian2"] == "2.5.1", "The stress oracle requires Brian2 2.5.1")
_check(versions["numpy"] == "1.22.3", "The stress oracle requires the recorded NumPy 1.22.3 promotion rules")
_save(out_dir / "case.json", {"spec": spec, "precision": precision,
"author_path": str(AUTHOR_PATH), "author_commit": AUTHOR_COMMIT, "author_sha256": source_hash,
"case_spec_sha256": _digest(case_dir / "spec.json"), "case_arrays_sha256": _digest(case_dir / "arrays.npz"),
"trace_queue": trace_queue, "validate_queue_proxy": validate_proxy,
"common_interface_preserved": True, "original_poi_interface_adopted": False,
"whole_brain": False, "learning": False})
result, coefficients = _run_once(spec, arrays, precision, out_dir / "instrumented",
source_bytes, trace_queue)
# Preserve the primary result even if the independent proxy check fails.
_save_arrays(out_dir / "state.npz", result)
_save(out_dir / "coefficients.json", coefficients)
validation = {"performed": False, "pass": None,
"reason": "Disabled by explicit case switch" if not validate_proxy else None}
if validate_proxy:
baseline, baseline_coefficients = _run_once(spec, arrays, precision, out_dir / "unwrapped-audit",
source_bytes, False)
checks = {key: bool(np.array_equal(value, baseline[key]))
for key, value in result.items() if not key.startswith("arrival_")}
coefficient_keys = ["gd", "c0", "cg", "cv", "threshold", "reset", "seed_volt"]
checks["coefficients"] = all(coefficients[key] == baseline_coefficients[key] for key in coefficient_keys)
validation = {"performed": True, "pass": bool(all(checks.values())),
"exact_arrays": checks,
"arrival_arrays_excluded": "An unwrapped queue is deliberately not observed"}
_save(out_dir / "proxy-validation.json", validation)
_check(not validate_proxy or validation["pass"], "Queue proxy changed an author result")
_check(_digest(AUTHOR_PATH) == source_hash, "Author source changed during the stress case")
summary = {"name": spec["name"], "precision": precision, "steps": int(spec["steps"]),
"nodes": len(arrays["nodes"]), "edges": len(arrays["edge_id"]),
"spikes": len(result["step"]), "requests": int(result["request_count"].sum()),
"applicable_requests": int(result["applicable_count"].sum()),
"rejected_requests": int(result["rejected_count"].sum()),
"due_arrivals": len(result["arrival_step"]),
"queue_trace_enabled": trace_queue, "queue_proxy_validated": validation["pass"],
"source_unchanged": True, "state_file": str(out_dir / "state.npz"),
"coefficient_file": str(out_dir / "coefficients.json"), "complete": True}
_save(out_dir / "summary.json", summary)
return {"arrays": result, "coefficients": coefficients, "summary": summary}
except Exception as error:
_save(out_dir / "error.json", {"type": type(error).__name__, "message": str(error)})
raise
@@ -0,0 +1,103 @@
"""Read-only selection and deterministic fixtures, before observing any dynamics."""
from datetime import datetime,timezone
import json
import subprocess
import numpy as np
import pandas as pd
import yaml
from context import MODEL,ROOT,EXP,AD,CFG,digest,save,sources,integrity
M=json.loads((AD/'protocol/preregistered.json').read_text())['config']['model']
def put(name,nodes,src,dst,raw,steps,description,edge_id=None,requests=None,external=None,v=None,g=None,metadata=None):
n=len(nodes);src=np.asarray(src,np.int32);dst=np.asarray(dst,np.int32);raw=np.asarray(raw,np.int64)
eid=np.arange(len(src),dtype=np.int64) if edge_id is None else np.asarray(edge_id,np.int64)
order=np.lexsort((eid,src));src,dst,raw,eid=src[order],dst[order],raw[order],eid[order]
assert len(src) and len(src)==len(dst)==len(raw) and np.all((src>=0)&(src<n)) and np.all((dst>=0)&(dst<n))
arrays=dict(nodes=np.asarray(nodes,np.int32),src=src,dst=dst,signed_raw=raw,edge_id=eid,
requests=np.zeros((steps,n),np.uint8) if requests is None else requests,
external_mv=np.zeros((steps,n),np.float64) if external is None else external,
initial_v_mv=np.full(n,M['resting_potential_mv'],np.float64) if v is None else v,
initial_g_mv=np.zeros(n,np.float64) if g is None else g)
folder=EXP/'cases'/name;folder.mkdir()
np.savez_compressed(folder/'arrays.npz',**arrays)
save(folder/'spec.json',dict(name=name,dt_ms=M['dt_ms'],steps=steps,duration_ms=steps*M['dt_ms'],model=M,
description=description,trace_queue=True,validate_queue_proxy=True,selection=metadata))
def synthetic():
n=33;t=160;raw=np.array([16777216,1,-16777216,1]+[(-1 if i%2 else 1)*(i+1) for i in range(28)],np.int64)
req=np.zeros((t,n),np.uint8);req[0,:32]=1;g=np.zeros(n);g[-1]=.123
put('ordered_cancellation',np.arange(-n,0),np.arange(32),np.full(32,32),raw,t,
'Pure synthetic numerical cancellation; huge counts are nonphysiological; nonzero initial g.',requests=req,g=g)
n=257;t=600;src=np.r_[np.arange(256),0,256];dst=np.full(258,256);raw=np.r_[np.where(np.arange(256)%2,-1,1)*(np.arange(256)%17+1)*13,19,11]
req=np.zeros((t,n),np.uint8)
for step in (0,21,22,40,62,101,150,300):req[step,:256]=1
req[0,256]=1;ext=np.zeros((t,n));ext[[0,21,22,23,300],256]=[1,1,1,-1,.5]
put('fanin_bursts',np.arange(-n,0),src,dst,raw,t,'256-source signed fan-in, labelled duplicate/loop; refractory requests and external writes.',requests=req,external=ext)
n=4;t=160;req=np.zeros((t,n),np.uint8);req[0,[0,1,2]]=1;req[22,2]=1;req[44,[0,1]]=1
ext=np.zeros((t,n));ext[[0,21,22,23,43,44,45],2]=[2,2,2,2,-1,1,1]
put('release_selfloop',np.arange(-n,0),[0,1,2,2,3],[2,2,2,3,2],[1,400,10,100,-1],t,'Release at22 with same-tick force/arrival, loops and ring wraps.',requests=req,external=ext)
dt=.0001;tm=.020;ts=.005;rest=-.052;th=-.045;cv=np.exp(-dt/tm);c0=rest-rest*cv
center=(th-c0)/cv;vals=[np.nextafter(center,-np.inf),center,np.nextafter(center,np.inf)]
center32=np.float32((np.float32(th)-np.float32(c0))/np.float32(cv));x=center32
for _ in range(3):x=np.nextafter(x,np.float32(-np.inf))
for _ in range(7):vals.append(float(x));x=np.nextafter(x,np.float32(np.inf))
n=len(vals);put('threshold_rounding',np.arange(-n,0),[0],[0],[0],100,
'Preintegration voltages bracket strict-threshold landing in both precisions; inert zero fixture edge.',v=np.asarray(vals)*1000)
rng=np.random.default_rng(CFG['seed']);n=128;t=400;src=np.repeat(np.arange(n),8);dst=rng.integers(0,n,len(src),dtype=np.int32)
raw=rng.integers(1,601,len(src),dtype=np.int64)*np.where(src%3==0,-1,1)
req=(rng.random((t,n))<.001).astype(np.uint8);req[200:]=0
g=rng.uniform(-1,1,n);ext=np.zeros((t,n));ext[::37]=rng.uniform(-.5,.5,(len(ext[::37]),n))
put('signed_random',np.arange(-n,0),src,dst,raw,t,'Fixed seed128-cell signed random fixture; forcing stops at20ms; signs assigned per synthetic source.',requests=req,external=ext,g=g)
def real():
can=ROOT/'datasets/male-cns-v1.0/processed/neuron-only';ptr=np.load(can/'gpu/csr-indptr.int64.npy',mmap_mode='r');dst=np.load(can/'gpu/csr-indices.int32.npy',mmap_mode='r');raw=np.load(can/'gpu/csr-weights.int64.npy',mmap_mode='r')
meta=pd.read_parquet(can/'neuron-index/unified-metadata.parquet').sort_values('neuron_idx').reset_index(drop=True)
assert np.array_equal(meta.neuron_idx,np.arange(166700)) and len(dst)==25582938 and int(raw.sum())==124177617
pol=yaml.safe_load((ROOT/'models/male-cns-lif/runtime-v1/phase4c.yaml').read_text())['sign_policy']
signs=[]
for nt in meta.neurotransmitter:
s='' if pd.isna(nt) else str(nt).strip().lower()
signs.append(pol['unknown_sign'] if not s else pol['unclear_sign'] if s=='unclear' else pol['neurotransmitter_sign'][s])
signs=np.asarray(signs,np.int8);src=np.repeat(np.arange(166700,dtype=np.int32),np.diff(ptr))
degree=np.bincount(dst,minlength=166700);pos=np.bincount(dst[signs[src]>0],minlength=166700);neg=np.bincount(dst[signs[src]<0],minlength=166700)
ranked=np.lexsort((np.arange(166700),-degree));targets=[int(ranked[0]),int(next(i for i in ranked if pos[i]>=64 and neg[i]>=63))]
for name,target in zip(CFG['real_cases'],targets):
incoming=np.flatnonzero((dst==target)&(src!=target));incoming=incoming[np.lexsort((src[incoming],-raw[incoming]))]
if name=='highest_indegree':chosen=incoming[:127]
else:
pp=incoming[signs[src[incoming]]>0][:64];nn=incoming[signs[src[incoming]]<0][:63]
assert len(pp)==64 and len(nn)==63,'Insufficient mixed partners for predeclared graph'
chosen=np.r_[pp,nn]
nodes=np.sort(np.unique(np.r_[src[chosen],target])).astype(np.int32);assert len(nodes)==128
lookup=np.full(166700,-1,np.int32);lookup[nodes]=np.arange(len(nodes),dtype=np.int32)
ids=np.concatenate([np.arange(ptr[s],ptr[s+1],dtype=np.int64) for s in nodes]);ids=ids[lookup[dst[ids]]>=0]
local_src=lookup[src[ids]];local_dst=lookup[dst[ids]];signed=raw[ids]*signs[src[ids]]
req=np.zeros((500,len(nodes)),np.uint8);stim=np.sort(src[chosen[:16]])
for tick in (0,25,50,100,150,200):req[tick,lookup[stim]]=1
put(name,nodes,local_src,local_dst,signed,500,'Deterministically selected real induced graph; canonical edge IDs/counts and frozen signs.',edge_id=ids,requests=req,
metadata=dict(target_neuron_idx=target,target_body_id=int(meta.body_id.iloc[target]),global_in_degree=int(degree[target]),positive_in_edges=int(pos[target]),negative_in_edges=int(neg[target]),stimulated_neuron_idx=stim.tolist()))
meta.iloc[nodes][['neuron_idx','body_id','annotation_type','neurotransmitter']].to_csv(EXP/'cases'/name/'identities.csv',index=False)
pd.DataFrame(dict(edge_id=ids,source_neuron_idx=src[ids],target_neuron_idx=dst[ids],raw_count=raw[ids],source_sign=signs[src[ids]],signed_raw=signed)).to_csv(EXP/'cases'/name/'anatomical-edges.csv',index=False)
ac=ROOT/'experiments/phase4ac-lif-semantics'
for name in CFG['prior_motifs']:
e=pd.read_csv(ac/'motifs'/f'{name}-edges.csv');nodes=np.load(ac/'motifs'/f'{name}-state.npz')['nodes'];look={int(v):i for i,v in enumerate(nodes)}
req=np.zeros((1000,len(nodes)),np.uint8)
if name!='no_input':req[0,look[1381]]=1
put('motif_'+name,nodes,[look[int(x)] for x in e.source_neuron_idx],[look[int(x)] for x in e.target_neuron_idx],e.raw_count.to_numpy()*signs[e.source_neuron_idx],1000,
'Same selected4AC/4AD motif,100ms prefix; no new edge or weight.',edge_id=e.edge_id.to_numpy(),requests=req)
def main():
assert not (EXP/'protocol').exists()
for d in ('protocol','cases','oracle','gpu','regression','audit','performance'):(EXP/d).mkdir()
can=ROOT/'datasets/male-cns-v1.0/processed/neuron-only'
inputs=[can/'gpu/csr-indptr.int64.npy',can/'gpu/csr-indices.int32.npy',can/'gpu/csr-weights.int64.npy',can/'neuron-index/unified-metadata.parquet',ROOT/'models/male-cns-lif/runtime-v1/phase4c.yaml',AD/'external-resources/model.py',AD/'manifest.json']
save(EXP/'protocol/preregistered.json',dict(timestamp_utc=datetime.now(timezone.utc).isoformat(),config=CFG,sources=sources(),
inputs=[dict(path=str(p.relative_to(ROOT)),sha256=digest(p)) for p in inputs],before_model_execution=True,
parent_commit=subprocess.check_output(['git','-C',str(ROOT),'rev-parse','HEAD'],text=True).strip()))
save(EXP/'regression/old-integrity-before.json',integrity())
synthetic();real()
save(EXP/'protocol/case-hashes.json',[dict(path=str(p.relative_to(EXP)),sha256=digest(p)) for p in sorted((EXP/'cases').rglob('*')) if p.is_file()])
print('PREPARED',len(list((EXP/'cases').iterdir())),'fixed cases; old integrity PASS',flush=True)
if __name__=='__main__':main()
@@ -0,0 +1,87 @@
"""Fixed Brian/GPU stress suite, with source-checked immutable inputs."""
import argparse
import json
import os
from pathlib import Path
import sys
import time
import numpy as np
from context import ROOT,MODEL,EXP,CFG,save,check
def load(path):
with np.load(path,allow_pickle=False) as f:return {k:f[k].copy() for k in f.files}
def compare(a,b,case,precision):
checks={k:bool(np.array_equal(a[k],b[k])) for k in ('step','neuron_idx','counts','ref')}
ah=np.asarray(a['history'],np.float64);bh=np.asarray(b['history'],np.float64)
assert ah.shape==bh.shape
delta=np.abs(ah[:,:2]-bh[:,:2]);limit=CFG['max_abs_volt'][precision]
checks['history_ref_exact']=bool(np.array_equal(ah[:,2],bh[:,2]))
checks['finite_states']=bool(np.isfinite(ah).all() and np.isfinite(bh).all())
checks['state_tolerance']=bool(delta.max()<=limit)
checks['weights_exact']=bool(np.array_equal(a['weights'],b['weights']))
# Distinguish logical event requests from realised requests after refractory.
for key in ('request_attempted','request_accepted','request_rejected'):
checks[key+'_exact']=bool(key in a and key in b and np.array_equal(a[key],b[key]))
data=load(case/'arrays.npz');targets=data['dst']
def arrival_sequence(result):
tick=result['arrival_step'];eid=result['arrival_edge_index'];guard=result['arrival_accepted']
assert tick.shape==eid.shape==guard.shape
# Different targets can interleave; order of writes to one target cannot.
order=np.lexsort((np.arange(len(eid)),targets[eid],tick))
return np.stack((tick[order],eid[order],guard[order]),axis=1)
aa=arrival_sequence(a);bb=arrival_sequence(b)
checks['arrival_order_and_guard_exact']=bool(np.array_equal(aa,bb))
first=None
for t in range(len(ah)):
if not np.array_equal(a['neuron_idx'][a['step']==t],b['neuron_idx'][b['step']==t]):first=t;break
bits_equal=a['history'].dtype==b['history'].dtype and np.array_equal(np.ascontiguousarray(a['history']).view(np.uint8),np.ascontiguousarray(b['history']).view(np.uint8))
return dict(pass_all=all(checks.values()),checks=checks,bitwise_history=bool(bits_equal),
max_abs_volt=float(delta.max()),mean_abs_volt=float(delta.mean()),fixed_limit_volt=limit,
first_event_divergence_step=first,oracle_events=int(len(a['step'])),gpu_events=int(len(b['step'])),
due_edges=int(len(aa)),accepted_edges=int(a['arrival_accepted'].sum()),
attempted_requests=int(a['request_attempted'].sum()),accepted_requests=int(a['request_accepted'].sum()))
def oracle():
from oracle_stress import run_case
t0=time.perf_counter();rows=[]
for case in sorted((EXP/'cases').iterdir()):
for precision in CFG['precisions']:
output=EXP/'oracle'/(case.name+'-'+precision)
run_case(case,output,precision)
r=load(output/'state.npz');assert np.isfinite(r['history']).all()
rows.append(dict(case=case.name,precision=precision,events=len(r['step'])))
print(json.dumps(rows[-1]),flush=True)
save(EXP/'regression/oracle-complete.json',dict(cases=rows,wall_seconds=time.perf_counter()-t0))
def gpu():
assert (EXP/'regression/oracle-complete.json').exists(),'No complete oracle suite'
os.environ['CUDA_VISIBLE_DEVICES']=CFG['gpu_uuid'];os.environ['CUDA_HOME']='/usr/local/cuda'
os.environ['PATH']=str(ROOT/'models/male-cns-lif/.venv/bin')+':/usr/local/cuda/bin:'+os.environ['PATH']
os.environ['MAX_JOBS']='2';os.environ['TORCH_CUDA_ARCH_LIST']='8.9'
sys.path.insert(0,str(ROOT/'models/male-cns-lif/terminal-scope-diagnostic-v1'))
from diagnostic import queues
queues()
from gpu_backend import run_gpu
t0=time.perf_counter();rows=[]
for case in sorted((EXP/'cases').iterdir()):
for precision in CFG['precisions']:
queues();refdir=EXP/'oracle'/(case.name+'-'+precision)
reference=load(refdir/'state.npz');coeff=json.loads((refdir/'coefficients.json').read_text());first=None
for repeat in range(CFG['gpu_repeats']):
out=EXP/'gpu'/(case.name+'-'+precision+f'-r{repeat}')
run_gpu(case,out,precision,coeff,EXP/'extension-build',weights=reference['weights'])
result=load(out/'state.npz');cmp=compare(reference,result,case,precision)
if first is None:first=result
rep={k:bool(np.array_equal(v,result[k])) for k,v in first.items()}
row=dict(case=case.name,precision=precision,repeat=repeat,comparison=cmp,repeat_arrays=rep,
pass_all=cmp['pass_all'] and all(rep.values()))
save(out/'comparison.json',row);rows.append(row)
print(json.dumps(dict(case=case.name,precision=precision,repeat=repeat,pass_all=row['pass_all'],
error_volt=cmp['max_abs_volt'],events=cmp['oracle_events'],bitwise_history=cmp['bitwise_history'])),flush=True)
save(EXP/'regression/gpu-complete.json',dict(pass_all=all(x['pass_all'] for x in rows),trials=rows,
wall_seconds=time.perf_counter()-t0,whole_brain_tested=False,learning_achieved=False))
if __name__=='__main__':
p=argparse.ArgumentParser();p.add_argument('mode',choices=['oracle','gpu']);args=p.parse_args()
check();globals()[args.mode]()
@@ -0,0 +1,33 @@
"""Run old frozen trials into new output only; no author-model equivalence claim."""
import os
import sys
import time
import numpy as np
import pandas as pd
import yaml
from context import ROOT,EXP,CFG,save,check
def main():
check();assert not (EXP/'regression/frozen-five.json').exists()
rtpath=ROOT/'models/male-cns-lif/deterministic-runtime-v2';pc=yaml.safe_load((rtpath/'phase4p.yaml').read_text())
os.environ['CUDA_VISIBLE_DEVICES']=CFG['gpu_uuid'];os.environ['CUDA_HOME']='/usr/local/cuda'
os.environ['PATH']=str(ROOT/'models/male-cns-lif/.venv/bin')+':/usr/local/cuda/bin:'+os.environ['PATH']
os.environ['MAX_JOBS']='2';os.environ['TORCH_CUDA_ARCH_LIST']='8.9'
sys.path.insert(0,str(ROOT/'models/male-cns-lif/terminal-scope-diagnostic-v1'))
from diagnostic import queues
queues();sys.path.insert(0,str(ROOT/'models/male-cns-lif'));sys.path.insert(0,str(rtpath))
import torch,phase4p
assert torch.cuda.mem_get_info()[0]>2*1024**3
torch.cuda.reset_peak_memory_stats();t0=time.perf_counter();rt=phase4p.CanonicalRuntime(pc,rtpath,EXP);rows=[]
for spec in pc['reference_trials']:
queues();base=ROOT/'experiments/phase4p-deterministic-runtime/phase4c-v2/reference'/spec['name']
ev=pd.read_parquet(base/'spike-events.parquet')
ref=dict(step=ev.step.to_numpy(),neuron_idx=ev.neuron_idx.to_numpy(),counts=np.load(base/'per-neuron-spike-counts.int32.npy'),v=np.load(base/'final-voltage.float32.npy'),g=np.load(base/'final-g.float32.npy'))
r=rt.run(np.flatnonzero(rt.pops[spec['population']]),spec['rate_hz'],spec['seed'],spec['duration_ms'])
row=dict(trial=spec['name'],exact={k:bool(np.array_equal(v,r[k])) for k,v in ref.items()},wall_seconds=r['metrics']['wall_seconds'])
rows.append(row);print(row,flush=True)
save(EXP/'regression/frozen-five.json',dict(pass_all=all(all(r['exact'].values()) for r in rows),trials=rows,
wall_seconds=time.perf_counter()-t0,peak_allocated_vram=int(torch.cuda.max_memory_allocated()),oracle='old frozen deterministic-v2 only'))
assert all(all(r['exact'].values()) for r in rows)
if __name__=='__main__':main()