diff --git a/CHANGELOG.md b/CHANGELOG.md index ae89a7f9..fcfb5bb5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,34 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/). clones), and the `ndi.fun.epoch` and `ndi.setup` entries are re-examined against NDI-matlab main. +- **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. + +- **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 235897fe..6f58a7b5 100644 --- a/src/ndi/daq/fun/ingest_epoch_documents.py +++ b/src/ndi/daq/fun/ingest_epoch_documents.py @@ -19,11 +19,12 @@ 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 typing import Any from pydantic import ConfigDict, validate_call @@ -37,7 +38,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 +68,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. @@ -81,27 +85,44 @@ 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. + 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: 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) + # 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 @@ -112,27 +133,26 @@ def tick() -> None: if number_of_workers > 0: reader, mreaders = payload - ctx = multiprocessing.get_context("spawn") - 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): k - for k, (ef, eid) in enumerate(zip(epochFiles, epochIds)) - } - for future in as_completed(futures): - docs_per_epoch[futures[future]] = future.result() - tick() - except BaseException: - executor.shutdown(wait=True, cancel_futures=True) - raise - 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 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 + ) + update_progress(k, 1) return docs_per_epoch, False @@ -140,20 +160,159 @@ def tick() -> 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 + daqreader: Any, + daqmetadatareaders: list[Any], + ef: list[str], + epoch_id: str, + epoch_progress: Callable[[float], Any] | None = None, ) -> list[Any]: """Create the documents for a single epoch (no database access).""" docs: list[Any] = [] if daqreader is not None: - new_d = daqreader.ingest_epochfiles(ef, epoch_id) - docs.extend(new_d if isinstance(new_d, list) else [new_d]) - for mreader in daqmetadatareaders: - new_d = mreader.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]) + docs.extend(_metadata_docs(daqmetadatareaders, ef, epoch_id)) return docs +def _reader_plans_segments(daqreader: Any) -> bool: + """Does DAQREADER use ndi_daq_reader_mfdaq's ingest_epochfiles? + + 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 + 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..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: "8feb9166e" + matlab_last_sync_hash: "569fdab67" 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,24 @@ 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. + + 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/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..7058e6ec 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,366 @@ 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. + + 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 + 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 + """ + 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 + + def ingest_epochfiles_plan( + self, + epochfiles: list[str], + epoch_id: str, + ) -> tuple[Any, list[dict[str, Any]]]: + """ + 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 + + # 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: 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"] + + 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] + + # 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] + 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 + ) + # 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}." + ) + + 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, + ) + 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] + ) + ).ravel()[0] + ) + + 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), + } + ) + return d, plan + + 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( self, epochfiles: list[str], @@ -1290,3 +1658,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..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: "770c1eb3" + 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 @@ -414,11 +414,129 @@ 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). + + 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" 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. + 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: - 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/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)) + ) 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