LFP filter output depends on available RAM at populate time

[edeno]

Problem

FirFilterParameters.filter_data_nwb filters each valid interval through one of
two branches, chosen at runtime by a free-memory check in
src/spyglass/common/common_filter.py (req_mem < MEM_USE_LIMIT * psutil.virtual_memory().available):

branch	signal passed to the filter	what the filter sees past the interval edge
fits in RAM	data_on_disk[start:stop] — a copy that ends at the interval boundary	nothing → implicit zero padding
does not fit	data_on_disk — the whole dataset	real neighboring samples

A mode='full' convolution with output_index_bounds=[filter_delay, filter_delay + n] needs input from -left_context to n - 1 + right_context
relative to the interval, so the branches genuinely see different input at the
edges.

Measured

20 000-sample interval, standard 0–400 Hz LFP filter (6401 taps, filter_delay
= 3200), unit-variance input:

• 6398 of 20 000 output samples differ
• max |Δ| = 0.160 = 96 % of the filtered signal's standard deviation, at the
  first sample
• mean |Δ| over the leading filter_delay samples = 2.8 % of the std
• the in-memory branch reproduces a zero-padded np.convolve reference exactly;
  the on-disk branch reproduces a real-context reference exactly

At 30 kHz, filter_delay = 3200 samples ≈ 106.7 ms per edge (213.3 ms
across both edges of an interval).

Consequence: the same raw file and the same LFPSelection populated on a
laptop versus a big-memory node produce different LFP. Which branch ran was
never recorded.

Relationship to PR #1635

The branch structure and the resulting edge behavior are inherited from the
ghostipy-backed code; #1635 did not change them. It did normalize two of the
input_index_bounds in filter_data_nwb to exclusive-stop semantics, verified
numerically inert (the block loop reads by block and clips at the array bounds,
so the stop only feeds an internal length check).

Which value-producing path is the outlier

Two code paths actually produce filtered values with real neighboring context —
filter_data_nwb's on-disk branch and filter_data — and one, the in-memory
branch of filter_data_nwb, zero-pads. (filter_data_nwb's sizing pass calls
fir.describe_output, which computes shapes only and reads nothing, so it is
not a third path.)

Decide the policy first

This is not a purely mechanical bug. Two questions have to be answered before
implementing anything:

1. Should interval edges draw on samples outside valid_times at all? The
   filter_delay samples beyond an interval boundary are, by construction,
   outside the valid interval — often across a gap, so from a different time.
   Using them is a scientific choice, not a neutral default.
2. If not, isolate how? Zero padding introduces a ramp transient at every
   boundary; reflection avoids the transient but is a new convention.

Whatever is chosen, both branches must implement it identically. Maintainers
should sign off on the policy explicitly — none of the options below is neutral.

Options

A. Standardize on real neighboring context

Extend the in-memory read so it covers the filter's full support, and pass
input_index_bounds relative to that window. Both branches then see identical
input.

Exact support — do not assume symmetry, it only holds for odd-length
filters and filter_coeff is caller-supplied. The bounds shift is the
alignment-critical part: the data window moves, but the requested output window
must not, so input_index_bounds has to be re-expressed relative to lo.

delay = (len(filter_coeff) - 1) // 2
left_context = len(filter_coeff) - 1 - delay
right_context = delay

lo = max(start - left_context, 0)
hi = min(stop + right_context, n_samples)

data = data_on_disk[lo:hi]                    # extended window
input_index_bounds = [start - lo, stop - lo]  # SAME output, re-based on lo

# timestamps stay on the ORIGINAL interval -- see the note below
timestamps = timestamps_on_disk[start:stop]

For the odd-length filters spyglass designs, left_context == right_context == delay. For an even-length custom filter they differ, and a symmetric pad would
under-read on the left.

• Changes fits-in-RAM populates; leaves memory-constrained (on-disk) ones
  unchanged.
• At the dataset's own start/end the window clips, and the on-disk branch clips
  identically, so no special case is needed.
• Extra RAM: (left_context + right_context) * n_electrodes — 410 KB for a
  6401-tap filter over 32 int16 channels.
• Verified: reproduces the on-disk result to ≤ 2.8e-16 across interior,
  dataset-start, dataset-end, whole-dataset, and shorter-than-filter intervals.
  Not bit-identical — the two paths land FFT block boundaries differently, so
  the residual is round-off.

Implementation notes:

• Keep timestamps = timestamps_on_disk[start:stop] unextended;
  extracted_ts = timestamps[0::decimation] must still start at the interval
  start. Extending the timestamp window too would silently shift every LFP
  timestamp — the one easy way to get this change wrong.

• Include the extended window in the req_mem estimate. It currently sizes
  from interval_samples; with the halo the data term must size from the
  extended window, while the timestamp term stays on the interval:

  req_mem = (
      interval_samples * timestamps_on_disk[0].itemsize
      + (hi - lo) * n_electrodes * data_on_disk[0][0].itemsize
  )

