From e4fe593a42484f788a99abe6fec1a9ecbcfa0b19 Mon Sep 17 00:00:00 2001
From: Christian Brodbeck <christianbrodbeck@nyu.edu>
Date: Fri, 7 Aug 2026 14:06:14 -0600
Subject: [PATCH 01/11] Test

---
 mne/preprocessing/tests/test_maxwell.py | 35 +++++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/mne/preprocessing/tests/test_maxwell.py b/mne/preprocessing/tests/test_maxwell.py
index f74195292c4..b02bc3ae43c 100644
--- a/mne/preprocessing/tests/test_maxwell.py
+++ b/mne/preprocessing/tests/test_maxwell.py
@@ -54,6 +54,7 @@
     _trans_sss_basis,
 )
 from mne.rank import _compute_rank_int, _get_rank_sss, compute_rank
+from mne.transforms import rot_to_quat
 from mne.utils import (
     _record_warnings,
     assert_meg_snr,
@@ -2056,6 +2057,40 @@ def test_find_bads_maxwell_flat():
     assert noisy == want_noisy
 
 
+@pytest.mark.slowtest
+@testing.requires_testing_data
+def test_find_bads_maxwell_head_pos():
+    """Test that each interval uses its own head positions."""
+    raw = read_raw_fif(raw_fname, allow_maxshield="yes")
+    raw.pick("meg", exclude=()).crop(0, 10)  # two 5 s intervals
+    raw.load_data()
+    # One head position per second, translating steadily so that using the positions
+    # from the wrong time window gives a different result.
+    trans = raw.info["dev_head_t"]["trans"]
+    head_pos = np.zeros((11, 10))
+    head_pos[:, 0] = raw._first_time + np.arange(11.0)
+    head_pos[:, 1:4] = rot_to_quat(trans[:3, :3])
+    head_pos[:, 4:7] = trans[:3, 3]
+    head_pos[:, 6] += np.arange(11) * 5e-3  # 5 mm/s
+    kwargs = dict(
+        origin=(0.0, 0.0, 0.04),
+        regularize=None,
+        bad_condition="ignore",
+        min_count=1,
+        return_scores=True,
+        h_freq=None,  # keep the two calls below operating on identical data
+    )
+    scores = find_bad_channels_maxwell(raw, head_pos=head_pos, **kwargs)[2]
+    assert scores["bins"].shape == (2, 2)
+    # Processing the second interval on its own must give the same scores as
+    # processing it as part of the whole recording.
+    raw_crop = raw.copy().crop(5.0)
+    pos_crop = head_pos[5:]
+    scores_crop = find_bad_channels_maxwell(raw_crop, head_pos=pos_crop, **kwargs)[2]
+    assert scores_crop["bins"].shape == (1, 2)
+    assert_allclose(scores_crop["scores_noisy"][:, 0], scores["scores_noisy"][:, 1])
+
+
 @pytest.mark.parametrize(
     "regularize, n, int_order",
     [

From 480c2a5ced06e1f4b03269bdcc4f123205d424e9 Mon Sep 17 00:00:00 2001
From: Christian Brodbeck <christianbrodbeck@nyu.edu>
Date: Fri, 7 Aug 2026 14:12:46 -0600
Subject: [PATCH 02/11] fix find_bad_channels_maxwell

---
 mne/_ola.py                  | 26 ++++++++++++++++++++------
 mne/preprocessing/maxwell.py | 21 +++++++++++++++++----
 mne/tests/test_ola.py        | 17 +++++++++++++++++
 3 files changed, 54 insertions(+), 10 deletions(-)

diff --git a/mne/_ola.py b/mne/_ola.py
index e182928f4a5..64c2eb72527 100644
--- a/mne/_ola.py
+++ b/mne/_ola.py
@@ -23,6 +23,11 @@ class _Interp2:
         arrays that must be interpolated.
     interp : str
         Can be 'zero', 'linear', 'hann', or 'cos2' (same as hann).
+    start : int
+        The control-point index at which feeding begins, for processing a segment
+        that does not start at the first control point. Feeding ``n`` points then
+        yields exactly what a continuous pass would yield for ``start`` to
+        ``start + n``, interpolation phase included.
 
     Notes
     -----
@@ -42,7 +47,9 @@ class _Interp2:
 
     """
 
-    def __init__(self, control_points, values, interp="hann", *, name="Interp2"):
+    def __init__(
+        self, control_points, values, interp="hann", *, name="Interp2", start=0
+    ):
         # set up interpolation
         self.control_points = np.array(control_points, int).ravel()
         if not np.array_equal(np.unique(self.control_points), self.control_points):
@@ -76,8 +83,11 @@ def val(pt):
             values = val
         self.values = values
         self.n_last = None
-        self._position = 0  # start at zero
-        self._left_idx = 0
+        self._position = start
+        # The last control point at or before start is the one in effect there, so
+        # feeding resumes mid-interval with the correct interpolation phase.
+        left_idx = np.searchsorted(self.control_points, start, "right") - 1
+        self._left_idx = max(left_idx, 0)
         self._left = self._right = self._use_interp = None
         self.name = name
         known_types = ("cos2", "linear", "zero", "hann")
@@ -94,9 +104,13 @@ def feed_generator(self, n_pts):
         logger.debug(f"    ~ {self.name} Feed {n_pts} ({self._position}-{stop})")
         used = np.zeros(n_pts, bool)
         if self._left is None:  # first one
-            logger.debug(f"    ~   {self.name} Eval @ 0 ({self.control_points[0]})")
-            self._left = self.values(self.control_points[0])
-            if len(self.control_points) == 1:
+            left_idx = self._left_idx
+            logger.debug(
+                f"    ~   {self.name} Eval @ {left_idx} "
+                f"({self.control_points[left_idx]})"
+            )
+            self._left = self.values(self.control_points[left_idx])
+            if left_idx == len(self.control_points) - 1:  # nothing to interpolate to
                 self._right = self._left
         n_used = 0
 
diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py
index 2f013961293..277c346dc40 100644
--- a/mne/preprocessing/maxwell.py
+++ b/mne/preprocessing/maxwell.py
@@ -741,6 +741,7 @@ def _run_maxwell_filter(
     st_fixed,
     st_overlap,
     mc,
+    mc_start=0,
 ):
     # Eventually find_bad_channels_maxwell could be sped up by moving this
     # outside the loop (e.g., in the prep function) but regularization depends
@@ -776,7 +777,7 @@ def _run_maxwell_filter(
 
     # This must be initialized inside _run_maxwell_filter because
     # find_bad_channels_maxwell modifies good_mask
-    mc.initialize(_get_this_decomp_trans, info["dev_head_t"], S_recon)
+    mc.initialize(_get_this_decomp_trans, info["dev_head_t"], S_recon, mc_start)
     update_kwargs.update(reg_moments=mc.reg_moments_0)
 
     # Process each valid block of data separately
@@ -865,6 +866,10 @@ class _MoveComp:
     """Perform movement compensation."""
 
     def __init__(self, pos, head_frame, raw, interp, reconstruct):
+        #   pos[0]: (n_pos, 4, 4): the dev_head_t transformation matrices
+        #   pos[1]: (n_pos,): sample indices into the recording, starting at 0
+        #   pos[2]: (n_pos, 9): rotation quaternion (:3), translation (3:6),
+        #       goodness of fit, error and velocity (6:9)
         self.pos = pos
         self.sfreq = raw.info["sfreq"]
         self.interp = interp
@@ -890,24 +895,32 @@ def get_decomp_by_offset(self, offset):
         op_resid -= np.dot(S_decomp[:, n_use_in:], pS_decomp[n_use_in:])
         return op_sss, op_in, op_resid
 
-    def initialize(self, get_decomp, dev_head_t, S_recon):
+    def initialize(self, get_decomp, dev_head_t, S_recon, start=0):
         """Secondary initialization."""
+        # Head positions are indexed relative to the start of the recording, so a
+        # segment that does not begin there (find_bad_channels_maxwell processes
+        # one interval at a time) must say where it starts, otherwise it would be
+        # compensated with the positions from the beginning of the recording.
+        self.start = start
         self.smooth = _Interp2(
             self.pos[1],
             self.get_decomp_by_offset,
             interp=self.interp,
             name="MC",
+            start=start,
         )
         _, _, pS_decomp, self.reg_moments_0, _ = get_decomp(dev_head_t, t=0.0)
         self.n_good = pS_decomp.shape[1]
         self.S_recon = S_recon
-        self.offset = 0
+        self.offset = start
         self.get_decomp = get_decomp
         # For the average passes
         self.last_avg_quat = np.nan * np.ones(6)
 
     def get_avg_op(self, *, start, stop):
         """Apply an average transformation over the next interval."""
+        # _COLA counts from the start of the segment, self.pos from the recording
+        start, stop = start + self.start, stop + self.start
         n_positions, avg_quat = _trans_lims(self.pos, start, stop)[1:]
         if not np.allclose(avg_quat, self.last_avg_quat, atol=1e-7):
             self.last_avg_quat = avg_quat
@@ -2941,7 +2954,7 @@ def find_bad_channels_maxwell(
             chunk_raw._data[:] = orig_data
             delta = chunk_raw.get_data(these_picks)
             with use_log_level(_verbose_safe_false()):
-                _run_maxwell_filter(chunk_raw, copy=False, **params)
+                _run_maxwell_filter(chunk_raw, copy=False, mc_start=start, **params)
 
             if n_iter == 1 and len(chunk_flats):
                 logger.info(
diff --git a/mne/tests/test_ola.py b/mne/tests/test_ola.py
index 0528b68f363..e1d7938f3a4 100644
--- a/mne/tests/test_ola.py
+++ b/mne/tests/test_ola.py
@@ -84,6 +84,23 @@ def test_interp_2pt():
         assert_allclose(out, expected)
 
 
+@pytest.mark.parametrize("interp", ("zero", "linear", "hann"))
+def test_interp_2pt_start(interp):
+    """Test that starting mid-stream matches a continuous pass."""
+    # control points deliberately uneven, so that most starts land mid-interval
+    control_points = [0, 100, 240, 390, 520, 680, 800]
+    values = np.array(control_points, float)
+    want = _Interp2(control_points, values, interp).feed(900)[0]
+    for start in range(0, 900, 7):
+        interper = _Interp2(control_points, values, interp, start=start)
+        assert_allclose(interper.feed(900 - start)[0], want[start:], atol=1e-12)
+    # state must carry across feeds of differing size, as when a segment is read
+    # in buffer-sized blocks
+    interper = _Interp2(control_points, values, interp, start=300)
+    out = np.concatenate([interper.feed(n)[0] for n in (70, 130, 100, 300)])
+    assert_allclose(out, want[300:900], atol=1e-12)
+
+
 @pytest.mark.parametrize("ndim", (1, 2, 3))
 def test_cola(ndim):
     """Test COLA processing."""

From 379c56428ec1caacd2f078430e7c6cc23dca316f Mon Sep 17 00:00:00 2001
From: Christian Brodbeck <christianbrodbeck@nyu.edu>
Date: Fri, 7 Aug 2026 14:19:09 -0600
Subject: [PATCH 03/11] changelog

---
 doc/changes/dev/14139.bugfix.rst | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 doc/changes/dev/14139.bugfix.rst

diff --git a/doc/changes/dev/14139.bugfix.rst b/doc/changes/dev/14139.bugfix.rst
new file mode 100644
index 00000000000..cbb72c37ae1
--- /dev/null
+++ b/doc/changes/dev/14139.bugfix.rst
@@ -0,0 +1 @@
+Fix bug in :func:`mne.preprocessing.find_bad_channels_maxwell` where, when ``head_pos`` was provided, every interval was processed using the head positions from the beginning of the recording rather than from its own time window, by `Christian Brodbeck`_.

From 3d0d311ff26532a5de9573b46982a06785036bb8 Mon Sep 17 00:00:00 2001
From: Christian Brodbeck <christianbrodbeck@nyu.edu>
Date: Fri, 7 Aug 2026 20:47:23 -0600
Subject: [PATCH 04/11] doc

---
 mne/preprocessing/maxwell.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py
index 277c346dc40..271335d8d4f 100644
--- a/mne/preprocessing/maxwell.py
+++ b/mne/preprocessing/maxwell.py
@@ -896,11 +896,11 @@ def get_decomp_by_offset(self, offset):
         return op_sss, op_in, op_resid
 
     def initialize(self, get_decomp, dev_head_t, S_recon, start=0):
-        """Secondary initialization."""
-        # Head positions are indexed relative to the start of the recording, so a
-        # segment that does not begin there (find_bad_channels_maxwell processes
-        # one interval at a time) must say where it starts, otherwise it would be
-        # compensated with the positions from the beginning of the recording.
+        """Secondary initialization.
+
+        :attr:`self.pos` is indexed relative to the start of the recording;
+        ``start`` adds an index offset when processing data in chunks.
+        """
         self.start = start
         self.smooth = _Interp2(
             self.pos[1],
@@ -919,7 +919,7 @@ def initialize(self, get_decomp, dev_head_t, S_recon, start=0):
 
     def get_avg_op(self, *, start, stop):
         """Apply an average transformation over the next interval."""
-        # _COLA counts from the start of the segment, self.pos from the recording
+        # Start and stop are relative to the start set at .initialize()
         start, stop = start + self.start, stop + self.start
         n_positions, avg_quat = _trans_lims(self.pos, start, stop)[1:]
         if not np.allclose(avg_quat, self.last_avg_quat, atol=1e-7):

From 8c22c0a8f0a619ccb179bf451a5ac423efb73065 Mon Sep 17 00:00:00 2001
From: Christian Brodbeck <christianbrodbeck@nyu.edu>
Date: Fri, 7 Aug 2026 21:03:43 -0600
Subject: [PATCH 05/11] fix changelog filename

---
 doc/changes/dev/{14139.bugfix.rst => 14142.bugfix.rst} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename doc/changes/dev/{14139.bugfix.rst => 14142.bugfix.rst} (100%)

diff --git a/doc/changes/dev/14139.bugfix.rst b/doc/changes/dev/14142.bugfix.rst
similarity index 100%
rename from doc/changes/dev/14139.bugfix.rst
rename to doc/changes/dev/14142.bugfix.rst

From 91dd47461410d49ad67918db40754acf0f2d200d Mon Sep 17 00:00:00 2001
From: Christian Brodbeck <christianbrodbeck@nyu.edu>
Date: Fri, 7 Aug 2026 21:10:58 -0600
Subject: [PATCH 06/11] fix error message

---
 mne/preprocessing/maxwell.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py
index 271335d8d4f..d69403dfb42 100644
--- a/mne/preprocessing/maxwell.py
+++ b/mne/preprocessing/maxwell.py
@@ -771,8 +771,8 @@ def _run_maxwell_filter(
     if not 0.0 < st_duration <= max_samps + 1.0:
         raise ValueError(
             f"st_duration ({st_duration / sfreq:0.1f}s) must be between 0 and the "
-            "longest contiguous duration of the data "
-            "({max_samps / sfreq:0.1f}s)."
+            f"longest contiguous duration of the data "
+            f"({max_samps / sfreq:0.1f}s)."
         )
 
     # This must be initialized inside _run_maxwell_filter because

From 6fb35e6f04a4aea55c3a0e6e05307622481489f7 Mon Sep 17 00:00:00 2001
From: Christian Brodbeck <christianbrodbeck@nyu.edu>
Date: Sat, 8 Aug 2026 07:02:05 -0600
Subject: [PATCH 07/11] test maxwell_filter

---
 mne/preprocessing/tests/test_maxwell.py | 62 +++++++++++++++++++++----
 1 file changed, 54 insertions(+), 8 deletions(-)

diff --git a/mne/preprocessing/tests/test_maxwell.py b/mne/preprocessing/tests/test_maxwell.py
index b02bc3ae43c..65195bc458f 100644
--- a/mne/preprocessing/tests/test_maxwell.py
+++ b/mne/preprocessing/tests/test_maxwell.py
@@ -209,6 +209,21 @@ def read_crop(fname, lims=(0, None)):
     return raw.copy().crop(*lims)
 
 
+def _linear_head_pos(raw, n_pos):
+    """Get head positions one second apart, translating steadily along z.
+
+    The steady motion is what makes reading the positions from the wrong time window
+    give a different result.
+    """
+    trans = raw.info["dev_head_t"]["trans"]
+    head_pos = np.zeros((n_pos, 10))
+    head_pos[:, 0] = raw._first_time + np.arange(float(n_pos))
+    head_pos[:, 1:4] = rot_to_quat(trans[:3, :3])
+    head_pos[:, 4:7] = trans[:3, 3]
+    head_pos[:, 6] += np.arange(n_pos) * 5e-3  # 5 mm/s
+    return head_pos
+
+
 # For backward compat and to be most like MaxFilter, we make "maxwell_filter"
 # the one that behaves like MaxFilter. _maxwell_filter is left to
 # be the advanced/better one.
@@ -1793,6 +1808,44 @@ def test_mf_skips():
     assert_allclose(data_sc, data_cs, atol=1e-20)
 
 
+@pytest.mark.slowtest
+@testing.requires_testing_data
+@pytest.mark.parametrize("st_duration", (None, 2.0))
+def test_mf_skips_head_pos(st_duration):
+    """Test that segments after a skip use their own head positions."""
+    raw = read_raw_fif(raw_fname, allow_maxshield="yes")
+    raw.pick("meg", exclude=()).crop(0, 16).load_data()
+    head_pos = _linear_head_pos(raw, 17)
+    # A 2 s skip, leaving segments of 3 s and 11 s. The latter must stay above the 10 s
+    # chunk size that st_duration=None falls back to.
+    raw.set_annotations(
+        mne.Annotations(
+            onset=[raw._first_time + 3.0],
+            duration=[2.0],
+            description=["bad_acq_skip"],
+            orig_time=raw.info["meas_date"],
+        )
+    )
+    kwargs = dict(
+        origin=(0.0, 0.0, 0.04),
+        regularize=None,
+        bad_condition="ignore",
+        st_duration=st_duration,
+    )
+    # use the default mc_interp="hann" rather than the "zero" of this module's
+    # maxwell_filter, so that each sample blends two head positions and picking the
+    # wrong one on either side of the interval shows up
+    raw_sss = _maxwell_filter_ola(raw, head_pos=head_pos, **kwargs)
+    # Processing the second segment on its own must give the same result as processing
+    # it as part of the whole recording.
+    raw_crop = raw.copy().crop(5.0).set_annotations(None)
+    raw_crop_sss = _maxwell_filter_ola(raw_crop, head_pos=head_pos[5:], **kwargs)
+    for picks in ("meg", "chpi"):  # chpi holds the head positions written back out
+        assert_allclose(
+            raw_sss.get_data(picks, tmin=5.0), raw_crop_sss.get_data(picks), atol=1e-20
+        )
+
+
 @pytest.mark.slowtest
 @testing.requires_testing_data
 @pytest.mark.parametrize(
@@ -2064,14 +2117,7 @@ def test_find_bads_maxwell_head_pos():
     raw = read_raw_fif(raw_fname, allow_maxshield="yes")
     raw.pick("meg", exclude=()).crop(0, 10)  # two 5 s intervals
     raw.load_data()
-    # One head position per second, translating steadily so that using the positions
-    # from the wrong time window gives a different result.
-    trans = raw.info["dev_head_t"]["trans"]
-    head_pos = np.zeros((11, 10))
-    head_pos[:, 0] = raw._first_time + np.arange(11.0)
-    head_pos[:, 1:4] = rot_to_quat(trans[:3, :3])
-    head_pos[:, 4:7] = trans[:3, 3]
-    head_pos[:, 6] += np.arange(11) * 5e-3  # 5 mm/s
+    head_pos = _linear_head_pos(raw, 11)
     kwargs = dict(
         origin=(0.0, 0.0, 0.04),
         regularize=None,

From 200600e400a2ebd3cee3362e6b8acf73a9583575 Mon Sep 17 00:00:00 2001
From: Christian Brodbeck <christianbrodbeck@nyu.edu>
Date: Sat, 8 Aug 2026 08:26:18 -0600
Subject: [PATCH 08/11] Add offset to COLA

---
 mne/_ola.py           | 14 +++++++++++++-
 mne/tests/test_ola.py | 28 +++++++++++++++++++++++++++-
 2 files changed, 40 insertions(+), 2 deletions(-)

diff --git a/mne/_ola.py b/mne/_ola.py
index 64c2eb72527..45f6ef05497 100644
--- a/mne/_ola.py
+++ b/mne/_ola.py
@@ -258,6 +258,11 @@ class _COLA:
         The window to use. Default is "hann".
     tol : float
         The tolerance for COLA checking.
+    offset : int
+        The index of the first sample that will be fed. ``offset`` is added to
+        the ``start`` and ``stop`` passed to ``process``.
+        Use it when processing a data segment that does not start at the
+        beginning of the recording that ``process`` is based on.
 
     Notes
     -----
@@ -292,8 +297,10 @@ def __init__(
         tol=1e-10,
         *,
         name="COLA",
+        offset=0,
         verbose=None,
     ):
+        self._offset = _ensure_int(offset, "offset")
         n_samples = _ensure_int(n_samples, "n_samples")
         n_overlap = _ensure_int(n_overlap, "n_overlap")
         n_total = _ensure_int(n_total, "n_total")
@@ -420,7 +427,12 @@ def feed(self, *datas, verbose=None, **kwargs):
                 raise RuntimeError("internal indexing error")
             start = self._store.idx
             stop = self._store.idx + this_len
-            outs = self._process(*this_proc, start=start, stop=stop, **kwargs)
+            outs = self._process(
+                *this_proc,
+                start=start + self._offset,
+                stop=stop + self._offset,
+                **kwargs,
+            )
             if self._out_buffers is None:
                 max_len = np.max(self.stops - self.starts)
                 self._out_buffers = [
diff --git a/mne/tests/test_ola.py b/mne/tests/test_ola.py
index e1d7938f3a4..e8a66ba224d 100644
--- a/mne/tests/test_ola.py
+++ b/mne/tests/test_ola.py
@@ -4,7 +4,7 @@
 
 import numpy as np
 import pytest
-from numpy.testing import assert_allclose
+from numpy.testing import assert_allclose, assert_array_equal
 
 from mne._ola import _COLA, _Interp2, _Storer
 
@@ -143,3 +143,29 @@ def processor(x, *, start, stop):
                             cola.feed(signal[..., n_input : n_input + next_len])
                             n_input += next_len
                         assert_allclose(out, signal / 2.0, atol=1e-7)
+
+
+def test_cola_offset():
+    """Test that COLA shifts the limits it hands to the processor."""
+    n_total, n_samples, n_overlap, offset = 1000, 100, 50, 4321
+    limits = list()
+
+    def processor(x, *, start, stop):
+        limits.append((start, stop))
+        return (x,)
+
+    signal = np.zeros(n_total)
+    runs = list()
+    for use_offset in (0, offset):
+        limits.clear()
+        _COLA(
+            processor,
+            np.zeros(n_total),
+            n_total,
+            n_samples,
+            n_overlap,
+            1000.0,
+            offset=use_offset,
+        ).feed(signal)
+        runs.append(np.array(limits))
+    assert_array_equal(runs[1] - offset, runs[0])

From 823ba12704e841de5ee5515edb40eb4e17e8d5e3 Mon Sep 17 00:00:00 2001
From: Christian Brodbeck <christianbrodbeck@nyu.edu>
Date: Sat, 8 Aug 2026 08:27:24 -0600
Subject: [PATCH 09/11] Fix maxwell_filter head position with
 skip_by_annotations

---
 doc/changes/dev/14142.bugfix.rst |  2 +-
 mne/preprocessing/maxwell.py     | 47 ++++++++++++++++++++------------
 2 files changed, 31 insertions(+), 18 deletions(-)

diff --git a/doc/changes/dev/14142.bugfix.rst b/doc/changes/dev/14142.bugfix.rst
index cbb72c37ae1..f3be8325cd3 100644
--- a/doc/changes/dev/14142.bugfix.rst
+++ b/doc/changes/dev/14142.bugfix.rst
@@ -1 +1 @@
-Fix bug in :func:`mne.preprocessing.find_bad_channels_maxwell` where, when ``head_pos`` was provided, every interval was processed using the head positions from the beginning of the recording rather than from its own time window, by `Christian Brodbeck`_.
+Fix bug in :func:`mne.preprocessing.find_bad_channels_maxwell` and :func:`mne.preprocessing.maxwell_filter`: data is processed in chunks, but when ``head_pos`` was provided, the head positions from the beginning of the recording was used for each chunk rather than from the chunk's time window. In ``find_bad_channels_maxwell`` this affected every interval, and in ``maxwell_filter`` every segment following a segment skipped due to ``skip_by_annotation``, by `Christian Brodbeck`_.
diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py
index d69403dfb42..a7f3bb67b23 100644
--- a/mne/preprocessing/maxwell.py
+++ b/mne/preprocessing/maxwell.py
@@ -741,7 +741,7 @@ def _run_maxwell_filter(
     st_fixed,
     st_overlap,
     mc,
-    mc_start=0,
+    raw_start=0,
 ):
     # Eventually find_bad_channels_maxwell could be sped up by moving this
     # outside the loop (e.g., in the prep function) but regularization depends
@@ -777,11 +777,15 @@ def _run_maxwell_filter(
 
     # This must be initialized inside _run_maxwell_filter because
     # find_bad_channels_maxwell modifies good_mask
-    mc.initialize(_get_this_decomp_trans, info["dev_head_t"], S_recon, mc_start)
+    mc.initialize(_get_this_decomp_trans, info["dev_head_t"], S_recon)
     update_kwargs.update(reg_moments=mc.reg_moments_0)
 
     # Process each valid block of data separately
     for onset, end in zip(onsets, ends):
+        # head positions are indexed relative to the recording, but onset and end are
+        # relative to raw, which can itself be a chunk of the recording
+        segment_start = raw_start + onset
+        mc.set_start(segment_start)
         n = end - onset
         assert n > 0
         tsss_valid = n >= st_duration
@@ -809,6 +813,7 @@ def _run_maxwell_filter(
             sfreq,
             window,
             name="tSSS-COLA",
+            offset=segment_start,
         )
 
         # Generate time points to break up data into equal-length windows
@@ -895,13 +900,26 @@ def get_decomp_by_offset(self, offset):
         op_resid -= np.dot(S_decomp[:, n_use_in:], pS_decomp[n_use_in:])
         return op_sss, op_in, op_resid
 
-    def initialize(self, get_decomp, dev_head_t, S_recon, start=0):
+    def initialize(self, get_decomp, dev_head_t, S_recon):
         """Secondary initialization.
 
-        :attr:`self.pos` is indexed relative to the start of the recording;
-        ``start`` adds an index offset when processing data in chunks.
+        Call :meth:`set_start` before feeding data.
         """
-        self.start = start
+        _, _, pS_decomp, self.reg_moments_0, _ = get_decomp(dev_head_t, t=0.0)
+        self.n_good = pS_decomp.shape[1]
+        self.S_recon = S_recon
+        self.get_decomp = get_decomp
+        # For the average passes
+        self.last_avg_quat = np.nan * np.ones(6)
+
+    def set_start(self, start):
+        """Position at a given sample of the recording, to process a segment from there.
+
+        :attr:`self.pos` is indexed relative to the start of the recording, so a segment
+        that does not begin there has to be told where it does, both to read the right
+        head positions and to resume interpolation with the right phase.
+        """
+        self.offset = start
         self.smooth = _Interp2(
             self.pos[1],
             self.get_decomp_by_offset,
@@ -909,18 +927,13 @@ def initialize(self, get_decomp, dev_head_t, S_recon, start=0):
             name="MC",
             start=start,
         )
-        _, _, pS_decomp, self.reg_moments_0, _ = get_decomp(dev_head_t, t=0.0)
-        self.n_good = pS_decomp.shape[1]
-        self.S_recon = S_recon
-        self.offset = start
-        self.get_decomp = get_decomp
-        # For the average passes
-        self.last_avg_quat = np.nan * np.ones(6)
 
     def get_avg_op(self, *, start, stop):
-        """Apply an average transformation over the next interval."""
-        # Start and stop are relative to the start set at .initialize()
-        start, stop = start + self.start, stop + self.start
+        """Apply an average transformation over the next interval.
+
+        ``start`` and ``stop`` are relative to the start of the recording, like
+        :attr:`self.offset`.
+        """
         n_positions, avg_quat = _trans_lims(self.pos, start, stop)[1:]
         if not np.allclose(avg_quat, self.last_avg_quat, atol=1e-7):
             self.last_avg_quat = avg_quat
@@ -2954,7 +2967,7 @@ def find_bad_channels_maxwell(
             chunk_raw._data[:] = orig_data
             delta = chunk_raw.get_data(these_picks)
             with use_log_level(_verbose_safe_false()):
-                _run_maxwell_filter(chunk_raw, copy=False, mc_start=start, **params)
+                _run_maxwell_filter(chunk_raw, copy=False, raw_start=start, **params)
 
             if n_iter == 1 and len(chunk_flats):
                 logger.info(

From a9f200d5681232a6c52e6f1ef812bfc59241251f Mon Sep 17 00:00:00 2001
From: Christian Brodbeck <christianbrodbeck@nyu.edu>
Date: Sat, 8 Aug 2026 09:24:43 -0600
Subject: [PATCH 10/11] Consistent variable names

---
 mne/_ola.py                  | 36 ++++++++++++++++++------------------
 mne/preprocessing/maxwell.py | 26 +++++++++++++-------------
 mne/tests/test_ola.py        | 12 ++++++------
 3 files changed, 37 insertions(+), 37 deletions(-)

diff --git a/mne/_ola.py b/mne/_ola.py
index 45f6ef05497..22c3bad91a5 100644
--- a/mne/_ola.py
+++ b/mne/_ola.py
@@ -23,11 +23,12 @@ class _Interp2:
         arrays that must be interpolated.
     interp : str
         Can be 'zero', 'linear', 'hann', or 'cos2' (same as hann).
-    start : int
-        The control-point index at which feeding begins, for processing a segment
-        that does not start at the first control point. Feeding ``n`` points then
-        yields exactly what a continuous pass would yield for ``start`` to
-        ``start + n``, interpolation phase included.
+    offset : int
+        The position of the first point that will be fed, in the same units as
+        ``control_points``. Use it to process a segment that does not start at the
+        beginning of the signal the control points refer to: feeding ``n`` points then
+        yields exactly what a continuous pass would yield for ``offset`` to
+        ``offset + n``, interpolation phase included.
 
     Notes
     -----
@@ -48,7 +49,7 @@ class _Interp2:
     """
 
     def __init__(
-        self, control_points, values, interp="hann", *, name="Interp2", start=0
+        self, control_points, values, interp="hann", *, name="Interp2", offset=0
     ):
         # set up interpolation
         self.control_points = np.array(control_points, int).ravel()
@@ -83,10 +84,10 @@ def val(pt):
             values = val
         self.values = values
         self.n_last = None
-        self._position = start
-        # The last control point at or before start is the one in effect there, so
+        self._position = offset
+        # The last control point at or before offset is the one in effect there, so
         # feeding resumes mid-interval with the correct interpolation phase.
-        left_idx = np.searchsorted(self.control_points, start, "right") - 1
+        left_idx = np.searchsorted(self.control_points, offset, "right") - 1
         self._left_idx = max(left_idx, 0)
         self._left = self._right = self._use_interp = None
         self.name = name
@@ -259,10 +260,9 @@ class _COLA:
     tol : float
         The tolerance for COLA checking.
     offset : int
-        The index of the first sample that will be fed. ``offset`` is added to
-        the ``start`` and ``stop`` passed to ``process``.
-        Use it when processing a data segment that does not start at the
-        beginning of the recording that ``process`` is based on.
+        The index of the first sample that will be fed. Use it to process a segment
+        that does not start at the beginning of the signal that ``process`` is based
+        on: ``offset`` is added to the ``start`` and ``stop`` handed to ``process``.
 
     Notes
     -----
@@ -410,12 +410,12 @@ def feed(self, *datas, verbose=None, **kwargs):
                 this_window = np.pad(
                     self._window, (0, this_len - len(this_window)), "constant"
                 )
-                for offset in range(self._step, len(this_window), self._step):
-                    n_use = len(this_window) - offset
-                    this_window[offset:] += self._window[:n_use]
+                for shift in range(self._step, len(this_window), self._step):
+                    n_use = len(this_window) - shift
+                    this_window[shift:] += self._window[:n_use]
             if self._idx == 0:
-                for offset in range(self._n_samples - self._step, 0, -self._step):
-                    this_window[:offset] += self._window[-offset:]
+                for n_use in range(self._n_samples - self._step, 0, -self._step):
+                    this_window[:n_use] += self._window[-n_use:]
             this_proc = [in_[..., :this_len].copy() for in_ in self._in_buffers]
             logger.debug(
                 f"    * {self.name}[:] Processing {start}:{stop} "
diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py
index a7f3bb67b23..1304c52dc25 100644
--- a/mne/preprocessing/maxwell.py
+++ b/mne/preprocessing/maxwell.py
@@ -741,7 +741,7 @@ def _run_maxwell_filter(
     st_fixed,
     st_overlap,
     mc,
-    raw_start=0,
+    raw_offset=0,
 ):
     # Eventually find_bad_channels_maxwell could be sped up by moving this
     # outside the loop (e.g., in the prep function) but regularization depends
@@ -784,8 +784,8 @@ def _run_maxwell_filter(
     for onset, end in zip(onsets, ends):
         # head positions are indexed relative to the recording, but onset and end are
         # relative to raw, which can itself be a chunk of the recording
-        segment_start = raw_start + onset
-        mc.set_start(segment_start)
+        segment_offset = raw_offset + onset
+        mc.set_offset(segment_offset)
         n = end - onset
         assert n > 0
         tsss_valid = n >= st_duration
@@ -813,7 +813,7 @@ def _run_maxwell_filter(
             sfreq,
             window,
             name="tSSS-COLA",
-            offset=segment_start,
+            offset=segment_offset,
         )
 
         # Generate time points to break up data into equal-length windows
@@ -903,7 +903,7 @@ def get_decomp_by_offset(self, offset):
     def initialize(self, get_decomp, dev_head_t, S_recon):
         """Secondary initialization.
 
-        Call :meth:`set_start` before feeding data.
+        Call :meth:`set_offset` before feeding data.
         """
         _, _, pS_decomp, self.reg_moments_0, _ = get_decomp(dev_head_t, t=0.0)
         self.n_good = pS_decomp.shape[1]
@@ -912,20 +912,20 @@ def initialize(self, get_decomp, dev_head_t, S_recon):
         # For the average passes
         self.last_avg_quat = np.nan * np.ones(6)
 
-    def set_start(self, start):
-        """Position at a given sample of the recording, to process a segment from there.
+    def set_offset(self, offset):
+        """Position at the given sample of the recording to process a segment there.
 
-        :attr:`self.pos` is indexed relative to the start of the recording, so a segment
-        that does not begin there has to be told where it does, both to read the right
-        head positions and to resume interpolation with the right phase.
+        ``pos`` is indexed relative to the start of the recording, so a segment that
+        does not begin there has to be told where it does, both to read the right head
+        positions and to resume interpolation with the right phase.
         """
-        self.offset = start
+        self.offset = offset
         self.smooth = _Interp2(
             self.pos[1],
             self.get_decomp_by_offset,
             interp=self.interp,
             name="MC",
-            start=start,
+            offset=offset,
         )
 
     def get_avg_op(self, *, start, stop):
@@ -2967,7 +2967,7 @@ def find_bad_channels_maxwell(
             chunk_raw._data[:] = orig_data
             delta = chunk_raw.get_data(these_picks)
             with use_log_level(_verbose_safe_false()):
-                _run_maxwell_filter(chunk_raw, copy=False, raw_start=start, **params)
+                _run_maxwell_filter(chunk_raw, copy=False, raw_offset=start, **params)
 
             if n_iter == 1 and len(chunk_flats):
                 logger.info(
diff --git a/mne/tests/test_ola.py b/mne/tests/test_ola.py
index e8a66ba224d..0e412abba0e 100644
--- a/mne/tests/test_ola.py
+++ b/mne/tests/test_ola.py
@@ -85,18 +85,18 @@ def test_interp_2pt():
 
 
 @pytest.mark.parametrize("interp", ("zero", "linear", "hann"))
-def test_interp_2pt_start(interp):
+def test_interp_2pt_offset(interp):
     """Test that starting mid-stream matches a continuous pass."""
-    # control points deliberately uneven, so that most starts land mid-interval
+    # control points deliberately uneven, so that most offsets land mid-interval
     control_points = [0, 100, 240, 390, 520, 680, 800]
     values = np.array(control_points, float)
     want = _Interp2(control_points, values, interp).feed(900)[0]
-    for start in range(0, 900, 7):
-        interper = _Interp2(control_points, values, interp, start=start)
-        assert_allclose(interper.feed(900 - start)[0], want[start:], atol=1e-12)
+    for offset in range(0, 900, 7):
+        interper = _Interp2(control_points, values, interp, offset=offset)
+        assert_allclose(interper.feed(900 - offset)[0], want[offset:], atol=1e-12)
     # state must carry across feeds of differing size, as when a segment is read
     # in buffer-sized blocks
-    interper = _Interp2(control_points, values, interp, start=300)
+    interper = _Interp2(control_points, values, interp, offset=300)
     out = np.concatenate([interper.feed(n)[0] for n in (70, 130, 100, 300)])
     assert_allclose(out, want[300:900], atol=1e-12)
 

From 454535e0be1b9b03d715266b5c6142cbba645f88 Mon Sep 17 00:00:00 2001
From: Christian Brodbeck <christianbrodbeck@nyu.edu>
Date: Sat, 8 Aug 2026 09:45:10 -0600
Subject: [PATCH 11/11] Small fixes

---
 doc/changes/dev/14142.bugfix.rst        | 2 +-
 mne/_ola.py                             | 2 +-
 mne/preprocessing/maxwell.py            | 6 ++++--
 mne/preprocessing/tests/test_maxwell.py | 2 +-
 4 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/doc/changes/dev/14142.bugfix.rst b/doc/changes/dev/14142.bugfix.rst
index f3be8325cd3..e71b445c182 100644
--- a/doc/changes/dev/14142.bugfix.rst
+++ b/doc/changes/dev/14142.bugfix.rst
@@ -1 +1 @@
-Fix bug in :func:`mne.preprocessing.find_bad_channels_maxwell` and :func:`mne.preprocessing.maxwell_filter`: data is processed in chunks, but when ``head_pos`` was provided, the head positions from the beginning of the recording was used for each chunk rather than from the chunk's time window. In ``find_bad_channels_maxwell`` this affected every interval, and in ``maxwell_filter`` every segment following a segment skipped due to ``skip_by_annotation``, by `Christian Brodbeck`_.
+Fix bug in :func:`mne.preprocessing.find_bad_channels_maxwell` and :func:`mne.preprocessing.maxwell_filter`: data is processed in chunks, but when ``head_pos`` was provided, the head positions from the beginning of the recording were used for each chunk rather than from the chunk's time window. In ``find_bad_channels_maxwell`` this affected every interval after the first, and in ``maxwell_filter`` every segment following a segment skipped due to ``skip_by_annotation``, by `Christian Brodbeck`_.
diff --git a/mne/_ola.py b/mne/_ola.py
index 22c3bad91a5..77176148825 100644
--- a/mne/_ola.py
+++ b/mne/_ola.py
@@ -84,7 +84,7 @@ def val(pt):
             values = val
         self.values = values
         self.n_last = None
-        self._position = offset
+        self._position = offset = _ensure_int(offset, "offset")
         # The last control point at or before offset is the one in effect there, so
         # feeding resumes mid-interval with the correct interpolation phase.
         left_idx = np.searchsorted(self.control_points, offset, "right") - 1
diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py
index 1304c52dc25..43de3e778f3 100644
--- a/mne/preprocessing/maxwell.py
+++ b/mne/preprocessing/maxwell.py
@@ -741,7 +741,7 @@ def _run_maxwell_filter(
     st_fixed,
     st_overlap,
     mc,
-    raw_offset=0,
+    raw_offset=0,  # time offset of ``raw`` relative to ``mc``
 ):
     # Eventually find_bad_channels_maxwell could be sped up by moving this
     # outside the loop (e.g., in the prep function) but regularization depends
@@ -911,6 +911,7 @@ def initialize(self, get_decomp, dev_head_t, S_recon):
         self.get_decomp = get_decomp
         # For the average passes
         self.last_avg_quat = np.nan * np.ones(6)
+        self.smooth = None  # set_offset positions us in the recording
 
     def set_offset(self, offset):
         """Position at the given sample of the recording to process a segment there.
@@ -932,7 +933,7 @@ def get_avg_op(self, *, start, stop):
         """Apply an average transformation over the next interval.
 
         ``start`` and ``stop`` are relative to the start of the recording, like
-        :attr:`self.offset`.
+        ``offset``.
         """
         n_positions, avg_quat = _trans_lims(self.pos, start, stop)[1:]
         if not np.allclose(avg_quat, self.last_avg_quat, atol=1e-7):
@@ -957,6 +958,7 @@ def get_avg_op(self, *, start, stop):
         return self.op_in_avg, self.op_resid_avg, n_positions
 
     def feed(self, data, good_mask, st_only):
+        assert self.smooth is not None  # set_offset must be called first
         n_samp = data.shape[1]
         pos_data, n_pos = _trans_lims(
             self.pos, self.offset, self.offset + data.shape[-1]
diff --git a/mne/preprocessing/tests/test_maxwell.py b/mne/preprocessing/tests/test_maxwell.py
index 65195bc458f..30983191766 100644
--- a/mne/preprocessing/tests/test_maxwell.py
+++ b/mne/preprocessing/tests/test_maxwell.py
@@ -212,7 +212,7 @@ def read_crop(fname, lims=(0, None)):
 def _linear_head_pos(raw, n_pos):
     """Get head positions one second apart, translating steadily along z.
 
-    The steady motion is what makes reading the positions from the wrong time window
+    The steady motion makes reading the positions from the wrong time window
     give a different result.
     """
     trans = raw.info["dev_head_t"]["trans"]