Files

151 lines
15 KiB
Python
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
"""Post-run descriptive audit only; no selection, simulation or model mutation."""
from pathlib import Path
import json,hashlib,importlib.util
import numpy as np
import pandas as pd
ROOT=Path('/data/digital-drosophila')
exp=ROOT/'experiments/phase4z-pn-feedback-audit'
assert not (exp/'checksums.sha256').exists(),'Sealed output is immutable'
(exp/'audit/descriptive-analysis-provenance.json').write_text(json.dumps({
'source':str(Path(__file__).resolve()),'sha256':hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
'post_run_descriptive_analysis':True,'not_part_of_simulation_preregistration':True,
'new_simulations_or_selection':False},indent=2)+'\n')
meta=pd.read_parquet(exp/'audit/incoming-source-metadata.parquet').set_index('neuron_idx')
watch=pd.read_parquet(exp/'audit/watched-metadata.parquet').set_index('neuron_idx')
d=pd.read_parquet(exp/'drive/edge-accounting.parquet');d=d[d.period=='late'].copy()
d['condition']=d.trial.str.split('-',n=2).str[-1]
rep=pd.read_parquet(exp/'replay/results.parquet');rep['condition']=rep.trial.str.split('-',n=2).str[-1]
path=pd.read_parquet(exp/'ancestry/event-supported-edges.parquet')
q=pd.read_parquet(exp/'input-audit/paired-summary.parquet')
current=json.loads((exp/'protocol/input-hashes.json').read_text())
prior=json.loads((ROOT/'experiments/phase4y-dpm-hypothesis/protocol/input-hashes.json').read_text())
old={x['odor']:x for x in prior if x['seed']==131101}
input_checks=[{'odor_index':x['odor_index'],
'forcing_hash_exact':x['forcing']==old[x['odor_index']]['forcing_hash'],
'targets_rates_hash_exact':x['targets_rates']==old[x['odor_index']]['targets_rates_hash']} for x in current]
assert all(x['forcing_hash_exact'] and x['targets_rates_hash_exact'] for x in input_checks)
(exp/'regression/input-hashes-vs-Y.json').write_text(json.dumps({'pass':True,'checks':input_checks},indent=2)+'\n')
rows=[]
for (trial,node),g in d.groupby(['trial','target_neuron_idx']):
den=g.positive_arrival_sum.sum()
for role,z in g.groupby('source_role'):
rows.append({'trial':trial,'target_neuron_idx':node,'source_role':role,
'positive_fraction':float(z.positive_arrival_sum.sum()/den) if den else None,
'positive_arrival_sum':float(z.positive_arrival_sum.sum()),'negative_arrival_sum':float(z.negative_arrival_sum.sum())})
pd.DataFrame(rows).to_parquet(exp/'drive/late-population-accounting.parquet',index=False)
top=[];timing=[];transitions=[]
for trial,g in d[d.target_neuron_idx==82312].groupby('trial'):
a=g[g.source_role=='ALPN'].sort_values(['positive_arrival_sum','edge_id'],ascending=[False,True]).iloc[0]
top.append({'trial':trial,'source_neuron_idx':int(a.source_neuron_idx),'body_id':int(meta.loc[a.source_neuron_idx,'body_id']),
'type':str(meta.loc[a.source_neuron_idx,'annotation_type']),'edge_id':int(a.edge_id),
'positive_fraction':float(a.positive_arrival_sum/g.positive_arrival_sum.sum()),'edge_event_count':int(a.edge_event_count)})
pd.DataFrame(top).to_parquet(exp/'drive/dominant-late-PN.parquet',index=False)
motifs=[]
for trial,g in path[path.endpoint=='first_late'].groupby('trial'):
g=g.sort_values('backward_hop')
pairs=g.groupby(['source_neuron_idx','target_neuron_idx']).size()
for (src,dst),n in pairs.items():
if (dst,src) in pairs.index:
motifs.append({'trial':trial,'source_neuron_idx':int(src),'target_neuron_idx':int(dst),
'forward_occurrences_on_ancestry':int(n),'reverse_occurrences_on_ancestry':int(pairs.loc[(dst,src)]),
'claim':'actual time-supported reciprocal motif on ranked path, not unique necessity'})
pd.DataFrame(motifs).to_parquet(exp/'ancestry/reciprocal-event-motifs.parquet',index=False)
loop=path[(path.endpoint=='first_late')&path.source_neuron_idx.isin([1381,1511])&path.target_neuron_idx.isin([1381,1511])]
for trial,g in loop.groupby('trial'):
timing.append({'trial':trial,'recurrent_edge_events_on_path':len(g),
'first_emission_ms':float(g.source_step.min()/10),'last_arrival_ms':float(g.arrival_step.max()/10),
'all_delays_exact':bool(((g.arrival_step-g.source_step)==18).all()),
'all_arrivals_precede_target':bool((g.arrival_step<=g.target_step).all())})
pd.DataFrame(timing).to_parquet(exp/'ancestry/loop-timing.parquet',index=False)
for r in q.itertuples():
z=pd.read_parquet(exp/'input-audit'/f'{r.odor_index}-{r.seed}-paired-requests.parquet')
diff=z[z.status_changed]
for (before,after),g in diff.groupby(['status_reference','status_off']):
transitions.append({'odor_index':r.odor_index,'reference_status':before,'off_status':after,'count':len(g)})
pd.DataFrame(transitions).to_parquet(exp/'input-audit/status-transitions.parquet',index=False)
isirows=[]
for f in sorted((exp/'states').glob('*-watched.npz')):
trial=f.name.removesuffix('-watched.npz');z=np.load(f);st=z['state'];ids=z['watch']
for j,node in enumerate(ids):
steps=np.flatnonzero(st[5000:,6,j])+5000;isi=np.diff(steps)
isirows.append({'trial':trial,'neuron_idx':int(node),'late_natural_spikes':len(steps),
'late_rate_hz':2*len(steps),'ISI_min_ms':float(isi.min()/10) if len(isi) else None,
'ISI_mean_ms':float(isi.mean()/10) if len(isi) else None,'ISI_max_ms':float(isi.max()/10) if len(isi) else None,
'ISI_equal_refractory_fraction':float(np.mean(isi==22)) if len(isi) else None,
'late_post_g_mean':float(st[5000:,5,j].mean(dtype=np.float64)),
'late_post_g_min':float(st[5000:,5,j].min()),'late_post_g_max':float(st[5000:,5,j].max())})
isidf=pd.DataFrame(isirows);isidf.to_parquet(exp/'drive/late-ISI-and-g.parquet',index=False)
# Read-only coverage audit of the frozen graded APL replacement, not an intervention.
canonical=ROOT/'datasets/male-cns-v1.0/processed/neuron-only'
fullmeta=pd.read_parquet(canonical/'neuron-index/unified-metadata.parquet').sort_values('neuron_idx').reset_index(drop=True)
assert np.array_equal(fullmeta.neuron_idx.to_numpy(),np.arange(166700))
ptr=np.load(canonical/'gpu/csr-indptr.int64.npy',mmap_mode='r')
dst=np.load(canonical/'gpu/csr-indices.int32.npy',mmap_mode='r')
raw=np.load(canonical/'gpu/csr-weights.int64.npy',mmap_mode='r')
aplids=np.flatnonzero(fullmeta.annotation_type.eq('APL'));assert len(aplids)==2
aplrows=[]
for src in aplids:
for eid in range(int(ptr[src]),int(ptr[src+1])):
node=int(dst[eid]);m=fullmeta.iloc[node]
aplrows.append({'edge_id':eid,'source_neuron_idx':int(src),'source_body_id':int(fullmeta.iloc[src].body_id),
'target_neuron_idx':node,'target_body_id':int(m.body_id),'target_class':m.annotation_class,
'target_type':m.annotation_type,'raw_synapse_count':int(raw[eid]),
'frozen_APL_fast_edge_zeroed':True,'graded_replacement_target':bool(m.annotation_class=='Kenyon_Cell')})
apldf=pd.DataFrame(aplrows);apldf.to_parquet(exp/'audit/APL-target-coverage.parquet',index=False)
omitted=apldf[~apldf.graded_replacement_target]
focus=omitted[omitted.target_neuron_idx.isin(watch[watch.annotation_class=='ALPN'].index)]
spec=importlib.util.find_spec('phase4p_revalidate');assert spec and spec.origin
srcfile=Path(spec.origin);assert srcfile.exists()
aplsummary={'scope':'static coverage of frozen DeterministicAPLSim.install and graded output matrix; not proof of functional effect',
'APL_neurons':len(aplids),'all_outgoing_edges':len(apldf),'all_outgoing_synapse_count':int(apldf.raw_synapse_count.sum()),
'non_KC_edges_without_graded_replacement':len(omitted),'non_KC_synapse_count_without_graded_replacement':int(omitted.raw_synapse_count.sum()),
'non_KC_target_neurons':int(omitted.target_neuron_idx.nunique()),'watched_PN_edges':len(focus),
'watched_PN_synapse_count':int(focus.raw_synapse_count.sum()),'source':str(srcfile),
'source_sha256':hashlib.sha256(srcfile.read_bytes()).hexdigest(),'baseline_not_modified':True,
'cannot_replace_terminal_inhibition_with_whole_PN_inhibition_without_new_hypothesis':True}
(exp/'audit/APL-target-coverage.json').write_text(json.dumps(aplsummary,indent=2)+'\n')
lines=['# Phase 4Z — mechanism interpretation','',
'Post-run descriptive analysis, not a new candidate gate. No tuning, learning, new stimulus or sign hypothesis.',
'The earlier generic label PN feedback is refined below using actual ALLN cell identities. Do not assume the residual loop is KC→MBON→PN feedback merely because that pathway exists in the literature.','',
'## Watched identities','', '| neuron_idx | body_id | class | type | metadata transmitter | raw ground_truth |','|---:|---:|---|---|---|---|']
for idx,r in watch.iterrows():lines.append(f'| {idx} | {r.body_id} | {r.annotation_class} | {r.annotation_type} | {r.neurotransmitter} | {r.nt_ground_truth} |')
lines+=['','## Representative first-late ancestry','',
'Odor 0 / seed 131101 / frozen reference. Ordered source spike → edge arrival → target spike (ms). These are ranked event-supported contributions; no unique predecessor assertion.',
'| source idx/type | target idx/type | source time | arrival time | target time | contacts | effective weight |','|---|---|---:|---:|---:|---:|---:|']
example=path[(path.trial=='0-131101-frozen_reference')&(path.endpoint=='first_late')].sort_values('backward_hop').head(5).iloc[::-1]
for r in example.itertuples():lines.append(f'| {r.source_neuron_idx}/{r.source_type} | {r.target_neuron_idx}/{r.target_type} | {r.source_step/10:.1f} | {r.arrival_step/10:.1f} | {r.target_step/10:.1f} | {r.raw_synapse_count} | {r.effective_weight:.6f} |')
lines+=['','## Dominant PN arrival contribution','', '| trial | PN idx/type | positive arrival fraction |','|---|---|---:|']
for r in top:lines.append(f'| {r["trial"]} | {r["source_neuron_idx"]}/{r["type"]} | {100*r["positive_fraction"]:.3f}% |')
lines+=['','Fractions are sums of positive synaptic arrivals, not percentages of causal responsibility. Negative input, recorded APL and membrane dynamics are retained in exact CPU replay.',
f'The preidentified 1381↔1511 event motif occurs in {len(timing)}/12 first-late ancestry paths. Full temporal edges and alternative inputs are retained.',
'Conditional omission results are in replay/results.parquet and REPORT.md. They retain other source events; loss of spikes there is not evidence that a whole-brain ablation would produce the same result.',
'Critically, omitting ALLN arrivals to a watched ALLN does not silence it in conditional replay: recorded PN and other inputs can still sustain high firing. Thus the 1381↔1511 motif cannot be called the sole cause or an isolated autonomous oscillator.',
'Removing the existing VA6_adPN→KC82312 edge in conditional replay eliminates target late firing; removing ALLN input to VA6_adPN also eliminates its late firing. These are two conditional links, not a demonstrated whole-brain necessity chain.','',
'## Late firing of preidentified ALLNs','', '| neuron_idx | late rate range Hz | min/max ISI ms | fraction ISI = refractory, range |','|---:|---:|---:|---:|']
for node,g in isidf[isidf.neuron_idx.isin([1197,1381,1511])].groupby('neuron_idx'):
lines.append(f'| {node} | {g.late_rate_hz.min()}–{g.late_rate_hz.max()} | {g.ISI_min_ms.min():.2f}/{g.ISI_max_ms.max():.2f} | {g.ISI_equal_refractory_fraction.min():.5f}–{g.ISI_equal_refractory_fraction.max():.5f} |')
lines+=['','This is a measured property of the mathematical model, not a validated biological firing rate. Window-boundary counting can put spikes/window slightly above 1/refractory. No new biological runaway threshold is inferred.','',
'## Frozen APL replacement target coverage','',
f'The frozen implementation zeroes all {len(apldf)} APL fast outgoing edges and reintroduces graded inhibition only on KC targets. {len(omitted)} non-KC edges ({int(omitted.raw_synapse_count.sum())} anatomical contacts) have no direct graded replacement in this implementation.',
f'This includes {len(focus)} edges / {int(focus.raw_synapse_count.sum())} anatomical contacts onto the seven watched PNs. Raw edges and counts remain intact. These edges were not removed by this audit.',
'This is an implementation coverage gap relative to the anatomical graph, not evidence that restoring a whole-neuron negative current would be biologically correct. APL inhibition of PN boutons is spatial/terminal-specific and a point-neuron approximation may not capture it.',
'Primary study: Prisco et al., eLife 10:e74172 (2021), https://pubmed.ncbi.nlm.nih.gov/34964714/ . It reports reciprocal APL inhibition targeting PN boutons and KC dendrites; it does not provide these MaleCNS edge-specific effective parameters.',
'No frozen APL parameter, edge, sign or output target was modified. Functional necessity of this coverage gap remains untested.','',
'## What input matching means','',
f'All six paired request schedules match. Status differs for {int(q.status_changed_requests.sum())} requests in total; actual output at requested times differs for {int(q.delivered_output_differences.sum())} requests.',
'An accepted→preexisting transition changes the accepted-new count without removing the output spike. A blocked transition can change actual delivery. Thus accepted-new totals alone are insufficient to describe input matching.',
'The first differing request is state dependent; this is an endogenous response to an existing intervention, not evidence of different external RNG. No clamp, injected substitute events or refractory override was introduced.',
'The earlier closed-loop intervention remains interpretable as a total system effect under matched requests. It is not an estimate of a direct DPM effect holding actual PN output trains fixed.','',
'## Literature boundary','',
'Excitatory local antennal-lobe circuits and reciprocal local-neuron/PN excitation have experimental support. They must not automatically be relabeled inhibitory because the model has persistent activity.',
'Shang et al., Cell 2007: https://pubmed.ncbi.nlm.nih.gov/17289577/',
'Huang et al., Neuron 2010, functional connectivity of excitatory local interneurons: https://doi.org/10.1016/j.neuron.2010.08.025',
'These experiments do not validate uniform current weights, this MaleCNS individual, or its simulated persistent firing. No receptor, adaptation or electrical-synapse parameters are inferred or inserted.',
'Next priority is independent assessment of local AL-loop dynamics and input-interface causal controls, not changing dopamine learning rate. No official model replacement is justified.']
(exp/'audit/MECHANISM.md').write_text('\n'.join(lines)+'\n',encoding='utf-8')
print(json.dumps({'dominant_PNs':sorted({r['source_neuron_idx'] for r in top}),
'dominant_PN_positive_fraction_range':[min(r['positive_fraction'] for r in top),max(r['positive_fraction'] for r in top)],
'trials_with_preidentified_reciprocal_motif':len(timing),
'status_changed_requests':int(q.status_changed_requests.sum()),'delivered_output_differences':int(q.delivered_output_differences.sum())}),flush=True)