Files
2026-09-30 09:26:21 +08:00

432 lines
24 KiB
Python

"""Independent extraction/replay of the frozen legacy APL component.
Only the initializer prefix before ``self.base`` and two original APL update
statements are executed. No old simulator is instantiated, no learning code is
executed, and no frozen artifact is written. The Torch replay is a component
oracle, not an independent CPU whole-brain closed-loop oracle.
"""
from __future__ import annotations
import argparse
import ast
import copy
from functools import lru_cache
import hashlib
import json
import math
from pathlib import Path
import struct
import types
import numpy as np
import pyarrow.parquet as pq
import torch
ROOT = Path('/data/digital-drosophila')
CANONICAL = ROOT/'datasets/male-cns-v1.0/processed/neuron-only'
SOURCE = ROOT/'models/male-cns-lif/deterministic-runtime-v2/phase4p_revalidate.py'
SOURCE_SHA256 = 'c7e693a2da70acdce4f53af5c2853e52e8333acc93f88a46a5e8083eac639346'
N, A, K = 166700, 2, 4064
BUNDLE_ARRAYS = {
'input_matrix': (np.float32, (A, N)), 'output_matrix': (np.float32, (A, K)),
'input_raw': (np.float32, (A, N)), 'output_raw': (np.float32, (A, K)),
'kc_globals': (np.int32, (K,)), 'apl_globals': (np.int32, (A,)),
'pn_globals': (np.int32, None), 'apl_out_ids': (np.int64, (5845,)),
}
def require(condition, message):
if not condition:
raise ValueError(message)
def digest(path):
h = hashlib.sha256()
with Path(path).open('rb') as source:
for block in iter(lambda: source.read(8 << 20), b''):
h.update(block)
return h.hexdigest()
def array_hash(value):
return hashlib.sha256(np.ascontiguousarray(value).tobytes()).hexdigest()
def bit_equal(a, b):
a, b = np.asarray(a), np.asarray(b)
return bool(a.dtype == b.dtype and a.shape == b.shape and
np.array_equal(np.ascontiguousarray(a).view(np.uint8),
np.ascontiguousarray(b).view(np.uint8)))
def save(path, value):
path = Path(path)
path.parent.mkdir(parents=True, exist_ok=True)
with path.open('x', encoding='utf-8') as target:
json.dump(value, target, indent=2, sort_keys=True, allow_nan=False)
target.write('\n')
class ReadOnlyNamespace:
"""The extracted prefix can read rt, but cannot mutate or add rt fields."""
def __init__(self, **values):
object.__setattr__(self, '_values', types.MappingProxyType(values))
def __getattr__(self, name):
return self._values[name]
def __setattr__(self, name, value):
raise AttributeError('Read-only canonical runtime view')
def _code(nodes, filename):
module = ast.Module(body=copy.deepcopy(nodes), type_ignores=[])
return compile(ast.fix_missing_locations(module), str(filename), 'exec')
def _self_attribute(target, field):
return isinstance(target, ast.Attribute) and isinstance(target.value, ast.Name) and \
target.value.id == 'self' and target.attr == field
def extract_statements():
data = SOURCE.read_bytes()
require(hashlib.sha256(data).hexdigest() == SOURCE_SHA256, 'Frozen legacy source SHA-256 mismatch')
text = data.decode('utf-8')
tree = ast.parse(text, filename=str(SOURCE))
classes = [node for node in tree.body if isinstance(node, ast.ClassDef) and node.name == 'DeterministicAPLSim']
require(len(classes) == 1, 'Expected exactly one frozen DeterministicAPLSim')
methods = {node.name: node for node in classes[0].body if isinstance(node, ast.FunctionDef)}
init, run = methods['__init__'], methods['run']
cuts = [i for i, node in enumerate(init.body) if isinstance(node, ast.Assign) and
any(_self_attribute(target, 'base') for target in node.targets)]
require(len(cuts) == 1 and cuts[0] > 0, 'Missing/ambiguous self.base cutoff')
prefix = init.body[:cuts[0]]
# Executing no import or file-writing call is allowed in the extracted prefix.
require(not any(isinstance(node, (ast.Import, ast.ImportFrom, ast.With, ast.AsyncWith))
for stmt in prefix for node in ast.walk(stmt)), 'Unexpected side-effect construct in initializer prefix')
expressions = [node for node in ast.walk(run) if isinstance(node, ast.Expr) and
isinstance(node.value, ast.Call) and isinstance(node.value.func, ast.Attribute) and
node.value.func.attr == 'add_' and isinstance(node.value.func.value, ast.Call) and
isinstance(node.value.func.value.func, ast.Attribute) and
isinstance(node.value.func.value.func.value, ast.Name) and
node.value.func.value.func.value.id == 'apl_state' and node.value.func.value.func.attr == 'mul_']
inhibition = [node for node in ast.walk(run) if isinstance(node, ast.Assign) and
any(isinstance(target, ast.Name) and target.id == 'inhibition' for target in node.targets)]
decay = [node for node in ast.walk(run) if isinstance(node, ast.Assign) and
any(isinstance(target, ast.Name) and target.id == 'apl_decay' for target in node.targets)]
require(len(expressions) == len(inhibition) == len(decay) == 1, 'Ambiguous frozen APL operator extraction')
require(expressions[0].lineno < inhibition[0].lineno, 'Unexpected frozen APL update order')
rhs = inhibition[0].value
require(isinstance(rhs, ast.BinOp) and isinstance(rhs.op, ast.Mult), 'Unexpected inhibition gain expression')
gain_expression = ast.Expression(copy.deepcopy(rhs.right))
provenance = {'source': str(SOURCE), 'sha256': SOURCE_SHA256,
'initializer_line_start': prefix[0].lineno, 'initializer_line_end': prefix[-1].end_lineno,
'excluded_self_base_line': init.body[cuts[0]].lineno,
'decay_line': decay[0].lineno, 'state_update_line': expressions[0].lineno,
'inhibition_line': inhibition[0].lineno,
'initializer_source': '\n'.join(ast.get_source_segment(text, node) for node in prefix),
'decay_source': ast.get_source_segment(text, decay[0]),
'state_update_source': ast.get_source_segment(text, expressions[0]),
'inhibition_source': ast.get_source_segment(text, inhibition[0]),
'learning_initializer_suffix_executed': False, 'original_run_loop_executed': False}
return (_code(prefix, SOURCE), _code(decay, SOURCE),
_code([expressions[0], inhibition[0]], SOURCE),
compile(ast.fix_missing_locations(gain_expression), str(SOURCE), 'eval'), provenance)
def _canonical_files():
return {name: CANONICAL/path for name, path in {
'metadata': 'neuron-index/unified-metadata.parquet',
'indptr': 'gpu/csr-indptr.int64.npy', 'indices': 'gpu/csr-indices.int32.npy',
'raw': 'gpu/csr-weights.int64.npy'}.items()}
@lru_cache(maxsize=2)
def reconstruct(canonical_hashes, device_name):
"""Cache only reconstructions tied to exact canonical bytes and source hash."""
files = _canonical_files()
require(dict(canonical_hashes) == {name: digest(path) for name, path in files.items()},
'Canonical bytes changed before matrix reconstruction')
arrays = {name: np.load(files[name], mmap_mode='r', allow_pickle=False)
for name in ('indptr', 'indices', 'raw')}
require(arrays['indptr'].dtype == np.int64 and arrays['indices'].dtype == np.int32 and
arrays['raw'].dtype == np.int64, 'Canonical CSR dtypes changed')
require(arrays['indptr'].shape == (N + 1,) and arrays['indptr'][-1] == 25582938 and
arrays['indices'].shape == arrays['raw'].shape == (25582938,), 'Canonical graph shape changed')
require(int(arrays['raw'].sum()) == 124177617 and bool(np.all(arrays['raw'] > 0)), 'Canonical raw weights changed')
table = pq.read_table(files['metadata'])
require(table.num_rows == N, 'Canonical metadata count changed')
rt = ReadOnlyNamespace(n=N, dt=.1, device=torch.device(device_name), table=table,
indptr_np=arrays['indptr'], indices_np=arrays['indices'], raw_np=arrays['raw'])
view = types.SimpleNamespace()
init_code, decay_code, update_code, gain_code, provenance = extract_statements()
scope = {'self': view, 'rt': rt, 'np': np, 'torch': torch}
exec(init_code, scope)
require(not hasattr(view, 'base') and not hasattr(view, 'identity'), 'Initializer crossed the learning cutoff')
expected = {'input_matrix': view.input_gpu.detach().cpu().numpy().copy(),
'output_matrix': view.output_gpu.detach().cpu().numpy().copy(),
'input_raw': scope['input_matrix'], 'output_raw': scope['output_matrix'],
'kc_globals': view.kc, 'apl_globals': view.apl, 'pn_globals': view.pn,
'apl_out_ids': view.apl_out_ids}
for value in expected.values():
value.flags.writeable = False
return view, expected, provenance, (decay_code, update_code, gain_code)
class TorchObservationProxy:
"""Observe the result of the actual extracted torch.mv call unchanged."""
def __init__(self):
self.drive = None
self.calls = 0
def __getattr__(self, name):
return getattr(torch, name)
def mv(self, *args, **kwargs):
value = torch.mv(*args, **kwargs)
self.drive = value
self.calls += 1
return value
def _up(value):
return np.nextafter(np.asarray(value, dtype=np.float64), np.inf)
def _add(a, b):
return _up(np.asarray(a, dtype=np.float64) + np.asarray(b, dtype=np.float64))
def _mul(a, b):
return _up(np.asarray(a, dtype=np.float64) * np.asarray(b, dtype=np.float64))
def _gamma(count, unit_roundoff):
require(count * unit_roundoff < 1, 'Gamma bound is undefined')
# Outward rounding on the division; the denominator is rounded downward.
numerator = _mul(count, unit_roundoff)
denominator = np.nextafter(1.0 - numerator, -np.inf)
return _up(numerator / denominator)
def _positive_dot_upper(state, matrix):
result = np.zeros(matrix.shape[1], np.float64)
for index in range(len(state)):
result = _add(result, _mul(state[index], matrix[index]))
return result
def _normal_or_zero(value, tiny):
value = np.asarray(value)
return bool(np.all((value == 0) | (value >= tiny)))
def _new_metric():
return {'max_abs_error': 0.0, 'sum_abs_error': 0.0, 'elements': 0,
'max_relative_error_nonzero_reference': 0.0, 'zero_reference_nonzero_actual': 0,
'max_error_bound': 0.0, 'max_bound_utilization': 0.0, 'bound_violations': 0}
def _accumulate(metric, actual, reference, bound):
error = np.abs(np.asarray(actual, np.float64) - reference)
error_upper = np.where(error == 0, 0.0, _up(error))
nonzero = reference != 0
ratio = np.divide(error, np.abs(reference), out=np.zeros_like(error), where=nonzero)
usage = np.divide(error_upper, bound, out=np.zeros_like(error), where=bound > 0)
metric['max_abs_error'] = max(metric['max_abs_error'], float(error.max(initial=0)))
metric['sum_abs_error'] += float(error.sum(dtype=np.float64))
metric['elements'] += error.size
metric['max_relative_error_nonzero_reference'] = max(metric['max_relative_error_nonzero_reference'], float(ratio.max(initial=0)))
metric['zero_reference_nonzero_actual'] += int(np.count_nonzero((~nonzero) & (error != 0)))
metric['max_error_bound'] = max(metric['max_error_bound'], float(np.max(bound, initial=0)))
metric['max_bound_utilization'] = max(metric['max_bound_utilization'], float(usage.max(initial=0)))
metric['bound_violations'] += int(np.count_nonzero(error_upper > bound))
def validate(gpu_case_dir, bundle_dir, out_json):
"""Validate one completed fixture; write one exclusive pure-JSON audit."""
gpu_case_dir, bundle_dir, out_json = Path(gpu_case_dir), Path(bundle_dir), Path(out_json)
require(not out_json.exists(), 'Refusing to overwrite an existing component audit')
report = {'pass_all': False, 'scope': 'Frozen APL component extraction and recorded-spike replay only',
'independent_cpu_closed_loop_oracle': False, 'brain_or_learning_parameters_changed': False,
'candidate_function_called': False, 'gpu_case_dir': str(gpu_case_dir), 'bundle_dir': str(bundle_dir)}
try:
require(digest(SOURCE) == SOURCE_SHA256, 'Frozen source changed before component validation')
require(torch.cuda.is_available() and torch.cuda.device_count() == 1 and
'4090' in torch.cuda.get_device_name(), 'Require the one declared RTX4090')
environment = {'torch': str(torch.__version__), 'numpy': np.__version__,
'gpu': torch.cuda.get_device_name(), 'cuda': torch.version.cuda,
'cuda_matmul_allow_tf32': torch.backends.cuda.matmul.allow_tf32,
'float32_matmul_precision': torch.get_float32_matmul_precision(),
'deterministic_algorithms': torch.are_deterministic_algorithms_enabled(),
'precision_flags_changed': False}
report['environment'] = environment
canonical_hashes = tuple(sorted((name, digest(path)) for name, path in _canonical_files().items()))
view, reconstructed, provenance, code = reconstruct(canonical_hashes, 'cuda:0')
bundle = {}
matrix_checks = {}
for name, (dtype, shape) in BUNDLE_ARRAYS.items():
value = np.load(bundle_dir/(name+'.npy'), mmap_mode='r', allow_pickle=False)
require(value.dtype == dtype and value.ndim == len(shape or (0,)) and
(shape is None or value.shape == shape), 'Bundle schema mismatch: '+name)
bundle[name] = value
matrix_checks[name] = bit_equal(value, reconstructed[name])
scales = json.loads((bundle_dir/'scales.json').read_text())
for name in ('input_scale', 'output_scale'):
matrix_checks[name] = struct.pack('<d', float(scales[name])) == struct.pack('<d', float(getattr(view, name)))
report['matrix_checks'] = matrix_checks
report['extracted_source'] = provenance
report['canonical_sha256'] = dict(canonical_hashes)
report['bundle_sha256'] = {name+'.npy': digest(bundle_dir/(name+'.npy')) for name in BUNDLE_ARRAYS}
report['bundle_sha256']['scales.json'] = digest(bundle_dir/'scales.json')
report['reconstructed_array_sha256'] = {name: array_hash(value) for name, value in reconstructed.items()}
bundle_provenance = json.loads((bundle_dir/'provenance.json').read_text())
require(Path(bundle_provenance['canonical_root']).resolve() == CANONICAL.resolve(),
'Bundle provenance names another canonical dataset')
report['bundle_sha256']['provenance.json'] = digest(bundle_dir/'provenance.json')
require(all(matrix_checks.values()), 'Bundle differs from the actual frozen initializer prefix')
traces = {name: np.load(gpu_case_dir/('apl-'+name+'.float32.npy'), mmap_mode='r', allow_pickle=False)
for name in ('drive', 'state', 'inhibition')}
steps = len(traces['state'])
require(steps > 0, 'Empty APL trajectory cannot pass')
for name, width in (('drive', A), ('state', A), ('inhibition', K)):
require(traces[name].dtype == np.float32 and traces[name].shape == (steps, width) and
bool(np.isfinite(traces[name]).all()) and bool(np.all(traces[name] >= 0)), 'Invalid APL trace: '+name)
summary = json.loads((gpu_case_dir/'summary.json').read_text())
require(summary.get('complete') is True and summary.get('pass_all') is True and
summary.get('steps') == steps and summary.get('completed_ticks') == steps,
'GPU fixture did not complete its internal execution checks')
with np.load(gpu_case_dir/'state.npz', allow_pickle=False) as archive:
event_step, event_node = archive['step'].copy(), archive['neuron_idx'].copy()
require(event_step.dtype == event_node.dtype == np.int32 and event_step.ndim == 1 and
event_step.shape == event_node.shape, 'Invalid sparse actual spike event schema')
require(bool(np.all((event_step >= 0) & (event_step < steps))) and
bool(np.all((event_node >= 0) & (event_node < N))), 'Sparse event outside its clock/identity range')
require(np.array_equal(np.lexsort((event_node, event_step)), np.arange(len(event_step))), 'Actual events not sorted')
keys = event_step.astype(np.int64)*N + event_node
require(len(keys) < 2 or bool(np.all(np.diff(keys) > 0)), 'Duplicate actual spikes')
boundaries = np.searchsorted(event_step, np.arange(steps + 1))
proxy = TorchObservationProxy()
# Reconstructed original matrices are already on the original GPU/dtype.
replay_rt = types.SimpleNamespace(spike_u8=torch.zeros(N, dtype=torch.uint8, device=view.device))
shadow = types.SimpleNamespace(input_gpu=view.input_gpu, output_gpu=view.output_gpu, rt=replay_rt, dt=.1)
scope = {'self': shadow, 'np': np, 'math': math, 'torch': proxy,
'apl_state': torch.zeros(A, dtype=torch.float32, device=view.device)}
decay_code, update_code, gain_code = code
exec(decay_code, scope)
decay32 = np.float32(scope['apl_decay'])
gain32 = np.float32(eval(gain_code, {'np': np}))
require(0 < decay32 < 1 and gain32 > 0, 'Invalid extracted positive APL constants')
exact = {name: True for name in traces}
first_difference = None
I64 = np.asarray(bundle['input_matrix'], np.float64)
O64 = np.asarray(bundle['output_matrix'], np.float64)
require(bool(np.all(I64 >= 0)) and bool(np.all(O64 >= 0)), 'Positive-sum theory requires nonnegative matrices')
d, q = float(decay32), float(gain32)
u32, u64 = 2.0**-24, 2.0**-53
gd32, gd64 = _gamma(N-1, u32), _gamma(N-1, u64)
g232, g264 = _gamma(2, u32), _gamma(2, u64)
gh32, gh64 = _gamma(2*A+1, u32), _gamma(2*A+1, u64)
state64 = np.zeros(A, np.float64)
e32 = np.zeros(A, np.float64)
e64 = np.zeros(A, np.float64)
previous32 = np.zeros(A, np.float64)
statistics = {name: _new_metric() for name in traces}
normal_arithmetic = _normal_or_zero(I64, np.finfo(np.float32).tiny) and _normal_or_zero(O64, np.finfo(np.float32).tiny)
max_active_terms = 0
with torch.no_grad():
for tick in range(steps):
fired = event_node[boundaries[tick]:boundaries[tick+1]]
replay_rt.spike_u8.zero_()
if len(fired):
replay_rt.spike_u8[torch.as_tensor(fired.astype(np.int64), device=view.device)] = 1
before_calls = proxy.calls
exec(update_code, scope)
require(proxy.calls == before_calls + 1, 'Original update did not issue exactly one mv')
actual = {'drive': proxy.drive.detach().cpu().numpy(),
'state': scope['apl_state'].detach().cpu().numpy(),
'inhibition': scope['inhibition'].detach().cpu().numpy()}
for name in traces:
same = bit_equal(actual[name], traces[name][tick])
exact[name] &= same
if not same and first_difference is None:
first_difference = {'tick': tick, 'field': name,
'max_abs_error': float(np.max(np.abs(actual[name].astype(np.float64)-traces[name][tick])))}
# Float64 reference is measurement only and never feeds either brain.
selected = I64[:, fired]
drive64 = selected.sum(axis=1, dtype=np.float64)
max_active_terms = max(max_active_terms, int(np.count_nonzero(selected, axis=1).max(initial=0)))
Dupper = _up(drive64 / np.nextafter(1.0-gd64, -np.inf))
drive_bound = _add(_mul(gd32, Dupper), _mul(gd64, Dupper))
state32 = traces['state'][tick].astype(np.float64)
normal_arithmetic &= _normal_or_zero(d*previous32, np.finfo(np.float32).tiny)
normal_arithmetic &= _normal_or_zero(O64*state32[:, None], np.finfo(np.float32).tiny)
e32 = _add(_mul(d, e32), _add(_mul(g232, _add(_mul(d, previous32), Dupper)),
_mul(1.0+u32, _mul(gd32, Dupper))))
e64 = _add(_mul(d, e64), _add(_mul(g264, _add(_mul(d, state64), Dupper)),
_mul(1.0+u64, _mul(gd64, Dupper))))
state64 = d*state64 + drive64
inhibition64 = np.matmul(state64, O64)*q
h32bound = _add(_mul(q, _positive_dot_upper(e32, O64)),
_mul(gh32, _mul(q, _positive_dot_upper(state32, O64))))
h64bound = _add(_mul(q, _positive_dot_upper(e64, O64)),
_mul(gh64, _mul(q, _positive_dot_upper(state64, O64))))
normal_arithmetic &= _normal_or_zero(q*np.matmul(state32, O64), np.finfo(np.float32).tiny)
_accumulate(statistics['drive'], traces['drive'][tick], drive64, drive_bound)
_accumulate(statistics['state'], traces['state'][tick], state64, _add(e32, e64))
_accumulate(statistics['inhibition'], traces['inhibition'][tick], inhibition64, _add(h32bound, h64bound))
previous32 = state32
for metric in statistics.values():
metric['mean_abs_error'] = metric.pop('sum_abs_error') / metric['elements']
theory_applicable = bool(normal_arithmetic and not environment['cuda_matmul_allow_tf32'] and
environment['float32_matmul_precision'] == 'highest')
report.update(steps=steps, actual_spike_events=len(event_step),
ast_replay_bitwise_exact=exact, first_ast_difference=first_difference,
value_units={'drive': 'legacy normalized-contact spike drive',
'state': 'legacy normalized-drive accumulator', 'inhibition': 'legacy mV'},
decay_f32_bits=f'{int(decay32.view(np.uint32)):08x}', gain_f32_bits=f'{int(gain32.view(np.uint32)):08x}',
float64_measurement=statistics,
theoretical_bound={'applicable': theory_applicable, 'normal_or_zero_arithmetic': bool(normal_arithmetic),
'gamma_definition': 'gamma_n=(n*u)/(1-n*u), upward-rounded; u32=2^-24,u64=2^-53',
'drive_n': N-1, 'max_observed_nonzero_drive_terms': max_active_terms,
'gamma_drive_f32': float(gd32), 'gamma_drive_f64': float(gd64),
'state_recurrence': 'E_p[t] <= d*E_p[t-1] + gamma_2,p*(d*abs(s_p[t-1])+D_upper) + (1+u_p)*gamma_(N-1),p*D_upper',
'output_bound': 'q*sum(B*E_p) + gamma_(2*A+1),p*q*sum(B*abs(s_p)); compare with E32+E64',
'assumptions': 'nonnegative f32 matrices, spikes exactly 0/1, fixed stored f32 decay/gain, normal finite RN binary32/64 operations; no reduced-precision matmul',
'bound_arithmetic': 'Each positive bound operation is rounded toward +infinity via nextafter; denominator toward -infinity',
'tolerance_fitted_to_observed_error': False,
'f64_values_used_to_drive_brain': False,
'no_cpu_gpu_bitwise_expectation_for_f64_measurement': True},
input_sha256={name: digest(gpu_case_dir/name) for name in
('state.npz', 'summary.json', 'apl-drive.float32.npy', 'apl-state.float32.npy', 'apl-inhibition.float32.npy')})
report['pass_all'] = bool(all(matrix_checks.values()) and all(exact.values()) and theory_applicable and
all(value['bound_violations'] == 0 for value in statistics.values()))
report['limitations'] = [
'AST shadow and candidate use the same original Torch/CUDA primitives; this validates preserved component semantics, not an independent CPU whole-brain oracle.',
'Float64 is a non-feedback mathematical measurement with explicit gamma error bounds, not replacement APL dynamics.',
'Replay uses the recorded actual spike train; it does not establish equivalence for untested closed-loop inputs.',
'The N-1 positive-sum bound is deliberately conservative; passing it does not prove physiological validity or learning.']
save(out_json, report)
return report
except Exception as error:
report['error'] = {'type': type(error).__name__, 'message': str(error)}
if not out_json.exists():
save(out_json, report)
raise
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('gpu_case_dir', type=Path)
parser.add_argument('bundle_dir', type=Path)
parser.add_argument('out_json', type=Path)
args = parser.parse_args()
result = validate(args.gpu_case_dir, args.bundle_dir, args.out_json)
print(json.dumps({'pass_all': result['pass_all'], 'steps': result['steps']}), flush=True)
if not result['pass_all']:
raise SystemExit(1)