Testing inference: multiple patterns in a long sequence¶
The trained readout from notebook 11 is loaded from the cache and probed in retrieval mode: all 16 motifs are concatenated into one long input stream, each flanked by N_iti_time = 4*N_pretime = 200 ms of silence (no spontaneous activity), and only its first num_delay ms is injected as a trigger. The network - with no further plasticity - must then re-emit each motif from the correct onset onward. This is a chain-retrieval stress test of the fixed weights (no retraining).
Notebook summary¶
Pipeline¶
Load
<datetag>_vanilla.pthviamnesis_chains.load()(StochasticSpikingPattern).Concatenate the 16 motifs (with long ITI pads) into a single
(N_neuron, T)stream; build the per-motif trigger masks; run one singleforward_passon the stream.Score every consecutive
window_size=60ms window withget_f1score(spikes_window, target_window)against the frozen motifs, cached as<datetag>_score.npyunder the standard.lockprotocol.Plot the F1 trajectory and an animated version (
retrieval.pdf/retrieval.mp4).Compare the failures with those predicted by the cascade model of notebook 14 (
get_cascade_prediction, no fitted parameter): recall of each motif, and onset of the gap ofnum_delaysilent bins at which its replay stops (a table, no figure).
Main outputs¶
cache
../cached_data/<datetag>_score.npy(+.lock);figures
retrieval.pdf/png/svgand (when figpath is active) theretrieval.mp4animation.
[1]:
# load() rebuilds the HD_SNN architecture, restores <datetag>_vanilla.pth and
# attaches a default StochasticSpikingPattern (see mnesis_chains.load()).
from mnesis_boilerplate import torch, np, pd, plt, data_cache, cache_lock, figpath, RECOMPUTE, Params
from mnesis_boilerplate import printfig, get_f1score, get_scores, get_cascade_prediction
from mnesis_chains import load
# RECOMPUTE = True
datetag = '2026-09-07'
Using device: mps
[2]:
opt = Params()
model_filename = data_cache / f"{opt.datetag}_vanilla.pth"
hd = load(opt, model_filename)
print(f"Model weights loaded from {model_filename} - the pattern generator is the {hd.pattern_object.desc}")
Model weights loaded from ../cached_data/2026-09-07_vanilla.pth - the pattern generator is the Stochastic spike pattern generator
Retrieval across a long, concatenated sequence with all motifs¶
Construction of the input stream: concat_motifs pads every motif with N_iti_time = 4*N_pretime = 200 ms of silence on both sides and concatenates the 16 motifs head-to-tail, giving (N_neuron, N_pattern * (N_time + 2*N_iti_time)) tensors:
target_full: the frozen motifs - the reference to score against;cue_full: a seeded noisy realisation (pattern_object(seed=seed));evoked: 1 over each motif’s full duration - the ideal retrieval envelope;mask: 1 on each motif’s firstnum_delaybins, 0 elsewhere - keeps only the trigger;input_spikes = cue_full * mask: the concatenated triggers (spontaneous activity is silenced here, unlikeget_input_spikes, to isolate pure trigger-driven retrieval).
Scanning: get_f1score is evaluated in consecutive window_size=60 ms windows starting at onsets, giving the per-position retrieval quality along the stream.
[3]:
# One-shot scan: cached F1 scores per 60-ms window along the concatenated stream.
window_size = 60
N_iti_time = opt.N_pretime * 4
T_full = opt.N_pattern * (opt.N_time + 2*N_iti_time)
def concat_motifs(x):
"""Pad each motif of x (N_pattern, N_neuron, N_time) with N_iti_time silent bins
on both sides and concatenate them in time -> (N_neuron, T_full)."""
x = torch.nn.functional.pad(x, (N_iti_time, N_iti_time))
return x.permute((1, 0, 2)).reshape((opt.N_neuron, T_full))
# the frozen motifs are the reference; a seeded noisy realisation provides the triggers
target_full = concat_motifs(hd.pattern_object.frozen_target)
cue_full = concat_motifs(hd.pattern_object(seed=opt.seed))
ones = torch.ones_like(hd.pattern_object.frozen_target)
evoked = concat_motifs(ones)
mask = ones.clone()
mask[:, :, opt.num_delay:] = 0
mask = concat_motifs(mask)
onsets = np.arange(0, T_full - window_size + 1, window_size)
npy_filename = data_cache / f"{hd.opt.datetag}_score.npy"
lock_filename = npy_filename.with_suffix('.lock')
# RECOMPUTE = True # the cache key ignores window_size: comment out only to reuse a score computed with the same window
if RECOMPUTE:
npy_filename.unlink(missing_ok=True)
lock_filename.unlink(missing_ok=True)
scores_np = None
try:
scores_np = np.load(npy_filename)
lock_filename.unlink(missing_ok=True) # in case the lock file was not unlinked
print(f"Scores loaded from {npy_filename}")
except FileNotFoundError:
with cache_lock(npy_filename) as acquired:
if acquired:
print(f"Score file not found: {npy_filename}, evaluating the score.")
with torch.no_grad():
_, _, spikes = hd.forward_pass((cue_full * mask).unsqueeze(0).to(opt.device))
spikes = spikes.squeeze(0)
scores_np = np.array([get_f1score(spikes[:, t:t+window_size], target_full[:, t:t+window_size]).item()
for t in onsets])
np.save(npy_filename, scores_np)
else:
print(f"Score file is locked: {lock_filename}, passing.")
# reference curves, sampled at the onset of each scoring window
evoked_np = evoked.mean(dim=0).cpu().numpy()[onsets]
mask_np = mask.mean(dim=0).cpu().numpy()[onsets]
Scores loaded from ../cached_data/2026-09-07_score.npy
F1 trajectory of the retrieval scan¶
Curves vs the onset time (ms): evoked_np (ideal 1 while a motif should be re-emitted - the reference envelope), mask_np (1 only on each trigger’s num_delay-ms window) and scores_np (the measured per-window F1). Retrieval should light up after every trigger.
Each motif is numbered above its envelope, from 1 (as the ordinals of the paper captions). The number of the motif shown in the rasters of notebook 11 (opt.i_pattern = 1, 0-based, so number 2; Figure fig:target of the paper) is in bold.
[9]:
def plot_retrieval(ax):
# colours follow the default cycle (C0, C1, C2), as described in the caption of fig:retrieval
ax.axhline(0, color='0.7', ls='--', lw=0.8) # grey dashed guides at the extremal F1 values
ax.axhline(1, color='0.7', ls='--', lw=0.8)
ax.plot(onsets, evoked_np, alpha=.5, label='ideal retrieval')
# F1 shown on the motifs only: in the ITIs the target is empty, so F1 is 0 whatever the network does
ax.plot(onsets, mask_np, alpha=1., lw=1.5, label='trigger')
line, = ax.plot(onsets, np.where(evoked_np > 0, scores_np, np.nan), lw=2, alpha=.5, label='MNESIS')
# area under the F1 curve, in the colour of the line (fill_between draws from its own colour cycle)
ax.fill_between(onsets, 0, scores_np, where=evoked_np > 0, color=line.get_color(), alpha=.2, lw=0)
# motif number (1-based, as the ordinals of the caption of fig:retrieval) above the centre of each motif,
# in bold for the motif shown in the rasters of notebook 11 (fig:target);
# clip_on keeps the labels outside the viewport of the animation hidden
for k in range(opt.N_pattern):
ax.text(k*(opt.N_time + 2*N_iti_time) + N_iti_time + opt.N_time/2, 1.03, f'{k+1}',
ha='center', va='bottom', fontsize='small', clip_on=True,
fontweight='bold' if k == opt.i_pattern else 'normal')
# headroom for the labels, with the ticks kept on the [0, 1] range of the F1-score
ax.set_ylim(-.05, 1.15)
ax.set_yticks(np.linspace(0, 1, 6))
ax.set_xlabel("Time (ms)")
ax.set_ylabel("F1-score")
# the curves span the whole [0, 1] range: the legend goes above the axes
ax.legend(loc='lower center', bbox_to_anchor=(.5, .95), ncols=3, frameon=False)
if scores_np is not None:
fig, ax = plt.subplots()
plot_retrieval(ax)
printfig(fig, 'retrieval', fig_width=opt.fig_width, fig_height=opt.fig_height, figpath=figpath)
Saving as ../figures/retrieval.pdf
Saving as ../figures/retrieval.png
Saving as ../figures/retrieval.svg
Failures predicted by the cascade model¶
The replay is a cascade (notebook 14): the spikes of a motif in a bin are emitted if at least a fraction \(\vartheta\) of the \(K\) spikes of that motif in the \(D\) preceding bins are present, and each emitted spike joins the support of the next ones. A gap of \(D\) silent bins in a motif leaves its next spike without any support (\(K = 0\)): the replay cannot cross it, whatever the trigger. Applied to the \(M = 16\) frozen motifs with no fitted parameter
(get_cascade_prediction, which simulates the cascade bin after bin on the actual motifs), the model predicts which motifs stop, where, and with which recall. The table compares, for each motif:
the predicted recall with the frozen trigger window and with the seeded noisy realisation that provides it (
cue), and the recall measured with the protocol of notebook 11 (same cue, spontaneous pad) - the model predicts the recall, not the precision, as it leaves out the crosstalk between motifs, so the measured precision is shown alongside;the predicted onset of the first gap (ms from the onset of the motif) and the breakdown measured in the stream above: the end of the last \(60\)-ms window with \(F_1 \geq 0.5\), when the replay does not reach the end of the motif (the resolution is that of the windows).
[5]:
# Cascade model (no fitted parameter) on the frozen motifs, against the replay measured with the
# protocol of notebook 11 and the breakdown measured in the stream above.
cue = hd.pattern_object(seed=opt.seed)
prediction = get_cascade_prediction(hd.pattern_object.frozen_target, opt.num_delay, opt.lif_threshold)
prediction_cue = get_cascade_prediction(hd.pattern_object.frozen_target, opt.num_delay, opt.lif_threshold, cue=cue)
with torch.no_grad():
_, _, spikes = hd.forward_pass(hd.get_input_spikes(target=cue, seed=opt.seed))
spikes_evoked = spikes[:, :, (opt.N_pretime+opt.num_delay):(opt.N_time+opt.N_pretime)]
target_evoked = hd.pattern_object.frozen_target[:, :, opt.num_delay:]
precision_motif, recall_motif, _ = [score.cpu().numpy() for score in get_scores(spikes_evoked, target_evoked, dim=(1, 2))]
precision_pooled, recall_pooled, _ = [score.item() for score in get_scores(spikes_evoked, target_evoked)]
# breakdown in the stream: end of the last window with F1 >= 0.5, when the replay does not reach the end of the motif
t_break = np.full(opt.N_pattern, np.nan)
if scores_np is not None:
for i_motif in range(opt.N_pattern):
start = i_motif * (opt.N_time + 2*N_iti_time) + N_iti_time # onset of the motif in the stream
inside = (onsets >= start + opt.num_delay) & (onsets + window_size <= start + opt.N_time)
good = scores_np[inside] >= .5
if not good[-1]:
t_break[i_motif] = onsets[inside][np.flatnonzero(good)[-1]] + window_size - start if good.any() else opt.num_delay
table = pd.DataFrame({'recall predicted': prediction['recall'], 'recall predicted (cue)': prediction_cue['recall'],
'recall measured': recall_motif, 'precision measured': precision_motif,
'gap onset predicted (ms)': prediction['t_gap'], 'breakdown measured (ms)': t_break},
index=pd.Index(np.arange(1, opt.N_pattern + 1), name='motif')) # numbered from 1, as in the figure
print(table.round(3).to_string())
# pooled recall predicted on the motifs (weighted by their number of spikes), as the measured one
n_spikes = target_evoked.sum(dim=(1, 2)).cpu().numpy()
for name, pred in [('frozen trigger', prediction), ('cue', prediction_cue)]:
print(f'predicted pooled recall ({name}) = {(pred["recall"] * n_spikes).sum() / n_spikes.sum():.3f}; ', end='')
print(f'measured: recall = {recall_pooled:.3f}, precision = {precision_pooled:.3f} (not predicted)')
# motif numbers from 1, as plain integers for printing
print(f'motifs predicted to stop at a gap: {(np.flatnonzero(prediction["recall"] < 1) + 1).tolist()}; '
f'motifs whose replay breaks down in the stream: {(np.flatnonzero(~np.isnan(t_break)) + 1).tolist()}; '
f'motifs predicted with a lower recall from the cue: {(np.flatnonzero(prediction_cue["recall"] < prediction["recall"]) + 1).tolist()}')
recall predicted recall predicted (cue) recall measured precision measured gap onset predicted (ms) breakdown measured (ms)
motif
1 1.000 1.000 1.000 0.546 NaN NaN
2 1.000 1.000 1.000 1.000 NaN NaN
3 1.000 1.000 1.000 1.000 NaN NaN
4 0.175 0.175 0.175 1.000 233.0 280.0
5 1.000 1.000 1.000 0.586 NaN NaN
6 1.000 1.000 1.000 1.000 NaN NaN
7 1.000 1.000 1.000 1.000 NaN NaN
8 1.000 1.000 1.000 1.000 NaN NaN
9 1.000 1.000 1.000 1.000 NaN NaN
10 0.221 0.221 0.221 1.000 258.0 280.0
11 1.000 1.000 1.000 0.653 NaN NaN
12 0.294 0.294 0.294 1.000 268.0 300.0
13 0.564 0.564 0.564 0.462 568.0 580.0
14 1.000 1.000 1.000 1.000 NaN NaN
15 1.000 1.000 1.000 1.000 NaN NaN
16 1.000 1.000 1.000 1.000 NaN NaN
predicted pooled recall (frozen trigger) = 0.825; predicted pooled recall (cue) = 0.825; measured: recall = 0.825, precision = 0.835 (not predicted)
motifs predicted to stop at a gap: [4, 10, 12, 13]; motifs whose replay breaks down in the stream: [4, 10, 12, 13]; motifs predicted with a lower recall from the cue: []
Animated version¶
Same figure seen through a viewport of n_view=100 scoring windows (6 s) swept along the stream by steps of 5 windows (ffmpeg, 10 fps -> retrieval.mp4); same colours as the static figure.
[6]:
if scores_np is not None and figpath is not None:
filename = figpath / 'retrieval.mp4'
# RECOMPUTE = True # comment out to reuse a previously computed video
if not filename.exists() or RECOMPUTE:
import matplotlib.animation as animation
cm = 1/2.54 # centimetres in inches (matplotlib figsize)
n_view, step = 100, 5 # viewport width and shift between frames (scoring windows)
# the curves are drawn once; each frame only moves the x limits
fig, ax = plt.subplots(figsize=(opt.fig_width*cm, opt.fig_height*cm))
plot_retrieval(ax)
def animate(frame):
ax.set_xlim(onsets[frame], onsets[frame + n_view - 1])
return []
anim = animation.FuncAnimation(fig, animate, frames=range(0, len(onsets) - n_view + 1, step))
anim.save(filename, writer='ffmpeg', fps=10)
plt.close(fig)