Merge branch 'main' into perf/mef-session-cache

Author: bruAristimunha

doc/changes/dev/14250.newfeature.rst

@@ -0,0 +1 @@
+Speed up :func:`mne.io.read_raw_egi` on simple-binary files by reading the event channels in blocks instead of one sample at a time, by `Bruno Aristimunha`_.

doc/changes/dev/14251.newfeature.rst

@@ -0,0 +1 @@
+Speed up :func:`mne.io.read_raw_persyst` and :func:`mne.io.read_raw_nihon` by decoding data in cache-sized blocks rather than materializing the whole request, by `Bruno Aristimunha`_.

doc/changes/dev/14252.newfeature.rst

@@ -0,0 +1 @@
+Speed up :func:`mne.io.read_raw_fif` on integer-format files by normalizing the byte order of each data buffer before calibrating it, by `Bruno Aristimunha`_.

mne/_fiff/utils.py

@@ -88,7 +88,19 @@ def _mult_cal_one(data_view, one, idx, cals, mult):
             # (was three passes plus a full float64 temporary).
             # Benchmark (128 ch x 1024 samples): ~85 -> ~30 us per call
             # on BrainVision/FIF window reads.
-            np.multiply(one[idx], cals.reshape(-1, 1), out=data_view, casting="unsafe")
+            one = one[idx]
+            swapped_ints = not one.dtype.isnative and one.dtype.kind in "iu"
+            if swapped_ints and not one.flags.c_contiguous:
+                # FIFF stores sample-major, so `one` is the transpose of a
+                # big-endian integer tag and consecutive samples of a channel
+                # sit a row apart. np.multiply has no fast loop for input that
+                # is byte-swapped *and* strided, so it swaps element by element
+                # while it also casts and transposes; swapping up front in one
+                # pass is cheaper. Floats are excluded because there the swap
+                # costs about what it saves, and contiguous input already has a
+                # fast loop.
+                one = one.astype(one.dtype.newbyteorder("="))
+            np.multiply(one, cals.reshape(-1, 1), out=data_view, casting="unsafe")
         else:
             one = np.asarray(one, dtype=data_view.dtype)
             np.take(one, idx, axis=0, out=data_view)

mne/io/egi/egi.py

@@ -79,15 +79,27 @@ def my_fread(*x, **y):
     return info
 
 
+# read whole frames a few MB at a time rather than the whole recording at once
+_EVENT_BLOCK_BYTES = 4 * 1024**2
+
+
 def _read_events(fid, info):
     """Read events."""
-    events = np.zeros([info["n_events"], info["n_segments"] * info["n_samples"]])
+    n_samples = info["n_samples"]
+    # Each sample is one frame of n_channels data values followed by n_events
+    # event values, so read whole frames and keep the event rows. Seeking past
+    # the data channels instead costs a seek and a read per sample, which
+    # dominates the time to open a long recording.
+    n_rows = info["n_channels"] + info["n_events"]
+    events = np.zeros([info["n_events"], info["n_segments"] * n_samples])
     fid.seek(36 + info["n_events"] * 4, 0)  # skip header