B. Delete the in-memory branch

One code path, no psutil check, no MEM_USE_LIMIT; the divergence becomes
structurally impossible. Same edge policy as Option A (real context), ~25 lines
lighter.

Measured locally (2 min × 32 ch @ 30 kHz, decimation_factor=15, chunked HDF5
on SSD): always-on-disk was 12 % faster, not slower.

Caveat: that is local SSD. Analysis files commonly live on NFS, where many
small chunk reads can behave far worse than one bulk read. Repeat the
measurement on the deployment filesystem before adopting.

C. Isolate each interval

Standardize on padding instead of real context. Zero padding is what the
in-memory branch already does, so it would change only the on-disk branch;
reflection is a new convention and would change both. Avoids pulling samples
across gaps between valid intervals, at the cost of an edge artifact (a ramp for
zeros, a mirror for reflection).

Recommendation

Answer the policy question first. If real context is intended, implement
Option A with the asymmetric support above, as its own PR; revisit Option B
once NFS numbers exist. If isolation is intended, Option C — and then Option A's
mechanics do not apply.

Either way, record an edge-context algorithm version alongside the result so the
convention in force is recoverable from the data.

Verification plan

1. The regression test that would have caught this, against
   filter_data_nwb directly. This is the one that matters and the one that
   does not exist: filter_data_nwb has no direct test at all
   (test_filter_data in tests/common/test_filter.py is
   @pytest.mark.skip), and the timestamp-alignment test added in Vendor FIR filter code, drop ghostipy dependency #1635
   exercises filter_data, which has only one branch. Force each memory branch
   on the same small NWB fixture and assert the two outputs agree. Make the
   branch choice injectable (a _fits_in_memory(...) helper or a max_memory
   argument) rather than monkeypatching psutil, so the test states its
   intent. It must fail before the fix -- that is the point.
2. Numeric test against a np.convolve reference over the extended window, for
   interior, dataset-start, dataset-end, and shorter-than-filter intervals.
3. Timestamp test against filter_data_nwb: stored LFP timestamps unchanged
   by the extension, and still aligned with the stored data (an impulse's peak
   lands on the impulse's own timestamp). Guards the trap above, which is the
   easiest way to get this change wrong.
4. Under Option A, on a large-memory machine the pre-fix result should
   change and then match what the forced-low-memory (on-disk) path produces.
   Force each branch explicitly rather than relying on ambient free memory.
5. Cover both time-axis layouts ((time, electrode) and (electrode, time)),
   an even-length filter, decimation, multiple intervals, and an h5py signal.

---

Found while reviewing #1635. Deliberately not fixed there: it changes results
and needs the policy decision above.

[edeno]

Additional context from a multi-agent review of #1635 — two findings that sharpen the scope here.

filter_data and filter_data_nwb already disagree, independent of RAM. filter_data (which LFPBand uses) always passes the whole array with absolute bounds, i.e. the real-context convention. So the isolate-each-interval option (C) would change LFPBand as well as LFP, not just the on-disk branch — worth folding into the policy decision rather than discovering during implementation.

The input_index_bounds docstring was actively misleading, and is now fixed in #1635. It read as [start, stop) of the input, which sounds like a promise that only that range is read. It is not — the filter draws support from neighbouring samples outside the window, which is precisely the mechanism here. That wording is a large part of why the divergence is easy to write and hard to notice, so the corrected docstring is worth reading before implementing any option.

Measurement corroborated independently. Two reviewers reproduced the divergence from scratch at different filter lengths: 200 of 800 samples differing (101-tap, max abs 0.164) and 100 of 800 (51-tap, max abs 0.158) — in both cases exactly 2 * filter_delay samples per lane, split across the two edges, matching the left_context/right_context formula in the plan above. One also confirmed the in-memory branch reproduces an isolated-segment convolution exactly while the on-disk branch reproduces a continuous-stream convolution exactly, so the two conventions are each internally consistent — this is a choice between them, not a bug in either.

[edeno]

Plan revised at the top of this issue. Substantive changes:

• The alignment-critical step is now spelled out. The data window moves but the requested output window must not, so input_index_bounds has to be re-based on lo:

  data = data_on_disk[lo:hi]
  input_index_bounds = [start - lo, stop - lo]
  timestamps = timestamps_on_disk[start:stop]   # NOT extended

• The memory estimate needs two different lengths. Timestamps stay on the interval; only the data term grows to (hi - lo).

• Verification step 1 is now explicitly against filter_data_nwb. That method has no direct test at all (test_filter_data is @pytest.mark.skip), and the timestamp-alignment test added in Vendor FIR filter code, drop ghostipy dependency #1635 exercises filter_data, which has only one branch — so it cannot catch this. The branch-parity test has to force each memory branch on an NWB fixture, and must fail before the fix.

• Stale references to describe_dims=True replaced with fir.describe_output (renamed in Vendor FIR filter code, drop ghostipy dependency #1635), and stale line numbers dropped in favour of naming the code.