diff --git a/dascore/io/sintela/protobuf_utils.py b/dascore/io/sintela/protobuf_utils.py index 2a5205a3..68a4cbf8 100644 --- a/dascore/io/sintela/protobuf_utils.py +++ b/dascore/io/sintela/protobuf_utils.py @@ -78,6 +78,34 @@ } # sample_count is a uint32 field; contiguity checks wrap at this bound. _SAMPLE_COUNT_MODULUS = 2**32 +# Bytes of a payload pulled by metadata-only paths, where protobuf puts the +# `header` submessage; ample for any header, a fraction of a typical packet. +_HEADER_PREFIX_SIZE = 1 << 16 +# Protobuf wire key for field 1, wire type 2 (length-delimited): `header`. +_HEADER_FIELD_KEY = 0x0A +# A seek skips a payload's transfer but costs a platter rotation (~18 ms +# measured on a USB HDD), more than sequentially reading the ~1.6 MB it saves. +# Smaller remainders are read and dropped, keeping readahead engaged. +_SEEK_SKIP_THRESHOLD = 1 << 21 +# Window for the backwards search from EOF for the final record's framing: +# the first packet's size plus slack for jitter, and a wider second attempt +# used only when that finds nothing (a recording opening with a short warm-up +# packet does not predict its own final packet's size). +_TAIL_SEARCH_SLACK = 1 << 12 +_MAX_TAIL_SEARCH = 1 << 23 +# A base-128 varint encodes at most a 64-bit value, so it never exceeds this. +_MAX_VARINT_BYTES = 10 +# How far endpoint timestamps may drift from the span their sample counts +# imply. Recorders resynchronize their clock against the counter, so valid +# files do NOT agree exactly: across a 10,308 file archive a quarter drifted, +# worst 56 ms over a 60 s span (~0.09%). Set an order of magnitude above that, +# since a false rejection silently costs a full read while what is worth +# catching is wrong by seconds to hours. The per-packet term is one packet, +# not several: a declared length is corruption-controlled, so scaling by it +# would let a malformed file widen its own budget. +_ENDPOINT_TIME_TOLERANCE = 1 / 100 +_ENDPOINT_TIME_PACKETS = 1 +_ENDPOINT_TIME_FLOOR_NS = 10_000_000 DIMS_TS = ("time", "distance") DIMS_BAND = ("time", "distance", "band") DIMS_FFT = ("time", "distance", "frequency") @@ -154,8 +182,67 @@ def _timestamp_to_dt64(timestamp) -> np.datetime64 | None: return np.datetime64(seconds, "s") + np.timedelta64(nanos, "ns") -def _iter_envelope_records(resource, *, strict: bool) -> Iterator[EnvelopeRecord]: - """Read all MTLV envelope records from a binary stream.""" +def _read_varint(buf: bytes, pos: int) -> tuple[int | None, int]: + """ + Return the varint starting at ``pos`` and the position after it. + + Gives up after ``_MAX_VARINT_BYTES``: a valid varint never runs longer, and + an unterminated run of continuation bytes would otherwise build an integer + as wide as the buffer, one 7-bit shift at a time. + """ + value = shift = 0 + end = min(len(buf), pos + _MAX_VARINT_BYTES) + while pos < end: + byte = buf[pos] + pos += 1 + value |= (byte & 0x7F) << shift + if not byte & 0x80: + return value, pos + shift += 7 + return None, pos + + +def _leading_header_bytes(payload_prefix: bytes) -> bytes | None: + """ + Return the serialized ``header`` submessage from the front of a payload. + + Every Sintela data packet carries ``header`` in field 1. The wire format + permits any field order, but implementations emit ascending field numbers, + so in practice the header sits at the front of the payload and can be + recovered from a short prefix -- which is what lets the metadata paths + ignore the samples behind it, ~99.99% of a recording. + + Returns None when the prefix does not open with field 1 or does not hold + the whole submessage, so a differently ordered payload is still handled + correctly, just read in full. + """ + if not payload_prefix or payload_prefix[0] != _HEADER_FIELD_KEY: + return None + length, pos = _read_varint(payload_prefix, 1) + if length is None or pos + length > len(payload_prefix): + return None + return payload_prefix[: pos + length] + + +def _stream_size(resource) -> int: + """Return the total length of a seekable binary stream.""" + resource.seek(0, 2) + return resource.tell() + + +def _iter_envelope_records( + resource, *, strict: bool, headers_only: bool = False +) -> Iterator[EnvelopeRecord]: + """ + Read all MTLV envelope records from a binary stream. + + With ``headers_only`` a data packet's payload is truncated to its leading + ``header`` submessage and the samples behind it skipped, so a metadata-only + pass never holds them. A payload whose header cannot be recovered from the + prefix is re-read whole, so a yielded payload always parses on its own. + META and unrecognized records are never truncated: field 1 of a META + payload is an ordinary string, not a header. + """ def _stop(message): """Raise in strict mode; otherwise signal a clean stop to the caller.""" @@ -163,8 +250,40 @@ def _stop(message): raise InvalidFiberFileError(message) return True + # Bytes of the current payload still to be stepped over. Deferred to the + # top of the next iteration so a consumer that stops early -- format + # detection returns on the first data tag -- never pays for a record it + # did not ask for. + pending = 0 + + def _advance(): + """Step over the previous payload's tail, seeking only when it pays.""" + nonlocal pending + if pending >= _SEEK_SKIP_THRESHOLD: + resource.seek(pending, 1) + elif pending > 0: + # Streaming a small remainder beats a seek on spinning media and + # keeps the kernel's sequential readahead engaged. + resource.read(pending) + pending = 0 + + def _read_header_payload(offset, size): + """Return just the header submessage of the payload at `offset`.""" + nonlocal pending + prefix = resource.read(min(size, _HEADER_PREFIX_SIZE)) + header = _leading_header_bytes(prefix) + if header is None: + resource.seek(offset) + return resource.read(size) + pending = size - len(prefix) + return header + + # Payload truncation is judged against the stream's length rather than a + # short read, because a skipped tail never produces one. + file_size = _stream_size(resource) resource.seek(0) while True: + _advance() # Each record opens with a 4-byte magic word; an empty read here is a # clean end-of-file, while a short or wrong magic is a malformed record. magic = resource.read(4) @@ -185,15 +304,21 @@ def _stop(message): size = struct.unpack(" str | None: """Return the first supported data tag in a file without using protobuf.""" - for record in _iter_envelope_records(resource, strict=False): + for record in _iter_envelope_records(resource, strict=False, headers_only=True): if record.tag == META_TAG: continue if record.tag in _TAG_TO_PACKET: @@ -204,6 +329,77 @@ def get_supported_family_tag(resource) -> str | None: return None +def _parse_frame(buf: bytes, index: int = 0) -> tuple[str, int] | None: + """Return the (tag, payload size) of the MTLV framing at ``buf[index:]``.""" + if index + 12 > len(buf): + return None + if struct.unpack(" file_size: + return None + return (*frame, buf[12:]) + + +def _find_first_data_record(resource, file_size: int): + """ + Walk from the start of the file to its first data record. + + Returns the META metadata picked up on the way plus that record as + (tag, size, payload prefix), or None when the framing is unreadable or the + file holds no data records at all. + """ + meta = ParsedMeta() + offset = 0 + while offset < file_size: + record = _read_record_head(resource, offset, file_size) + if record is None: + return None + tag, size, _prefix = record + if tag in _TAG_TO_PACKET: + return meta, record + if tag == META_TAG: + resource.seek(offset + 12) + meta = _parse_meta(resource.read(size)) + offset += 12 + size + return None + + +def _find_last_record(resource, file_size: int, window: int): + """ + Locate the file's final record by searching backwards from its end. + + The real last record is the one whose framing lands exactly on EOF. A + stray magic word in sample data would need a length hitting that same + byte, so the first match walking backwards is the record itself. + """ + start = max(0, file_size - window) + resource.seek(start) + buf = resource.read(file_size - start) + index = len(buf) + while (index := buf.rfind(struct.pack(" np.datetime64 | None: ) +def _parse_packet(tag: str, payload: bytes, messages, decode_error): + """Decode one data record's payload into its packet message.""" + msg = messages[_TAG_TO_PACKET[tag]]() + try: + msg.ParseFromString(payload) + except decode_error as exc: + out = f"Failed to parse Sintela protobuf {tag} payload: {exc}" + raise InvalidFiberFileError(out) from exc + return msg + + def _parse_records( records: Iterable[EnvelopeRecord], *, scan_mode: bool = False ) -> tuple[list[Any], ParsedMeta]: @@ -477,16 +684,10 @@ def _parse_records( if tag == META_TAG: meta = _parse_meta(record.payload) continue - packet_name = _TAG_TO_PACKET.get(tag) - if packet_name is None: + if tag not in _TAG_TO_PACKET: first_unsupported_tag = first_unsupported_tag or tag continue - msg = messages[packet_name]() - try: - msg.ParseFromString(record.payload) - except decode_error as exc: - out = f"Failed to parse Sintela protobuf {tag} payload: {exc}" - raise InvalidFiberFileError(out) from exc + msg = _parse_packet(tag, record.payload, messages, decode_error) if scan_mode: # Omitting the sample fields from the scan descriptor stops them # being *decoded*, but protobuf still retains their raw bytes as @@ -713,8 +914,21 @@ class TimeseriesMetadata(_PacketMetadata): channel_step: int @classmethod - def from_parsed(cls, parsed: list[tuple[str, Any]], meta: ParsedMeta): - """Validate timeseries headers and build shared attrs/coords.""" + def from_parsed( + cls, + parsed: list[tuple[str, Any]], + meta: ParsedMeta, + *, + total_samples: int | None = None, + ): + """ + Validate timeseries headers and build shared attrs/coords. + + ``total_samples`` lets a caller holding only the file's first and last + packets supply the length it derived from their sample counts. Without + it ``parsed`` is taken to be every packet in the file and their + contiguity is checked here. + """ headers = [msg.header for _tag, msg in parsed] common_headers = [h.common_header for h in headers] fields = [ @@ -750,22 +964,24 @@ def from_parsed(cls, parsed: list[tuple[str, Any]], meta: ParsedMeta): raise InvalidFiberFileError( "Dropped samples in Sintela protobuf stream." ) - sample_counts = [int(h.sample_count) for h in headers] - num_samples_per_packet = [f.num_samples for f in fields] - for current, nxt, count in zip( - sample_counts, - sample_counts[1:], - num_samples_per_packet[:-1], - strict=False, - ): - # sample_count is a uint32 on the wire, so a long acquisition (or a - # recorder-wide counter) can wrap to zero mid-file. Compare modulo - # 2**32 so a wrapped-but-contiguous packet is not read as a gap. - if (current + count) % _SAMPLE_COUNT_MODULUS != nxt: - raise InvalidFiberFileError( - "Non-contiguous Sintela protobuf sample counts." - ) - total_samples = sum(num_samples_per_packet) + if total_samples is None: + sample_counts = [int(h.sample_count) for h in headers] + num_samples_per_packet = [f.num_samples for f in fields] + for current, nxt, count in zip( + sample_counts, + sample_counts[1:], + num_samples_per_packet[:-1], + strict=False, + ): + # sample_count is a uint32 on the wire, so a long acquisition + # (or a recorder-wide counter) can wrap to zero mid-file. + # Compare modulo 2**32 so a wrapped-but-contiguous packet is + # not read as a gap. + if (current + count) % _SAMPLE_COUNT_MODULUS != nxt: + raise InvalidFiberFileError( + "Non-contiguous Sintela protobuf sample counts." + ) + total_samples = sum(num_samples_per_packet) first_time = _common_header_time(common_headers[0]) if first_time is None: raise InvalidFiberFileError("Missing Sintela protobuf start time.") @@ -801,31 +1017,96 @@ def from_parsed(cls, parsed: list[tuple[str, Any]], meta: ParsedMeta): attrs=attrs, ) + def _fill_packet(self, data, index: int, tag: str, msg) -> int: + """Copy one packet's samples into ``data`` at ``index``, return its rows.""" + packet = np.asarray(msg.samples, dtype=np.float32) + rows = int(msg.header.num_samples) + expected = rows * self.num_channels + if not packet.size and msg.raw_frames: + # Timeseries packets may carry samples in the packed `raw_frames` + # blob instead of the repeated `samples` field. That encoding is + # undocumented here, so fail with a specific message rather than a + # confusing payload-size mismatch. + msg_ = ( + f"Sintela protobuf {tag} packets store samples in " + "raw_frames, which DASCore cannot yet decode." + ) + raise InvalidFiberFileError(msg_) + if packet.size != expected: + raise InvalidFiberFileError( + "Unexpected Sintela protobuf TS sample payload size." + ) + data[index : index + rows] = packet.reshape(rows, self.num_channels) + return rows + def decode(self, parsed: list[tuple[str, Any]]): """Decode timeseries packets into data, coords, and attrs.""" data = np.empty(self.shape, dtype=self.dtype) index = 0 for tag, msg in parsed: - packet = np.asarray(msg.samples, dtype=np.float32) - rows = int(msg.header.num_samples) - expected = rows * self.num_channels - if not packet.size and msg.raw_frames: - # Timeseries packets may carry samples in the packed - # `raw_frames` blob instead of the repeated `samples` field. - # That encoding is undocumented here, so fail with a specific - # message rather than a confusing payload-size mismatch. - msg_ = ( - f"Sintela protobuf {tag} packets store samples in " - "raw_frames, which DASCore cannot yet decode." + index += self._fill_packet(data, index, tag, msg) + return data, self.coords, self.attrs + + def decode_stream(self, resource, meta: ParsedMeta): + """ + Fill the patch array packet by packet straight from the stream. + + The endpoint shortcut already established the shape, so samples go + straight to their final home and each decoded packet is released at + once; holding every packet *and* the output array, as the generic path + must, costs roughly twice the patch. + + Headers are copied into fresh messages before the packet is dropped. + Clearing the samples in place would not do: protobuf releases an arena + only when the whole message dies, so a cleared packet still owns the + space its samples occupied. The copies go to ``from_parsed`` at the + end, so full cross-packet validation still runs and supplies the coords + and attrs returned. + """ + messages = _get_proto_messages(include_sample_fields=True) + header_messages = _get_proto_messages(include_sample_fields=False) + decode_error = _get_protobuf_decode_error() + data = np.empty(self.shape, dtype=self.dtype) + parsed: list[tuple[str, Any]] = [] + index = 0 + for record in _iter_envelope_records(resource, strict=True): + if record.tag == META_TAG: + # Every record is visited here, so META is picked up wherever + # it sits, matching the read-everything path. The endpoint + # scan only sees the ones before the first data packet. + meta = _parse_meta(record.payload) + continue + if record.tag not in _TAG_TO_PACKET: + continue + if record.tag not in TS_TAGS: + raise InvalidFiberFileError( + "Mixed Sintela protobuf packet families are unsupported." ) - raise InvalidFiberFileError(msg_) - if packet.size != expected: + msg = _parse_packet(record.tag, record.payload, messages, decode_error) + if index + int(msg.header.num_samples) > self.shape[0]: + # More samples than the endpoints implied: the packets in + # between are not the contiguous run this path assumed. raise InvalidFiberFileError( - "Unexpected Sintela protobuf TS sample payload size." + "Non-contiguous Sintela protobuf sample counts." ) - data[index : index + rows] = packet.reshape(rows, self.num_channels) - index += rows - return data, self.coords, self.attrs + index += self._fill_packet(data, index, record.tag, msg) + light = header_messages[_TAG_TO_PACKET[record.tag]]() + # Round-tripped rather than copied: the sample-bearing and + # header-only classes come from separate descriptor pools, so + # CopyFrom rejects them as different types. A header is ~100 bytes. + light.header.ParseFromString(msg.header.SerializeToString()) + parsed.append((record.tag, light)) + del msg + metadata = type(self).from_parsed(parsed, meta) + if index != self.shape[0]: + # Contiguity, validated just above, should force the packet lengths + # to sum to the endpoint-derived total. Checked unconditionally + # anyway: the alternative to raising is handing back the + # uninitialized tail of the output array. + raise InvalidFiberFileError( + "Sintela protobuf packets do not fill the expected sample count." + ) + return data, metadata.coords, metadata.attrs class BandMetadata(_PacketMetadata): @@ -1042,8 +1323,153 @@ def decode(self, parsed: list[tuple[str, Any]]): } +def _parse_packet_header(tag: str, payload_prefix: bytes): + """Parse just the header submessage of a data packet from a payload prefix.""" + header_bytes = _leading_header_bytes(payload_prefix) + if header_bytes is None: + return None + messages = _get_proto_messages(include_sample_fields=False) + try: + return _parse_packet(tag, header_bytes, messages, _get_protobuf_decode_error()) + except InvalidFiberFileError: + # Leave the diagnostic to the read-everything path. + return None + + +def _fits_in_file(total_samples: int, num_channels: int, file_size: int) -> bool: + """ + Return whether a file this size could hold that many samples. + + Guards the *over*-estimate direction only: a float32 sample costs at least + four bytes on the wire, so a length implying more bytes than exist cannot + be real -- the shape of a counter that ran backwards, which the modular + difference turns into billions of samples. Under-estimates pass trivially; + ``_endpoint_time_agrees`` catches those. + """ + if total_samples <= 0 or num_channels <= 0: + return False + return total_samples * num_channels * 4 <= file_size + + +def _endpoint_time_agrees(first_header, last_header, total_samples: int) -> bool: + """ + Check the endpoint timestamps against the span their sample counts imply. + + The counters cannot tell one contiguous recording from two concatenated + ones, or from a counter that reset partway; the timestamps are an + independent witness, since for a contiguous run the elapsed time between + the stamps should equal the samples between them over the sample rate. + + Coarse by design: recorder clocks resynchronize as a recording runs, so a + valid file's stamps drift by milliseconds and the tolerance must absorb + that, which leaves room for a small mid-file gap to hide. It catches the + gross disagreement -- seconds to hours -- left by concatenation, a reset + counter, a reconfiguration, or a misidentified final record. + + Returns False when either stamp is missing: an unverifiable shortcut is not + worth taking when the full read is always available. + """ + first_time = _common_header_time(first_header.common_header) + last_time = _common_header_time(last_header.common_header) + if first_time is None or last_time is None: + return False + rate = float(first_header.common_header.sample_rate) + if not np.isfinite(rate) or rate <= 0: + return False + # Samples strictly between the two stamps, i.e. excluding the last packet. + leading = total_samples - int(last_header.num_samples) + expected_ns = round(leading / rate * 1e9) + elapsed_ns = int((last_time - first_time).astype("timedelta64[ns]").astype(int)) + packet_ns = round(int(last_header.num_samples) / rate * 1e9) + tolerance = max( + _ENDPOINT_TIME_PACKETS * packet_ns, + expected_ns * _ENDPOINT_TIME_TOLERANCE, + _ENDPOINT_TIME_FLOOR_NS, + ) + return abs(elapsed_ns - expected_ns) <= tolerance + + +def _get_endpoint_metadata(resource): + """ + Summarize a timeseries recording from its first and last packets alone. + + A summary needs only header fields, and ``sample_count`` numbers the + samples preceding each packet, so the last packet's count plus its length + is the total. Scanning a half-gigabyte recording becomes two small reads. + + This assumes the packets in between are one contiguous, homogeneous run. + Three checks test that assumption: each endpoint's declared length must fit + its own record, the total must fit the file, and the endpoint timestamps + must match the span those samples imply. A concatenated file, a reset + counter, or a bogus tail match fails one of them. + + Returns ``(metadata, meta)``, or None when any check fails or the layout + does not suit the shortcut, leaving the caller to read the whole file. + Only META preceding the first data packet is seen here; ``read`` picks up + any that appear later. + """ + file_size = _stream_size(resource) + head = _find_first_data_record(resource, file_size) + if head is None: + return None + meta, (tag, size, prefix) = head + # Only the timeseries family has an evenly sampled time coord derivable + # from the endpoints; BAND and FFT time coords list every packet's stamp. + if tag not in TS_TAGS: + return None + # Capped as well as floored: the first packet's declared size comes off + # disk, and an implausible one would otherwise make the window swallow the + # whole file in a single buffer. + tail = _find_last_record( + resource, file_size, min(size + _TAIL_SEARCH_SLACK, _MAX_TAIL_SEARCH) + ) + if tail is None: + # Sizing the window from the first packet assumes the last one is + # about as big. A recording that opens with a short warm-up packet + # breaks that, so widen once rather than give up and read everything + # -- the retry costs one read, and only for files that need it. + tail = _find_last_record(resource, file_size, min(file_size, _MAX_TAIL_SEARCH)) + if tail is None or tail[0] != tag: + return None + first = _parse_packet_header(tag, prefix) + last = _parse_packet_header(tag, tail[2]) + if first is None or last is None: + return None + channels = int(first.header.common_header.num_channels) + # Each endpoint's declared length has to fit the record carrying it. The + # time check below cannot see the last packet's own length -- it cancels + # out of the elapsed-time comparison -- so without this a packet could + # claim millions of samples it has no room for and still be believed. + if not ( + _fits_in_file(int(first.header.num_samples), channels, size) + and _fits_in_file(int(last.header.num_samples), channels, tail[1]) + ): + return None + total_samples = ( + int(last.header.sample_count) - int(first.header.sample_count) + ) % _SAMPLE_COUNT_MODULUS + int(last.header.num_samples) + if not _fits_in_file(total_samples, channels, file_size): + return None + if not _endpoint_time_agrees(first.header, last.header, total_samples): + return None + try: + metadata = TimeseriesMetadata.from_parsed( + [(tag, first), (tag, last)], meta, total_samples=total_samples + ) + except InvalidFiberFileError: + # Endpoints that fail validation may be a genuinely bad file or a + # false tail match; either way the full path decides, reporting the + # same error for the former rather than trusting a doubtful header. + return None + return metadata, meta + + def read_payload(resource): """Decode a Sintela protobuf file into data, coords, and attrs.""" + endpoints = _get_endpoint_metadata(resource) + if endpoints is not None: + metadata, meta = endpoints + return metadata.decode_stream(resource, meta) records = _iter_envelope_records(resource, strict=True) parsed, meta = _parse_records(records, scan_mode=False) return _decode_family(parsed, meta) @@ -1051,8 +1477,13 @@ def read_payload(resource): def scan_payload(resource) -> list[ScanPayload]: """Decode a Sintela protobuf file and return FiberIO scan payloads.""" - records = _iter_envelope_records(resource, strict=True) - parsed, meta = _parse_records(records, scan_mode=True) - family_cls = _FAMILY_CLASSES[_validate_single_family(parsed)] - shape, coords, attrs, dtype = family_cls.from_parsed(parsed, meta).scan() + endpoints = _get_endpoint_metadata(resource) + if endpoints is not None: + metadata = endpoints[0] + else: + records = _iter_envelope_records(resource, strict=True, headers_only=True) + parsed, meta = _parse_records(records, scan_mode=True) + family_cls = _FAMILY_CLASSES[_validate_single_family(parsed)] + metadata = family_cls.from_parsed(parsed, meta) + shape, coords, attrs, dtype = metadata.scan() return [make_scan_payload(attrs=attrs, coords=coords, shape=shape, dtype=dtype)] diff --git a/docs/changelog.qmd b/docs/changelog.qmd index 8cdde9dd..9c657766 100644 --- a/docs/changelog.qmd +++ b/docs/changelog.qmd @@ -5,6 +5,7 @@ The [releases page](https://github.com/DASDAE/dascore/releases) tracks changes f ## Unreleased API Changes - Missing optional dependency messages now name the package to install and give the command to install it. The import name and the package name are not always the same, so `dc.scan` reporting `{'google.protobuf.descriptor_pb2': 10308}` left users guessing; it now reports `protobuf (10308 files)` along with ``Install with `pip install protobuf` or `uv pip install protobuf` ``. `MissingOptionalDependencyError` also carries an `install_name` attribute, and an `ImportError` raised inside an installed package is reported as a failed import rather than a missing install. +- **Scanning a Sintela *timeseries* protobuf recording no longer reads the whole file.** The scan read every packet's samples to reach its headers, so indexing a directory cost a full sequential read of every file in it. A timeseries scan now derives its summary from the first and last packets alone: `sample_count` numbers the samples preceding each packet, so the endpoints give the total directly. Measured over a directory of 10,308 recordings of 505 MB each, `dc.scan` went from 2.05 s to 0.02 s per file (5.9 hours to 3.3 minutes), reading 1.3 MB per file instead of 505 MB. The shortcut is taken only when two independent checks agree that the file is one contiguous run — the derived length must fit in the file's bytes, and the endpoint timestamps must match the span those samples imply — so a concatenated file, a reset sample counter, or a spurious match falls back to reading everything. The one behavioural consequence: a scan no longer detects a gap between the endpoints when the timestamps corroborate it (a paused and resumed acquisition), and reports the span they imply; `read` validates every packet and still raises `InvalidFiberFileError`. BAND and FFT scans need every packet's timestamp and so still visit every record; they are now handed each packet's header without its samples, which bounds their memory but not their I/O. Reading a timeseries recording also streams packets into the output array instead of decoding all of them first, cutting peak memory from ~2.0x the patch to ~1.02x. - The fiber IO format readers now share two helpers instead of each carrying its own copy of the same scaffolding: [`make_scan_payload`](`dascore.io.make_scan_payload`) builds one `FiberIO.scan` payload (taking `dims` and `shape` from the coords unless given), and `dascore.io.utils.build_patches` performs the common `read` tail of trim, drop-if-empty, attach attrs. Two side effects for readers: an already empty source now yields no patch from the `APSensing` and `HDAS` readers rather than a zero-size one (matching every other format), and the `GDR_DAS` and `Neubrex` readers declare `time`/`distance` explicitly rather than absorbing them from `**kwargs`. The unused, never-populated `ProdMLPatchAttrs` classes are removed from `dascore.io.prodml.core` and `dascore.io.dashdf5.core`; `ProdMLRawPatchAttrs` in `dascore.io.prodml.utils` is the one the reader uses. Relatedly, evenly sampled coordinates are now built from a sample count rather than a hand-computed stop, so a source declaring zero samples yields an empty coordinate instead of raising a validation error. - Fixed chunking a coordinate whose units are not the canonical SI unit (e.g. a distance in feet). Plan trims carry canonical SI magnitudes, which were applied to the patch coordinate as bare numbers, so each piece covered the wrong physical interval and samples were silently dropped — chunking a 300 channel patch in feet returned 61 channels. The trim now converts to the coordinate's own units at load, the same conversion `Spool.select` already performed. - **`Spool.chunk` accepts quantities as the chunk length.** A quantity in the coordinate's own units works (`chunk(time=10 * dc.units.s)`, `chunk(distance=100 * dc.units.ft)`), where any quantity previously raised `NotImplementedError`. A quantity of information (`chunk(time=25 * dc.units.megabytes)`) chunks so each patch's *data array* is at most the requested size; the sample count is floored, and a partition mixing element types is sized against the dtype assembly upcasts to. `overlap` accepts both forms, and `Spool.chunk_plan(...).params["size"]` reports what a size resolved to. The index schema records each patch's element dtype to make this possible, so its version is bumped and existing indexes must be deleted and rebuilt. An attr named `dtype` is now reserved and stays unindexed. diff --git a/tests/test_io/test_sintela/test_protobuf.py b/tests/test_io/test_sintela/test_protobuf.py index ad3b8208..dad39123 100644 --- a/tests/test_io/test_sintela/test_protobuf.py +++ b/tests/test_io/test_sintela/test_protobuf.py @@ -4,6 +4,7 @@ from __future__ import annotations +import gc import struct import warnings from functools import cache @@ -106,15 +107,27 @@ def _get_num_samples(payload: bytes) -> int: return int(msg.header.num_samples) -def _build_ts_payloads(): - """Create two contiguous timeseries packets.""" +def _build_ts_payloads(n_packets: int = 2): + """ + Create a run of contiguous timeseries packets. + + Packet stamps advance by the packet's own duration (3 samples at 2 Hz = + 1.5 s), so the timestamps and the sample counters tell the same story -- + which is what a real recording does, and what the scan shortcut checks. + """ packet_cls = _get_test_proto_messages()["TimeseriesPacket"] packets = [] - for offset, sample_count in enumerate((0, 3)): + for offset in range(n_packets): + sample_count = offset * 3 + elapsed_ns = offset * 3 * 1_000_000_000 // 2 msg = packet_cls() hdr = msg.header common = hdr.common_header - _set_timestamp(common.time, 1_700_000_000 + offset) + _set_timestamp( + common.time, + 1_700_000_000 + elapsed_ns // 1_000_000_000, + elapsed_ns % 1_000_000_000, + ) common.num_channels = 2 common.sample_rate = 2.0 common.channel_spacing = 10.0 @@ -163,6 +176,11 @@ def _build_band_payloads(): return packets +def _build_real_fft_payloads(): + """Create two real-valued FFT packets.""" + return _build_fft_payloads(complex_data=False) + + def _build_fft_payloads(*, complex_data: bool): """Create two FFT packets.""" packet_cls = _get_test_proto_messages()["FFTPacket"] @@ -216,6 +234,74 @@ def _mutate_all(records, message_type: str, mutator): return out +class _BytesReader: + """A minimal seekable binary reader over a bytes object.""" + + def __init__(self, data): + self._data = data + self._pos = 0 + + def seek(self, pos, whence=0): + """Seek from the start, or from the end when whence is 2.""" + self._pos = len(self._data) + pos if whence == 2 else pos + + def tell(self): + """Return the current position.""" + return self._pos + + def read(self, size=-1): + """Read up to size bytes, or the remainder when size is negative.""" + end = len(self._data) if size < 0 else self._pos + size + out = self._data[self._pos : end] + self._pos = min(end, len(self._data)) + return out + + +class _CountingReader: + """ + A binary handle that records how many bytes were actually read. + + Wrapping the handle (rather than sampling process memory) makes the cost + of a scan directly observable: protobuf holds sample data in C++ memory + that Python's allocation tracing never sees, but every byte still has to + come through a ``read`` call first. + """ + + def __init__(self, handle): + self._handle = handle + self.bytes_read = 0 + + def read(self, size=-1): + """Read from the wrapped handle, accumulating the byte count.""" + out = self._handle.read(size) + self.bytes_read += len(out) + return out + + def __getattr__(self, name): + """Delegate seek/tell and friends to the wrapped handle.""" + return getattr(self._handle, name) + + +def _bytes_read_by(func, path) -> int: + """Return how many bytes ``func`` pulls off disk for ``path``.""" + with path.open("rb") as handle: + reader = _CountingReader(handle) + func(reader) + return reader.bytes_read + + +def del_samples_beyond(msg, keep: int): + """Trim a packet's sample array down to ``keep`` values.""" + del msg.samples[keep:] + + +def _pad_samples(records, message_type: str, pad: int): + """Return records whose packets each carry ``pad`` extra samples.""" + return _mutate_all( + records, message_type, lambda msg: msg.samples.extend([0.0] * pad) + ) + + @pytest.fixture() def ts_records(): """Return baseline timeseries records.""" @@ -363,28 +449,6 @@ def test_complex_fft_not_labeled_power_spectral_density( # real packets keep the power-spectral-density label assert fiber_io.read(real_path)[0].attrs.data_type == "power_spectral_density" - def test_scan_memory_independent_of_sample_count(self, ts_records): - """ - Metadata-only scans must not hold the sample bytes in memory. - - Protobuf keeps undeclared wire fields as unknown fields, so omitting - the sample declarations is not by itself enough. The invariant is - that scan-retained size does not grow with the sample payload. - """ - - def _retained(pad_samples: int) -> int: - packet_cls = _get_test_proto_messages()["TimeseriesPacket"] - records = [] - for tag, payload in ts_records: - msg = packet_cls() - msg.ParseFromString(payload) - msg.samples.extend([0.0] * pad_samples) - records.append(_envelope((tag, msg.SerializeToString()))) - parsed, _ = sintela_utils._parse_records(records, scan_mode=True) - return sum(msg.ByteSize() for _tag, msg in parsed) - - assert _retained(0) == _retained(5_000) - def test_raw_frame_only_tags_are_not_detected( self, fiber_io, write_sintela_file, ts_records ): @@ -493,19 +557,33 @@ def test_mixed_families_raise( with pytest.raises(InvalidFiberFileError, match="Mixed Sintela protobuf"): fiber_io.scan(path) - def test_non_contiguous_timeseries_raises( + def test_non_contiguous_timeseries_caught_on_read_not_scan( self, fiber_io, write_sintela_file, ts_records ): - """Timeseries packets with gaps or reordering should fail.""" + """ + A self-consistent gap is reported by read, not by scan. + + A paused and resumed acquisition looks, from the endpoints alone, like + a longer continuous one: the timestamps corroborate the counters. The + summary spans the gap and read raises. A gap whose timestamps do *not* + corroborate is caught earlier; see the concatenation test. + """ records = _mutate_record( ts_records, 1, "TimeseriesPacket", - lambda msg: setattr(msg.header, "sample_count", 10), + # 10 samples precede this packet rather than 3, and its stamp moves + # out to match: 10 samples at 2 Hz is 5 s after the first. + lambda msg: ( + setattr(msg.header, "sample_count", 10), + _set_timestamp(msg.header.common_header.time, 1_700_000_005), + ), ) path = write_sintela_file("non_contiguous.pb", records) + # The last packet reports 10 samples before it and adds 3 of its own. + assert _payload_to_summary(fiber_io.scan(path)[0]).shape[0] == 13 with pytest.raises(InvalidFiberFileError, match="Non-contiguous"): - fiber_io.scan(path) + fiber_io.read(path) def test_bad_magic_returns_false(self, fiber_io, tmp_path): """Invalid magic bytes should not identify as the format.""" @@ -607,6 +685,9 @@ def test_truncated_payload_raises(self, fiber_io, tmp_path): assert not fiber_io.get_format(path) with pytest.raises(InvalidFiberFileError, match="Truncated"): fiber_io.scan(path) + # read takes the full-payload path, which detects the short tail too + with pytest.raises(InvalidFiberFileError, match="Truncated"): + fiber_io.read(path) def test_bad_protobuf_payloads_raise_invalid_fiber_file_error( self, fiber_io, write_sintela_file @@ -627,6 +708,119 @@ def test_bad_protobuf_payloads_raise_invalid_fiber_file_error( ): fiber_io.scan(data_path) + def test_read_rejects_family_change_between_endpoints( + self, fiber_io, write_sintela_file, ts_records, band_records + ): + """ + A foreign packet between matching endpoints is caught on read. + + Both endpoints are timeseries, so the scan shortcut accepts the file; + the streaming read visits every packet and rejects the BAND one. + """ + records = [ts_records[0], band_records[0], ts_records[1]] + path = write_sintela_file("ts_band_ts.pb", records) + with pytest.raises(InvalidFiberFileError, match="Mixed Sintela protobuf"): + fiber_io.read(path) + + def test_read_rejects_more_samples_than_endpoints_imply( + self, fiber_io, write_sintela_file + ): + """ + Extra samples between the endpoints are caught before they overrun. + + The output array is sized from the endpoints, so a middle packet that + pushes past that total has to fail rather than write out of bounds. + """ + records = _build_ts_payloads(3) + # Renumber the last packet so the endpoints imply two packets' worth + # of samples while three packets are actually present. + records = _mutate_record( + records, + 2, + "TimeseriesPacket", + lambda msg: setattr(msg.header, "sample_count", 3), + ) + path = write_sintela_file("ts_overrun.pb", records) + with pytest.raises(InvalidFiberFileError, match="Non-contiguous"): + fiber_io.read(path) + + def test_short_final_packet_is_measured_from_the_last_endpoint( + self, fiber_io, write_sintela_file + ): + """ + A recording ending in a partial packet reports its true length. + + Recorders close a file whenever the acquisition stops, so the final + packet is routinely shorter than the rest. The endpoint total has to + add *that* packet's length to the counter difference; adding the first + packet's length instead happens to give the same answer for a uniform + file, which is what every other fixture here is. + """ + records = _build_ts_payloads(3) + # 3 + 3 + 1 samples: counts stay 0, 3, 6 and the tail is short. + records = _mutate_record( + records, + 2, + "TimeseriesPacket", + lambda msg: ( + setattr(msg.header, "num_samples", 1), + del_samples_beyond(msg, msg.header.common_header.num_channels), + ), + ) + path = write_sintela_file("ts_short_tail.pb", records) + # Asserted on the shortcut itself, not just the final answer: adding + # the wrong packet's length makes the endpoint total disagree with the + # timestamps, and the fallback then returns the right shape anyway. + with path.open("rb") as handle: + endpoints = sintela_utils._get_endpoint_metadata(handle) + assert endpoints is not None + assert endpoints[0].shape[0] == 7 + assert _payload_to_summary(fiber_io.scan(path)[0]).shape[0] == 7 + assert fiber_io.read(path)[0].shape[0] == 7 + + def test_reset_sample_counter_falls_back_instead_of_trusting_endpoints( + self, fiber_io, write_sintela_file + ): + """ + A counter that restarts mid-file must not be read as billions of samples. + + `sample_count` is a uint32, so the endpoint difference is taken modulo + 2**32 to survive a wrap. A counter that *resets* instead looks like an + enormous backwards jump, which would size the patch from a number the + file cannot possibly hold. The shortcut declines and the full path + reports the gap. + """ + records = _build_ts_payloads(3) + # Start high so the counter appears to restart: the last packet's + # count is far below the first, and the modular difference is enormous. + records = _mutate_record( + records, + 0, + "TimeseriesPacket", + lambda msg: setattr(msg.header, "sample_count", 1000), + ) + path = write_sintela_file("ts_counter_reset.pb", records) + with pytest.raises(InvalidFiberFileError, match="Non-contiguous"): + fiber_io.scan(path) + with pytest.raises(InvalidFiberFileError, match="Non-contiguous"): + fiber_io.read(path) + + def test_read_picks_up_meta_after_the_first_data_packet( + self, fiber_io, write_sintela_file, ts_records + ): + """ + Read honours a META record wherever it appears, as it always has. + + The endpoint scan only walks as far as the first data packet, so a + trailing META is invisible to it; the streaming read visits every + record and must not lose metadata the old full-decode path collected. + """ + records = [ts_records[0], ("META", _build_meta_payload()), ts_records[1]] + path = write_sintela_file("ts_trailing_meta.pb", records) + patch = fiber_io.read(path)[0] + assert patch.attrs.instrument_manufacturer == "Sintela" + assert patch.attrs.serial_number == "SN123" + def test_timeseries_read_rejects_dropped_samples( self, fiber_io, write_sintela_file, ts_records ): @@ -949,6 +1143,562 @@ def _raise(*args, **kwargs): fiber_io.read(path) +_FAMILY_BUILDERS = [ + (_build_ts_payloads, "TimeseriesPacket"), + (_build_band_payloads, "BandPacket"), + (_build_real_fft_payloads, "FFTPacket"), +] + + +class TestSintelaProtobufScanCost: + """Scanning must not pay for the sample data it does not report.""" + + def test_scan_reads_do_not_grow_with_recording_length(self, write_sintela_file): + """ + A longer timeseries recording must not cost a longer scan. + + This is what makes indexing a directory of large recordings tractable: + the summary is derived from the first and last packets, so a file with + sixteen times as many packets is read no more than a short one. The + packets are padded past the header prefix so the reads being compared + are the bounded ones, not an artifact of tiny fixtures. + """ + + def _write(name, n_packets): + records = _pad_samples( + _build_ts_payloads(n_packets), "TimeseriesPacket", 20_000 + ) + return write_sintela_file(name, records) + + short = _write("ts_short.pb", 4) + long = _write("ts_long.pb", 64) + assert long.stat().st_size > 8 * short.stat().st_size + + short_bytes = _bytes_read_by(sintela_utils.scan_payload, short) + long_bytes = _bytes_read_by(sintela_utils.scan_payload, long) + assert short_bytes == long_bytes + assert long_bytes < long.stat().st_size / 4 + + def test_read_and_scan_take_the_endpoint_shortcut( + self, monkeypatch, write_sintela_file, ts_records, sintela_protobuf_path + ): + """ + Both entry points reach the fast path, for synthetic and real files. + + Every behavioural test here passes whether or not the shortcut engages, + because the full path reports the same errors and returns the same + data. Only the cost differs, so the call sites are checked directly: + without this, `read` could quietly revert to decoding every packet up + front and nothing would notice. + """ + path = write_sintela_file("ts_shortcut.pb", ts_records) + for target in (path, Path(sintela_protobuf_path)): + with target.open("rb") as handle: + assert sintela_utils._get_endpoint_metadata(handle) is not None + + used = [] + original = sintela_utils.TimeseriesMetadata.decode_stream + + def _spy(self, resource, meta): + used.append(True) + return original(self, resource, meta) + + monkeypatch.setattr( + sintela_utils.TimeseriesMetadata, "decode_stream", _spy, raising=True + ) + with path.open("rb") as handle: + sintela_utils.read_payload(handle) + assert used, "read_payload did not take the streaming endpoint path" + + def test_sample_counter_rollover_keeps_the_shortcut( + self, write_sintela_file, ts_records + ): + """ + A counter that wraps mid-file is still summarized from its endpoints. + + Dropping the modulus would make the difference negative, which the + file-size bound rejects: the answer stays right but arrives by the slow + route, so only a cost assertion catches it. + """ + packet_cls = _get_test_proto_messages()["TimeseriesPacket"] + modulus = sintela_utils._SAMPLE_COUNT_MODULUS + records = [] + for index, (tag, payload) in enumerate(ts_records): + msg = packet_cls() + msg.ParseFromString(payload) + msg.header.sample_count = ( + modulus - msg.header.num_samples if index == 0 else 0 + ) + records.append((tag, msg.SerializeToString())) + path = write_sintela_file("ts_rollover_shortcut.pb", records) + with path.open("rb") as handle: + endpoints = sintela_utils._get_endpoint_metadata(handle) + assert endpoints is not None + assert endpoints[0].shape[0] == 6 + + def test_find_last_record_ignores_a_decoy_frame_inside_the_payload(self): + """ + Only the frame whose extent ends on EOF is the final record. + + Sample bytes can spell out a magic word and a tag; anchoring on EOF is + what distinguishes the real framing from that noise, and the backwards + search reaches the decoy first. + """ + magic = struct.pack("= 6, "the streaming path did not decode packet by packet" + assert max(live) == baseline, ( + "a previously decoded packet was still alive when the next was " + f"decoded (live counts: {live}, baseline: {baseline})" + ) + held = [cls for group in kept for cls in group] + assert held, "nothing was handed to from_parsed" + assert not any( + any(field.name == "samples" for field in cls.DESCRIPTOR.fields) + for cls in held + ), "the read retains messages that are able to hold samples" + + @pytest.mark.parametrize(("builder", "message_type"), _FAMILY_BUILDERS) + def test_header_only_records_exclude_samples( + self, write_sintela_file, builder, message_type + ): + """ + Every family's scan path is handed headers, never sample payloads. + + The families that cannot use the endpoint shortcut still have to visit + each packet; they must do it without pulling the samples along, so the + payload handed to the parser stays the same size as the recording's + samples grow. + """ + + def _payload_bytes(pad: int) -> int: + records = _pad_samples(builder(), message_type, pad) + path = write_sintela_file(f"headers_{message_type}_{pad}.pb", records) + with path.open("rb") as handle: + return sum( + len(record.payload) + for record in sintela_utils._iter_envelope_records( + handle, strict=True, headers_only=True + ) + ) + + assert _payload_bytes(0) == _payload_bytes(50_000) + + def test_large_payloads_are_skipped_not_streamed(self, write_sintela_file): + """ + A payload large enough to be worth a seek is stepped over, not read. + + Small remainders are streamed and discarded because a seek costs a + platter rotation on spinning media, which buys back more than the + transfer it saves. Past the threshold the seek wins, and this pins that + the large-payload branch really does skip the bytes. + """ + # Fixed size, not derived from the threshold: a fixture that scales + # with the constant would keep passing however the constant is retuned. + n_floats = 1_000_000 + records = _pad_samples(_build_band_payloads(), "BandPacket", n_floats) + path = write_sintela_file("band_large.pb", records) + size = path.stat().st_size + # Guard on the payload bytes the fixture actually produces, so this + # cannot silently turn into a permanent skip. + payload_bytes = size / len(records) + if payload_bytes <= sintela_utils._SEEK_SKIP_THRESHOLD: + pytest.skip("seek threshold retuned above this fixture's payload size") + + def _walk(handle): + for _record in sintela_utils._iter_envelope_records( + handle, strict=True, headers_only=True + ): + pass + + assert _bytes_read_by(_walk, path) < size / 4 + + @pytest.mark.parametrize(("builder", "message_type"), _FAMILY_BUILDERS) + def test_scan_memory_independent_of_sample_count(self, builder, message_type): + """ + Metadata-only scans must not retain the sample bytes. + + Protobuf keeps undeclared wire fields as unknown fields, so omitting + the sample declarations is not by itself enough. The invariant is that + scan-retained size does not grow with the sample payload. + """ + + def _retained(pad: int) -> int: + records = [ + _envelope(record) + for record in _pad_samples(builder(), message_type, pad) + ] + parsed, _ = sintela_utils._parse_records(records, scan_mode=True) + return sum(msg.ByteSize() for _tag, msg in parsed) + + assert _retained(0) == _retained(5_000) + + +class TestSintelaProtobufWireHelpers: + """ + Tests for the hand-rolled wire-format helpers. + + These read protobuf framing directly so a packet's header can be taken + from a short prefix. Each returns None rather than raising when a payload + does not fit that shape, which sends the caller to the read-everything + path; the fallbacks are exercised here since a well-formed file never + reaches them. + """ + + def test_read_varint_spans_multiple_bytes(self): + """Varints longer than one byte decode, and a truncated one is None.""" + # 300 encodes as two bytes; the continuation bit drives the loop. + assert sintela_utils._read_varint(b"\xac\x02", 0) == (300, 2) + assert sintela_utils._read_varint(b"\x05", 0) == (5, 1) + # A varint whose continuation bit never terminates is unusable. + assert sintela_utils._read_varint(b"\xac", 0) == (None, 1) + + def test_leading_header_bytes_rejects_unusable_prefixes(self): + """Only a complete field-1 submessage at the front is usable.""" + header = b"\x0a\x03abc" + assert sintela_utils._leading_header_bytes(header) == header + # Empty, or not starting with field 1 (here field 2). + assert sintela_utils._leading_header_bytes(b"") is None + assert sintela_utils._leading_header_bytes(b"\x12\x03abc") is None + # Field 1 declared longer than the prefix actually holds. + assert sintela_utils._leading_header_bytes(b"\x0a\x7fab") is None + + def test_parse_frame_rejects_short_and_wrong_magic(self): + """Framing needs twelve bytes opening with the magic word.""" + frame = ( + struct.pack("