Digital Drosophila Archive
330 lines
19 KiB
Python
330 lines
19 KiB
Python
"""Non-feedback float64 attribution of the frozen, recorded graded APL trace.
|
|
|
|
No files, Torch, GPU, neuronal state, or counterfactual simulation are used.
|
|
The reference is mathematical source accounting, NOT bitwise CUDA replay.
|
|
Full bundle output is generated even for graph fixtures; the caller must not
|
|
interpret inhibition for absent fixture KCs as inhibition actually applied.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import math
|
|
import numpy as np
|
|
|
|
|
|
N = 166700
|
|
GROUPS = ("design_ALPN", "non_design_ALPN", "KC")
|
|
DECAY32 = np.float32(math.exp(-.1 / 20.))
|
|
GAIN32 = np.float32(.01)
|
|
|
|
|
|
def _require(ok, message):
|
|
if not ok:
|
|
raise ValueError(message)
|
|
|
|
|
|
def _ids(value, name):
|
|
value = np.asarray(value)
|
|
_require(value.dtype == np.int32 and value.ndim == 1, name + " must be an int32 vector")
|
|
_require(np.all((value >= 0) & (value < N)) and np.all(value[1:] > value[:-1]),
|
|
name + " must contain sorted unique global IDs")
|
|
return value
|
|
|
|
|
|
def _positive_array(value, shape, name):
|
|
value = np.asarray(value)
|
|
_require(value.dtype == np.float32 and value.shape == shape, name + " has the wrong shape/dtype")
|
|
_require(np.isfinite(value).all() and np.all(value >= 0), name + " must be finite and nonnegative")
|
|
return value
|
|
|
|
|
|
def _up(value):
|
|
return np.nextafter(np.asarray(value, np.float64), np.inf)
|
|
|
|
|
|
def _add(a, b):
|
|
return _up(np.asarray(a, np.float64) + np.asarray(b, np.float64))
|
|
|
|
|
|
def _mul(a, b):
|
|
return _up(np.asarray(a, np.float64) * np.asarray(b, np.float64))
|
|
|
|
|
|
def _gamma(count, unit_roundoff):
|
|
_require(count * unit_roundoff < 1, "Undefined gamma bound")
|
|
numerator = _mul(count, unit_roundoff)
|
|
denominator = np.nextafter(1. - numerator, -np.inf)
|
|
return _up(numerator / denominator)
|
|
|
|
|
|
def _sum_upper(values):
|
|
"""Outward-rounded sum on axis zero; inputs are nonnegative."""
|
|
result = np.zeros(values.shape[1:], np.float64)
|
|
for value in values:
|
|
result = _add(result, value)
|
|
return result
|
|
|
|
|
|
def _sum_groups(values):
|
|
# Explicit separate additions; never rely on a backend reduction over G.
|
|
return np.add(np.add(values[0], values[1]), values[2])
|
|
|
|
|
|
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(value, tiny):
|
|
return bool(np.isfinite(value).all() and np.all((value == 0) | (value >= tiny)))
|
|
|
|
|
|
def _metric():
|
|
return dict(elements=0, sum_abs_error=0., max_abs_error=0., max_error_bound=0.,
|
|
max_bound_utilization=0., bound_violations=0, zero_reference_nonzero_actual=0)
|
|
|
|
|
|
def _measure(metric, actual, reference, bound):
|
|
_require(np.isfinite(reference).all() and np.isfinite(bound).all(), "Nonfinite reference/bound")
|
|
error = np.abs(actual.astype(np.float64) - reference)
|
|
upper = np.where(error == 0, 0., _up(error))
|
|
impossible_zero = (reference == 0) & (actual != 0)
|
|
bad = (upper > bound) | impossible_zero
|
|
usage = np.divide(upper, bound, out=np.zeros_like(upper), where=bound > 0)
|
|
metric['elements'] += error.size
|
|
metric['sum_abs_error'] += float(error.sum(dtype=np.float64))
|
|
metric['max_abs_error'] = max(metric['max_abs_error'], float(error.max(initial=0)))
|
|
metric['max_error_bound'] = max(metric['max_error_bound'], float(bound.max(initial=0)))
|
|
metric['max_bound_utilization'] = max(metric['max_bound_utilization'], float(usage.max(initial=0)))
|
|
metric['bound_violations'] += int(bad.sum())
|
|
metric['zero_reference_nonzero_actual'] += int(impossible_zero.sum())
|
|
|
|
|
|
def decompose_apl_sources(I32, O32, pn_globals, kc_globals, design_pn_globals,
|
|
event_step, event_node, drive32, state32, h32, *,
|
|
decay32, gain32, standard_fp32, chunk_size=512, callback=None,
|
|
first_kc_tick=None):
|
|
"""Attribute APL input/output to three exact anatomical source groups.
|
|
|
|
I32 is the full (2,166700) normalized float32 input matrix; O32 is (2,K),
|
|
with columns in sorted kc_globals order. Traces are drive/state float32
|
|
(T,2), inhibition float32 (T,K), including ALL bundle KC output columns.
|
|
Actual event vectors are int32 and sorted/unique by tick/global neuron.
|
|
Design IDs must be a subset of pn_globals; fixtures' design definition is
|
|
supplied explicitly by the caller, never inferred from response outcome.
|
|
|
|
standard_fp32 must be True only after historical metadata verifies TF32 is
|
|
disabled and matmul precision is 'highest'. This function cannot inspect a
|
|
GPU environment. d/q must have the exact frozen float32 constant bits.
|
|
|
|
Return summary (pure JSON values), tick_rows (T*3*2 JSON dictionaries),
|
|
first_kc_tick (int32[K], -1 for never fired), and per_kc cumulative arrays.
|
|
Optional first_kc_tick is checked against actual events, never overrides
|
|
them. No T*3*K history is retained; source cumulative storage is O(3*K).
|
|
callback(tick, dict) optionally receives copies of drive/state group arrays
|
|
(3,2), plus their totals, never a full inhibition history.
|
|
|
|
Early means strictly tick < observed first KC spike; same-tick APL output
|
|
follows threshold and cannot cause that spike. Never-fired KC output and
|
|
output strictly before the first spike of ANY KC are reported separately.
|
|
"""
|
|
_require(standard_fp32 is True, "Historical standard FP32 arithmetic must be verified by caller")
|
|
_require(isinstance(decay32, np.float32) and isinstance(gain32, np.float32) and
|
|
decay32.view(np.uint32) == DECAY32.view(np.uint32) and
|
|
gain32.view(np.uint32) == GAIN32.view(np.uint32), "Frozen float32 decay/gain bits changed")
|
|
_require(isinstance(chunk_size, (int, np.integer)) and not isinstance(chunk_size, (bool, np.bool_)) and
|
|
chunk_size > 0, "chunk_size must be a positive integer")
|
|
_require(callback is None or callable(callback), "callback must be callable")
|
|
pn = _ids(pn_globals, 'pn_globals')
|
|
kc = _ids(kc_globals, 'kc_globals')
|
|
design = _ids(design_pn_globals, 'design_pn_globals')
|
|
_require(len(kc) > 0 and not np.intersect1d(pn, kc).size, "PN/KC identities overlap or KC set is empty")
|
|
_require(np.isin(design, pn).all(), "Design PN is outside the anatomical PN population")
|
|
k = len(kc)
|
|
I = _positive_array(I32, (2, N), 'I32').astype(np.float64)
|
|
O = _positive_array(O32, (2, k), 'O32').astype(np.float64)
|
|
state32 = np.asarray(state32)
|
|
_require(state32.ndim == 2 and len(state32) > 0, "Empty/invalid state trace")
|
|
steps = len(state32)
|
|
state32 = _positive_array(state32, (steps, 2), 'state32')
|
|
drive32 = _positive_array(drive32, (steps, 2), 'drive32')
|
|
h32 = np.asarray(h32)
|
|
_require(h32.dtype == np.float32 and h32.shape == (steps, k), "h32 has wrong shape/dtype")
|
|
# Validate without allocating a dense T*K temporary boolean array.
|
|
for start in range(0, k, chunk_size):
|
|
value = h32[:, start:start + chunk_size]
|
|
_require(np.isfinite(value).all() and np.all(value >= 0), "h32 must be finite and nonnegative")
|
|
labels = np.full(N, -1, np.int8)
|
|
labels[pn], labels[design], labels[kc] = 1, 0, 2
|
|
support = np.any(I != 0, axis=0)
|
|
unknown = np.flatnonzero(support & (labels < 0))
|
|
_require(not len(unknown), "Nonzero APL input outside design/non-design PN/KC: " + str(unknown[:10].tolist()))
|
|
ticks, nodes = np.asarray(event_step), np.asarray(event_node)
|
|
_require(ticks.dtype == nodes.dtype == np.int32 and ticks.ndim == nodes.ndim == 1 and
|
|
ticks.shape == nodes.shape, "Actual events must be matching int32 vectors")
|
|
_require(np.all((ticks >= 0) & (ticks < steps)) and np.all((nodes >= 0) & (nodes < N)), "Event out of range")
|
|
_require(np.all((ticks[1:] > ticks[:-1]) | ((ticks[1:] == ticks[:-1]) & (nodes[1:] > nodes[:-1]))),
|
|
"Actual events must be sorted and unique")
|
|
ptr = np.searchsorted(ticks, np.arange(steps + 1))
|
|
lookup = np.full(N, -1, np.int32)
|
|
lookup[kc] = np.arange(k, dtype=np.int32)
|
|
first = np.full(k, steps, np.int32)
|
|
is_kc = lookup[nodes] >= 0
|
|
np.minimum.at(first, lookup[nodes[is_kc]], ticks[is_kc])
|
|
ever = first < steps
|
|
first_any = int(first.min()) if ever.any() else None
|
|
first[~ever] = -1
|
|
if first_kc_tick is not None:
|
|
supplied = np.asarray(first_kc_tick)
|
|
_require(supplied.dtype == np.int32 and supplied.shape == (k,) and np.array_equal(supplied, first),
|
|
"Supplied first_kc_tick disagrees with observed actual KC events")
|
|
d, q = float(decay32), float(gain32)
|
|
u32, u64 = 2.**-24, 2.**-53
|
|
gd32, gd64 = _gamma(N - 1, u32), _gamma(N - 1, u64)
|
|
gs32, gs64 = _gamma(2, u32), _gamma(2, u64)
|
|
gh32, gh64, gg64 = _gamma(5, u32), _gamma(5, u64), _gamma(2, u64)
|
|
group_state = np.zeros((3, 2), np.float64)
|
|
e64 = np.zeros_like(group_state)
|
|
e32 = np.zeros(2, np.float64)
|
|
previous32 = np.zeros(2, np.float64)
|
|
statistics = {name: _metric() for name in ('drive', 'state', 'inhibition')}
|
|
tiny32, tiny64 = np.finfo(np.float32).tiny, np.finfo(np.float64).tiny
|
|
normal = _normal(I, tiny32) and _normal(O, tiny32)
|
|
sums = {name: dict(actual_spikes=0, input_supported_spikes=0, first_actual_spike_tick=None,
|
|
first_nonzero_drive_tick=None, cumulative_drive_by_apl=np.zeros(2),
|
|
cumulative_state_by_apl=np.zeros(2), h_sum_mV=0., h_before_first_KC_sum_mV=0.,
|
|
h_never_fired_KC_sum_mV=0., h_before_first_any_KC_sum_mV=0.) for name in GROUPS}
|
|
rows = []
|
|
per_kc = dict(neuron_idx=kc.copy(), first_spike_tick=first.copy(), ever_fired=ever.copy())
|
|
for name in GROUPS:
|
|
for suffix in ('pre_first_h_sum_mV', 'never_fired_h_sum_mV', 'h_sum_mV'):
|
|
per_kc[name + '_' + suffix] = np.zeros(k, np.float64)
|
|
O_total, O_never = O.sum(axis=1), O[:, ~ever].sum(axis=1)
|
|
for tick in range(steps):
|
|
fired = nodes[ptr[tick]:ptr[tick + 1]]
|
|
group_drive = np.zeros_like(group_state)
|
|
spike_counts = np.zeros(3, np.int64)
|
|
for group, name in enumerate(GROUPS):
|
|
selected = fired[labels[fired] == group]
|
|
spike_counts[group] = len(selected)
|
|
group_drive[group] = I[:, selected].sum(axis=1, dtype=np.float64)
|
|
acc = sums[name]
|
|
acc['actual_spikes'] += len(selected)
|
|
acc['input_supported_spikes'] += int(support[selected].sum())
|
|
if len(selected) and acc['first_actual_spike_tick'] is None:
|
|
acc['first_actual_spike_tick'] = tick
|
|
if np.any(group_drive[group] > 0) and acc['first_nonzero_drive_tick'] is None:
|
|
acc['first_nonzero_drive_tick'] = tick
|
|
# Each source-group positive reduction independently bounds its exact
|
|
# real sum. Outward sum across groups bounds the common exact drive.
|
|
Dupper = _up(group_drive / np.nextafter(1. - gd64, -np.inf))
|
|
total_Dupper = _sum_upper(Dupper)
|
|
drive_ref = _sum_groups(group_drive)
|
|
drive_bound = _add(_mul(gd32, total_Dupper),
|
|
_add(_sum_upper(_mul(gd64, Dupper)), _mul(gg64, _sum_upper(group_drive))))
|
|
e32 = _add(_mul(d, e32), _add(_mul(gs32, _add(_mul(d, previous32), total_Dupper)),
|
|
_mul(1. + u32, _mul(gd32, total_Dupper))))
|
|
e64 = _add(_mul(d, e64), _add(_mul(gs64, _add(_mul(d, group_state), Dupper)),
|
|
_mul(1. + u64, _mul(gd64, Dupper))))
|
|
normal &= _normal(d * previous32, tiny32) and _normal(d * group_state, tiny64)
|
|
group_state = np.add(np.multiply(d, group_state), group_drive)
|
|
state_ref = _sum_groups(group_state)
|
|
state_bound = _add(e32, _add(_sum_upper(e64), _mul(gg64, _sum_upper(group_state))))
|
|
observed_state = state32[tick].astype(np.float64)
|
|
normal &= _normal(drive32[tick], tiny32) and _normal(observed_state, tiny32)
|
|
_measure(statistics['drive'], drive32[tick], drive_ref, drive_bound)
|
|
_measure(statistics['state'], state32[tick], state_ref, state_bound)
|
|
for start in range(0, k, chunk_size):
|
|
stop = min(k, start + chunk_size)
|
|
matrix = O[:, start:stop]
|
|
h_group = np.empty((3, stop - start), np.float64)
|
|
h_error = np.empty_like(h_group)
|
|
for group in range(3):
|
|
products = group_state[group, :, None] * matrix
|
|
dot = np.add(products[0], products[1])
|
|
h_group[group] = np.multiply(q, dot)
|
|
normal &= _normal(products, tiny64) and _normal(h_group[group], tiny64)
|
|
h_error[group] = _add(_mul(q, _positive_dot_upper(e64[group], matrix)),
|
|
_mul(gh64, _mul(q, _positive_dot_upper(group_state[group], matrix))))
|
|
name = GROUPS[group]
|
|
per_kc[name + '_h_sum_mV'][start:stop] += h_group[group]
|
|
early_mask = ever[start:stop] & (tick < first[start:stop])
|
|
per_kc[name + '_pre_first_h_sum_mV'][start:stop] += h_group[group] * early_mask
|
|
per_kc[name + '_never_fired_h_sum_mV'][start:stop] += h_group[group] * ~ever[start:stop]
|
|
h_ref = _sum_groups(h_group)
|
|
h32_error = _add(_mul(q, _positive_dot_upper(e32, matrix)),
|
|
_mul(gh32, _mul(q, _positive_dot_upper(observed_state, matrix))))
|
|
bound = _add(h32_error, _add(_sum_upper(h_error), _mul(gg64, _sum_upper(h_group))))
|
|
actual_products = observed_state[:, None] * matrix
|
|
actual_dot = np.add(actual_products[0], actual_products[1])
|
|
normal &= _normal(actual_products, tiny32) and _normal(q * actual_dot, tiny32)
|
|
normal &= _normal(h32[tick, start:stop], tiny32)
|
|
_measure(statistics['inhibition'], h32[tick, start:stop], h_ref, bound)
|
|
# Aggregate output per source and per APL, never T*G*K. This is linear
|
|
# f64 measurement; grouping/reduction is not claimed bitwise identical.
|
|
early = ever & (tick < first)
|
|
O_early = O[:, early].sum(axis=1)
|
|
before_any = first_any is not None and tick < first_any
|
|
for group, name in enumerate(GROUPS):
|
|
acc = sums[name]
|
|
acc['cumulative_drive_by_apl'] += group_drive[group]
|
|
acc['cumulative_state_by_apl'] += group_state[group]
|
|
for apl in range(2):
|
|
amplitude = q * group_state[group, apl]
|
|
h_sum = float(amplitude * O_total[apl])
|
|
h_early = float(amplitude * O_early[apl])
|
|
h_never = float(amplitude * O_never[apl])
|
|
h_any = h_sum if before_any else 0.
|
|
row = dict(tick=tick, group=name, apl_index=apl, source_spikes=int(spike_counts[group]),
|
|
drive64=float(group_drive[group, apl]), state64=float(group_state[group, apl]),
|
|
h_sum_mV=h_sum, h_before_first_KC_sum_mV=h_early,
|
|
h_never_fired_KC_sum_mV=h_never, h_before_first_any_KC_sum_mV=h_any)
|
|
rows.append(row)
|
|
for key in ('h_sum_mV', 'h_before_first_KC_sum_mV', 'h_never_fired_KC_sum_mV', 'h_before_first_any_KC_sum_mV'):
|
|
acc[key] += row[key]
|
|
if callback is not None:
|
|
callback(tick, dict(drive_by_group=group_drive.copy(), state_by_group=group_state.copy(),
|
|
drive_total=drive_ref.copy(), state_total=state_ref.copy()))
|
|
previous32 = observed_state
|
|
total_h = sum(value['h_sum_mV'] for value in sums.values())
|
|
total_drive = sum(float(value['cumulative_drive_by_apl'].sum()) for value in sums.values())
|
|
for acc in sums.values():
|
|
acc['generated_h_fraction'] = None if total_h == 0 else acc['h_sum_mV'] / total_h
|
|
acc['drive_fraction'] = None if total_drive == 0 else float(acc['cumulative_drive_by_apl'].sum()) / total_drive
|
|
for key in ('cumulative_drive_by_apl', 'cumulative_state_by_apl'):
|
|
acc[key] = acc[key].tolist()
|
|
for metric in statistics.values():
|
|
metric['mean_abs_error'] = metric.pop('sum_abs_error') / metric['elements']
|
|
finite_metrics = all(np.isfinite([row[key] for row in rows]).all() for key in
|
|
('drive64', 'state64', 'h_sum_mV', 'h_before_first_KC_sum_mV', 'h_never_fired_KC_sum_mV', 'h_before_first_any_KC_sum_mV'))
|
|
_require(finite_metrics, "Nonfinite aggregate accounting")
|
|
_require(all(np.isfinite(value).all() for value in per_kc.values()), "Nonfinite per-KC accounting")
|
|
passed = bool(normal and all(metric['bound_violations'] == 0 for metric in statistics.values()))
|
|
summary = dict(pass_all=passed, steps=steps, neurons=N, apl=2, bundle_KC=k,
|
|
actual_events=len(nodes), groups=sums, group_order=list(GROUPS),
|
|
input_support_partition_exact=True, unknown_nonzero_input_columns=0,
|
|
source_population_counts={name: int(np.count_nonzero(labels == index)) for index, name in enumerate(GROUPS)},
|
|
input_supported_population_counts={name: int(np.count_nonzero(support & (labels == index))) for index, name in enumerate(GROUPS)},
|
|
first_any_KC_spike_tick=first_any, ever_fired_KC=int(ever.sum()), never_fired_KC=int((~ever).sum()),
|
|
decay_f32_bits=f'{int(decay32.view(np.uint32)):08x}', gain_f32_bits=f'{int(gain32.view(np.uint32)):08x}',
|
|
comparison=statistics, theoretical_bound=dict(applicable=bool(normal),
|
|
normal_or_zero_arithmetic=bool(normal), historical_standard_fp32_verified_by_caller=True,
|
|
gamma_definition='gamma_n=n*u/(1-n*u); u32=2^-24,u64=2^-53; outward rounded',
|
|
drive_terms_bound=N-1, state_operation_bound=2, output_operation_bound=5, group_sum_operation_bound=2,
|
|
state_recurrence='E_p<=d*E_p_prev+gamma_2,p*(d*abs(s_p_prev)+D_upper)+(1+u_p)*gamma_(N-1),p*D_upper',
|
|
grouped_f64_bound='Sum source-group error bounds, plus gamma_2,64*upward_sum(group_values) for per-tick compared drive/state/h totals only',
|
|
output_bound='q*upward_dot(E_p,O)+gamma_5,p*q*upward_dot(s_p,O); f32 error plus all f64 group/output/sum errors',
|
|
tolerance_fitted_to_observed_error=False,
|
|
scope='Per-tick summed-source drive/state/h trace comparisons only; not cumulative source fractions'),
|
|
early_definition='strict tick<each observed first KC spike; never-fired columns separate; same-tick end APL excluded',
|
|
value_units=dict(drive64='legacy normalized-contact spike drive', state64='legacy normalized-drive accumulator',
|
|
h='legacy mV per state-update step; sums span output KC and recorded ticks, not instantaneous conductance'),
|
|
generated_full_bundle_trace=True, fixture_absent_KC_output_is_not_claimed_applied=True,
|
|
row_source_spikes_repeated_for_each_APL=True, nonfeedback=True, CUDA_float32_bitwise_replay=False,
|
|
cumulative_and_fraction_metrics_descriptive_float64=True,
|
|
cumulative_source_fraction_error_bounds_verified=False,
|
|
source_majority_robustness_proven=False,
|
|
independent_brain_oracle=False, lesion_or_counterfactual=False, parameters_changed=False)
|
|
return dict(summary=summary, tick_rows=rows, first_kc_tick=first, per_kc=per_kc)
|