-    for si in range(info["n_samples"]):
-        # skip data channels
-        fid.seek(info["n_channels"] * info["dtype"].itemsize, 1)
-        # read event channels
-        events[:, si] = np.fromfile(fid, info["dtype"], info["n_events"])
+    n_block = max(1, _EVENT_BLOCK_BYTES // (n_rows * info["dtype"].itemsize))
+    for start in range(0, n_samples, n_block):
+        n_read = min(n_block, n_samples - start)
+        frames = np.fromfile(fid, info["dtype"], n_read * n_rows)
+        events[:, start : start + n_read] = frames.reshape(n_read, n_rows).T[
+            info["n_channels"] :
+        ]
     return events
 
 

mne/io/nihon/nihon.py

@@ -418,6 +418,10 @@ def _map_ch_to_specs(ch_name, chan_labels_upper):
     return out
 
 
+# decode in cache-sized blocks rather than one huge one (1.8x on a 106 MB file)
+_BLOCK_BYTES = 1024**2
+
+
 @fill_doc
 class RawNihon(BaseRaw):
     """Raw object from a Nihon Kohden EEG file.
@@ -566,13 +570,20 @@ def _read_segment_file(self, data, idx, fi, start, stop, cals, mult):
                 rel_start = start - ends[start_block - 1]
             start_offset = datastart + rel_start * n_channels * 2
 
+            # Decode a few MB at a time: each step below builds a temporary the
+            # size of the block, so reading the whole request at once pushes
+            # them all out of cache.
+            n_times = stop - start
+            n_block = max(1, _BLOCK_BYTES // 2 // n_channels)
             with open(self.filenames[fi], "rb") as fid:
-                to_read = (stop - start) * n_channels
                 fid.seek(start_offset)
-                block_data = np.fromfile(fid, "<u2", to_read) + 0x8000
-                block_data = block_data.astype(np.int16)
-                block_data = block_data.reshape(n_channels, -1, order="F")
-                block_data = block_data[:-1] * cal  # cast to float64
-                block_data += offsets
-                block_data *= gains
-                _mult_cal_one(data, block_data, idx, cals, mult)
+                for sample_start in range(0, n_times, n_block):
+                    n_read = min(n_block, n_times - sample_start)
+                    block_data = np.fromfile(fid, "<u2", n_read * n_channels) + 0x8000
+                    block_data = block_data.astype(np.int16)
+                    block_data = block_data.reshape(n_channels, -1, order="F")
+                    block_data = block_data[:-1] * cal  # cast to float64
+                    block_data += offsets
+                    block_data *= gains
+                    data_view = data[:, sample_start : sample_start + n_read]
+                    _mult_cal_one(data_view, block_data, idx, cals, mult)

mne/io/persyst/persyst.py

@@ -13,7 +13,7 @@
 
 from ..._fiff.constants import FIFF
 from ..._fiff.meas_info import create_info
-from ..._fiff.utils import _mult_cal_one
+from ..._fiff.utils import _read_segments_file
 from ...annotations import Annotations
 from ...utils import _check_fname, fill_doc, logger, verbose, warn
 from ..base import BaseRaw
@@ -55,6 +55,11 @@ def read_raw_persyst(
     return RawPersyst(fname, preload, verbose)
 
 
+# read in cache-sized blocks rather than one huge one (2.3x on a 107 MB file);
+# see _read_segments_file() for why a smaller block is faster
+_BLOCK_BYTES = 16 * 1024**2
+
+
 @fill_doc
 class RawPersyst(BaseRaw):
     """Raw object from a Persyst file.
@@ -267,31 +272,19 @@ def _read_segment_file(self, data, idx, fi, start, stop, cals, mult):
         binary files. In addition, it stores the calibration to convert
         data to uV in the lay file.
         """
-        dtype = self._raw_extras[fi]["dtype"]
-        n_chs = self._raw_extras[fi]["n_chs"]
-        dat_fname = self.filenames[fi]
-
-        # compute samples count based on start and stop
-        time_length_samps = stop - start
-
-        # read data from .dat file into array of correct size, then calibrate
-        # records = recnum rows x inf columns
-        count = time_length_samps * n_chs
-
-        # seek the dat file
-        with open(dat_fname, "rb") as dat_file_ID:
-            # allow offset to occur
-            dat_file_ID.seek(n_chs * dtype.itemsize * start, 1)
-
-            # read in the actual record starting at possibly offset
-            record = np.fromfile(dat_file_ID, dtype=dtype, count=count)
-
-        # chs * rows
-        # cast as float32; more than enough precision
-        record = np.reshape(record, (n_chs, -1), order="F").astype(np.float32)
-
-        # calibrate to convert to V and handle mult
-        _mult_cal_one(data, record, idx, cals, mult)
+        _read_segments_file(
+            self,
+            data,
+            idx,
+            fi,
+            start,
+            stop,
+            cals,
+            mult,
+            dtype=self._raw_extras[fi]["dtype"],
+            n_channels=self._raw_extras[fi]["n_chs"],
+            max_block_bytes=_BLOCK_BYTES,
+        )
 
 
 def _get_subjectinfo(patient_dict):