From da0b2eb0d813649354c500273e6103078838fb5e Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 28 Oct 2026 18:24:29 +0000 Subject: [PATCH 1/2] Ingest: write the epoch's data, as NDI-matlab does (+ within-epoch progress) Ingestion made each epoch's documents without any data, so nothing could be read from an ingested epoch ("channel_list.bin not found"). - ndi_daq_reader_mfdaq.ingest_epochfiles (NDI-matlab 6ff0264e5, #1089): the channel list (channel_list.bin, MATLAB's saveStructArray format) and the data in compressed segments -- analog in the reader's underlying type, digital, time, event/marker/text -- plus ProgressFcn, the fraction of the epoch done after each segment. - ingestEpochDocuments (b192ae605): fractional ProgressFcn while an epoch is underway, also from worker processes (Manager queue). - Intan reader: underlying_datatype (uint16, as MATLAB); aux inputs read at their own sample rate; digital channels as NDR numbers them (16 were requested and NDR refused); reads past the end are not padded. - Metadata readers store readmetadata's parameters (data.bin) and readmetadata_ingested reads them back. - mfdaq_epoch_channel.writeToFile writes MATLAB's format; readFromFile still reads the older JSON files. Needs ndi-compress's compress_ephys format options (Waltham-Data-Science/NDI-compress-python#4). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_011v7R5rnKYwncJkrNZRXhtw --- CHANGELOG.md | 18 + src/ndi/daq/fun/ingest_epoch_documents.py | 113 ++++++- src/ndi/daq/fun/ndi_matlab_python_bridge.yaml | 13 +- src/ndi/daq/metadatareader/__init__.py | 42 ++- .../ndi_matlab_python_bridge.yaml | 12 +- src/ndi/daq/mfdaq.py | 313 ++++++++++++++++++ src/ndi/daq/ndi_matlab_python_bridge.yaml | 48 ++- src/ndi/daq/reader/mfdaq/intan.py | 93 ++++-- .../mfdaq/ndi_matlab_python_bridge.yaml | 36 +- src/ndi/file/type/mfdaq_epoch_channel.py | 27 +- .../file/type/ndi_matlab_python_bridge.yaml | 5 + tests/test_ingest_data_writer.py | 273 +++++++++++++++ 12 files changed, 934 insertions(+), 59 deletions(-) create mode 100644 tests/test_ingest_data_writer.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 3c2fe34f..7247647a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,24 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/). ### Added +- **Ingestion stores the data** (mirroring NDI-matlab's + `ndi.daq.reader.mfdaq/ingest_epochfiles`, with VH-Lab/NDI-matlab#1089's + within-epoch progress). Ingesting an epoch used to make its documents + without any data, so an ingested epoch could not be read (`channel_list.bin + not found`). `ndi_daq_reader_mfdaq.ingest_epochfiles` now writes the channel + list and the epoch's data in compressed segments as MATLAB does -- analog + data in the reader's underlying type (Intan: uint16, via the new + `ndi_daq_reader_mfdaq_intan.underlying_datatype`), digital, time, and + event/marker/text channels -- and metadata readers store their parameters + (`data.bin`), which `readmetadata_ingested` now reads back. The channel list + is written in MATLAB's `saveStructArray` format so NDI-matlab can read it. + `ingest_epochfiles` takes a `ProgressFcn`, and `ingestEpochDocuments` + reports fractional progress as segments are written, in worker processes + too. The Intan reader now reads auxiliary inputs at their own sample rate + and its digital channels as NDR numbers them (both as MATLAB), and no + longer pads reads that run past the end of the file. Needs ndi-compress with + Waltham-Data-Science/NDI-compress-python#4. + - **Ingestion progress bars and parallel per-epoch work** (mirroring VH-Lab/NDI-matlab#1079). `ndi_session.ingest()` now shows progress through `ndi_gui_component_ProgressBarWindow` -- one bar for the session and one per diff --git a/src/ndi/daq/fun/ingest_epoch_documents.py b/src/ndi/daq/fun/ingest_epoch_documents.py index 235897fe..12678191 100644 --- a/src/ndi/daq/fun/ingest_epoch_documents.py +++ b/src/ndi/daq/fun/ingest_epoch_documents.py @@ -19,11 +19,13 @@ from __future__ import annotations import copy +import functools import multiprocessing import pickle import warnings from collections.abc import Callable -from concurrent.futures import ProcessPoolExecutor, as_completed +from concurrent.futures import FIRST_COMPLETED, ProcessPoolExecutor, wait +from queue import Empty from typing import Any from pydantic import ConfigDict, validate_call @@ -37,7 +39,7 @@ def ingestEpochDocuments( # noqa: N802 (MATLAB mirror) daqmetadatareaders: list[Any], epochFiles: list[list[str]], epochIds: list[str], - ProgressFcn: Callable[[int, int], Any] | None = None, + ProgressFcn: Callable[[float, int], Any] | None = None, NumberOfWorkers: int | None = None, ) -> tuple[list[list[Any]], bool]: """Build the ingested-data documents for a set of epochs. @@ -67,7 +69,10 @@ def ingestEpochDocuments( # noqa: N802 (MATLAB mirror) epochFiles: One list of file names per epoch. epochIds: One epoch id per epoch. ProgressFcn: Called as ``ProgressFcn(numDone, numTotal)`` in the calling - process after each epoch completes. + process as epochs progress. numDone counts epochs and is fractional + while an epoch is underway: readers whose ingest_epochfiles is + ndi_daq_reader_mfdaq's report each segment they write. It is called + only when numDone changes, so each whole number is reported once. NumberOfWorkers: Overrides the preferences for this call (passed to ndi.fun.parallelWorkers); 0 forces serial. @@ -83,23 +88,36 @@ def ingestEpochDocuments( # noqa: N802 (MATLAB mirror) be pickled falls back to serial with a warning rather than failing. Completed epochs are reported as they finish, so ProgressFcn may see epochs out of order, but the returned list is always in epoch order. + Within-epoch progress travels from the workers on a + ``multiprocessing.Manager`` queue (MATLAB uses a DataQueue) that the + caller drains while it waits. """ num_epochs = len(epochFiles) if len(epochIds) != num_epochs: raise ValueError("EPOCHFILES and EPOCHIDS must have the same number of entries.") docs_per_epoch: list[list[Any]] = [[] for _ in range(num_epochs)] - done = 0 - - def tick() -> None: - nonlocal done - done += 1 + # per-epoch fractions done + frac = [0.0] * num_epochs + + def update_progress(k: int, fraction: float) -> None: + # Fractions only move forward, and ProgressFcn is called only when the + # total changes, so a reader's own final report and the epoch's + # completion do not count the epoch twice. + new_fraction = max(frac[k], min(max(float(fraction), 0.0), 1.0)) + if new_fraction == frac[k]: + return + frac[k] = new_fraction if ProgressFcn is not None: try: - ProgressFcn(done, num_epochs) + ProgressFcn(sum(frac), num_epochs) except Exception as exc: # noqa: BLE001 - progress must never break an ingest warnings.warn(f"Progress reporting failed: {exc}", stacklevel=3) + # only readers using ndi_daq_reader_mfdaq's ingest_epochfiles take a + # ProgressFcn; any other reader reports per epoch only + within_epoch = ProgressFcn is not None and _reader_reports_progress(daqreader) + number_of_workers = 0 if num_epochs > 1: number_of_workers = _resolve_workers(NumberOfWorkers) @@ -113,26 +131,50 @@ def tick() -> None: if number_of_workers > 0: reader, mreaders = payload ctx = multiprocessing.get_context("spawn") + manager = ctx.Manager() if within_epoch else None + queue = manager.Queue() if manager is not None else None + + def drain() -> None: + while queue is not None: + try: + k, fraction = queue.get_nowait() + except Empty: + return + update_progress(k, fraction) + executor = ProcessPoolExecutor( max_workers=min(number_of_workers, num_epochs), mp_context=ctx ) try: futures = { - executor.submit(_ingestOneEpoch, reader, mreaders, ef, eid): k + executor.submit(_ingestOneEpoch, reader, mreaders, ef, eid, queue, k): k for k, (ef, eid) in enumerate(zip(epochFiles, epochIds)) } - for future in as_completed(futures): - docs_per_epoch[futures[future]] = future.result() - tick() + pending = set(futures) + while pending: + finished, pending = wait(pending, timeout=0.25, return_when=FIRST_COMPLETED) + drain() + for future in finished: + docs_per_epoch[futures[future]] = future.result() + update_progress(futures[future], 1) + drain() except BaseException: executor.shutdown(wait=True, cancel_futures=True) raise + finally: + if manager is not None: + manager.shutdown() executor.shutdown(wait=True) return docs_per_epoch, True for k, (ef, eid) in enumerate(zip(epochFiles, epochIds)): - docs_per_epoch[k] = _ingestOneEpoch(daqreader, daqmetadatareaders, ef, eid) - tick() + epoch_progress = None + if within_epoch: + epoch_progress = functools.partial(update_progress, k) + docs_per_epoch[k] = _ingestOneEpoch( + daqreader, daqmetadatareaders, ef, eid, epoch_progress=epoch_progress + ) + update_progress(k, 1) return docs_per_epoch, False @@ -141,12 +183,27 @@ def tick() -> None: def _ingestOneEpoch( # noqa: N802 (mirrors MATLAB's local function) - daqreader: Any, daqmetadatareaders: list[Any], ef: list[str], epoch_id: str + daqreader: Any, + daqmetadatareaders: list[Any], + ef: list[str], + epoch_id: str, + queue: Any = None, + k: int = 0, + epoch_progress: Callable[[float], Any] | None = None, ) -> list[Any]: - """Create the documents for a single epoch (no database access).""" + """Create the documents for a single epoch (no database access). + + In a worker process, progress is sent as ``(k, fraction)`` on QUEUE; in + the calling process it goes to EPOCH_PROGRESS. + """ + if queue is not None: + epoch_progress = functools.partial(_send_progress, queue, k) docs: list[Any] = [] if daqreader is not None: - new_d = daqreader.ingest_epochfiles(ef, epoch_id) + if epoch_progress is None: + new_d = daqreader.ingest_epochfiles(ef, epoch_id) + else: + new_d = daqreader.ingest_epochfiles(ef, epoch_id, ProgressFcn=epoch_progress) docs.extend(new_d if isinstance(new_d, list) else [new_d]) for mreader in daqmetadatareaders: new_d = mreader.ingest_epochfiles(ef, epoch_id) @@ -154,6 +211,26 @@ def _ingestOneEpoch( # noqa: N802 (mirrors MATLAB's local function) return docs +def _send_progress(queue: Any, k: int, fraction: float) -> None: + """Send epoch K's progress from a worker process to the caller.""" + queue.put((k, fraction)) + + +def _reader_reports_progress(daqreader: Any) -> bool: + """Does DAQREADER's ingest_epochfiles take a ProgressFcn? + + True when the method is ndi_daq_reader_mfdaq's own; a subclass that + overrides ingest_epochfiles keeps its own signature and is called without + one. + """ + if daqreader is None: + return False + from ..mfdaq import ndi_daq_reader_mfdaq + + method = getattr(type(daqreader), "ingest_epochfiles", None) + return method is ndi_daq_reader_mfdaq.ingest_epochfiles + + def _resolve_workers(override: int | None) -> int: """Worker count for this call, following NDI's parallel preferences.""" from ndi.fun.utils import parallelWorkers diff --git a/src/ndi/daq/fun/ndi_matlab_python_bridge.yaml b/src/ndi/daq/fun/ndi_matlab_python_bridge.yaml index 3617093c..674d5000 100644 --- a/src/ndi/daq/fun/ndi_matlab_python_bridge.yaml +++ b/src/ndi/daq/fun/ndi_matlab_python_bridge.yaml @@ -14,7 +14,7 @@ functions: - name: ingestEpochDocuments type: function matlab_path: "+ndi/+daq/+fun/ingestEpochDocuments.m" - matlab_last_sync_hash: "8feb9166e" + matlab_last_sync_hash: "b192ae605" python_path: "ndi/daq/fun/ingest_epoch_documents.py" input_arguments: - name: daqreader @@ -31,7 +31,7 @@ functions: type_python: "list[str]" - name: ProgressFcn type_matlab: "function handle or []" - type_python: "Callable[[int, int], Any] | None" + type_python: "Callable[[float, int], Any] | None" default: "None" - name: NumberOfWorkers type_matlab: "double ([]: use the preferences)" @@ -55,3 +55,12 @@ functions: falls back to serial with a warning. MATLAB's thread-pool refusal has no counterpart (Python workers here are always processes). A snake_case alias ingest_epoch_documents is bound to the same function. + + NDI-matlab b192ae605 (VH-Lab/NDI-matlab#1089): ProgressFcn's numDone is + fractional while an epoch is underway. A reader whose ingest_epochfiles + is ndi_daq_reader_mfdaq's own is passed a ProgressFcn and reports each + segment it writes; the total is reported only when it changes, so each + whole number is reported once. Ported. MATLAB's DataQueue sends from + the workers become a multiprocessing.Manager queue that the caller + drains while it waits on the futures; the DefiningClass check becomes + an identity check on the ingest_epochfiles function. diff --git a/src/ndi/daq/metadatareader/__init__.py b/src/ndi/daq/metadatareader/__init__.py index fc47145a..39dc2886 100644 --- a/src/ndi/daq/metadatareader/__init__.py +++ b/src/ndi/daq/metadatareader/__init__.py @@ -187,10 +187,21 @@ def readmetadata_ingested( if doc is None: return [] - # Extract metadata from document - # This would typically involve reading compressed data - # For now, return empty list - return [] + import os + + import ndicompress + + from ...database_fun import copydocfile2temp + + tname, tname_without_ext = copydocfile2temp(doc, session, "data.bin", ".nbf.tgz") + try: + parameters = ndicompress.expand_metadata(tname_without_ext) + finally: + os.remove(tname) + # MATLAB's jsondecode returns a struct array for a list of objects + if isinstance(parameters, dict): + parameters = [parameters] + return parameters if parameters is not None else [] def get_ingested_document( self, @@ -232,6 +243,12 @@ def ingest_epochfiles( """ Create a document for ingested epoch metadata. + MATLAB equivalent: ndi.daq.metadatareader/ingest_epochfiles + + The parameters returned by :meth:`readmetadata` are stored, compressed + with ``ndicompress.compress_metadata``, as the document's ``data.bin``. + The document is not added to any database. + Args: epochfiles: List of file paths epoch_id: ndi_epoch_epoch identifier @@ -249,13 +266,16 @@ def ingest_epochfiles( ) doc.set_dependency_value("daqmetadatareader_id", self.id) - # Read and store metadata - try: - self.readmetadata(epochfiles) - # In a full implementation, we would compress and store - # the parameters as binary data - except (ValueError, FileNotFoundError): - pass + import os + import tempfile + import uuid + + import ndicompress + + parameters = self.readmetadata(epochfiles) + metadatafile = os.path.join(tempfile.gettempdir(), uuid.uuid4().hex) + ndicompress.compress_metadata(parameters, metadatafile) + doc.add_file("data.bin", metadatafile + ".nbf.tgz") return doc diff --git a/src/ndi/daq/metadatareader/ndi_matlab_python_bridge.yaml b/src/ndi/daq/metadatareader/ndi_matlab_python_bridge.yaml index 72692447..05d2781a 100644 --- a/src/ndi/daq/metadatareader/ndi_matlab_python_bridge.yaml +++ b/src/ndi/daq/metadatareader/ndi_matlab_python_bridge.yaml @@ -84,7 +84,11 @@ classes: output_arguments: - name: parameters type_python: "list[dict[str, Any]]" - decision_log: "Exact match." + decision_log: > + Ported: expands the ingested document's data.bin with + ndicompress.expand_metadata (it returned an empty list). A single + JSON object comes back as a one-element list, as MATLAB's + jsondecode makes a struct array. - name: get_ingested_document input_arguments: @@ -110,7 +114,11 @@ classes: output_arguments: - name: doc type_python: "Any" - decision_log: "Exact match." + decision_log: > + Ported: readmetadata's parameters are stored with + ndicompress.compress_metadata as the document's data.bin (it made + the document without them and ignored readmetadata errors, which + now propagate as in MATLAB). - name: newdocument input_arguments: [] diff --git a/src/ndi/daq/mfdaq.py b/src/ndi/daq/mfdaq.py index 80933653..d057245f 100644 --- a/src/ndi/daq/mfdaq.py +++ b/src/ndi/daq/mfdaq.py @@ -7,7 +7,13 @@ from __future__ import annotations +import logging +import os +import tempfile +import uuid +import warnings from abc import abstractmethod +from collections.abc import Callable from dataclasses import dataclass from enum import Enum from typing import Any @@ -18,6 +24,8 @@ from .daqsystemstring import ndi_daq_daqsystemstring from .reader_base import ndi_daq_reader +logger = logging.getLogger(__name__) + class ChannelType(Enum): """Channel types for multi-function DAQ systems.""" @@ -583,6 +591,236 @@ def epochtimes2samples( return s + #: Samples per stored segment of analog/time data, as MATLAB's ingest_epochfiles. + SAMPLE_ANALOG_SEGMENT = 1_000_000 + #: Samples per stored segment of digital data, as MATLAB's ingest_epochfiles. + SAMPLE_DIGITAL_SEGMENT = 10_000_000 + + #: File-record prefix for each channel type's segments. + INGEST_FILE_PREFIX = { + "analog_in": "ai", + "analog_out": "ao", + "auxiliary_in": "ax", + "digital_in": "di", + "digital_out": "do", + "eventmarktext": "evmktx", + "time": "ti", + } + + def ingest_epochfiles( + self, + epochfiles: list[str], + epoch_id: str, + ProgressFcn: Callable[[float], Any] | None = None, # noqa: N803 (MATLAB name) + ) -> Any: + """ + Create a document holding the data read by this reader for one epoch. + + MATLAB equivalent: ndi.daq.reader.mfdaq/ingest_epochfiles + + Creates an ndi_document of type ``daqreader_mfdaq_epochdata_ingested`` + with the epoch's channel list (``channel_list.bin``) and its data, + compressed with ``ndicompress`` into segments of + ``SAMPLE_ANALOG_SEGMENT`` (analog, auxiliary, time) or + ``SAMPLE_DIGITAL_SEGMENT`` (digital) samples per channel group, plus + one file per group of event/marker/text channels. Analog data are + stored in the type given by :meth:`underlying_datatype` (e.g. Intan's + uint16). The document is not added to any database; its files are + temporary files that the database takes over when it is added. + + Args: + epochfiles: Files for this epoch + epoch_id: The epoch id + ProgressFcn: Called as ``ProgressFcn(F)`` with F, the fraction of + this epoch done (0..1), after each segment is written. Each + channel group counts in proportion to its number of channels. + Errors in it are turned into warnings. + + Returns: + ndi_document of type daqreader_mfdaq_epochdata_ingested + + Python-specific Notes: + Sample numbers passed to :meth:`readchannels_epochsamples` are + 1-based, as that method takes them. Digital-output segments are + stored as ``do_group*`` (the document's file list and the ingested + reader use ``do``); MATLAB writes them as ``di_group*``, which + collides with digital-input segments of the same group. + """ + import ndicompress + + from ..document import ndi_document + from ..file.type.mfdaq_epoch_channel import ndi_file_type_mfdaq__epoch__channel + from ..fun.utils import channelname2prefixnumber + + # the epoch table (clock and t0_t1), as the base reader writes it + base_doc = super().ingest_epochfiles(epochfiles, epoch_id) + epochtable = base_doc.document_properties["daqreader_epochdata_ingested"]["epochtable"] + + d = ndi_document( + "ingestion/daqreader_mfdaq_epochdata_ingested", + daqreader_mfdaq_epochdata_ingested={ + "parameters": { + "sample_analog_segment": self.SAMPLE_ANALOG_SEGMENT, + "sample_digital_segment": self.SAMPLE_DIGITAL_SEGMENT, + } + }, + daqreader_epochdata_ingested={"epochtable": epochtable}, + epochid={"epochid": epoch_id}, + ) + d.set_dependency_value("daqreader_id", self.id) + + # Step 1: the channel list, with each channel's sample rate, offset and scale + ch = self.getchannelsepoch(epochfiles) + names = [c.name if hasattr(c, "name") else c["name"] for c in ch] + types = [c.type if hasattr(c, "type") else c["type"] for c in ch] + numbers = [channelname2prefixnumber(n)[1] for n in names] + sample_rates = [float("nan")] * len(ch) + for rate_type in ( + "analog_in", + "analog_out", + "auxiliary_in", + "digital_in", + "digital_out", + "time", + ): + idx = [i for i, t in enumerate(types) if t == rate_type] + if idx: + sr = np.atleast_1d( + self.samplerate(epochfiles, rate_type, [numbers[i] for i in idx]) + ) + for j, i in enumerate(idx): + sample_rates[i] = float(sr[j]) + records = [] + for i, c in enumerate(ch): + _, poly, _ = self.underlying_datatype(epochfiles, types[i], numbers[i]) + poly = np.atleast_2d(poly) + time_channel = c.time_channel if hasattr(c, "time_channel") else c.get("time_channel") + records.append( + { + "name": names[i], + "type": types[i], + "time_channel": float("nan") if time_channel is None else time_channel, + "sample_rate": sample_rates[i], + "offset": float(poly[0, 0]), + "scale": float(poly[0, 1]), + } + ) + + mec = ndi_file_type_mfdaq__epoch__channel().create_properties(records) + channel_file_name = _ingest_temp_name() + b, errmsg = mec.writeToFile(channel_file_name) + if not b: + raise RuntimeError(errmsg) + d.add_file("channel_list.bin", channel_file_name) + + ci = mec.channel_information + + # Step 2: the data, one channel type and group at a time; event, marker + # and text channels are stored together + all_types, _ = self.channel_types() + emt = ("event", "marker", "text") + store_types = [t for t in all_types if t not in emt] + ["eventmarktext"] + + def entries_of(type_name: str) -> list[Any]: + if type_name == "eventmarktext": + return [c for c in ci if c.type in emt] + return [c for c in ci if c.type == type_name] + + weight_total = sum(len(entries_of(t)) for t in store_types) + weight_done = 0 + _report_ingest_progress(ProgressFcn, 0) + + for type_name in store_types: + entries = entries_of(type_name) + prefix = self.INGEST_FILE_PREFIX[type_name] + for group in sorted({c.group for c in entries}): + members = [c for c in entries if c.group == group] + weight_here = len(members) + channels_here = [c.number for c in members] + underlying_format, mypoly, datasize = self.underlying_datatype( + epochfiles, type_name, channels_here + ) + # event/marker/text channels have no scalar sample rate + if type_name == "eventmarktext": + rates = [float("nan")] + else: + rates = np.atleast_1d( + self.samplerate(epochfiles, [c.type for c in members], channels_here) + ).tolist() + unique_rates = {r for r in rates if not np.isnan(r)} + if len(unique_rates) > 1: + raise ValueError( + f"Sample rates are not all identical for {type_name} group {group}." + ) + + dataclass_here = members[0].dataclass + if dataclass_here in ("ephys", "digital", "time"): + if not unique_rates: + raise ValueError("Analog records have a NaN sample rate.") + t0t1 = self.t0_t1(epochfiles) + # epochtimes2samples is 0-based; S1 is the 1-based last sample + S1 = ( + int( + np.asarray( + self.epochtimes2samples( + type_name, channels_here, epochfiles, t0t1[0][-1] + ) + ).ravel()[0] + ) + + 1 + ) + segment = ( + self.SAMPLE_DIGITAL_SEGMENT + if dataclass_here == "digital" + else self.SAMPLE_ANALOG_SEGMENT + ) + s_starts = list(range(1, S1 + 1, segment)) + for s, s0 in enumerate(s_starts, start=1): + logger.info( + "Working on %s ingestion segment %d of %d.", + {"ephys": "ephys/analog"}.get(dataclass_here, dataclass_here), + s, + len(s_starts), + ) + s1 = min(s0 + segment - 1, S1) + data = self.readchannels_epochsamples( + type_name, channels_here, epochfiles, s0, s1 + ) + filename_here = _ingest_temp_name() + if dataclass_here == "ephys": + data = _scaled2underlying(data, mypoly) + _compress_ephys_underlying( + data, filename_here, underlying_format, datasize + ) + elif dataclass_here == "digital": + data = _scaled2underlying(data, mypoly) + ndicompress.compress_digital(data, filename_here) + else: + ndicompress.compress_time(np.asarray(data, dtype=float), filename_here) + d.add_file(f"{prefix}_group{group}_seg.nbf_{s}", filename_here + ".nbf.tgz") + _report_ingest_progress( + ProgressFcn, + (weight_done + weight_here * s / len(s_starts)) / weight_total, + ) + elif dataclass_here == "eventmarktext": + logger.info("Working on event/marker/text ingestion segment 1 of 1.") + channeltype = [c.type for c in members] + T, D = self.readevents_epochsamples_native( + channeltype, channels_here, epochfiles, -np.inf, np.inf + ) + filename_here = _ingest_temp_name() + ndicompress.compress_eventmarktext( + channeltype, channels_here, T, D, filename_here + ) + d.add_file(f"evmktx_group{group}_seg.nbf_1", filename_here + ".nbf.tgz") + else: + raise ValueError(f"Unknown channel type {type_name}.") + weight_done += weight_here + _report_ingest_progress(ProgressFcn, weight_done / max(weight_total, 1)) + + _report_ingest_progress(ProgressFcn, 1) + return d + def underlying_datatype( self, epochfiles: list[str], @@ -1290,3 +1528,78 @@ def epochtimes2samples_ingested( s[pos_inf] = s_end return s + + +def _ingest_temp_name() -> str: + """A new temporary file name (no extension), as MATLAB's ndi.file.temp_name.""" + return os.path.join(tempfile.gettempdir(), uuid.uuid4().hex) + + +def _report_ingest_progress(progress_fcn: Callable[[float], Any] | None, fraction: float) -> None: + """Call PROGRESS_FCN(FRACTION); progress reporting must never break an ingest. + + MATLAB equivalent: ndi.daq.reader.mfdaq.reportIngestProgress + """ + if progress_fcn is None: + return + try: + progress_fcn(min(max(float(fraction), 0.0), 1.0)) + except Exception as exc: # noqa: BLE001 - progress must never break an ingest + warnings.warn(f"Progress reporting failed: {exc}", stacklevel=2) + + +def _scaled2underlying(data: Any, poly: Any) -> np.ndarray: + """Scaled data to underlying values: D ./ scale + offset, per column. + + MATLAB equivalent: ndi.compress.scaled2underlying + """ + poly = np.atleast_2d(np.asarray(poly, dtype=float)) + data = np.asarray(data, dtype=float) + if data.ndim == 1: + data = data[:, np.newaxis] + return data / poly[:, 1][np.newaxis, :] + poly[:, 0][np.newaxis, :] + + +_FLOAT_TYPES = { + "float64": np.float64, + "double": np.float64, + "float32": np.float32, + "single": np.float32, +} + + +def _compress_ephys_underlying( + data: np.ndarray, filename: str, underlying_format: str, datasize: int +) -> None: + """Compress underlying-valued DATA as MATLAB's ingest does. + + MATLAB calls ndi.compress.compress_ephys(DATA, DATASIZE, FILENAME, + 'data_format', ndi.compress.datatypestring(FORMAT), 'data_is_unsigned', ...), + whose integer path rounds the samples to integers of the underlying type. + The C encoder takes the samples in that type, so they are rounded and + clipped to its range here. + """ + import ndicompress + + fmt = underlying_format.lower() + if fmt in _FLOAT_TYPES: + ndicompress.compress_ephys( + np.asarray(data, dtype=_FLOAT_TYPES[fmt]), + filename, + output_size_bits=int(datasize), + data_is_unsigned=False, + data_format="float", + ) + return + dtype = np.dtype(fmt) + if not np.issubdtype(dtype, np.integer): + raise ValueError(f"Unknown underlying data type {underlying_format}.") + info = np.iinfo(dtype) + values = np.clip(np.round(data), info.min, info.max).astype(dtype) + ndicompress.compress_ephys( + values, + filename, + output_size_bits=int(datasize), + data_is_unsigned=fmt.startswith("u"), + data_format="integer", + ) diff --git a/src/ndi/daq/ndi_matlab_python_bridge.yaml b/src/ndi/daq/ndi_matlab_python_bridge.yaml index f50dfe3e..35fe83f5 100644 --- a/src/ndi/daq/ndi_matlab_python_bridge.yaml +++ b/src/ndi/daq/ndi_matlab_python_bridge.yaml @@ -396,7 +396,7 @@ classes: - name: mfdaq_reader type: class matlab_path: "+ndi/+daq/+reader/mfdaq.m" - matlab_last_sync_hash: "770c1eb3" + matlab_last_sync_hash: "6ff0264e5" decision_log: > NDI-matlab 770c1eb3 clamps readchannels_epochsamples_ingested's window to what the epoch has: a finite s0 below the first sample @@ -414,11 +414,57 @@ classes: port: this port's ingest (ndi.daq.reader_base.ingest_epochfiles) never calls samplerate at all. samplerate here is reached only from times2samples and samples2times, which event channels do not use. + + NDI-matlab 6ff0264e5 (VH-Lab/NDI-matlab#1089) adds ingest_epochfiles' + ProgressFcn. With it, this port gains ingest_epochfiles itself: until + now ingestion here made the epoch documents without any data, so an + ingested epoch could not be read. It now writes channel_list.bin and + the compressed segments as MATLAB does (see the method's entry). python_path: "ndi/daq/mfdaq.py" python_class: "ndi_daq_reader_mfdaq" inherits: "ndi.daq.reader" methods: + - name: ingest_epochfiles + input_arguments: + - name: epochfiles + type_matlab: "cell array of char" + type_python: "list[str]" + - name: epoch_id + type_matlab: "char" + type_python: "str" + - name: ProgressFcn + type_matlab: "function handle or []" + type_python: "Callable[[float], Any] | None" + default: "None" + output_arguments: + - name: d + type_python: "ndi_document" + decision_log: > + Ported from NDI-matlab 6ff0264e5. Same document, file records, + segment sizes (1e6 analog/time, 1e7 digital), channel list + (sample rate, offset and scale per channel, written in MATLAB's + saveStructArray format), per-group sample-rate check, progress + weights and priority of log messages (logger.info). Analog data are + converted to the reader's underlying_datatype and compressed with + ndicompress.compress_ephys's data_format/data_is_unsigned/ + output_size_bits, which mirror MATLAB's compress_ephys (they need + ndi-compress with Waltham-Data-Science/NDI-compress#14). + Differences: (1) readchannels_epochsamples takes 1-based samples, + and epochtimes2samples here returns 0-based ones, so S1 is the + latter plus 1. (2) Digital-output segments are written as + do_group*, which is what the document's file list and + readchannels_epochsamples_ingested use; MATLAB writes di_group*, + colliding with digital inputs of the same group. (3) The time + channel group computes its own sample range; MATLAB reuses the + previous group's S0/S1/channels (the same values whenever the time + channel shares the rate of the channel type before it). + (4) MATLAB maps an underlying 'float64' through + ndi.compress.datatypestring, which accepts only MATLAB class names + ('double'), so a reader that keeps analog data in doubles cannot + ingest there; here 'float64'/'float32' (and 'double'/'single') are + stored as floats. + - name: epochclock input_arguments: - name: epochfiles diff --git a/src/ndi/daq/reader/mfdaq/intan.py b/src/ndi/daq/reader/mfdaq/intan.py index d21fab76..c7e40fda 100644 --- a/src/ndi/daq/reader/mfdaq/intan.py +++ b/src/ndi/daq/reader/mfdaq/intan.py @@ -184,6 +184,7 @@ def readchannels_epochsamples( n_samples = s1 - s0 + 1 result = np.zeros((n_samples, len(channel))) + n_read = n_samples # Cache NDR reads by channel type to avoid re-reading the file ndr_cache: dict[str, np.ndarray] = {} @@ -196,6 +197,7 @@ def readchannels_epochsamples( ) data = ndr_cache["time"] n = min(n_samples, data.shape[0]) + n_read = min(n_read, n) result[:n, col] = data[:n, 0] elif ct == "analog_in": @@ -211,6 +213,7 @@ def readchannels_epochsamples( ) data = ndr_cache["amp"] n = min(n_samples, data.shape[0]) + n_read = min(n_read, n) result[:n, col] = data[:n, ch_num - 1] else: # Board ADC channel @@ -233,62 +236,114 @@ def readchannels_epochsamples( adc_idx = ch_num - n_amp - 1 data = ndr_cache["adc"] n = min(n_samples, data.shape[0]) + n_read = min(n_read, n) result[:n, col] = data[:n, adc_idx] elif ct == "auxiliary_in": + # Aux inputs are read at their own sample rate, as MATLAB + # does: S0 and S1 count aux samples (see samplerate). if "aux" not in ndr_cache: n_aux = len(header.get("aux_input_channels", [])) if n_aux > 0: - # Read all aux data to interpolate correctly + aux_sr = header["frequency_parameters"].get("aux_input_sample_rate", sr / 4) ndr_cache["aux"] = read_Intan_RHD2000_datafile( filepath, channeltype="aux", channel=list(range(1, n_aux + 1)), - t0=0.0, - t1=float("inf"), + t0=(s0 - 1) / aux_sr, + t1=(s1 - 1) / aux_sr, ) if "aux" in ndr_cache: - ch_idx = ch_num - 1 - aux_data = ndr_cache["aux"] - # Aux is sampled at 1/4 rate; interpolate to main rate - i0 = s0 - 1 - i1 = s1 # exclusive end - aux_i0 = i0 // 4 - aux_i1 = i1 // 4 - if aux_i1 > aux_i0 and aux_i1 <= aux_data.shape[0]: - result[:, col] = np.interp( - np.arange(i0, i1), - np.arange(aux_i0 * 4, aux_i1 * 4, 4), - aux_data[aux_i0:aux_i1, ch_idx], - ) + data = ndr_cache["aux"] + n = min(n_samples, data.shape[0]) + n_read = min(n_read, n) + result[:n, col] = data[:n, ch_num - 1] elif ct == "digital_in": if "din" not in ndr_cache: + # NDR numbers the enabled digital channels 1..N and + # maps each to its bit (native_order) ndr_cache["din"] = read_Intan_RHD2000_datafile( filepath, channeltype="din", - channel=list(range(1, 17)), + channel=list(range(1, len(header.get("board_dig_in_channels", [])) + 1)), t0=t0_sec, t1=t1_sec, ) data = ndr_cache["din"] n = min(n_samples, data.shape[0]) + n_read = min(n_read, n) result[:n, col] = data[:n, ch_num - 1] elif ct == "digital_out": if "dout" not in ndr_cache: + # NDR numbers the enabled digital channels 1..N and + # maps each to its bit (native_order) ndr_cache["dout"] = read_Intan_RHD2000_datafile( filepath, channeltype="dout", - channel=list(range(1, 17)), + channel=list(range(1, len(header.get("board_dig_out_channels", [])) + 1)), t0=t0_sec, t1=t1_sec, ) data = ndr_cache["dout"] n = min(n_samples, data.shape[0]) + n_read = min(n_read, n) result[:n, col] = data[:n, ch_num - 1] - return result + # As MATLAB, return only the samples the file has: a request that runs + # past the end of the recording is not padded with zeros. + return result[:n_read] + + def underlying_datatype( + self, + epochfiles: list[str], + channeltype: str, + channel: int | list[int], + ) -> tuple[str, np.ndarray, int]: + """ + Get the underlying data type for channels. + + MATLAB equivalent: ndi.daq.reader.mfdaq.intan/underlying_datatype + + Amplifier and auxiliary channels are stored by Intan as uint16 + samples, so ingestion stores them that way: amplifier samples convert + as (raw - 32768) * 0.195 (microvolts), auxiliary samples as + raw * 3.74e-5 (volts). + + Args: + epochfiles: Files for this epoch + channeltype: Channel type + channel: Channel number(s) + + Returns: + Tuple of (datatype, polynomial, datasize), with one polynomial row + per channel. + + Python-specific Notes: + The Python reader also returns the board ADC inputs as analog_in + channels numbered after the amplifier channels (MATLAB's reader + does not list them). They are uint16 too, converting as + (raw - 32768) * 0.0003125 (volts), so their rows differ. + """ + if isinstance(channel, (int, np.integer)): + channel = [channel] + channel = [int(c) for c in channel] + channeltype = standardize_channel_type(channeltype) + + if channeltype in ("analog_in", "analog_out"): + header = self._get_header(epochfiles) or {} + n_amp = len(header.get("amplifier_channels", [])) + poly = np.array( + [ + [32768.0, 0.195] if (c <= n_amp or not n_amp) else [32768.0, 0.0003125] + for c in channel + ] + ).reshape(len(channel), 2) + return "uint16", poly, 16 + if channeltype == "auxiliary_in": + return "uint16", np.tile([0.0, 3.74e-05], (len(channel), 1)), 16 + return super().underlying_datatype(epochfiles, channeltype, channel) def t0_t1(self, epochfiles: list[str]) -> list[tuple[float, float]]: header = self._get_header(epochfiles) diff --git a/src/ndi/daq/reader/mfdaq/ndi_matlab_python_bridge.yaml b/src/ndi/daq/reader/mfdaq/ndi_matlab_python_bridge.yaml index 3b6a6656..ce8ebed9 100644 --- a/src/ndi/daq/reader/mfdaq/ndi_matlab_python_bridge.yaml +++ b/src/ndi/daq/reader/mfdaq/ndi_matlab_python_bridge.yaml @@ -61,6 +61,33 @@ classes: type_python: "list[ChannelInfo]" decision_log: "Exact match. Delegates to ndi_daq_reader_SpikeInterfaceReader." + - name: underlying_datatype + input_arguments: + - name: epochfiles + type_matlab: "cell array of char" + type_python: "list[str]" + - name: channeltype + type_matlab: "char" + type_python: "str" + - name: channel + type_matlab: "integer | array" + type_python: "int | list[int]" + output_arguments: + - name: datatype + type_python: "str" + - name: p + type_python: "np.ndarray" + - name: datasize + type_python: "int" + decision_log: > + Ported: uint16/16 bits with [32768 0.195] for analog channels and + [0 3.74e-05] for auxiliary inputs; time, digital and event types + as the base class. One polynomial row per channel. The Python + reader also lists the board ADC inputs as analog_in channels after + the amplifier channels (MATLAB's does not); their rows are + [32768 0.0003125], matching how readchannels_epochsamples scales + them. + - name: readchannels_epochsamples input_arguments: - name: channeltype @@ -83,7 +110,14 @@ classes: type_python: "np.ndarray" decision_log: > Semantic Parity: channel and s0/s1 are 1-indexed. - Delegates to ndi_daq_reader_SpikeInterfaceReader. + Reads through ndr.format.intan. As MATLAB: auxiliary inputs are + read at their own sample rate (s0/s1 count aux samples, as + samplerate reports them; they were interpolated to the amplifier + rate, which disagreed with samplerate); digital channels are the + enabled channels 1..N, each mapped by NDR to its bit (16 were + requested, which NDR rejects); and a request past the end of the + recording returns only the samples the file has rather than zero + padding. - name: samplerate input_arguments: diff --git a/src/ndi/file/type/mfdaq_epoch_channel.py b/src/ndi/file/type/mfdaq_epoch_channel.py index 0620366d..33da9df7 100644 --- a/src/ndi/file/type/mfdaq_epoch_channel.py +++ b/src/ndi/file/type/mfdaq_epoch_channel.py @@ -16,6 +16,7 @@ from __future__ import annotations import json +import math from dataclasses import asdict, dataclass, field from pathlib import Path from typing import Any @@ -54,6 +55,12 @@ } +def _int_or_none(value: Any) -> int | None: + """An integer field read from a file; NaN (a missing value) becomes None.""" + value = float(value) + return None if math.isnan(value) else int(value) + + @dataclass class ChannelInfo: """Information about a single channel in an MFDAQ recording.""" @@ -262,7 +269,7 @@ def readFromFile(self, filename: str) -> ndi_file_type_mfdaq__epoch__channel: ChannelInfo( name=str(rec.get("name", "")), type=str(rec.get("type", "")), - time_channel=int(rec.get("time_channel", 1)), + time_channel=_int_or_none(rec.get("time_channel", 1)), sample_rate=float(rec.get("sample_rate", 0.0)), offset=float(rec.get("offset", 0.0)), scale=float(rec.get("scale", 1.0)), @@ -278,10 +285,14 @@ def readFromFile(self, filename: str) -> ndi_file_type_mfdaq__epoch__channel: def writeToFile(self, filename: str) -> tuple[bool, str]: """ - Write channel information to a JSON file. + Write channel information to a file. MATLAB equivalent: ndi.file.type.mfdaq_epoch_channel/writeToFile + Writes the tab-delimited format of ``vlt.file.saveStructArray``, as + MATLAB does, so NDI-matlab can read the channel list of an epoch + ingested by NDI-python. A missing ``time_channel`` is written as NaN. + Args: filename: Path to write @@ -289,10 +300,16 @@ def writeToFile(self, filename: str) -> tuple[bool, str]: Tuple of (success, error_message) """ try: - data = {"channel_information": [ch.to_dict() for ch in self.channel_information]} + from vlt.file import saveStructArray + + records = [] + for ch in self.channel_information: + rec = ch.to_dict() + if rec["time_channel"] is None: + rec["time_channel"] = float("nan") + records.append(rec) Path(filename).parent.mkdir(parents=True, exist_ok=True) - with open(filename, "w") as f: - json.dump(data, f, indent=2) + saveStructArray(filename, records) return True, "" except Exception as e: return False, str(e) diff --git a/src/ndi/file/type/ndi_matlab_python_bridge.yaml b/src/ndi/file/type/ndi_matlab_python_bridge.yaml index 3ecfaf6b..0899d9df 100644 --- a/src/ndi/file/type/ndi_matlab_python_bridge.yaml +++ b/src/ndi/file/type/ndi_matlab_python_bridge.yaml @@ -119,6 +119,11 @@ classes: (snake_case) and returned None. Renamed to 'writeToFile' and updated to return 2-tuple (success_bool, error_message) matching MATLAB. Synchronized 2026-03-13. + Now writes MATLAB's format (vlt.file.saveStructArray, tab + delimited) rather than JSON, so NDI-matlab can read the channel + list of an epoch ingested by NDI-python; a missing time_channel is + written as NaN. readFromFile still reads the older JSON files and + reads a NaN time_channel back as None. # --- Channel Queries --- - name: channels_of_type diff --git a/tests/test_ingest_data_writer.py b/tests/test_ingest_data_writer.py new file mode 100644 index 00000000..eac40595 --- /dev/null +++ b/tests/test_ingest_data_writer.py @@ -0,0 +1,273 @@ +"""Ingestion writes the epoch's data, as NDI-matlab does. + +MATLAB counterpart: ndi.daq.reader.mfdaq/ingest_epochfiles (with its +ProgressFcn, VH-Lab/NDI-matlab#1089), ndi.daq.metadatareader/ingest_epochfiles, +ndi.file.type.mfdaq_epoch_channel/writeToFile and +ndi.daq.reader.mfdaq.intan/underlying_datatype. + +Ingestion used to create the epoch documents without any data, so nothing +could be read from an ingested epoch. These tests ingest a copy of the example +Intan session and read every channel type back from the ingested files, +comparing it with what the reader returns from the raw file. +""" + +from __future__ import annotations + +import json + +import numpy as np +import pytest + +from ndi.daq.fun import ingestEpochDocuments +from ndi.daq.mfdaq import ndi_daq_reader_mfdaq +from ndi.file.type.mfdaq_epoch_channel import ChannelInfo, ndi_file_type_mfdaq__epoch__channel +from ndi.gui.component.ProgressBarWindow import ndi_gui_component_ProgressBarWindow as PBW +from tests.test_ingest_progress_parallel import _have_intan, _intan_session, _prefs + +ndicompress = pytest.importorskip("ndicompress") + +#: channel type, channel numbers, largest difference allowed (time channels +#: are stored by a fit, so they round-trip to ~1e-12 s, not exactly) +_CHANNELS = [ + ("analog_in", list(range(1, 33)), 0.0), + ("auxiliary_in", [1, 2, 3], 0.0), + ("digital_in", [1], 0.0), + ("time", [1], 1e-9), +] + + +@pytest.fixture(autouse=True) +def _silent(): + previous = PBW.silentModeDefault() + PBW.silentModeDefault(True) + yield + PBW.silentModeDefault(previous) + + +def _daq(session): + d = session.daqsystem_load(name="intan1") + return d[0] if isinstance(d, list) else d + + +def _epochs(session): + """(raw epoch files, 'epochid://' files) for each epoch.""" + return [ + (e["underlying_epochs"]["underlying"], ["epochid://" + e["epoch_id"]]) + for e in _daq(session).epochtable() + ] + + +def _ingested(session, monkeypatch): + _prefs(monkeypatch, OpenPoolIfNone=False) + ok, msg = session.ingest() + assert ok, msg + return session + + +def _assert_reads_match(session, windows=None): + reader = _daq(session).daqreader + for files, ingested in _epochs(session): + for channeltype, numbers, tolerance in _CHANNELS: + if windows is None: + last = int( + reader.epochtimes2samples( + channeltype, numbers, files, reader.t0_t1(files)[0][1] + ) + ) + spans = [(0, last)] + else: + spans = windows + for s0, s1 in spans: + raw = reader.readchannels_epochsamples(channeltype, numbers, files, s0 + 1, s1 + 1) + stored = reader.readchannels_epochsamples_ingested( + channeltype, numbers, ingested, s0, s1, session + ) + assert stored.shape == raw.shape, (channeltype, s0, s1) + assert np.max(np.abs(stored - raw)) <= tolerance, (channeltype, s0, s1) + + +@pytest.mark.skipif(not _have_intan(), reason="Intan example data (.rhd) not found") +class TestIntanIngestion: + def test_every_channel_type_reads_back(self, tmp_path, monkeypatch): + session = _ingested(_intan_session(tmp_path / "exp1"), monkeypatch) + _assert_reads_match(session) + + def test_document_files_and_channel_list(self, tmp_path, monkeypatch): + session = _ingested(_intan_session(tmp_path / "exp1"), monkeypatch) + reader = _daq(session).daqreader + files, ingested = _epochs(session)[0] + doc = reader.getingesteddocument(ingested, session) + names = sorted(f["name"] for f in doc.document_properties["files"]["file_info"]) + assert names == [ + "ai_group1_seg.nbf_1", + "ax_group1_seg.nbf_1", + "channel_list.bin", + "di_group1_seg.nbf_1", + "ti_group1_seg.nbf_1", + ] + params = doc.document_properties["daqreader_mfdaq_epochdata_ingested"]["parameters"] + assert params == {"sample_analog_segment": 1_000_000, "sample_digital_segment": 10_000_000} + + # channel_list.bin is MATLAB's tab-delimited saveStructArray format + fobj = session.database_openbinarydoc(doc, "channel_list.bin") + text = fobj.read() + fobj.close() + text = text.decode() if isinstance(text, bytes) else text + assert text.splitlines()[0].split("\t") == [ + "name", + "type", + "time_channel", + "sample_rate", + "offset", + "scale", + "number", + "group", + "dataclass", + ] + channels = reader.getchannelsepoch_ingested(ingested, session) + assert len(channels) == len(reader.getchannelsepoch(files)) + ai1 = next(c for c in channels if c.type == "analog_in" and c.number == 1) + assert (ai1.sample_rate, ai1.offset, ai1.scale) == (20000.0, 32768.0, 0.195) + assert (ai1.group, ai1.dataclass) == (1, "ephys") + aux1 = next(c for c in channels if c.type == "auxiliary_in" and c.number == 1) + assert (aux1.sample_rate, aux1.offset, aux1.scale) == (5000.0, 0.0, 3.74e-05) + + def test_analog_is_stored_as_uint16(self, tmp_path, monkeypatch): + session = _ingested(_intan_session(tmp_path / "exp1"), monkeypatch) + reader = _daq(session).daqreader + _, ingested = _epochs(session)[0] + doc = reader.getingesteddocument(ingested, session) + import os + import tarfile + import tempfile + + from ndicompress.header import read_ndi_header + + from ndi.database_fun import copydocfile2temp + + for name in ("ai_group1_seg.nbf_1", "ax_group1_seg.nbf_1"): + tname, _ = copydocfile2temp(doc, session, name, ".nbf.tgz") + with tempfile.TemporaryDirectory() as d, tarfile.open(tname, "r:gz") as tar: + member = next(m for m in tar.getmembers() if m.name.endswith(".nbh")) + tar.extract(member, path=d, filter="data") + header = read_ndi_header(os.path.join(d, member.name)) + os.remove(tname) + assert header["original_isunsigned"] == 1, name + assert header["original_bits_per_sample"] == 16, name + + def test_epochs_span_several_segments(self, tmp_path, monkeypatch): + # small segments, so each epoch is split across files and reads cross + # segment boundaries + monkeypatch.setattr(ndi_daq_reader_mfdaq, "SAMPLE_ANALOG_SEGMENT", 7001) + monkeypatch.setattr(ndi_daq_reader_mfdaq, "SAMPLE_DIGITAL_SEGMENT", 13003) + session = _ingested(_intan_session(tmp_path / "exp1"), monkeypatch) + reader = _daq(session).daqreader + _, ingested = _epochs(session)[0] + doc = reader.getingesteddocument(ingested, session) + names = {f["name"] for f in doc.document_properties["files"]["file_info"]} + assert "ai_group1_seg.nbf_9" in names # 60000 samples / 7001 + assert "di_group1_seg.nbf_5" in names # 60000 / 13003 + _assert_reads_match(session) + _assert_reads_match(session, windows=[(6990, 7015), (0, 5), (13000, 14010), (14990, 14999)]) + + def test_reader_reports_within_epoch_progress(self, tmp_path, monkeypatch): + monkeypatch.setattr(ndi_daq_reader_mfdaq, "SAMPLE_ANALOG_SEGMENT", 20000) + session = _intan_session(tmp_path / "exp1") + reader = _daq(session).daqreader + files, _ = _epochs(session)[0] + fractions = [] + reader.ingest_epochfiles(files, "epoch_x", ProgressFcn=fractions.append) + assert fractions[0] == 0 + assert fractions[-1] == 1 + assert all(b >= a for a, b in zip(fractions, fractions[1:])) + assert any(0 < f < 1 for f in fractions) + + def test_ingest_epoch_documents_reports_fractions(self, tmp_path, monkeypatch): + _prefs(monkeypatch, OpenPoolIfNone=False) + monkeypatch.setattr(ndi_daq_reader_mfdaq, "SAMPLE_ANALOG_SEGMENT", 20000) + session = _intan_session(tmp_path / "exp1") + dev = _daq(session) + epochs = _epochs(session) + calls = [] + docs, used_parallel = ingestEpochDocuments( + dev.daqreader, + [], + [files for files, _ in epochs], + [f"epoch_{k}" for k in range(len(epochs))], + ProgressFcn=lambda done, total: calls.append((done, total)), + ) + assert not used_parallel + done = [c[0] for c in calls] + assert all(b > a for a, b in zip(done, done[1:])), "reported only when it changes" + assert done[-1] == len(epochs) + assert any(d != int(d) for d in done), "fractional progress within an epoch" + for whole in range(1, len(epochs) + 1): + assert done.count(whole) == 1 + assert len(docs) == len(epochs) + + +class TestIntanUnderlyingDatatype: + def test_matches_matlab(self): + from ndi.daq.reader.mfdaq.intan import ndi_daq_reader_mfdaq_intan + + r = ndi_daq_reader_mfdaq_intan() + dtype, poly, size = r.underlying_datatype([], "analog_in", [1, 2]) + assert (dtype, size) == ("uint16", 16) + np.testing.assert_array_equal(poly, [[32768, 0.195], [32768, 0.195]]) + dtype, poly, size = r.underlying_datatype([], "auxiliary_in", 1) + assert (dtype, size) == ("uint16", 16) + np.testing.assert_array_equal(poly, [[0, 3.74e-05]]) + assert r.underlying_datatype([], "time", 1)[0] == "float64" + assert r.underlying_datatype([], "digital_in", 1)[2] == 8 + + +class TestChannelListFile: + def test_round_trip_in_matlab_format(self, tmp_path): + mec = ndi_file_type_mfdaq__epoch__channel().create_properties( + [ + {"name": "ai1", "type": "analog_in", "time_channel": 1, "sample_rate": 20000.0}, + {"name": "t1", "type": "time", "time_channel": None, "sample_rate": 20000.0}, + {"name": "mk1", "type": "marker", "time_channel": 1, "sample_rate": np.nan}, + ] + ) + fname = str(tmp_path / "channel_list.bin") + ok, msg = mec.writeToFile(fname) + assert ok, msg + with open(fname) as f: + assert f.readline().startswith("name\ttype\ttime_channel") + back = ndi_file_type_mfdaq__epoch__channel().readFromFile(fname).channel_information + assert [(c.name, c.type, c.number, c.group, c.dataclass) for c in back] == [ + (c.name, c.type, c.number, c.group, c.dataclass) for c in mec.channel_information + ] + t1 = next(c for c in back if c.name == "t1") + assert t1.time_channel is None + mk1 = next(c for c in back if c.name == "mk1") + assert np.isnan(mk1.sample_rate) + + def test_reads_older_json_files(self, tmp_path): + fname = tmp_path / "channel_list.bin" + ch = ChannelInfo(name="ai1", type="analog_in", sample_rate=1000.0, number=1, group=1) + fname.write_text(json.dumps({"channel_information": [ch.to_dict()]})) + back = ndi_file_type_mfdaq__epoch__channel().readFromFile(str(fname)) + assert back.channel_information == [ch] + + +class TestMetadataIngestion: + def test_parameters_are_stored_and_read_back(self, tmp_path): + from ndi.daq.metadatareader import ndi_daq_metadatareader + from ndi.session.dir import ndi_session_dir + + (tmp_path / "exp").mkdir() + tsv = tmp_path / "exp" / "stims.tsv" + tsv.write_text("stimid\tangle\n1\t0\n2\t90\n") + session = ndi_session_dir("exp", tmp_path / "exp") + reader = ndi_daq_metadatareader(r"stims\.tsv", session=session) + expected = reader.readmetadata([str(tsv)]) + assert len(expected) == 2 + + doc = reader.ingest_epochfiles([str(tsv)], "epoch_1") + names = [f["name"] for f in doc.document_properties["files"]["file_info"]] + assert names == ["data.bin"] + session.database_add(reader.newdocument()) + session.database_add(doc) + assert reader.readmetadata_ingested(["epochid://epoch_1"], session) == expected From 30435c0bccbfb83c7b60b5126465c55805588cbb Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 31 Oct 2026 13:59:09 +0000 Subject: [PATCH 2/2] Port NDI-matlab #1090 (ingest by segment) and #1091 (ingest compare) ndi_daq_reader_mfdaq.ingest_epochfiles is split into ingest_epochfiles_plan, ingest_segment and the static ingest_attach_segments, as in MATLAB, and ingest_epochfiles runs them in sequence. ingestEpochDocuments now plans each epoch in the caller and hands single segments (and each epoch's metadata readers) to the worker processes, so the pool stays busy however unequal the epochs are; progress is counted in the caller from the plan weights, which retires the Manager queue. New ndi.test.ingest.compare / mfdaq_compare compare a raw session with an ingested copy: epochs, channel lists, every sample of every continuous type (RelativeTolerance / TimeTolerance) and every event, with a report of every check. tests/test_ingest_compare.py mirrors ingestCompareTest, plus a parallel many-segment ingest that must compare equal, a corrupted sample that must be reported, and MATLAB's plan/segment/attach equivalence test. Bridge: mfdaq_reader -> 5c5784c6a with the three new methods, ingestEpochDocuments -> 569fdab67; fun/epoch (9821a12e6) and setup (fd7ddbf06) bumped with notes, as nothing there needs porting. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_011v7R5rnKYwncJkrNZRXhtw --- CHANGELOG.md | 10 + src/ndi/daq/fun/ingest_epoch_documents.py | 216 +++++++--- src/ndi/daq/fun/ndi_matlab_python_bridge.yaml | 14 +- src/ndi/daq/mfdaq.py | 278 +++++++++---- src/ndi/daq/ndi_matlab_python_bridge.yaml | 74 +++- src/ndi/fun/ndi_matlab_python_bridge.yaml | 16 +- src/ndi/ndi_matlab_python_bridge.yaml | 9 +- src/ndi/test/__init__.py | 5 + src/ndi/test/ingest/__init__.py | 9 + src/ndi/test/ingest/compare.py | 146 +++++++ src/ndi/test/ingest/mfdaq_compare.py | 390 ++++++++++++++++++ tests/test_ingest_compare.py | 204 +++++++++ 12 files changed, 1225 insertions(+), 146 deletions(-) create mode 100644 src/ndi/test/__init__.py create mode 100644 src/ndi/test/ingest/__init__.py create mode 100644 src/ndi/test/ingest/compare.py create mode 100644 src/ndi/test/ingest/mfdaq_compare.py create mode 100644 tests/test_ingest_compare.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 7247647a..04b078bc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,6 +26,16 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/). longer pads reads that run past the end of the file. Needs ndi-compress with Waltham-Data-Science/NDI-compress-python#4. +- **Parallel ingestion by segment, and an ingest comparison tool** (mirroring + VH-Lab/NDI-matlab#1090 and #1091). `ingest_epochfiles` is split into + `ingest_epochfiles_plan`, `ingest_segment` and the static + `ingest_attach_segments`, and `ingestEpochDocuments` now hands single + segments to the worker processes, so every worker stays busy even when one + epoch is much longer than the others. The new `ndi.test.ingest.compare` + and `mfdaq_compare` compare a raw session with an ingested copy of it -- + every epoch, channel list, sample and event -- and return a report of + every check, with no stopping or prompting on a mismatch. + - **Ingestion progress bars and parallel per-epoch work** (mirroring VH-Lab/NDI-matlab#1079). `ndi_session.ingest()` now shows progress through `ndi_gui_component_ProgressBarWindow` -- one bar for the session and one per diff --git a/src/ndi/daq/fun/ingest_epoch_documents.py b/src/ndi/daq/fun/ingest_epoch_documents.py index 12678191..6f58a7b5 100644 --- a/src/ndi/daq/fun/ingest_epoch_documents.py +++ b/src/ndi/daq/fun/ingest_epoch_documents.py @@ -25,7 +25,6 @@ import warnings from collections.abc import Callable from concurrent.futures import FIRST_COMPLETED, ProcessPoolExecutor, wait -from queue import Empty from typing import Any from pydantic import ConfigDict, validate_call @@ -86,11 +85,14 @@ def ingestEpochDocuments( # noqa: N802 (MATLAB mirror) workers -- and live only for this call. Readers are sent with their ``_session`` cleared (see the module docstring). A reader that cannot be pickled falls back to serial with a warning rather than failing. - Completed epochs are reported as they finish, so ProgressFcn may see - epochs out of order, but the returned list is always in epoch order. - Within-epoch progress travels from the workers on a - ``multiprocessing.Manager`` queue (MATLAB uses a DataQueue) that the - caller drains while it waits. + On workers, a reader whose ingest_epochfiles is ndi_daq_reader_mfdaq's + is ingested one segment per task (MATLAB's ``parfeval`` scheduler): + every epoch is planned in the calling process, the segments of all + epochs (and each epoch's metadata readers) are handed to the workers + as they free up, and the files are attached in plan order afterwards. + Progress is reported as each task finishes, weighted by segment size, + so ProgressFcn may see epochs out of order, but the returned list is + always in epoch order. The pool is used even for a single epoch. """ num_epochs = len(epochFiles) if len(epochIds) != num_epochs: @@ -114,12 +116,13 @@ def update_progress(k: int, fraction: float) -> None: except Exception as exc: # noqa: BLE001 - progress must never break an ingest warnings.warn(f"Progress reporting failed: {exc}", stacklevel=3) - # only readers using ndi_daq_reader_mfdaq's ingest_epochfiles take a - # ProgressFcn; any other reader reports per epoch only - within_epoch = ProgressFcn is not None and _reader_reports_progress(daqreader) + # readers using ndi_daq_reader_mfdaq's ingest_epochfiles take a ProgressFcn + # and can be ingested segment by segment; any other reader works (and + # reports) per epoch + plans_segments = _reader_plans_segments(daqreader) number_of_workers = 0 - if num_epochs > 1: + if num_epochs > 0: number_of_workers = _resolve_workers(NumberOfWorkers) payload = None @@ -130,46 +133,21 @@ def update_progress(k: int, fraction: float) -> None: if number_of_workers > 0: reader, mreaders = payload - ctx = multiprocessing.get_context("spawn") - manager = ctx.Manager() if within_epoch else None - queue = manager.Queue() if manager is not None else None - - def drain() -> None: - while queue is not None: - try: - k, fraction = queue.get_nowait() - except Empty: - return - update_progress(k, fraction) - - executor = ProcessPoolExecutor( - max_workers=min(number_of_workers, num_epochs), mp_context=ctx + docs_per_epoch = _ingest_on_pool( + daqreader, + reader, + mreaders, + epochFiles, + epochIds, + number_of_workers, + plans_segments, + update_progress, ) - try: - futures = { - executor.submit(_ingestOneEpoch, reader, mreaders, ef, eid, queue, k): k - for k, (ef, eid) in enumerate(zip(epochFiles, epochIds)) - } - pending = set(futures) - while pending: - finished, pending = wait(pending, timeout=0.25, return_when=FIRST_COMPLETED) - drain() - for future in finished: - docs_per_epoch[futures[future]] = future.result() - update_progress(futures[future], 1) - drain() - except BaseException: - executor.shutdown(wait=True, cancel_futures=True) - raise - finally: - if manager is not None: - manager.shutdown() - executor.shutdown(wait=True) return docs_per_epoch, True for k, (ef, eid) in enumerate(zip(epochFiles, epochIds)): epoch_progress = None - if within_epoch: + if plans_segments and ProgressFcn is not None: epoch_progress = functools.partial(update_progress, k) docs_per_epoch[k] = _ingestOneEpoch( daqreader, daqmetadatareaders, ef, eid, epoch_progress=epoch_progress @@ -182,22 +160,131 @@ def drain() -> None: ingest_epoch_documents = ingestEpochDocuments +def _ingest_on_pool( + daqreader: Any, + reader: Any, + mreaders: list[Any], + epoch_files: list[list[str]], + epoch_ids: list[str], + number_of_workers: int, + plans_segments: bool, + update_progress: Callable[[int, float], None], +) -> list[list[Any]]: + """Ingest the epochs on worker processes, one segment per task. + + MATLAB counterpart: ``ingestOnPool``. Every epoch is planned here + (ndi_daq_reader_mfdaq.ingest_epochfiles_plan reads only file headers and + writes the channel list); each planned segment, and each epoch's + metadata readers, is a task. Readers that do not plan segments run one + task per epoch. Tasks are handed out as workers free up, so the pool stays + busy however unequal the epochs are, and the segment files are attached + in plan order afterwards, so the documents match a serial ingest. + """ + from ..mfdaq import ndi_daq_reader_mfdaq + + num_epochs = len(epoch_files) + epoch_docs: list[Any] = [None] * num_epochs + plans: list[list[dict[str, Any]]] = [[] for _ in range(num_epochs)] + segment_files: list[list[Any]] = [[] for _ in range(num_epochs)] + meta_docs: list[list[Any]] = [[] for _ in range(num_epochs)] + docs_per_epoch: list[list[Any]] = [[] for _ in range(num_epochs)] + + # (epoch, kind, index into the plan, weight) + tasks: list[tuple[int, str, int, float]] = [] + for k in range(num_epochs): + if plans_segments: + epoch_docs[k], plans[k] = daqreader.ingest_epochfiles_plan(epoch_files[k], epoch_ids[k]) + segment_files[k] = [None] * len(plans[k]) + tasks.extend((k, "segment", t, e["weight"]) for t, e in enumerate(plans[k])) + if mreaders: + tasks.append((k, "metadata", 0, 0.0)) + else: + tasks.append((k, "epoch", 0, 1.0)) + remaining = [0] * num_epochs + for k, *_ in tasks: + remaining[k] += 1 + weight_done = [0.0] * num_epochs + + ctx = multiprocessing.get_context("spawn") + executor = ProcessPoolExecutor(max_workers=number_of_workers, mp_context=ctx) + try: + futures = {} + for index, (k, kind, t, _) in enumerate(tasks): + entry = plans[k][t] if kind == "segment" else None + future = executor.submit( + _run_task, + reader, + mreaders if kind != "segment" else [], + epoch_files[k], + epoch_ids[k], + kind, + entry, + ) + futures[future] = index + pending = set(futures) + while pending: + finished, pending = wait(pending, return_when=FIRST_COMPLETED) + for future in finished: + k, kind, t, weight = tasks[futures[future]] + out = future.result() + if kind == "segment": + segment_files[k][t] = out + elif kind == "metadata": + meta_docs[k] = out + else: + docs_per_epoch[k] = out + remaining[k] -= 1 + weight_done[k] += weight + update_progress(k, 1 if remaining[k] == 0 else min(weight_done[k], 0.999)) + except BaseException: + executor.shutdown(wait=True, cancel_futures=True) + raise + executor.shutdown(wait=True) + for k in range(num_epochs): + update_progress(k, 1) # e.g., an epoch with nothing to write + + if plans_segments: + for k in range(num_epochs): + d = ndi_daq_reader_mfdaq.ingest_attach_segments( + epoch_docs[k], plans[k], segment_files[k] + ) + docs_per_epoch[k] = [d, *meta_docs[k]] + return docs_per_epoch + + +def _run_task( + daqreader: Any, + daqmetadatareaders: list[Any], + ef: list[str], + epoch_id: str, + kind: str, + plan_entry: dict[str, Any] | None, +) -> Any: + """One worker task: a segment, an epoch's metadata, or a whole epoch.""" + if kind == "segment": + return daqreader.ingest_segment(ef, plan_entry) + if kind == "metadata": + return _metadata_docs(daqmetadatareaders, ef, epoch_id) + return _ingestOneEpoch(daqreader, daqmetadatareaders, ef, epoch_id) + + +def _metadata_docs(daqmetadatareaders: list[Any], ef: list[str], epoch_id: str) -> list[Any]: + """The metadata readers' documents for one epoch.""" + docs: list[Any] = [] + for mreader in daqmetadatareaders: + new_d = mreader.ingest_epochfiles(ef, epoch_id) + docs.extend(new_d if isinstance(new_d, list) else [new_d]) + return docs + + def _ingestOneEpoch( # noqa: N802 (mirrors MATLAB's local function) daqreader: Any, daqmetadatareaders: list[Any], ef: list[str], epoch_id: str, - queue: Any = None, - k: int = 0, epoch_progress: Callable[[float], Any] | None = None, ) -> list[Any]: - """Create the documents for a single epoch (no database access). - - In a worker process, progress is sent as ``(k, fraction)`` on QUEUE; in - the calling process it goes to EPOCH_PROGRESS. - """ - if queue is not None: - epoch_progress = functools.partial(_send_progress, queue, k) + """Create the documents for a single epoch (no database access).""" docs: list[Any] = [] if daqreader is not None: if epoch_progress is None: @@ -205,23 +292,18 @@ def _ingestOneEpoch( # noqa: N802 (mirrors MATLAB's local function) else: new_d = daqreader.ingest_epochfiles(ef, epoch_id, ProgressFcn=epoch_progress) docs.extend(new_d if isinstance(new_d, list) else [new_d]) - for mreader in daqmetadatareaders: - new_d = mreader.ingest_epochfiles(ef, epoch_id) - docs.extend(new_d if isinstance(new_d, list) else [new_d]) + docs.extend(_metadata_docs(daqmetadatareaders, ef, epoch_id)) return docs -def _send_progress(queue: Any, k: int, fraction: float) -> None: - """Send epoch K's progress from a worker process to the caller.""" - queue.put((k, fraction)) - - -def _reader_reports_progress(daqreader: Any) -> bool: - """Does DAQREADER's ingest_epochfiles take a ProgressFcn? +def _reader_plans_segments(daqreader: Any) -> bool: + """Does DAQREADER use ndi_daq_reader_mfdaq's ingest_epochfiles? - True when the method is ndi_daq_reader_mfdaq's own; a subclass that - overrides ingest_epochfiles keeps its own signature and is called without - one. + MATLAB counterpart: ``readerPlansSegments``. Such a reader takes a + ProgressFcn and can be ingested segment by segment + (ingest_epochfiles_plan / ingest_segment). A subclass that overrides + ingest_epochfiles keeps its own signature and is called once per epoch, + without one. """ if daqreader is None: return False diff --git a/src/ndi/daq/fun/ndi_matlab_python_bridge.yaml b/src/ndi/daq/fun/ndi_matlab_python_bridge.yaml index 674d5000..780a0a82 100644 --- a/src/ndi/daq/fun/ndi_matlab_python_bridge.yaml +++ b/src/ndi/daq/fun/ndi_matlab_python_bridge.yaml @@ -14,7 +14,7 @@ functions: - name: ingestEpochDocuments type: function matlab_path: "+ndi/+daq/+fun/ingestEpochDocuments.m" - matlab_last_sync_hash: "b192ae605" + matlab_last_sync_hash: "569fdab67" python_path: "ndi/daq/fun/ingest_epoch_documents.py" input_arguments: - name: daqreader @@ -64,3 +64,15 @@ functions: the workers become a multiprocessing.Manager queue that the caller drains while it waits on the futures; the DefiningClass check becomes an identity check on the ingest_epochfiles function. + + NDI-matlab 569fdab67 (VH-Lab/NDI-matlab#1090) schedules the parallel + work by segment instead of by epoch, so every worker stays busy even + with one long epoch: each epoch whose reader plans segments (its + ingest_epochfiles is ndi_daq_reader_mfdaq's own) is planned in the + caller, and each segment, each metadata reader's work and each other + epoch become one task; segments are attached in plan order as each + epoch's tasks finish. The pool is used whenever there is an epoch, as + in MATLAB. Ported: parfeval + fetchNext becomes a spawn + ProcessPoolExecutor + as_completed, and progress is counted in the + caller from the plan weights, so the Manager queue of #1089's port is + gone. diff --git a/src/ndi/daq/mfdaq.py b/src/ndi/daq/mfdaq.py index d057245f..7058e6ec 100644 --- a/src/ndi/daq/mfdaq.py +++ b/src/ndi/daq/mfdaq.py @@ -628,6 +628,11 @@ def ingest_epochfiles( uint16). The document is not added to any database; its files are temporary files that the database takes over when it is added. + Runs :meth:`ingest_epochfiles_plan`, :meth:`ingest_segment` for each + planned segment, then :meth:`ingest_attach_segments`, so the segments + can also be written one at a time elsewhere (see + ``ndi.daq.fun.ingestEpochDocuments``). + Args: epochfiles: Files for this epoch epoch_id: The epoch id @@ -638,16 +643,63 @@ def ingest_epochfiles( Returns: ndi_document of type daqreader_mfdaq_epochdata_ingested + """ + d, plan = self.ingest_epochfiles_plan(epochfiles, epoch_id) + _report_ingest_progress(ProgressFcn, 0) + weight_total = sum(entry["weight"] for entry in plan) + weight_done = 0.0 + segment_files = [] + for entry in plan: + segment_files.append(self.ingest_segment(epochfiles, entry)) + weight_done += entry["weight"] + _report_ingest_progress(ProgressFcn, weight_done / max(weight_total, 1e-300)) + d = self.ingest_attach_segments(d, plan, segment_files) + _report_ingest_progress(ProgressFcn, 1) + return d - Python-specific Notes: - Sample numbers passed to :meth:`readchannels_epochsamples` are - 1-based, as that method takes them. Digital-output segments are - stored as ``do_group*`` (the document's file list and the ingested - reader use ``do``); MATLAB writes them as ``di_group*``, which - collides with digital-input segments of the same group. + def ingest_epochfiles_plan( + self, + epochfiles: list[str], + epoch_id: str, + ) -> tuple[Any, list[dict[str, Any]]]: """ - import ndicompress + Plan the ingestion of an epoch, segment by segment. + + MATLAB equivalent: ndi.daq.reader.mfdaq/ingest_epochfiles_plan + + Creates the ``daqreader_mfdaq_epochdata_ingested`` document for the + epoch, with its channel list already attached, and the plan: one dict + per data segment still to be written, in the order the segment files + are added to the document. Each entry can be computed independently + (and in another process) with :meth:`ingest_segment`; + :meth:`ingest_attach_segments` then adds the resulting files. + + Each entry has MATLAB's fields: ``kind`` (``"ephys"``, ``"digital"``, + ``"eventmarktext"`` or ``"time"``), ``channeltype``, ``channels``, + ``eventchanneltypes`` (eventmarktext: the type of each channel), + ``s0`` and ``s1`` (first and last sample, 1-based), ``poly``, + ``datasize``, ``underlying_format``, ``docfilename`` (e.g. + ``"ai_group1_seg.nbf_3"``), ``segment`` and ``numSegments`` (this + segment's number within its group), and ``weight`` (its share of the + epoch's work; the weights sum to 1). + + Args: + epochfiles: Files for this epoch + epoch_id: The epoch id + Returns: + ``(d, plan)``. The document is not added to any database. + + Python-specific Notes: + Sample numbers are 1-based, as :meth:`readchannels_epochsamples` + takes them; ``epochtimes2samples`` here is 0-based, so the last + sample is its answer plus 1. Digital-output segments are stored as + ``do_group*`` (the document's file list and the ingested reader use + ``do``); MATLAB writes them as ``di_group*``, which collides with + digital-input segments of the same group. The time group computes + its own sample range, where MATLAB reuses the group planned before + it (the same whenever the time channel shares that group's rate). + """ from ..document import ndi_document from ..file.type.mfdaq_epoch_channel import ndi_file_type_mfdaq__epoch__channel from ..fun.utils import channelname2prefixnumber @@ -715,8 +767,8 @@ def ingest_epochfiles( ci = mec.channel_information - # Step 2: the data, one channel type and group at a time; event, marker - # and text channels are stored together + # Step 2: plan the data, one channel type and group at a time; event, + # marker and text channels are stored together all_types, _ = self.channel_types() emt = ("event", "marker", "text") store_types = [t for t in all_types if t not in emt] + ["eventmarktext"] @@ -726,16 +778,17 @@ def entries_of(type_name: str) -> list[Any]: return [c for c in ci if c.type in emt] return [c for c in ci if c.type == type_name] - weight_total = sum(len(entries_of(t)) for t in store_types) - weight_done = 0 - _report_ingest_progress(ProgressFcn, 0) + # each channel group counts in proportion to its number of channels, + # which tracks the work of reading and compressing it + weight_total = max(sum(len(entries_of(t)) for t in store_types), 1) + plan: list[dict[str, Any]] = [] for type_name in store_types: entries = entries_of(type_name) prefix = self.INGEST_FILE_PREFIX[type_name] for group in sorted({c.group for c in entries}): members = [c for c in entries if c.group == group] - weight_here = len(members) + group_weight = len(members) / weight_total channels_here = [c.number for c in members] underlying_format, mypoly, datasize = self.underlying_datatype( epochfiles, type_name, channels_here @@ -753,72 +806,149 @@ def entries_of(type_name: str) -> list[Any]: f"Sample rates are not all identical for {type_name} group {group}." ) - dataclass_here = members[0].dataclass - if dataclass_here in ("ephys", "digital", "time"): - if not unique_rates: - raise ValueError("Analog records have a NaN sample rate.") - t0t1 = self.t0_t1(epochfiles) - # epochtimes2samples is 0-based; S1 is the 1-based last sample - S1 = ( - int( - np.asarray( - self.epochtimes2samples( - type_name, channels_here, epochfiles, t0t1[0][-1] - ) - ).ravel()[0] - ) - + 1 - ) - segment = ( - self.SAMPLE_DIGITAL_SEGMENT - if dataclass_here == "digital" - else self.SAMPLE_ANALOG_SEGMENT + entry: dict[str, Any] = { + "kind": "", + "channeltype": type_name, + "channels": channels_here, + "eventchanneltypes": [], + "s0": None, + "s1": None, + "poly": np.atleast_2d(mypoly), + "datasize": datasize, + "underlying_format": underlying_format, + "docfilename": "", + "segment": None, + "numSegments": None, + "weight": None, + } + kind = members[0].dataclass + if kind == "eventmarktext": + entry.update( + kind="eventmarktext", + eventchanneltypes=[c.type for c in members], + docfilename=f"evmktx_group{group}_seg.nbf_1", + segment=1, + numSegments=1, + weight=group_weight, ) - s_starts = list(range(1, S1 + 1, segment)) - for s, s0 in enumerate(s_starts, start=1): - logger.info( - "Working on %s ingestion segment %d of %d.", - {"ephys": "ephys/analog"}.get(dataclass_here, dataclass_here), - s, - len(s_starts), - ) - s1 = min(s0 + segment - 1, S1) - data = self.readchannels_epochsamples( - type_name, channels_here, epochfiles, s0, s1 - ) - filename_here = _ingest_temp_name() - if dataclass_here == "ephys": - data = _scaled2underlying(data, mypoly) - _compress_ephys_underlying( - data, filename_here, underlying_format, datasize + plan.append(entry) + continue + if kind not in ("ephys", "digital", "time"): + raise ValueError(f"Unknown channel type {type_name}.") + if not unique_rates: + raise ValueError("Analog records have a NaN sample rate.") + t0t1 = self.t0_t1(epochfiles) + # epochtimes2samples is 0-based; S1 is the 1-based last sample + S1 = ( + int( + np.asarray( + self.epochtimes2samples( + type_name, channels_here, epochfiles, t0t1[0][-1] ) - elif dataclass_here == "digital": - data = _scaled2underlying(data, mypoly) - ndicompress.compress_digital(data, filename_here) - else: - ndicompress.compress_time(np.asarray(data, dtype=float), filename_here) - d.add_file(f"{prefix}_group{group}_seg.nbf_{s}", filename_here + ".nbf.tgz") - _report_ingest_progress( - ProgressFcn, - (weight_done + weight_here * s / len(s_starts)) / weight_total, - ) - elif dataclass_here == "eventmarktext": - logger.info("Working on event/marker/text ingestion segment 1 of 1.") - channeltype = [c.type for c in members] - T, D = self.readevents_epochsamples_native( - channeltype, channels_here, epochfiles, -np.inf, np.inf + ).ravel()[0] ) - filename_here = _ingest_temp_name() - ndicompress.compress_eventmarktext( - channeltype, channels_here, T, D, filename_here + + 1 + ) + segment_length = ( + self.SAMPLE_DIGITAL_SEGMENT if kind == "digital" else self.SAMPLE_ANALOG_SEGMENT + ) + s_starts = list(range(1, S1 + 1, segment_length)) + for s, s0 in enumerate(s_starts, start=1): + plan.append( + { + **entry, + "kind": kind, + "s0": s0, + "s1": min(s0 + segment_length - 1, S1), + "docfilename": f"{prefix}_group{group}_seg.nbf_{s}", + "segment": s, + "numSegments": len(s_starts), + "weight": group_weight / len(s_starts), + } ) - d.add_file(f"evmktx_group{group}_seg.nbf_1", filename_here + ".nbf.tgz") - else: - raise ValueError(f"Unknown channel type {type_name}.") - weight_done += weight_here - _report_ingest_progress(ProgressFcn, weight_done / max(weight_total, 1)) + return d, plan - _report_ingest_progress(ProgressFcn, 1) + def ingest_segment(self, epochfiles: list[str], planEntry: dict[str, Any]) -> str: # noqa: N803 + """ + Read and compress one planned segment of an epoch. + + MATLAB equivalent: ndi.daq.reader.mfdaq/ingest_segment + + Reads the data described by PLANENTRY (one entry of the plan from + :meth:`ingest_epochfiles_plan`) and writes it, compressed, to a new + temporary file. Touches no database, so it can run in a worker + process. + + Args: + epochfiles: Files for this epoch + planEntry: One entry of the plan + + Returns: + The full path of the ``.nbf.tgz`` file written. + """ + import ndicompress + + kind = planEntry["kind"] + filename_here = _ingest_temp_name() + if kind in ("ephys", "digital", "time"): + logger.info( + "Working on %s ingestion segment %d of %d.", + {"ephys": "ephys/analog"}.get(kind, kind), + planEntry["segment"], + planEntry["numSegments"], + ) + data = self.readchannels_epochsamples( + planEntry["channeltype"], + planEntry["channels"], + epochfiles, + planEntry["s0"], + planEntry["s1"], + ) + if kind == "ephys": + _compress_ephys_underlying( + _scaled2underlying(data, planEntry["poly"]), + filename_here, + planEntry["underlying_format"], + planEntry["datasize"], + ) + elif kind == "digital": + ndicompress.compress_digital( + _scaled2underlying(data, planEntry["poly"]), filename_here + ) + else: + ndicompress.compress_time(np.asarray(data, dtype=float), filename_here) + elif kind == "eventmarktext": + logger.info("Working on event/marker/text ingestion segment 1 of 1.") + T, D = self.readevents_epochsamples_native( + planEntry["eventchanneltypes"], planEntry["channels"], epochfiles, -np.inf, np.inf + ) + ndicompress.compress_eventmarktext( + planEntry["eventchanneltypes"], planEntry["channels"], T, D, filename_here + ) + else: + raise ValueError(f"Unknown segment kind {kind}.") + return filename_here + ".nbf.tgz" + + @staticmethod + def ingest_attach_segments( + d: Any, plan: list[dict[str, Any]], segmentFiles: list[str] # noqa: N803 (MATLAB name) + ) -> Any: + """ + Add planned segment files to an ingested-epoch document, in plan order. + + MATLAB equivalent: ndi.daq.reader.mfdaq.ingest_attach_segments + + Args: + d: The document from :meth:`ingest_epochfiles_plan` + plan: Its plan + segmentFiles: ``segmentFiles[t]``, the file :meth:`ingest_segment` + wrote for ``plan[t]`` + + Returns: + The document, with the files added under each entry's ``docfilename``. + """ + for entry, filename in zip(plan, segmentFiles): + d.add_file(entry["docfilename"], filename) return d def underlying_datatype( diff --git a/src/ndi/daq/ndi_matlab_python_bridge.yaml b/src/ndi/daq/ndi_matlab_python_bridge.yaml index 35fe83f5..fed6f8e4 100644 --- a/src/ndi/daq/ndi_matlab_python_bridge.yaml +++ b/src/ndi/daq/ndi_matlab_python_bridge.yaml @@ -396,7 +396,7 @@ classes: - name: mfdaq_reader type: class matlab_path: "+ndi/+daq/+reader/mfdaq.m" - matlab_last_sync_hash: "6ff0264e5" + matlab_last_sync_hash: "5c5784c6a" decision_log: > NDI-matlab 770c1eb3 clamps readchannels_epochsamples_ingested's window to what the epoch has: a finite s0 below the first sample @@ -420,6 +420,13 @@ classes: now ingestion here made the epoch documents without any data, so an ingested epoch could not be read. It now writes channel_list.bin and the compressed segments as MATLAB does (see the method's entry). + + NDI-matlab 569fdab67 and 5c5784c6a (VH-Lab/NDI-matlab#1090) split + ingest_epochfiles into ingest_epochfiles_plan, ingest_segment and the + static ingest_attach_segments, so ingestEpochDocuments can spread one + epoch's segments over the workers. Ported, with the same plan fields + and weights; ingest_epochfiles is now plan, each segment, attach, as + in MATLAB. 5c5784c6a only marks a growing array for Code Analyzer. python_path: "ndi/daq/mfdaq.py" python_class: "ndi_daq_reader_mfdaq" inherits: "ndi.daq.reader" @@ -464,6 +471,71 @@ classes: ('double'), so a reader that keeps analog data in doubles cannot ingest there; here 'float64'/'float32' (and 'double'/'single') are stored as floats. + Since NDI-matlab 569fdab67 (VH-Lab/NDI-matlab#1090) it runs + ingest_epochfiles_plan, ingest_segment for each plan entry + (reporting ProgressFcn by the entries' weights) and + ingest_attach_segments, as MATLAB does. + + - name: ingest_epochfiles_plan + input_arguments: + - name: epochfiles + type_matlab: "cell array of char" + type_python: "list[str]" + - name: epoch_id + type_matlab: "char" + type_python: "str" + output_arguments: + - name: d + type_python: "ndi_document" + - name: plan + type_python: "list[dict[str, Any]]" + decision_log: > + Ported from NDI-matlab 569fdab67 (VH-Lab/NDI-matlab#1090). The + document (with channel_list.bin attached) and one plan entry per + segment, with MATLAB's fields (kind, channeltype, channels, + eventchanneltypes, s0, s1, poly, datasize, underlying_format, + docfilename, segment, numSegments, weight); MATLAB's struct array + becomes a list of dicts. s0/s1 stay 1-based, as MATLAB's, since + readchannels_epochsamples here takes 1-based samples. The + differences listed under ingest_epochfiles (do_group* file names, + the time group's own sample range, float underlying types) are + made here. + + - name: ingest_segment + input_arguments: + - name: epochfiles + type_matlab: "cell array of char" + type_python: "list[str]" + - name: planEntry + type_matlab: "struct (one plan entry)" + type_python: "dict[str, Any]" + output_arguments: + - name: filename + type_python: "str" + decision_log: > + Ported from NDI-matlab 569fdab67 (VH-Lab/NDI-matlab#1090). Reads and + compresses one plan entry to a temporary .nbf.tgz file and returns + its name; touches no database, so it runs in worker processes. + + - name: ingest_attach_segments + input_arguments: + - name: d + type_matlab: "ndi.document" + type_python: "ndi_document" + - name: plan + type_matlab: "struct array" + type_python: "list[dict[str, Any]]" + - name: segmentFiles + type_matlab: "cell array of char" + type_python: "list[str]" + output_arguments: + - name: d + type_python: "ndi_document" + decision_log: > + Ported from NDI-matlab 569fdab67 (VH-Lab/NDI-matlab#1090). Static in + both. Adds segmentFiles[t] to the document as plan[t]'s docfilename, + in plan order, so the document is the same whatever order the + segments were written in. - name: epochclock input_arguments: diff --git a/src/ndi/fun/ndi_matlab_python_bridge.yaml b/src/ndi/fun/ndi_matlab_python_bridge.yaml index 274827db..aab3db12 100644 --- a/src/ndi/fun/ndi_matlab_python_bridge.yaml +++ b/src/ndi/fun/ndi_matlab_python_bridge.yaml @@ -1161,7 +1161,7 @@ functions: - name: epochid2element type: function matlab_path: "+ndi/+fun/+epoch/epochid2element.m" - matlab_last_sync_hash: "fe64a9f53" + matlab_last_sync_hash: "9821a12e6" python_path: "ndi/fun/epoch.py" input_arguments: - name: session @@ -1209,10 +1209,16 @@ functions: input where this returns a dict keyed by epoch id, and MATLAB drives a ProgressBarWindow that has no counterpart in a headless port. + NDI-matlab 774a01f71 and 9821a12e6 change one line: when the input + names the same epoch id more than once, MATLAB now fills every one of + those entries (it used to index a cell with several trues and fail). + Nothing to port: the result here is a dict keyed by epoch id, so a + repeated epoch id is one key holding its elements, and nothing fails. + - name: filename2epochid type: function matlab_path: "+ndi/+fun/+epoch/filename2epochid.m" - matlab_last_sync_hash: "fe64a9f53" + matlab_last_sync_hash: "9821a12e6" python_path: "ndi/fun/epoch.py" input_arguments: - name: session @@ -1238,6 +1244,12 @@ functions: stored daq_system documents -- the same table when the documents are current, and stale when they are not. + NDI-matlab 774a01f71 and 9821a12e6 change one line: when the input + names the same filename more than once, MATLAB now fills every one of + those entries (it used to index a cell with several trues and fail). + Nothing to port: the result here is a dict keyed by filename, so a + repeated filename is one key holding its epoch ids, and nothing fails. + - name: t0_t1cell2array type: function matlab_path: "+ndi/+fun/+doc/t0_t1cell2array.m" diff --git a/src/ndi/ndi_matlab_python_bridge.yaml b/src/ndi/ndi_matlab_python_bridge.yaml index d7db2917..81274d43 100644 --- a/src/ndi/ndi_matlab_python_bridge.yaml +++ b/src/ndi/ndi_matlab_python_bridge.yaml @@ -1320,7 +1320,7 @@ not_applicable: # ========================================================================= - name: setup matlab_path: "+ndi/+setup/" - matlab_last_sync_hash: "03c79fa48" + matlab_last_sync_hash: "fd7ddbf06" python_path: "ndi/setup/" status: ported decision_log: > @@ -1354,6 +1354,13 @@ not_applicable: said so. Changing the signature now would break every caller for a cosmetic gain. + NDI-matlab 4f0e7234d..fd7ddbf06 change only +ndi/+setup/+conv/+vhlab: + the vhlab (Stacy) data importer, its approach and treatment tables, and + its validation and study drivers. Nothing to port: +setup/+conv, the + per-lab importers for one-off dataset conversions, has no counterpart + in this port (src/ndi/setup holds the lab setup functions and the DAQ, + epoch, file and sync helpers), and none of these commits touches those. + # ========================================================================= # ndi.docs namespace — MATLAB documentation generation # ========================================================================= diff --git a/src/ndi/test/__init__.py b/src/ndi/test/__init__.py new file mode 100644 index 00000000..a35fd7ed --- /dev/null +++ b/src/ndi/test/__init__.py @@ -0,0 +1,5 @@ +"""Tools for testing NDI itself on real data. + +MATLAB equivalent: the shipped helpers of +ndi/+test (MATLAB's own unit tests +live in tests/, as Python's do). +""" diff --git a/src/ndi/test/ingest/__init__.py b/src/ndi/test/ingest/__init__.py new file mode 100644 index 00000000..00715b25 --- /dev/null +++ b/src/ndi/test/ingest/__init__.py @@ -0,0 +1,9 @@ +"""Check that ingestion stored a session's data faithfully. + +MATLAB equivalent: +ndi/+test/+ingest +""" + +from .compare import compare +from .mfdaq_compare import mfdaq_compare + +__all__ = ["compare", "mfdaq_compare"] diff --git a/src/ndi/test/ingest/compare.py b/src/ndi/test/ingest/compare.py new file mode 100644 index 00000000..cae1a1a6 --- /dev/null +++ b/src/ndi/test/ingest/compare.py @@ -0,0 +1,146 @@ +"""Compare a raw session with an ingested copy of it, channel by channel. + +MATLAB equivalent: ndi.test.ingest.compare +""" + +from __future__ import annotations + +import math +import warnings +from pathlib import Path +from typing import Any + +from .mfdaq_compare import mfdaq_compare + + +def compare( + rawSession: Any, # noqa: N803 - MATLAB's name + ingestedSession: Any, # noqa: N803 + *, + Mode: str = "full", + ChunkSeconds: float = 60, + SpotSeconds: float = 300, + RelativeTolerance: float = 1e-12, + TimeTolerance: float = 1e-6, + DaqSystemNames: list[str] | tuple[str, ...] | str = (), + Verbose: bool = True, + ShowProgress: bool = True, +) -> tuple[bool, list[dict[str, Any]]]: + """Compare a raw session with an ingested copy of it. + + MATLAB equivalent: ndi.test.ingest.compare + + Compares the data of *rawSession*, a session that has NOT been ingested, + with *ingestedSession*, a copy of it that HAS been ingested (e.g., copy + the session folder and call ``S.ingest()`` on the copy). Each may be an + ndi session or the path of a session folder. Reads from the ingested copy + come only from its ingested documents, so a match shows that ingestion + stored the data faithfully. + + Every multifunction DAQ system of *rawSession* is matched by name to one + in *ingestedSession* and compared with :func:`mfdaq_compare`. A DAQ system + missing from *ingestedSession* is reported (check 'daqsystem'). + + Returns ``(b, report)``: *b* is True if everything compared equal; + *report* is the list described in :func:`mfdaq_compare`, for all DAQ + systems. Nothing stops or prompts on a mismatch. + + Args: + Mode: 'full' compares every sample; 'spot' one random window per + channel type per epoch. + ChunkSeconds: Length of each read in 'full' mode. + SpotSeconds: Length of the window in 'spot' mode. + RelativeTolerance: Largest relative difference allowed between raw and + ingested samples. + TimeTolerance: Largest difference allowed for time channels. + DaqSystemNames: Compare only these DAQ systems (default: all). + Verbose: Print progress and each failed check. + ShowProgress: Show a progress bar per DAQ system. + + Example:: + + b, report = compare(s_raw, s_ingested) + [r for r in report if not r["ok"]] # the checks that failed + """ + from ...daq.mfdaq import ndi_daq_reader_mfdaq + + s1 = _open_session(rawSession) + s2 = _open_session(ingestedSession) + daqs1 = _as_list(s1.daqsystem_load(name="(.*)")) + daqs2 = {str(d.name): d for d in _as_list(s2.daqsystem_load(name="(.*)"))} + wanted = [DaqSystemNames] if isinstance(DaqSystemNames, str) else list(DaqSystemNames) + + progress_bar = None + if ShowProgress: + try: + from ...gui.component.ProgressBarWindow import ndi_gui_component_ProgressBarWindow + + progress_bar = ndi_gui_component_ProgressBarWindow( + "NDI Ingest Compare", Overwrite=False + ) + except Exception as exc: # noqa: BLE001 - no window is not a failed comparison + warnings.warn(f"Could not open a progress window: {exc}", stacklevel=2) + + report: list[dict[str, Any]] = [] + for daq in daqs1: + name = str(daq.name) + if wanted and name not in wanted: + continue + if not isinstance(getattr(daq, "daqreader", None), ndi_daq_reader_mfdaq): + if Verbose: + print(f"Skipping {name}: not an ndi_daq_system_mfdaq.") + continue + other = daqs2.get(name) + if other is None: + msg = "DAQ system is missing from the ingested session." + report.append( + { + "daqsystem": name, + "epoch_id": "", + "check": "daqsystem", + "ok": False, + "message": msg, + "t0_t1": [], + "samples": 0, + "maxAbsDiff": math.nan, + } + ) + if Verbose: + print(f"MISMATCH {name}: {msg}") + continue + _, _, here = mfdaq_compare( + daq, + other, + Mode=Mode, + ChunkSeconds=ChunkSeconds, + SpotSeconds=SpotSeconds, + RelativeTolerance=RelativeTolerance, + TimeTolerance=TimeTolerance, + Verbose=Verbose, + ProgressBar=progress_bar, + ) + report.extend(here) + + b = all(r["ok"] for r in report) + if Verbose: + print(f"{len(report)} checks, {sum(not r['ok'] for r in report)} failed.") + return b, report + + +def _open_session(s: Any) -> Any: + """An ndi session, or a session opened from a folder path.""" + if isinstance(s, (str, Path)): + from ...session.dir import ndi_session_dir + + return ndi_session_dir(str(s)) + if hasattr(s, "daqsystem_load"): + return s + raise TypeError("Each session must be an ndi session or the path of a session folder.") + + +def _as_list(x: Any) -> list[Any]: + if x is None: + return [] + if isinstance(x, (list, tuple)): + return list(x) + return [x] diff --git a/src/ndi/test/ingest/mfdaq_compare.py b/src/ndi/test/ingest/mfdaq_compare.py new file mode 100644 index 00000000..180f878e --- /dev/null +++ b/src/ndi/test/ingest/mfdaq_compare.py @@ -0,0 +1,390 @@ +"""Compare raw and ingested data from two multifunction DAQ systems. + +MATLAB equivalent: ndi.test.ingest.mfdaq_compare +""" + +from __future__ import annotations + +import dataclasses +import math +import random +from typing import Any + +import numpy as np + +#: The continuous channel types compared, sample by sample. +CONTINUOUS_TYPES = ("analog_in", "analog_out", "auxiliary_in", "digital_in", "digital_out", "time") +#: The event channel types compared, over the whole epoch. +EVENT_TYPES = ("event", "marker", "text") + + +def mfdaq_compare( + daq1: Any, + daq2: Any, + *, + Mode: str = "full", + ChunkSeconds: float = 60, + SpotSeconds: float = 300, + RelativeTolerance: float = 1e-12, + TimeTolerance: float = 1e-6, + Verbose: bool = True, + ProgressBar: Any = None, +) -> tuple[bool, list[str], list[dict[str, Any]]]: + """Compare raw and ingested data from two ``ndi_daq_system_mfdaq`` objects. + + MATLAB equivalent: ndi.test.ingest.mfdaq_compare + + Compares the data read through *daq1*, a DAQ system of a session that has + NOT been ingested, with the data read through *daq2*, the same DAQ system + in a copy of that session that HAS been ingested. Reads through *daq2* + come only from the ingested documents, so a match shows that ingestion + stored the data faithfully. + + For every epoch of *daq1*, matched to *daq2*'s epoch by epoch id, it + checks: + + - that the epoch exists in *daq2* and is ingested there; + - the channel list (every channel, matched by name; every field the raw + reader sets, e.g., type, time channel, sample rate); + - every continuous channel type present (analog_in, analog_out, + auxiliary_in, digital_in, digital_out, time): the sample range and the + samples themselves (to *RelativeTolerance*; time channels to + *TimeTolerance*); + - every event channel type present (event, marker, text): the timestamps + and values over the whole epoch. + + Returns ``(b, errmsg, report)``. *b* is True if every check passed; + *errmsg* has one message per failed check; *report* has one dict per + check, with keys ``daqsystem``, ``epoch_id``, ``check`` ('epoch', + 'channel_list' or the channel type), ``ok``, ``message`` ('' when ok), + ``t0_t1``, ``samples`` (the number compared, continuous types) and + ``maxAbsDiff`` (the largest absolute difference, continuous types). + Nothing stops on a mismatch; every difference is reported and the + comparison continues with the next check. + + Args: + Mode: 'full' compares every sample of every epoch, *ChunkSeconds* at a + time; 'spot' compares one randomly placed window of *SpotSeconds* + per channel type per epoch. + ChunkSeconds: Length of each read in 'full' mode (bounds memory use). + SpotSeconds: Length of the window in 'spot' mode. + RelativeTolerance: Largest difference allowed between raw and ingested + samples, relative to the larger of the two values. Scaled channels + can differ in the last bit because the two paths round the scaling + differently; 0 demands exact equality. + TimeTolerance: Largest difference allowed for time channels. + Verbose: Print each epoch as it is compared, and each failed check. + ProgressBar: A ProgressBarWindow to add a bar to (None: no bar). + + Sample numbers here, as everywhere in this port, are 0-based: a raw + reader's ``readchannels_epochsamples`` takes 1-based samples and an + ingested read 0-based ones, and both are called accordingly. + """ + from ...daq.mfdaq import ndi_daq_reader_mfdaq + from ...daq.system import ndi_daq_system + from ...file.navigator import ndi_file_navigator + + if Mode not in ("full", "spot"): + raise ValueError("Mode must be 'full' or 'spot'.") + if ChunkSeconds <= 0 or SpotSeconds <= 0: + raise ValueError("ChunkSeconds and SpotSeconds must be positive.") + if RelativeTolerance < 0 or TimeTolerance < 0: + raise ValueError("RelativeTolerance and TimeTolerance must be nonnegative.") + for daq in (daq1, daq2): + if not isinstance(getattr(daq, "daqreader", None), ndi_daq_reader_mfdaq): + raise TypeError("Each DAQ system must be an ndi_daq_system_mfdaq.") + + options = { + "Mode": Mode, + "ChunkSeconds": ChunkSeconds, + "SpotSeconds": SpotSeconds, + "RelativeTolerance": RelativeTolerance, + "TimeTolerance": TimeTolerance, + "Verbose": Verbose, + } + report: list[dict[str, Any]] = [] + name = str(daq1.name) + + for daq in (daq1, daq2): + cache = getattr(daq.session, "cache", None) + if cache is not None: + cache.clear() + + et1 = daq1.epochtable() + et2 = daq2.epochtable() + by_id2 = {e["epoch_id"]: e for e in et2} + reader1 = daq1.daqreader + reader2 = daq2.daqreader + session2 = daq2.session + + progress_bar, bar_tag = None, "" + if ProgressBar is not None and et1: + progress_bar, bar_tag = ndi_daq_system.addIngestProgressBar( + ProgressBar, f"Comparing {name} ({len(et1)} epochs)" + ) + + def add(epoch_id, check, ok, msg, t0_t1, n=0, mx=math.nan): + report.append( + { + "daqsystem": name, + "epoch_id": epoch_id, + "check": check, + "ok": bool(ok), + "message": msg, + "t0_t1": t0_t1, + "samples": n, + "maxAbsDiff": mx, + } + ) + if not ok and Verbose: + print(f" MISMATCH {epoch_id} {check}: {msg}") + + for i, e1 in enumerate(et1, start=1): + epoch_id = e1["epoch_id"] + t0_t1 = list(e1["t0_t1"][0]) + files1 = e1["underlying_epochs"]["underlying"] + if Verbose: + print(f"Comparing {name} epoch {epoch_id} ({i} of {len(et1)})...") + try: + e2 = by_id2.get(epoch_id) + if e2 is None: + add(epoch_id, "epoch", False, "Epoch is missing from the ingested session.", t0_t1) + continue + files2 = e2["underlying_epochs"]["underlying"] + if not ndi_file_navigator.isingested(files2): + add(epoch_id, "epoch", False, "Epoch is not ingested in the second session.", t0_t1) + continue + add(epoch_id, "epoch", True, "", t0_t1) + + c1 = reader1.getchannelsepoch(files1) + c2 = reader2.getchannelsepoch_ingested(files2, session2) + ok, msg = _compare_channel_lists(c1, c2) + add(epoch_id, "channel_list", ok, msg, t0_t1) + + # continuous channels, grouped by sample rate so each read is uniform + for channeltype in CONTINUOUS_TYPES: + channels = [c.number for c in c1 if c.type == channeltype] + if not channels: + continue + sr = np.ravel( + np.asarray(reader1.samplerate(files1, channeltype, channels), dtype=float) + ) + for rate in sorted({r for r in sr if not np.isnan(r)}): + here = [ch for ch, r in zip(channels, sr) if r == rate] + ok, msg, n, mx = _compare_continuous( + reader1, + reader2, + session2, + files1, + files2, + channeltype, + here, + rate, + t0_t1, + options, + ) + add(epoch_id, channeltype, ok, msg, t0_t1, n, mx) + + # event channels, over the whole epoch + for channeltype in EVENT_TYPES: + channels = [c.number for c in c1 if c.type == channeltype] + if not channels: + continue + ok, msg = _compare_events( + reader1, reader2, session2, files1, files2, channeltype, channels, t0_t1 + ) + add(epoch_id, channeltype, ok, msg, t0_t1) + finally: + ndi_daq_system.updateIngestProgressBar(progress_bar, bar_tag, i / len(et1)) + + errmsg = [ + f"{r['daqsystem']} {r['epoch_id']} {r['check']}: {r['message']}" + for r in report + if not r["ok"] + ] + return all(r["ok"] for r in report), errmsg, report + + +def _compare_channel_lists(c1: list[Any], c2: list[Any]) -> tuple[bool, str]: + """Same channels (by name) and the same values in every shared field. + + MATLAB compares the fields the two structure arrays share; a raw reader's + list there often has only name, type and time_channel. Python's + ChannelInfo always has every field, holding its default where the reader + set nothing, so a field counts as the raw list's when some raw channel + holds a value other than the default.""" + d1 = {c.name: _as_dict(c) for c in c1} + d2 = {c.name: _as_dict(c) for c in c2} + if sorted(d1) != sorted(d2): + only1 = ", ".join(sorted(set(d1) - set(d2))) + only2 = ", ".join(sorted(set(d2) - set(d1))) + return False, f"Channel names differ. Only raw: {only1}; only ingested: {only2}." + defaults = {} + if c1 and dataclasses.is_dataclass(c1[0]): + defaults = { + f.name: f.default + for f in dataclasses.fields(c1[0]) + if f.default is not dataclasses.MISSING + } + shared = set.intersection(*(set(d) for d in [*d1.values(), *d2.values()])) + fields = sorted( + f + for f in shared + if f not in defaults or any(not _same(d[f], defaults[f]) for d in d1.values()) + ) + bad = [f for f in fields if any(not _same(d1[n][f], d2[n][f]) for n in d1)] + if bad: + return False, f"Channel list fields differ: {', '.join(bad)}." + return True, "" + + +def _as_dict(c: Any) -> dict[str, Any]: + if hasattr(c, "to_dict"): + return c.to_dict() + if dataclasses.is_dataclass(c): + return dataclasses.asdict(c) + return dict(c) if isinstance(c, dict) else dict(vars(c)) + + +def _same(a: Any, b: Any) -> bool: + """Equality in which NaN equals NaN (MATLAB's isequaln) and 1 == 1.0.""" + if isinstance(a, float) and isinstance(b, float) and math.isnan(a) and math.isnan(b): + return True + return a == b + + +def _compare_continuous( + reader1, reader2, session2, files1, files2, channeltype, channels, rate, t0_t1, options +) -> tuple[bool, str, int, float]: + """Compare the samples of one channel type at one sample rate.""" + n_compared = 0 + max_abs_diff = 0.0 + is_time = channeltype == "time" + + s1 = np.ravel(reader1.epochtimes2samples(channeltype, channels, files1, np.asarray(t0_t1))) + s2 = np.ravel( + reader2.epochtimes2samples_ingested( + channeltype, channels, files2, np.asarray(t0_t1), session2 + ) + ) + s1 = [int(round(float(x))) for x in s1] + s2 = [int(round(float(x))) for x in s2] + if s1 != s2: + return ( + False, + f"Sample range differs: raw {s1[0]}..{s1[-1]}, ingested {s2[0]}..{s2[-1]}.", + 0, + max_abs_diff, + ) + s_first, s_last = s1[0], s1[-1] + + if options["Mode"] == "full": + chunk = max(1, round(options["ChunkSeconds"] * rate)) + starts = range(s_first, s_last + 1, chunk) + else: + chunk = max(1, round(options["SpotSeconds"] * rate)) + span = max(0, s_last - s_first + 1 - chunk) + starts = [s_first + random.randint(0, span)] + + for s0 in starts: + s_end = min(s0 + chunk - 1, s_last) + try: + # a raw reader's samples are 1-based + d1 = np.asarray( + reader1.readchannels_epochsamples(channeltype, channels, files1, s0 + 1, s_end + 1) + ) + except Exception as exc: # noqa: BLE001 - report it and go on + return ( + False, + f"Reading the raw data failed at samples {s0}..{s_end}: {exc}", + n_compared, + max_abs_diff, + ) + try: + d2 = np.asarray( + reader2.readchannels_epochsamples_ingested( + channeltype, channels, files2, s0, s_end, session2 + ) + ) + except Exception as exc: # noqa: BLE001 - report it and go on + return ( + False, + f"Reading the ingested data failed at samples {s0}..{s_end}: {exc}", + n_compared, + max_abs_diff, + ) + if d1.shape != d2.shape: + return ( + False, + f"Sizes differ at samples {s0}..{s_end}: raw {list(d1.shape)}, " + f"ingested {list(d2.shape)}.", + n_compared, + max_abs_diff, + ) + n_compared += d1.shape[0] if d1.ndim else 0 + if d1.size == 0: + continue + a = d1.astype(float).reshape(d1.shape[0], -1) + b = d2.astype(float).reshape(d2.shape[0], -1) + d = np.abs(a - b) + d[np.isnan(a) & np.isnan(b)] = 0 + d[np.isnan(a) ^ np.isnan(b)] = np.inf + max_abs_diff = max(max_abs_diff, float(np.max(d))) + if is_time: + allowed = options["TimeTolerance"] + else: + allowed = options["RelativeTolerance"] * np.maximum(np.abs(a), np.abs(b)) + too_large = d > allowed + if np.any(too_large): + excess = np.where(too_large, d, -np.inf) + row, col = np.unravel_index(int(np.argmax(excess)), excess.shape) + return ( + False, + f"Data differ (max |difference| {max_abs_diff:g}) starting by sample " + f"{s0 + int(row)}, channel {channels[int(col)]} (read window {s0}..{s_end}).", + n_compared, + max_abs_diff, + ) + return True, "", n_compared, max_abs_diff + + +def _compare_events( + reader1, reader2, session2, files1, files2, channeltype, channels, t0_t1 +) -> tuple[bool, str]: + """Compare the timestamps and values of event-type channels.""" + try: + t1, d1 = reader1.readevents_epochsamples(channeltype, channels, files1, t0_t1[0], t0_t1[1]) + t2, d2 = reader2.readevents_epochsamples_ingested( + channeltype, channels, files2, t0_t1[0], t0_t1[1], session2 + ) + except Exception as exc: # noqa: BLE001 - report it and go on + return False, f"Reading events failed: {exc}" + if not _events_equal(t1, t2): + return False, "Event timestamps differ." + if not _events_equal(d1, d2): + return False, "Event values (codes/text) differ." + return True, "" + + +def _events_equal(a: Any, b: Any) -> bool: + """isequaln for event arrays: per channel when given as lists; every empty + array equal to every other (MATLAB's normalizeEmpty); NaN equal to NaN.""" + if isinstance(a, (list, tuple)) and isinstance(b, (list, tuple)): + a_is_channels = all(isinstance(x, (list, tuple, np.ndarray)) for x in a) + b_is_channels = all(isinstance(x, (list, tuple, np.ndarray)) for x in b) + if a_is_channels and b_is_channels: + return len(a) == len(b) and all(_events_equal(x, y) for x, y in zip(a, b)) + x = np.asarray(a, dtype=object if _is_text(a) else None) + y = np.asarray(b, dtype=object if _is_text(b) else None) + if x.size == 0 and y.size == 0: + return True + if x.dtype == object or y.dtype == object: + return [str(v) for v in np.ravel(x)] == [str(v) for v in np.ravel(y)] + if np.squeeze(x).shape != np.squeeze(y).shape: + return False + return bool(np.array_equal(np.squeeze(x), np.squeeze(y), equal_nan=True)) + + +def _is_text(a: Any) -> bool: + flat = a if isinstance(a, (list, tuple)) else np.ravel(np.asarray(a, dtype=object)) + return any(isinstance(v, str) for v in flat) diff --git a/tests/test_ingest_compare.py b/tests/test_ingest_compare.py new file mode 100644 index 00000000..17a47b0c --- /dev/null +++ b/tests/test_ingest_compare.py @@ -0,0 +1,204 @@ +"""ndi.test.ingest.compare / mfdaq_compare, and the plan/segment/attach split. + +MATLAB counterparts: tests/+ndi/+unittest/+ingestion/ingestCompareTest.m +(VH-Lab/NDI-matlab#1091) and ingestionProgressParallel.m's +testPlanAndSegmentsMatchIngestEpochfiles (VH-Lab/NDI-matlab#1090). + +The compare tests build the example Intan session, copy it, ingest the copy, +and compare the two: every check should pass in 'full' and 'spot' modes. Then +an epoch is added to the raw session only, which the comparison must report as +missing from the ingested copy. +""" + +from __future__ import annotations + +import shutil + +import numpy as np +import pytest + +from ndi.daq.mfdaq import ndi_daq_reader_mfdaq +from ndi.gui.component.ProgressBarWindow import ndi_gui_component_ProgressBarWindow as PBW +from ndi.test.ingest import compare, mfdaq_compare +from tests.test_ingest_progress_parallel import _have_intan, _intan_session, _prefs + +ndicompress = pytest.importorskip("ndicompress") + +pytestmark = pytest.mark.skipif(not _have_intan(), reason="Intan example data (.rhd) not found") + + +@pytest.fixture(autouse=True) +def _silent(): + previous = PBW.silentModeDefault() + PBW.silentModeDefault(True) + yield + PBW.silentModeDefault(previous) + + +def _daq(session): + d = session.daqsystem_load(name="intan1") + return d[0] if isinstance(d, list) else d + + +def _raw_and_ingested(tmp_path, monkeypatch, **prefs): + """The test session and an ingested copy of it.""" + from ndi.session.dir import ndi_session_dir + + _prefs(monkeypatch, **{"OpenPoolIfNone": False, **prefs}) + raw = _intan_session(tmp_path / "raw") + # build the epoch table first, so the epoch ids are written to disk and + # the copy shares them + _daq(raw).epochtable() + shutil.copytree(tmp_path / "raw", tmp_path / "ingested") + raw = ndi_session_dir("exp1", tmp_path / "raw") + ing = ndi_session_dir("exp1", tmp_path / "ingested") + ok, msg = ing.ingest() + assert ok, msg + return raw, ing + + +def _failures(report): + return "\n".join( + f"{r['daqsystem']} {r['epoch_id']} {r['check']}: {r['message']}" + for r in report + if not r["ok"] + ) + + +class TestIngestCompare: + def test_ingested_copy_matches_raw(self, tmp_path, monkeypatch): + raw, ing = _raw_and_ingested(tmp_path, monkeypatch) + b, report = compare(raw, ing, ShowProgress=False, Verbose=False) + assert b, _failures(report) + + checks = [r["check"] for r in report] + assert "channel_list" in checks + analog = [r for r in report if r["check"] == "analog_in"] + assert analog, "analog_in should be compared" + reader = _daq(raw).daqreader + files = { + e["epoch_id"]: e["underlying_epochs"]["underlying"] for e in _daq(raw).epochtable() + } + for r in analog: + # 'full' covers every sample of the epoch + s = np.ravel( + reader.epochtimes2samples("analog_in", [1], files[r["epoch_id"]], r["t0_t1"]) + ) + assert r["samples"] == int(s[-1]) - int(s[0]) + 1 + assert r["maxAbsDiff"] == 0 + + def test_spot_mode_and_paths(self, tmp_path, monkeypatch): + # sessions may also be given as folder paths; 'spot' mode passes too + raw, ing = _raw_and_ingested(tmp_path, monkeypatch) + b, report = compare( + raw.path, ing.path, Mode="spot", SpotSeconds=1, ShowProgress=False, Verbose=False + ) + assert b, _failures(report) + assert any(r["check"] == "analog_in" for r in report) + + def test_mfdaq_compare_outputs(self, tmp_path, monkeypatch): + raw, ing = _raw_and_ingested(tmp_path, monkeypatch) + b, errmsg, report = mfdaq_compare(_daq(raw), _daq(ing), Verbose=False) + assert b, errmsg + assert errmsg == [] + assert report + + def test_a_changed_sample_is_reported(self, tmp_path, monkeypatch): + raw, ing = _raw_and_ingested(tmp_path, monkeypatch) + real = ndi_daq_reader_mfdaq.readchannels_epochsamples_ingested + + def corrupted(self, channeltype, channel, *args, **kwargs): + data = np.array(real(self, channeltype, channel, *args, **kwargs), dtype=float) + if channeltype == "analog_in" and data.size: + data[0, 0] += 1 + return data + + monkeypatch.setattr(ndi_daq_reader_mfdaq, "readchannels_epochsamples_ingested", corrupted) + b, errmsg, report = mfdaq_compare(_daq(raw), _daq(ing), Verbose=False) + assert not b + bad = [r for r in report if not r["ok"]] + assert {r["check"] for r in bad} == {"analog_in"} + assert "Data differ" in bad[0]["message"] + + def test_missing_epoch_is_reported(self, tmp_path, monkeypatch): + raw, ing = _raw_and_ingested(tmp_path, monkeypatch) + # an epoch that only the raw session has (one epoch per directory in + # this port; see _intan_session). Only the data files are copied: the + # hidden .epochid.ndi file would give the new epoch the old one's id. + shutil.copytree( + tmp_path / "raw" / "t00002", + tmp_path / "raw" / "t00003", + ignore=shutil.ignore_patterns(".*"), + ) + raw.cache.clear() + b, report = compare(raw, ing, ShowProgress=False, Verbose=False) + assert not b, "An epoch missing from the ingested copy should fail the comparison." + missing = [r for r in report if r["check"] == "epoch" and not r["ok"]] + assert len(missing) == 1 + assert "missing" in missing[0]["message"] + + def test_parallel_ingest_by_segment_matches_raw(self, tmp_path, monkeypatch): + # small segments on two workers: the segments of both epochs are + # written out of order across the pool and attached afterwards + monkeypatch.setattr(ndi_daq_reader_mfdaq, "SAMPLE_ANALOG_SEGMENT", 7001) + monkeypatch.setattr(ndi_daq_reader_mfdaq, "SAMPLE_DIGITAL_SEGMENT", 13003) + raw, ing = _raw_and_ingested( + tmp_path, + monkeypatch, + OpenPoolIfNone=True, + UseDefaultNumberOfWorkers=False, + NumberOfWorkers=2, + ) + b, report = compare(raw, ing, ChunkSeconds=0.5, ShowProgress=False, Verbose=False) + assert b, _failures(report) + + def test_daq_system_names_filter(self, tmp_path, monkeypatch): + raw, ing = _raw_and_ingested(tmp_path, monkeypatch) + b, report = compare( + raw, ing, DaqSystemNames=["nosuchdaq"], ShowProgress=False, Verbose=False + ) + assert b + assert report == [] + + +def _ai_data(filename): + """The analog data in one segment file (.nbf.tgz).""" + data, _ = ndicompress.expand_ephys(filename) + return np.asarray(data) + + +class TestPlanAndSegments: + def test_plan_and_segments_match_ingest_epochfiles(self, tmp_path, monkeypatch): + """Planning, then ingesting the segments in any order, then attaching + them gives the same document as ingest_epochfiles (the parallel + scheduler relies on this).""" + monkeypatch.setattr(ndi_daq_reader_mfdaq, "SAMPLE_ANALOG_SEGMENT", 20000) + session = _intan_session(tmp_path / "exp1") + reader = _daq(session).daqreader + files = _daq(session).epochtable()[0]["underlying_epochs"]["underlying"] + + whole = reader.ingest_epochfiles(files, "epoch_x") + + d, plan = reader.ingest_epochfiles_plan(files, "epoch_x") + assert len(plan) > 1 + assert abs(sum(p["weight"] for p in plan) - 1) < 1e-12 + segment_files = [None] * len(plan) + for k in reversed(range(len(plan))): + segment_files[k] = reader.ingest_segment(files, plan[k]) + split = ndi_daq_reader_mfdaq.ingest_attach_segments(d, plan, segment_files) + + def names(doc): + return [f["name"] for f in doc.document_properties["files"]["file_info"]] + + def location(doc, name): + info = next( + f for f in doc.document_properties["files"]["file_info"] if f["name"] == name + ) + return info["locations"][0]["location"] + + assert names(split) == names(whole) + for name in names(whole): + if name.startswith("ai_"): + np.testing.assert_array_equal( + _ai_data(location(split, name)), _ai_data(location(whole, name)) + )