Digital Drosophila Archive
198 lines
12 KiB
Python
198 lines
12 KiB
Python
"""Non-feedback drive ledger. None of these summaries affect the replay."""
|
|
import numpy as np
|
|
|
|
|
|
def bit_equal(a, b):
|
|
a, b = np.asarray(a), np.asarray(b)
|
|
return 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 fraction(numerator, denominator):
|
|
return None if denominator == 0 else float(numerator / denominator)
|
|
|
|
|
|
class DriveAccounting:
|
|
def __init__(self, target_globals, incoming, raw_contacts, source_groups, group_names,
|
|
expected_events, coeff, expected_counts, histories=None):
|
|
self.targets = target_globals
|
|
self.incoming = incoming
|
|
self.raw = raw_contacts
|
|
self.groups = source_groups
|
|
self.names = group_names
|
|
self.coeff = coeff
|
|
self.expected_counts = expected_counts
|
|
self.active = expected_counts > 0
|
|
self.histories = histories or {}
|
|
self.expected = {}
|
|
for tick in np.unique(expected_events['step']):
|
|
self.expected[int(tick)] = expected_events['neuron_idx'][expected_events['step'] == tick]
|
|
n = len(target_globals)
|
|
self.n = n
|
|
self.first = np.full(n, -1, np.int32)
|
|
self.tick_rows = []
|
|
self.group_stats = {name: {key: np.zeros(n, np.float64 if 'SI' in key else np.int64)
|
|
for key in ('due_edges', 'accepted_edges', 'rejected_edges', 'due_contacts',
|
|
'accepted_contacts', 'rejected_contacts', 'accepted_exc_SI',
|
|
'accepted_inh_abs_SI', 'rejected_exc_SI', 'rejected_inh_abs_SI')}
|
|
for name in group_names}
|
|
self.metrics = {key: np.zeros(n, np.float64) for key in (
|
|
'APL_write_SI', 'APL_write_ref_SI', 'APL_write_available_SI', 'APL_pre_first_spike_SI',
|
|
'APL_same_tick_first_spike_SI', 'negative_g_ticks', 'refractory_ticks',
|
|
'g_end_time_integral_V_ms', 'V_end_time_integral_V_ms', 'accepted_exc_before_first_SI',
|
|
'accepted_inh_before_first_abs_SI', 'APL_zeroed_due_edges', 'APL_zeroed_due_contacts')}
|
|
self.min_g = np.full(n, np.inf)
|
|
self.max_v = np.full(n, -np.inf)
|
|
self.ledger = {key: np.zeros(n, np.longdouble) for key in ('decay', 'reset', 'arrivals', 'APL')}
|
|
self.ledger_abs = np.zeros(n, np.longdouble)
|
|
self.arrival_rounding = np.zeros(n, np.longdouble)
|
|
# Source/sign lanes are a conditional accounting, never fed to V/g.
|
|
self.shadow_names = [name + '_' + sign for name in group_names for sign in ('exc', 'inh')] + ['graded_APL']
|
|
self.shadow = np.zeros((len(self.shadow_names), n), np.float64)
|
|
self.shadow_integral = np.zeros_like(self.shadow)
|
|
self.initial_shadow_integral = np.zeros(n, np.float64)
|
|
self.shadow_max_residual = 0.0
|
|
self.callback_count = 0
|
|
self.initial_g = None
|
|
|
|
def __call__(self, tick, s):
|
|
fired_nodes = self.targets[s['fired']]
|
|
expected = self.expected.get(tick, np.empty(0, np.int32))
|
|
if not bit_equal(fired_nodes, expected):
|
|
raise ValueError(f'KC spike divergence at tick {tick}: computed={fired_nodes.tolist()} expected={expected.tolist()}')
|
|
if self.initial_g is None:
|
|
self.initial_g = s['old_g'].copy()
|
|
# Fixture initial g is a separate lane, not attributed to any input class.
|
|
self.initial_shadow = s['old_g'].copy()
|
|
for name, history in self.histories.items():
|
|
key = {'pre_v': 'v_pre_apl', 'pre_g': 'g_pre_apl', 'pre_ref': 'ref_pre_apl',
|
|
'v': 'v', 'g': 'g', 'ref': 'ref'}[name]
|
|
if not bit_equal(history[tick], s[key]):
|
|
raise ValueError(f'Fixture dense {name} divergence at tick {tick}')
|
|
before_first = self.first < 0
|
|
first_now = before_first & s['fired']
|
|
self.first[first_now] = tick
|
|
self.max_v = np.maximum(self.max_v, s['v_threshold'])
|
|
self.min_g = np.minimum(self.min_g, s['g'])
|
|
ref_mask = s['ref'] > 0
|
|
delta = s['delta_si']
|
|
self.metrics['APL_write_SI'] += delta
|
|
self.metrics['APL_write_ref_SI'] += delta * ref_mask
|
|
self.metrics['APL_write_available_SI'] += delta * ~ref_mask
|
|
# Strictly prior ticks: arrival/APL at first-spike tick cannot cause that spike.
|
|
self.metrics['APL_pre_first_spike_SI'] += delta * (before_first & ~s['fired'])
|
|
self.metrics['APL_same_tick_first_spike_SI'] += delta * first_now
|
|
self.metrics['negative_g_ticks'] += s['g'] < 0
|
|
self.metrics['refractory_ticks'] += ref_mask
|
|
self.metrics['g_end_time_integral_V_ms'] += s['g'] * .1
|
|
self.metrics['V_end_time_integral_V_ms'] += s['v'] * .1
|
|
stages = {'decay': ('g_decayed', 'old_g'), 'reset': ('g_after_reset', 'g_decayed'),
|
|
'arrivals': ('g_pre_apl', 'g_after_reset'), 'APL': ('g', 'g_pre_apl')}
|
|
for label, (after, before) in stages.items():
|
|
diff = s[after].astype(np.longdouble) - s[before].astype(np.longdouble)
|
|
self.ledger[label] += diff
|
|
self.ledger_abs += np.abs(diff)
|
|
due = s['due_edge_indices']
|
|
dst = self.incoming.dst_local[due]
|
|
weights = self.incoming.weights[due]
|
|
groups = self.groups[due]
|
|
accepted = s['arrival_accepted']
|
|
positive, negative = weights > 0, weights < 0
|
|
nominal = np.zeros(self.n, np.longdouble)
|
|
np.add.at(nominal, dst[accepted], weights[accepted].astype(np.longdouble))
|
|
actual = s['g_pre_apl'].astype(np.longdouble) - s['g_after_reset'].astype(np.longdouble)
|
|
self.arrival_rounding += actual - nominal
|
|
self.shadow[:, s['available']] *= self.coeff['gd']
|
|
self.initial_shadow[s['available']] *= self.coeff['gd']
|
|
self.shadow[:, s['fired']] = 0.0
|
|
self.initial_shadow[s['fired']] = 0.0
|
|
row = dict(tick=tick, time_ms=tick * .1, KC_spikes=int(s['fired'].sum()),
|
|
KC_negative_g=int(np.count_nonzero(s['g'] < 0)), KC_refractory=int(ref_mask.sum()),
|
|
APL_write_SI_sum=float(delta.sum()), APL_ref_write_SI_sum=float(delta[ref_mask].sum()))
|
|
for gi, name in enumerate(self.names):
|
|
mask = groups == gi
|
|
acc, rej = mask & accepted, mask & ~accepted
|
|
stat = self.group_stats[name]
|
|
for stem, selected in (('due', mask), ('accepted', acc), ('rejected', rej)):
|
|
np.add.at(stat[stem + '_edges'], dst[selected], 1)
|
|
np.add.at(stat[stem + '_contacts'], dst[selected], self.raw[due[selected]])
|
|
for stem, selected in (('accepted', acc), ('rejected', rej)):
|
|
for suffix, signmask, sign in (('exc_SI', positive, 1), ('inh_abs_SI', negative, -1)):
|
|
pick = selected & signmask
|
|
np.add.at(stat[stem + '_' + suffix], dst[pick], sign * weights[pick])
|
|
row[name + '_' + stem + '_' + suffix] = float((sign * weights[pick]).sum())
|
|
for sign_index, signmask in ((0, positive), (1, negative)):
|
|
selected = acc & signmask
|
|
np.add.at(self.shadow[2 * gi + sign_index], dst[selected], weights[selected])
|
|
self.shadow[-1] -= delta
|
|
self.shadow_integral += self.shadow * .1
|
|
self.initial_shadow_integral += self.initial_shadow * .1
|
|
self.shadow_max_residual = max(self.shadow_max_residual,
|
|
float(np.max(np.abs(self.shadow.sum(axis=0) + self.initial_shadow - s['g']))))
|
|
no_first_yet = before_first[dst] & ~s['fired'][dst]
|
|
for metric, selected, sign in (
|
|
('accepted_exc_before_first_SI', accepted & positive & no_first_yet, 1),
|
|
('accepted_inh_before_first_abs_SI', accepted & negative & no_first_yet, -1)):
|
|
np.add.at(self.metrics[metric], dst[selected], sign * weights[selected])
|
|
apl_id = self.names.index('APL')
|
|
zeroed = (groups == apl_id) & (weights == 0)
|
|
np.add.at(self.metrics['APL_zeroed_due_edges'], dst[zeroed], 1)
|
|
np.add.at(self.metrics['APL_zeroed_due_contacts'], dst[zeroed], self.raw[due[zeroed]])
|
|
for stratum, mask in (('ever_active', self.active), ('never_active', ~self.active)):
|
|
row[stratum + '_mean_g_V'] = None if not mask.any() else float(s['g'][mask].mean())
|
|
row[stratum + '_mean_V_V'] = None if not mask.any() else float(s['v'][mask].mean())
|
|
self.tick_rows.append(row)
|
|
self.callback_count += 1
|
|
|
|
def finish(self, final_g):
|
|
reconstructed = self.initial_g.astype(np.longdouble) + sum(self.ledger.values())
|
|
residual = reconstructed - final_g.astype(np.longdouble)
|
|
bound = (64 * np.finfo(np.longdouble).eps * (self.callback_count + 1) *
|
|
(self.ledger_abs + np.abs(self.initial_g.astype(np.longdouble)) + np.abs(final_g)))
|
|
assert np.all(np.abs(residual) <= bound), 'g ledger residual exceeds stated arithmetic bound'
|
|
data = dict(neuron_idx=self.targets, spike_count=self.expected_counts, first_spike_tick=self.first,
|
|
minimum_end_g_V=self.min_g, maximum_prethreshold_V=self.max_v,
|
|
best_threshold_margin_mV=(self.coeff['threshold'] - self.max_v) * 1000,
|
|
g_ledger_residual_V=residual.astype(np.float64),
|
|
g_ledger_bound_V=bound.astype(np.float64), arrival_rounding_V=self.arrival_rounding.astype(np.float64),
|
|
initial_g_state_area_V_ms=self.initial_shadow_integral,
|
|
source_g_area_residual_V_ms=(self.shadow_integral.sum(axis=0) + self.initial_shadow_integral -
|
|
self.metrics['g_end_time_integral_V_ms']),
|
|
**self.metrics)
|
|
for name, values in self.group_stats.items():
|
|
for metric, values2 in values.items():
|
|
data[name + '_' + metric] = values2
|
|
for name, value in self.ledger.items():
|
|
data['actual_g_delta_' + name + '_V'] = value.astype(np.float64)
|
|
for index, name in enumerate(self.shadow_names):
|
|
data['source_g_area_' + name + '_V_ms'] = self.shadow_integral[index]
|
|
for key, value in data.items():
|
|
assert np.isfinite(value).all(), 'Nonfinite accounting: ' + key
|
|
aggregates = {}
|
|
for label, mask in (('all', np.ones(self.n, bool)), ('ever_active', self.active), ('never_active', ~self.active)):
|
|
groups = {name: {key: float(value[mask].sum()) if 'SI' in key else int(value[mask].sum())
|
|
for key, value in metrics.items()} for name, metrics in self.group_stats.items()}
|
|
exc = sum((self.group_stats[n]['accepted_exc_SI'] for n in self.names), np.zeros(self.n))
|
|
alp = self.group_stats['ALPN']['accepted_exc_SI']
|
|
aggregates[label] = dict(neurons=int(mask.sum()), source_groups=groups,
|
|
no_accepted_exc_neurons=int(np.count_nonzero(mask & (exc == 0))),
|
|
no_accepted_ALPN_exc_neurons=int(np.count_nonzero(mask & (alp == 0))),
|
|
has_exc_but_never_active=int(np.count_nonzero(mask & (exc > 0) & ~self.active)),
|
|
metrics_sum={k: float(v[mask].sum()) for k, v in self.metrics.items()},
|
|
source_g_area_V_ms={name: float(self.shadow_integral[i, mask].sum()) for i, name in enumerate(self.shadow_names)},
|
|
initial_g_state_area_V_ms=float(self.initial_shadow_integral[mask].sum()),
|
|
source_g_area_residual_V_ms=float(data['source_g_area_residual_V_ms'][mask].sum()),
|
|
maximum_prethreshold_V=None if not mask.any() else float(self.max_v[mask].max()),
|
|
minimum_end_g_V=None if not mask.any() else float(self.min_g[mask].min()))
|
|
return data, dict(strata=aggregates, ledger_max_abs_error_V=float(np.max(np.abs(residual))),
|
|
ledger_max_bound_V=float(bound.max()), ledger_roundoff_bound_pass=True,
|
|
ledger_accumulator_mantissa_bits=int(np.finfo(np.longdouble).nmant),
|
|
ledger_bound_formula='64 * eps(longdouble) * (ticks+1) * (sum_abs_actual_stage_deltas + abs(initial_g) + abs(final_g)) per KC',
|
|
source_shadow_max_residual_V=self.shadow_max_residual,
|
|
source_area_max_residual_V_ms=float(np.max(np.abs(data['source_g_area_residual_V_ms']))),
|
|
source_shadow_meaning='Signed source/sign g bookkeeping under observed availability/reset only; non-feedback, not counterfactual',
|
|
cumulative_increment_units='SI volts; sums of individual increments, not instantaneous conductance',
|
|
source_state_area_units='V*ms; post-APL source-state left samples times dt .1ms',
|
|
activity_strata_outcome_defined=True, APL_cause_of_CAL_failure_proven=False,
|
|
full_CAL_all_tick_subset_state_independently_verified=False)
|