From 9ed9f23cc4ba7135e6298991db8e2345b7c31422 Mon Sep 17 00:00:00 2001 From: slazyverse <188746735+slazyverse@users.noreply.github.com> Date: Sat, 12 Sep 2026 03:51:31 +0530 Subject: [PATCH] =?UTF-8?q?feat(m3):=20M3/E5/#18=20=E2=80=94=20HEL1OS=20pa?= =?UTF-8?q?rsers?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Implements `SPEC-parsers` §2.5–§2.9 and §4 for HEL1OS: band light curves, PHA spectra, housekeeping, per-detector GTI, an event reader, and orbit identity with version precedence. No merge API, no coverage map, no aggregation and no dual write — those belong to the write path (#19). GOVERNANCE FIRST. Verifying against all 391 orbits proved three contract checks undecidable at float64 resolution: R-1 H3 rejected 647 of 1,564 spectra products and §2.8's header-span check rejected 94 of 391 housekeeping products, every one of them by one or two representable steps. `SPEC-parsers` r7 therefore defines `col_span`, adds the §5.1 time-representation allowance ε_t = 1 ms — the value CONTRADICTION-006 Defect A approved in code and never recorded in the contract — and gives the §2.5 and §2.8 time checks the F-06 id they lacked. Recorded as CONTRA-007, CLOSED by r7. §2.5's `mjd` non-decreasing rule is falsified by every event HDU in the archive (1,564 of 1,564; 37.9 M backward steps; largest 1,156 s). It is NOT amended: the parser enforces it and fails closed, and the falsification is recorded OPEN as CONTRA-008 for the maintainer to rule on. D6: the missing §2.5 row checks (span, V-EVT-2, `ener>0`), the declared `ener` and `CTR` units, CZT `pix`/`offsetchn`, and every §2.8 decisive column with its declared unit and archive dtype. Only `czt1temp` and `czt2temp` are required finite, as §2.8 states. D3: `_fits.py` moves from `parsers/solexs/` to `parsers/` so both instruments share the no-default FITS accessors without either importing the other, and gains `declared_unit()`. SoLEXS behaviour is unchanged and regression-tested. Verified with the real 391-orbit corpus: links, contracts, imports, architecture 563, unit 583, integration 146. Clean export from the staged tree: unit 581 (2 skipped), integration 100 (22 skipped), every real-data test skipping rather than failing. Archive-wide, 389 of 391 orbits now parse end to end; the 2 that do not are the known duplicate-HK archive defects, still terminating at F-16. Co-Authored-By: Claude Opus 5 --- contexts/ingest/parsers/__init__.py | 4 +- contexts/ingest/parsers/{solexs => }/_fits.py | 32 +- contexts/ingest/parsers/hel1os/__init__.py | 27 + contexts/ingest/parsers/hel1os/events.py | 317 ++++++ contexts/ingest/parsers/hel1os/gti.py | 142 +++ contexts/ingest/parsers/hel1os/hk.py | 308 ++++++ contexts/ingest/parsers/hel1os/lc.py | 328 +++++++ contexts/ingest/parsers/hel1os/orbit.py | 233 +++++ contexts/ingest/parsers/hel1os/spectra.py | 326 +++++++ contexts/ingest/parsers/solexs/gti.py | 2 +- contexts/ingest/parsers/solexs/lc.py | 2 +- contexts/ingest/parsers/solexs/pi.py | 2 +- contexts/ingest/tests/hel1os_fixtures.py | 264 +++++ contexts/ingest/tests/test_hel1os_parsers.py | 903 ++++++++++++++++++ specs/contradictions/CONTRA-007.md | 191 ++++ specs/contradictions/CONTRA-008.md | 103 ++ specs/contradictions/index.md | 2 + specs/parsers/SPEC-parsers.md | 40 +- specs/parsers/index.md | 4 +- .../test_hel1os_parse_with_platform.py | 506 ++++++++++ 20 files changed, 3721 insertions(+), 15 deletions(-) rename contexts/ingest/parsers/{solexs => }/_fits.py (69%) create mode 100644 contexts/ingest/parsers/hel1os/__init__.py create mode 100644 contexts/ingest/parsers/hel1os/events.py create mode 100644 contexts/ingest/parsers/hel1os/gti.py create mode 100644 contexts/ingest/parsers/hel1os/hk.py create mode 100644 contexts/ingest/parsers/hel1os/lc.py create mode 100644 contexts/ingest/parsers/hel1os/orbit.py create mode 100644 contexts/ingest/parsers/hel1os/spectra.py create mode 100644 contexts/ingest/tests/hel1os_fixtures.py create mode 100644 contexts/ingest/tests/test_hel1os_parsers.py create mode 100644 specs/contradictions/CONTRA-007.md create mode 100644 specs/contradictions/CONTRA-008.md create mode 100644 tests/integration/test_hel1os_parse_with_platform.py diff --git a/contexts/ingest/parsers/__init__.py b/contexts/ingest/parsers/__init__.py index e0904af..945f0ed 100644 --- a/contexts/ingest/parsers/__init__.py +++ b/contexts/ingest/parsers/__init__.py @@ -1,6 +1,8 @@ """Instrument parsers — canonicalise an acquired product into domain objects. -`parsers/solexs` is delivered by M3/E5/#17. `parsers/hel1os` is #18 and does not exist yet. +`parsers/solexs` is delivered by M3/E5/#17 and `parsers/hel1os` by M3/E5/#18. Both use the +strict FITS accessors in `_fits`, which lives here rather than inside either instrument's +package so that neither instrument depends on the other's internals. Acquisition and canonicalisation are one context with two internal module groups (ADR-0026), with the seam preserved as directory boundaries because it is free under ADR-0025. There is diff --git a/contexts/ingest/parsers/solexs/_fits.py b/contexts/ingest/parsers/_fits.py similarity index 69% rename from contexts/ingest/parsers/solexs/_fits.py rename to contexts/ingest/parsers/_fits.py index 7d0d5be..dc649eb 100644 --- a/contexts/ingest/parsers/solexs/_fits.py +++ b/contexts/ingest/parsers/_fits.py @@ -13,8 +13,16 @@ tolerated absence must say so at the call site, in one place, with a reason — as `.hk` absence is tolerated for v1.0 archives (§2.4) and an empty GTI is tolerated for SDD1 (F-12). -Three parsers use this (`lc`, `pi`, `gti`), which satisfies STD-11's two-instance rule; it -exists to centralise the ban rather than to anticipate a fourth. +RELOCATED BY M3/E5/#18. It was written for the three SoLEXS parsers and lived beside them; +the six HEL1OS parsers need the same ban, and importing it out of a sibling instrument's +package would make one instrument's internals another's dependency. Nine call sites now, +across two instruments — STD-11's two-instance rule satisfied twice over. + +The move changes no behaviour and no boundary: this module still imports only `domain.errors`, +it knows no instrument's conventions, SoLEXS's `lc`, `pi` and `gti` call it exactly as before, +and neither instrument package imports the other (asserted by the HEL1OS unit tests). #18 adds +one accessor, `declared_unit`, so that a declared unit is read with the same no-default +discipline as a keyword or a column. """ from __future__ import annotations @@ -93,4 +101,22 @@ def column(table: Any, name: str, *, source: str, hdu_name: str) -> Any: ) -__all__ = ["column", "expect", "fail", "hdu", "keyword"] +def declared_unit(table: Any, name: str, *, source: str, hdu_name: str) -> str | None: + """The unit a column declares (`TUNITn`), or `None` when the product declares none. + + `None` is returned as data, not replaced: an undeclared unit is a fact about the product + that a caller must handle at the call site, as §2.8 does for `czt2enth` (§8 A-4). The column + itself must exist — absence is F-04, exactly as in `column`. + """ + wanted = name.upper() + for candidate in table.columns: + if candidate.name.upper() == wanted: + return None if candidate.unit is None else str(candidate.unit) + fail( + "F-04", + f"/{source}#{hdu_name}/{name}", + f"column {name!r} absent; present: {list(table.columns.names)}", + ) + + +__all__ = ["column", "declared_unit", "expect", "fail", "hdu", "keyword"] diff --git a/contexts/ingest/parsers/hel1os/__init__.py b/contexts/ingest/parsers/hel1os/__init__.py new file mode 100644 index 0000000..10299e5 --- /dev/null +++ b/contexts/ingest/parsers/hel1os/__init__.py @@ -0,0 +1,27 @@ +"""HEL1OS parsers — `lc`, `spectra`, `hk`, `gti`, `events`, and orbit precedence (M3/E5/#18). + +Each implements its section of `SPEC-parsers@r7`: + + orbit §4 identity and version precedence. No merge API, no coverage map (#19). + lc §2.6 per-detector band rates; band edges parsed from EXTNAME; CTR unit read. + spectra §2.7 PHA spectra; 341 for CZT, 511 for CdTe; epoch resolution R-1 (r7 bound). + hk §2.8 housekeeping in ARCHIVE ORDER; every decisive column with its unit. + gti §2.9 per-detector good time intervals; lowercase column names. + events §2.5 photon events, row rules per r7. Exposed, and deliberately not ingestible. + +Where a comparison involves an MJD-derived time it uses §5.1's time-representation allowance +`ε_t` (CONTRA-007). §2.5's non-decreasing event rule is enforced as written although the archive +falsifies it; that falsification is recorded OPEN (CONTRA-008). + +TWO INSTRUMENTS, NO SHARED CONVENTIONS +-------------------------------------- +Nothing here imports `parsers.solexs`, and nothing there imports this. The two instruments +disagree on almost every convention — counts vs rates, undeclared vs declared units, Unix +seconds vs MJD, PI(340) vs PHA(341)/PHA(511), uppercase vs lowercase column names — and +F-07 and F-11 exist because assuming a shared convention across them is the specific error +that has already been made once. +""" + +from contexts.ingest.parsers.hel1os import events, gti, hk, lc, orbit, spectra + +__all__ = ["events", "gti", "hk", "lc", "orbit", "spectra"] diff --git a/contexts/ingest/parsers/hel1os/events.py b/contexts/ingest/parsers/hel1os/events.py new file mode 100644 index 0000000..560f3ce --- /dev/null +++ b/contexts/ingest/parsers/hel1os/events.py @@ -0,0 +1,317 @@ +"""HEL1OS `events/evt.fits` — photon event lists. §2.5 as amended at r7. + +THE SCOPE RULE IS PART OF THE SPECIFICATION +-------------------------------------------- +§2.5, binding: *"**Not required for the canonical tables** — retained for Phase 1a pile-up / +gain work. The 0.5.2 parser MUST expose an event reader but MUST NOT ingest events into the +canonical minute tables."* + +So this module reads events and produces no Observation. There is no `observations()` here, +and `test_the_event_reader_emits_no_observations` asserts that none appears — the canonical +tables are #19's, and an event stream that could flow into them would violate §2.5 by +default rather than by decision. + +VOLUME +------ +~85.7 GB across the corpus, and a single detector HDU can hold tens of millions of rows. +Everything here streams: `parse` reads headers only, and `events()` validates and yields one row +at a time from a freshly opened file. + +FOUR DETECTOR HDUs, OR F-03 +--------------------------- +§2.5: HDU1–4 are `CDTE1-EVENTS`, `CDTE2-EVENTS`, `CZT1-EVENTS`, `CZT2-EVENTS`, each carrying +`DETNAM`. F-03 fires if any is absent, and its rationale is *"silent detector loss"* — three +detectors' worth of events is not a smaller version of four, it is a different measurement. +The columns §2.5 lists are required in every HDU (F-04), and CZT HDUs additionally carry `pix` +and `offsetchn`, which the reader exposes rather than drops. + +`ener` IS ALREADY IN keV — AND THE DECLARATION IS READ +------------------------------------------------------ +§2.5: *"`ener` is **already energy-calibrated in keV** (unlike SoLEXS)."* `parse` reads each +HDU's declared `ener` unit and refuses anything but `keV` (F-07): stating keV over a different +declaration is exactly the cross-instrument assumption F-07 exists to prevent. + +THE ROW RULES — §2.5 AS AMENDED r7 +----------------------------------- +The rows are never held, so the row-level rules run inside `events()` as each row is read: + + mjd within [TSTART − ε_t, TSTOP + ε_t] F-06 §2.5 r7, §5.1 + V-EVT-2 |instant(mjd) − instant(utc-isot)| ≤ ½·r_isot + ε_t F-06 §2.5 r7 + mjd non-decreasing F-16 §2.5 + ener > 0 F-19 §2.5 + +`r_isot` is read from the string itself: the observed `YYYY-MM-DDTHH:MM:SS.sss` carries three +fractional digits, so it fixes the instant to 1 ms and a faithful string lies within 0.5 ms of +the `mjd` instant. `ε_t` (1 ms) absorbs float64 representation and nothing else. + +NON-DECREASING IS ENFORCED AS WRITTEN — AND THE ARCHIVE FALSIFIES IT +-------------------------------------------------------------------- +Real event HDUs step backward in `mjd` (CONTRA-008, OPEN). §2.5 has not been amended, so the +rule is enforced: a real stream terminates with F-16 at its first decrease, having yielded only +rows that passed every rule. Relaxing the rule here would be an amendment made in code, which +the contract forbids. Until the owner rules, the reader is exposed and every row rule is +exercised on real data, but a real stream does not run to completion. +""" + +from __future__ import annotations + +from collections.abc import Iterator +from dataclasses import dataclass +from datetime import datetime, timedelta, timezone +from pathlib import Path + +from contexts.ingest.parsers import _fits +from contexts.ingest.parsers.hel1os.lc import ( + MJD_UNIX_EPOCH, + SECONDS_PER_DAY, + TIME_REPRESENTATION_ALLOWANCE_S, + mjd_to_timestamp, +) +from contexts.ingest.parsers.hel1os.orbit import DETECTORS, family_of +from domain.values import Digest, Identifier, Timestamp + +#: §2.5 `OBSERVED`. Keyed by detector so a caller asks by detector, not by HDU index. +EVENT_HDU = { + "cdte1": "CDTE1-EVENTS", "cdte2": "CDTE2-EVENTS", + "czt1": "CZT1-EVENTS", "czt2": "CZT2-EVENTS", +} + +#: `DETNAM` as the archive spells it, per HDU. Case differs from the detector key, so the +#: comparison is case-insensitive — but presence and identity are still required. +DETNAM = {"cdte1": "CdTe1", "cdte2": "CdTe2", "czt1": "CZT1", "czt2": "CZT2"} + +ENERGY_UNIT = "keV" +ISOT_COLUMN = "utc-isot" + +#: §2.5's columns, required in every detector HDU; CZT additionally carries `pix`/`offsetchn`. +REQUIRED_COLUMNS = ("mjd", "hlsobt", "currtemp", "chn", "ener", "recnum", ISOT_COLUMN) +CZT_COLUMNS = ("pix", "offsetchn") + + +@dataclass(frozen=True) +class Event: + """One photon event. Energy in keV, as the archive declares it.""" + + mjd: float + hlsobt: float + energy_kev: float + channel: int + detector_temp_c: float + recnum: int + utc_isot: str + pixel: int | None = None + offset_channel: int | None = None + + @property + def valid_time(self) -> Timestamp: + return mjd_to_timestamp(self.mjd) + + +@dataclass(frozen=True) +class DetectorEvents: + """Header facts for one detector's event HDU. Rows are streamed, never held.""" + + detector: str + extname: str + detnam: str + rows: int + tstart_mjd: float + tstop_mjd: float + + @property + def instrument_id(self) -> Identifier: + return Identifier(f"hel1os-{self.detector}") + + +@dataclass(frozen=True) +class EventList: + """A validated `evt.fits`: all four detectors present, none read into memory.""" + + source_path: Path + source_digest: Digest + detectors: tuple[DetectorEvents, ...] + + def detector(self, name: str) -> DetectorEvents: + for candidate in self.detectors: + if candidate.detector == name.lower(): + return candidate + _fits.fail("F-03", f"/{self.source_path.name}", + f"no events for detector {name!r}; present: " + f"{[d.detector for d in self.detectors]}") + + def events(self, name: str) -> Iterator[Event]: + """Validate and yield one detector's events lazily, row by row (§2.5 r7).""" + from astropy.io import fits + + target = self.detector(name.lower()) + source = self.source_path.name + where = f"/{source}#{target.extname}" + allowance = TIME_REPRESENTATION_ALLOWANCE_S / SECONDS_PER_DAY + earliest, latest = target.tstart_mjd - allowance, target.tstop_mjd + allowance + czt = family_of(target.detector) == "czt" + + with fits.open(self.source_path) as hdul: + table = _fits.hdu(hdul, target.extname, source=source) + mjd = _fits.column(table, "mjd", source=source, hdu_name=target.extname) + hlsobt = _fits.column(table, "hlsobt", source=source, hdu_name=target.extname) + temp = _fits.column(table, "currtemp", source=source, hdu_name=target.extname) + chn = _fits.column(table, "chn", source=source, hdu_name=target.extname) + ener = _fits.column(table, "ener", source=source, hdu_name=target.extname) + recnum = _fits.column(table, "recnum", source=source, hdu_name=target.extname) + isot = _fits.column(table, ISOT_COLUMN, source=source, hdu_name=target.extname) + pix = (_fits.column(table, "pix", source=source, hdu_name=target.extname) + if czt else None) + offset = (_fits.column(table, "offsetchn", source=source, hdu_name=target.extname) + if czt else None) + + previous: float | None = None + for index in range(len(mjd)): + when = float(mjd[index]) + if not earliest <= when <= latest: + _fits.fail( + "F-06", f"{where}/mjd[{index}]", + f"mjd[{index}] = {when!r} is outside the header span " + f"[{target.tstart_mjd!r}, {target.tstop_mjd!r}] widened by ε_t = " + f"{TIME_REPRESENTATION_ALLOWANCE_S} s (§2.5 r7, §5.1)", + ) + + stamp = str(isot[index]) + residual, bound = isot_disagreement_s( + when, stamp, pointer=f"{where}/{ISOT_COLUMN}[{index}]") + if residual > bound: + _fits.fail( + "F-06", f"{where}/{ISOT_COLUMN}[{index}]", + f"V-EVT-2: {ISOT_COLUMN} {stamp!r} and mjd {when!r} disagree by " + f"{residual!r} s, beyond ½·r_isot + ε_t = {bound!r} s (§2.5 r7). Two " + f"representations of one instant that disagree are an ambiguous time.", + ) + + if previous is not None and when < previous: + _fits.fail( + "F-16", f"{where}/mjd[{index}]", + f"mjd decreases at row {index}: {previous!r} then {when!r}. §2.5 " + f"requires non-decreasing mjd; the archive falsifies the rule and the " + f"falsification is recorded OPEN in CONTRA-008, so it is enforced as " + f"written.", + ) + previous = when + + energy = float(ener[index]) + if not energy > 0: + _fits.fail( + "F-19", f"{where}/ener[{index}]", + f"ener[{index}] is {energy!r}; §2.5 validates ener > 0", + ) + + yield Event( + mjd=when, + hlsobt=float(hlsobt[index]), + energy_kev=energy, + channel=int(chn[index]), + detector_temp_c=float(temp[index]), + recnum=int(recnum[index]), + utc_isot=stamp, + pixel=int(pix[index]) if pix is not None else None, + offset_channel=int(offset[index]) if offset is not None else None, + ) + + +def isot_disagreement_s(mjd: float, isot: str, *, pointer: str) -> tuple[float, float]: + """V-EVT-2: how far apart the two representations are, and the bound they must meet. + + The bound is `½·r_isot + ε_t` (§2.5 r7). `r_isot` is the resolution the string carries — + 10⁻ⁿ s for n fractional digits — so it is read from the value, not chosen. The column is + named `utc-isot`; a string that declares a non-UTC offset contradicts that name and is + refused rather than converted. + """ + text = isot.strip() + try: + parsed = datetime.fromisoformat(text) + except ValueError: + _fits.fail("F-06", pointer, + f"{text!r} is not an ISO-8601 instant, so V-EVT-2 cannot be evaluated; " + f"§5 admits no ambiguous time") + if parsed.tzinfo is None: + parsed = parsed.replace(tzinfo=timezone.utc) + elif parsed.utcoffset() != timedelta(0): + _fits.fail("F-06", pointer, + f"{text!r} declares a non-UTC offset in a column named {ISOT_COLUMN!r}") + + fraction = text.partition(".")[2] + digits = len(fraction) - len(fraction.lstrip("0123456789")) + resolution = 10.0 ** -digits + residual = abs((mjd - MJD_UNIX_EPOCH) * SECONDS_PER_DAY - parsed.timestamp()) + return residual, 0.5 * resolution + TIME_REPRESENTATION_ALLOWANCE_S + + +def parse(path: Path, digest: Digest) -> EventList: + """Validate `evt.fits` headers: four detector HDUs, DETNAM, columns, `ener` unit, span.""" + from astropy.io import fits + + source = path.name + try: + opened = fits.open(path) + except OSError as exc: + _fits.fail("F-01", f"/{source}", f"not readable as FITS: {exc}") + + found: list[DetectorEvents] = [] + with opened as hdul: + present = {hdu.name.upper() for hdu in hdul} + missing = [ + EVENT_HDU[d] for d in DETECTORS if EVENT_HDU[d].upper() not in present + ] + if missing: + _fits.fail( + "F-03", f"/{source}", + f"event HDU(s) {missing} absent. §2.5 requires all four detectors; three " + f"detectors' events are a different measurement, not a smaller one.", + ) + + for detector in DETECTORS: + extname = EVENT_HDU[detector] + table = _fits.hdu(hdul, extname, source=source) + header = table.header + + detnam = str(_fits.keyword(header, "DETNAM", source=source, hdu_name=extname)) + if detnam.upper() != DETNAM[detector].upper(): + _fits.fail( + "F-03", f"/{source}#{extname}/DETNAM", + f"DETNAM is {detnam!r} in HDU {extname!r}, expected " + f"{DETNAM[detector]!r}. A mismatch means the HDU name and the detector " + f"it holds disagree.", + ) + + required = REQUIRED_COLUMNS + (CZT_COLUMNS if family_of(detector) == "czt" else ()) + declared = { + name: _fits.declared_unit(table, name, source=source, hdu_name=extname) + for name in required + } + if declared["ener"] != ENERGY_UNIT: + _fits.fail( + "F-07", f"/{source}#{extname}/ener", + f"ener declares unit {declared['ener']!r}; §2.5 records calibrated " + f"{ENERGY_UNIT!r}. An energy is never stated over a different declaration.", + ) + + tstart = float(_fits.keyword(header, "TSTART", source=source, hdu_name=extname)) + tstop = float(_fits.keyword(header, "TSTOP", source=source, hdu_name=extname)) + if not tstart < tstop: + _fits.fail("F-05", f"/{source}#{extname}/TSTART", + f"TSTART {tstart!r} / TSTOP {tstop!r} do not bound an interval") + + found.append(DetectorEvents( + detector=detector, + extname=extname, + detnam=detnam, + rows=int(_fits.keyword(header, "NAXIS2", source=source, hdu_name=extname)), + tstart_mjd=tstart, + tstop_mjd=tstop, + )) + + return EventList(source_path=path, source_digest=digest, detectors=tuple(found)) + + +__all__ = [ + "CZT_COLUMNS", "DETNAM", "ENERGY_UNIT", "EVENT_HDU", "ISOT_COLUMN", "REQUIRED_COLUMNS", + "DetectorEvents", "Event", "EventList", "isot_disagreement_s", "parse", +] diff --git a/contexts/ingest/parsers/hel1os/gti.py b/contexts/ingest/parsers/hel1os/gti.py new file mode 100644 index 0000000..aa12893 --- /dev/null +++ b/contexts/ingest/parsers/hel1os/gti.py @@ -0,0 +1,142 @@ +"""HEL1OS `aux/gti{czt,cdte}{1,2}.fits` — per-detector good time intervals. §2.9. + +`OBSERVED`: HDU1 `GTI_`; columns **lowercase** `tstart`, `tstop`, units undeclared; +`NAXIS2 = 1` on the sample orbit. + +THE CASE DIFFERENCE IS THE POINT +-------------------------------- +§2.9: *"Column-name case differs from SoLEXS (`START`/`STOP` uppercase) → all column access +MUST be case-insensitive (§8 A-2)."* `_fits.column` compares upper-cased names, so the same +helper reads both instruments without either knowing the other's spelling — and +`test_columns_are_read_case_insensitively` asserts it, because a case-sensitive lookup would +fail here as F-04 and look like a missing column rather than a spelling difference. + +UNITS ARE UNDECLARED, SO THE UNIT IS NOT ASSUMED +------------------------------------------------- +§2.9 records the units as undeclared. The values `OBSERVED` on the reference orbit are +`61017.00009887` / `61017.49984764`, which are MJD — the same encoding the other HEL1OS +products use, and far outside any plausible seconds interpretation. That reading is +*checked* rather than assumed: a bound outside the plausible MJD range terminates rather than +being reinterpreted, because a GTI misread by an epoch would silently include or exclude the +whole orbit. + +No `Σ(STOP−START+1) == EXPOSURE` check appears here. That is SoLEXS §2.3's rule, resting on +1-second sampling and a declared `EXPOSURE` keyword; §2.9 declares neither. Importing it +would be applying one instrument's convention to another, which F-07 forbids. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass +from pathlib import Path + +from contexts.ingest.parsers import _fits +from contexts.ingest.parsers.hel1os.lc import mjd_to_timestamp +from contexts.ingest.parsers.hel1os.orbit import family_of +from domain.values import Digest, Identifier, Timestamp + +#: The archive's HEL1OS coverage begins 2025-12-07 (MJD 61016). A bound outside a generous +#: window around the mission is not an MJD, and reinterpreting it would be a guess. +PLAUSIBLE_MJD = (50_000.0, 80_000.0) + + +@dataclass(frozen=True) +class Interval: + """One good time interval, in MJD. + + No `+1` second. That is SoLEXS §2.3's inclusive-second-mark convention, verified for a + 1 s-sampled product with a declared EXPOSURE; §2.9 declares neither, so the duration here + is the plain difference and is labelled as such. + """ + + start_mjd: float + stop_mjd: float + + def __post_init__(self) -> None: + for name, value in (("start", self.start_mjd), ("stop", self.stop_mjd)): + if not math.isfinite(value): + _fits.fail("F-16", f"/GTI/{name}", f"{name} is not finite ({value!r})") + if not PLAUSIBLE_MJD[0] <= value <= PLAUSIBLE_MJD[1]: + _fits.fail( + "F-05", f"/GTI/{name}", + f"{name} {value!r} is outside the plausible MJD range {PLAUSIBLE_MJD}. " + f"§2.9 leaves the unit undeclared; a bound that is not an MJD is not " + f"reinterpreted as another epoch.", + ) + if self.start_mjd >= self.stop_mjd: + _fits.fail("F-09", "/GTI", + f"start {self.start_mjd!r} is not before stop {self.stop_mjd!r}") + + @property + def duration_s(self) -> float: + return (self.stop_mjd - self.start_mjd) * 86_400.0 + + @property + def start_utc(self) -> Timestamp: + return mjd_to_timestamp(self.start_mjd) + + @property + def stop_utc(self) -> Timestamp: + return mjd_to_timestamp(self.stop_mjd) + + def covers(self, mjd: float) -> bool: + return self.start_mjd <= mjd <= self.stop_mjd + + +@dataclass(frozen=True) +class GoodTimeIntervals: + """A parsed HEL1OS GTI product.""" + + detector: str + source_path: Path + source_digest: Digest + intervals: tuple[Interval, ...] + detector_active: bool + + @property + def instrument_id(self) -> Identifier: + return Identifier(f"hel1os-{self.detector.lower()}") + + @property + def live_time_s(self) -> float: + return sum(interval.duration_s for interval in self.intervals) + + +def parse(path: Path, digest: Digest, *, detector: str) -> GoodTimeIntervals: + """Parse one `gti.fits`. F-12 applies: zero rows is legal, not an error.""" + from astropy.io import fits + + source = path.name + family_of(detector) + expected = f"GTI_{detector.upper()}" + + try: + opened = fits.open(path) + except OSError as exc: + _fits.fail("F-01", f"/{source}", f"not readable as FITS: {exc}") + + with opened as hdul: + table = _fits.hdu(hdul, expected, source=source) + rows = int(_fits.keyword(table.header, "NAXIS2", source=source, hdu_name=expected)) + + if rows == 0: + # F-12, the single deliberate non-terminating rule: the detector was inactive. + return GoodTimeIntervals(detector, path, digest, (), detector_active=False) + + starts = _fits.column(table, "tstart", source=source, hdu_name=expected) + stops = _fits.column(table, "tstop", source=source, hdu_name=expected) + intervals = tuple(Interval(float(a), float(b)) for a, b in zip(starts, stops)) + + for index in range(1, len(intervals)): + if intervals[index].start_mjd <= intervals[index - 1].stop_mjd: + _fits.fail( + "F-09", f"/{source}#{expected}/tstart[{index}]", + f"interval {index} starts at or before the previous stop; overlapping " + f"intervals would double-count live time", + ) + + return GoodTimeIntervals(detector, path, digest, intervals, detector_active=True) + + +__all__ = ["GoodTimeIntervals", "Interval", "PLAUSIBLE_MJD", "parse"] diff --git a/contexts/ingest/parsers/hel1os/hk.py b/contexts/ingest/parsers/hel1os/hk.py new file mode 100644 index 0000000..651c182 --- /dev/null +++ b/contexts/ingest/parsers/hel1os/hk.py @@ -0,0 +1,308 @@ +"""HEL1OS `aux/hk.fits` — housekeeping. §2.8, Phase 1a's key asset. + +THE PARSER PRESERVES ARCHIVE ORDER EXACTLY. IT PERFORMS NO SORTING. +-------------------------------------------------------------------- +§2.8, binding at r4, and stated as a general v2 principle rather than a HEL1OS special case: + +> *"The parser MUST preserve archive order exactly. It performs no sorting. The parser is a +> lossless representation of the archive — reading and transforming are separate acts, and a +> parser that silently reorders is no longer a faithful reader."* + +`mjd` here is a **measurement written in telemetry-arrival order**, not a sorted index. The +r0 requirement that it be non-decreasing was **falsified by the archive** and removed at r4. +What replaced it, with the header-span comparison made precise at r7: + + mjd finite F-16 + mjd unique F-16 — a repeated timestamp is a defect, not jitter + header-span consistency TSTART − ε_t ≤ min(mjd) and max(mjd) ≤ TSTOP + ε_t, else F-06 + inversion statistics RECORDED n_out_of_order and max_backward_step_s, never thresholded + czt1temp/czt2temp finite F-16 — exactly these two; §2.8 requires no others + suninfov in {0, 1} F-07 + +**No jitter threshold is defined.** §2.8: the magnitude is *reported*, never compared against +an invented tolerance. So `InversionStats` carries the numbers and nothing here compares them +to anything. `ε_t` is not such a tolerance: it is §5.1's allowance for float64 representation in +a comparison against a header bound, and it applies to that comparison alone (CONTRA-007 +Defect B, which also records CONTRADICTION-006's earlier implementation-only ruling). + +`OBSERVED` (orbit `HLS_20251208_000008`): 9,514 rows; 424 of 9,513 steps decrease; max +backward step 892.4 ms; 0 duplicates; range within the header span. Every one of those figures +is re-measured by the real-corpus tests. + +WHAT IS CAPTURED — §2.8's DECISIVE COLUMNS, BY THE ARCHIVE'S OWN NAMES +---------------------------------------------------------------------- +§2.8 lists the columns decisive for Phase 1a by concern: pile-up, saturation, gain/HV, thermal, +detector health, rates, pointing and time. Its detector-health line abbreviates per-detector +columns (`hotpixcnt`, `hotpixthr`, `hotpixlgcstat`, `bunpxctr`); the archive carries them only +with a detector prefix — `czt1hotpixcnt`, `czt2bunpxctr` — which is how §3 T4 names them too. +Those prefixed names are read and nothing beyond that spelling is inferred (CONTRA-007 +Observation E). + +Each column keeps its archive values and dtype — an integer counter stays an integer — and its +declared unit. Where §2.8 states a unit (V for the HV monitors, keV for `czt1enth`, degC for the +four temperatures, c/s for the four rates) the declaration is checked and a contradiction is +F-07. Where §2.8 states none, the declared unit is carried as the archive gives it, `None` +included, and no meaning is supplied. + +`suninfov` IS A FIRST-CLASS QUALITY FLAG +----------------------------------------- +§2.8, binding: *"data outside Sun-in-FOV is not solar signal. It MUST propagate to the +canonical tables."* Those tables are #19's, so this parser's obligation is to carry the flag +faithfully per row and to refuse a value outside {0, 1} rather than coercing it to a bool. + +THE `czt2enth` UNIT INCONSISTENCY IS RECORDED, NOT REPAIRED +------------------------------------------------------------ +§2.8: *"`czt1enth` has unit `keV` but `czt2enth` has unit=None for the same physical +quantity — a metadata inconsistency; the parser applies the `czt1enth` unit to both and +records the assumption (§8 A-4)."* So when `czt2enth` declares no unit it takes `czt1enth`'s and +`assumptions` says so; when it declares the same unit nothing is assumed; when it declares a +different one, the two readings of one quantity disagree and the parse terminates (F-07). +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass +from pathlib import Path + +from contexts.ingest.parsers import _fits +from contexts.ingest.parsers.hel1os.lc import ( + SECONDS_PER_DAY, + TIME_REPRESENTATION_ALLOWANCE_S, + mjd_to_timestamp, +) +from domain.values import Digest, Identifier, Timestamp + +HDU_NAME = "HLSHK" + +#: §2.8's decisive columns, by concern, in the archive's spelling. +PILEUP = ("cdte1pilectr", "cdte2pilectr") +SATURATION = ("czt1satctr1", "czt2satctr1") +GAIN_HV = ("czthvmon", "cdtehvmon", "czt1enth", "czt2enth", "cdte1enerthr", "cdte2enerthr") +THERMAL = ("czt1temp", "czt2temp", "cdte1temp", "cdte2temp") +DETECTOR_HEALTH = ( + "czt1hotpix", "czt2hotpix", "czt1hotpixcnt", "czt2hotpixcnt", + "czt1hotpixthr", "czt2hotpixthr", "czt1hotpixlgcstat", "czt2hotpixlgcstat", + "czt1bunpxctr", "czt2bunpxctr", "fehkstat", +) +RATES = ("czt1ctr", "czt2ctr", "cdte1ctr", "cdte2ctr") +POINTING = ("sunradeg", "sundecdeg", "suninfov", "sun2yawdeg", "sun2rolldeg", "sun2pitchdeg") +TIME = ( + "l0dhobt", "l0utcyr", "l0utcmon", "l0utcdy", "l0utchr", "l0utcmin", "l0utcsc", "l0utcmsc", +) + +#: Every decisive column carried in `columns`. `mjd` and `suninfov` are validated and carried +#: in their own fields. +CAPTURED = ( + PILEUP + SATURATION + GAIN_HV + THERMAL + DETECTOR_HEALTH + RATES + + tuple(name for name in POINTING if name != "suninfov") + TIME +) + +#: The units §2.8 states. A declared unit that contradicts one of these is F-07. +STATED_UNITS: dict[str, str] = { + "czthvmon": "V", "cdtehvmon": "V", "czt1enth": "keV", + "czt1temp": "degC", "czt2temp": "degC", "cdte1temp": "degC", "cdte2temp": "degC", + "czt1ctr": "c/s", "czt2ctr": "c/s", "cdte1ctr": "c/s", "cdte2ctr": "c/s", +} + +#: §2.8 requires exactly these finite. The CdTe temperatures are carried, NaN included. +REQUIRED_FINITE = ("czt1temp", "czt2temp") + +#: §8 A-4, recorded on the product whenever it is applied. +A4_ASSUMPTION = ( + "czt2enth carries unit=None in the archive while czt1enth carries keV for the same " + "physical quantity; the czt1enth unit is applied to both and the assumption is recorded " + "(SPEC-parsers r7 §2.8, §8 A-4)" +) + + +@dataclass(frozen=True) +class InversionStats: + """Recorded, never thresholded (§2.8 r4). + + §2.8 is explicit that no jitter threshold is defined and the magnitude is reported. So + this type has no `is_acceptable` and no limit to compare against — adding one would be + inventing the tolerance the amendment refused to invent. + """ + + rows: int + steps: int + n_out_of_order: int + max_backward_step_s: float + + @property + def max_backward_step_ms(self) -> float: + return self.max_backward_step_s * 1000.0 + + +@dataclass(frozen=True) +class Housekeeping: + """A parsed `hk.fits`, in archive order. + + `columns` holds every §2.8 decisive column as the archive's values; `units` holds each + column's declared unit, after A-4 where it applied. + """ + + source_path: Path + source_digest: Digest + mjd: tuple[float, ...] + suninfov: tuple[int, ...] + columns: dict[str, tuple] + units: dict[str, str | None] + header_tstart: float + header_tstop: float + inversions: InversionStats + assumptions: tuple[str, ...] = () + + @property + def instrument_id(self) -> Identifier: + return Identifier("hel1os") + + def valid_time(self, index: int) -> Timestamp: + return mjd_to_timestamp(self.mjd[index]) + + def sun_in_fov(self, index: int) -> bool: + """The first-class quality flag §2.8 requires to propagate.""" + return self.suninfov[index] == 1 + + +def parse(path: Path, digest: Digest) -> Housekeeping: + """Parse one `hk.fits`. Archive order preserved exactly; nothing is sorted.""" + from astropy.io import fits + + source = path.name + try: + opened = fits.open(path) + except OSError as exc: + _fits.fail("F-01", f"/{source}", f"not readable as FITS: {exc}") + + with opened as hdul: + table = _fits.hdu(hdul, HDU_NAME, source=source) + header = table.header + + mjd = tuple( + float(v) for v in _fits.column(table, "mjd", source=source, hdu_name=HDU_NAME) + ) + if not mjd: + _fits.fail("F-17", f"/{source}#{HDU_NAME}", + "HLSHK has no rows; §2.8 describes an orbit-resolved series") + for index, value in enumerate(mjd): + if not math.isfinite(value): + _fits.fail("F-16", f"/{source}#{HDU_NAME}/mjd[{index}]", + f"mjd[{index}] is not finite ({value!r})") + + # Duplicates remain F-16 (§2.8 r4): a repeated timestamp is a genuine defect, where + # an inversion is packet-arrival jitter. + if len(set(mjd)) != len(mjd): + _fits.fail( + "F-16", f"/{source}#{HDU_NAME}/mjd", + f"{len(mjd) - len(set(mjd))} duplicate timestamp(s). §2.8 r4 removed the " + f"non-decreasing requirement but kept this one: an inversion is jitter, a " + f"duplicate is a defect.", + ) + + header_tstart = float(_fits.keyword(header, "TSTART", source=source, + hdu_name=HDU_NAME)) + header_tstop = float(_fits.keyword(header, "TSTOP", source=source, + hdu_name=HDU_NAME)) + _check_header_span(mjd, header_tstart, header_tstop, source) + + suninfov_raw = _fits.column(table, "suninfov", source=source, hdu_name=HDU_NAME) + suninfov: list[int] = [] + for index, value in enumerate(suninfov_raw): + flag = int(value) + if flag not in (0, 1): + _fits.fail( + "F-07", f"/{source}#{HDU_NAME}/suninfov[{index}]", + f"suninfov[{index}] is {value!r}; §2.8 declares it a flag in {{0, 1}} " + f"and it is a first-class quality flag — coercing it would decide " + f"whether data is solar signal on the parser's authority", + ) + suninfov.append(flag) + + captured: dict[str, tuple] = {} + units: dict[str, str | None] = {} + for name in CAPTURED: + captured[name] = tuple( + _fits.column(table, name, source=source, hdu_name=HDU_NAME).tolist() + ) + units[name] = _fits.declared_unit(table, name, source=source, hdu_name=HDU_NAME) + units["suninfov"] = _fits.declared_unit(table, "suninfov", source=source, + hdu_name=HDU_NAME) + + for name, stated in STATED_UNITS.items(): + if units[name] != stated: + _fits.fail( + "F-07", f"/{source}#{HDU_NAME}/{name}", + f"{name} declares unit {units[name]!r}; §2.8 states {stated!r}. Reading it " + f"under the stated unit would assume a convention the product contradicts.", + ) + + assumptions: tuple[str, ...] = () + if units["czt2enth"] is None: + units["czt2enth"] = units["czt1enth"] + assumptions = (A4_ASSUMPTION,) + elif units["czt2enth"] != units["czt1enth"]: + _fits.fail( + "F-07", f"/{source}#{HDU_NAME}/czt2enth", + f"czt2enth declares {units['czt2enth']!r} and czt1enth declares " + f"{units['czt1enth']!r} for the same physical quantity; §8 A-4 covers an " + f"undeclared czt2enth unit, not a contradictory one", + ) + + for name in REQUIRED_FINITE: + for index, value in enumerate(captured[name]): + if not math.isfinite(value): + _fits.fail("F-16", f"/{source}#{HDU_NAME}/{name}[{index}]", + f"{name}[{index}] is not finite ({value!r}); §2.8 requires it") + + return Housekeeping( + source_path=path, + source_digest=digest, + mjd=mjd, + suninfov=tuple(suninfov), + columns=captured, + units=units, + header_tstart=header_tstart, + header_tstop=header_tstop, + inversions=inversion_stats(mjd), + assumptions=assumptions, + ) + + +def _check_header_span(mjd: tuple[float, ...], header_tstart: float, header_tstop: float, + source: str) -> None: + """§2.8 header-span consistency as amended at r7: bounds widened by `ε_t` only (§5.1).""" + allowance = TIME_REPRESENTATION_ALLOWANCE_S / SECONDS_PER_DAY + earliest, latest = min(mjd), max(mjd) + if earliest < header_tstart - allowance or latest > header_tstop + allowance: + _fits.fail( + "F-06", f"/{source}#{HDU_NAME}/mjd", + f"the mjd range [{earliest!r}, {latest!r}] is not inside the header span " + f"[{header_tstart!r}, {header_tstop!r}] widened by ε_t = " + f"{TIME_REPRESENTATION_ALLOWANCE_S} s (§2.8 r7, §5.1)", + ) + + +def inversion_stats(mjd: tuple[float, ...]) -> InversionStats: + """Count backward steps and measure the largest. Reported, never compared (§2.8 r4).""" + backward = 0 + largest = 0.0 + for index in range(1, len(mjd)): + step = mjd[index] - mjd[index - 1] + if step < 0: + backward += 1 + largest = max(largest, -step * SECONDS_PER_DAY) + return InversionStats( + rows=len(mjd), + steps=max(len(mjd) - 1, 0), + n_out_of_order=backward, + max_backward_step_s=largest, + ) + + +__all__ = [ + "A4_ASSUMPTION", "CAPTURED", "DETECTOR_HEALTH", "GAIN_HV", "HDU_NAME", "Housekeeping", + "InversionStats", "PILEUP", "POINTING", "RATES", "REQUIRED_FINITE", "SATURATION", + "STATED_UNITS", "THERMAL", "TIME", "inversion_stats", "parse", +] diff --git a/contexts/ingest/parsers/hel1os/lc.py b/contexts/ingest/parsers/hel1os/lc.py new file mode 100644 index 0000000..f38492e --- /dev/null +++ b/contexts/ingest/parsers/hel1os/lc.py @@ -0,0 +1,328 @@ +"""HEL1OS `lightcurve_{czt,cdte}{1,2}.fits` — per-detector band light curves. §2.6. + +THE OPPOSITE UNIT CONVENTION FROM SoLEXS +----------------------------------------- +§2.6: *"`CTR` is a **rate** with units **declared** — the opposite convention from SoLEXS +`.lc` (undeclared counts). The parser MUST NOT assume a shared convention across +instruments (F-07)."* + +So this module reads `cts/sec` and emits a rate, while `parsers/solexs/lc.py` reads +undeclared values that `HDUCLAS3` declares to be counts and emits counts. Two instruments, +two conventions, and the only thing stopping one from being read as the other is that each +parser reads its own product's declaration. Neither imports the other. + +BANDS ARE PARSED FROM `EXTNAME`, NEVER FROM POSITION +----------------------------------------------------- +§2.6, design rule: *"band edges are **parsed from `EXTNAME`** and validated against the +allowlist, never hardcoded by position — HDU order is not a contract."* + +`OBSERVED` in the archive: `CZT1_LC_BAND_20.00KEV_TO_40.00KEV` and its four siblings. An +`EXTNAME` outside the allowlist is **F-10** — *"never silently accept an unknown band"* — +because a band this code has not been told about would be ingested as though it were one it +had, and a rate would be attributed to the wrong energy range. + +`CTR`'s declared unit is read and must be `cts/sec` (F-07): the Observation's `cts/s` is stated +because the product declares a rate, never because this module expects one. + +EACH BAND KEEPS ITS OWN TIME AXIS +--------------------------------- +§2.6 records `NAXIS2 ≈ 43171` without saying whether the five bands of one detector share +timestamps. Archive-wide they often do not — every CdTe file and a few CZT files carry bands of +different lengths — so each `Band` keeps its own `MJD` column and nothing aligns them +(CONTRA-007 Observation D). + +`TIME_REPRESENTATION_ALLOWANCE_S` is defined here, beside the MJD conversion every HEL1OS module +already imports, because it belongs to comparisons of MJD-derived times: `SPEC-parsers@r7` §5.1. +""" + +from __future__ import annotations + +import math +import re +from collections.abc import Iterator +from dataclasses import dataclass +from pathlib import Path + +from contexts.ingest.parsers import _fits +from contexts.ingest.parsers.hel1os.orbit import family_of +from domain.entities import Observation +from domain.values import Digest, Identifier, Timestamp + +#: `CZT1_LC_BAND_20.00KEV_TO_40.00KEV` +BAND_EXTNAME = re.compile( + r"^(?P[A-Z0-9]+)_LC_BAND_(?P[0-9.]+)KEV_TO_(?P[0-9.]+)KEV$" +) + +#: §2.6's allowlist, per family. Four science bands plus one total, per detector. +#: An unlisted band terminates via F-10; the set is not extended by observation. +BANDS_KEV: dict[str, tuple[tuple[float, float], ...]] = { + "czt": ((20.0, 40.0), (40.0, 60.0), (60.0, 80.0), (80.0, 150.0), (18.0, 160.0)), + "cdte": ((5.0, 20.0), (20.0, 30.0), (30.0, 40.0), (40.0, 60.0), (1.8, 90.0)), +} + +#: The widest band of each family is the total, and it is the one that spans the others. +TOTAL_BAND: dict[str, tuple[float, float]] = {"czt": (18.0, 160.0), "cdte": (1.8, 90.0)} + +EXPECTED_BAND_COUNT = 5 + +#: MJD of the Unix epoch. HEL1OS states times in MJD (§2.5, §2.6), where SoLEXS states them +#: in Unix seconds — another convention that must not be shared across instruments. +MJD_UNIX_EPOCH = 40587.0 +SECONDS_PER_DAY = 86_400.0 + +#: `SPEC-parsers@r7` §5.1 — the time-representation allowance `ε_t`. Applied only where the +#: contract compares an MJD-derived time with a bound or with another representation of the +#: same instant (§2.5 span and V-EVT-2, §2.7 R-1 H3, §2.8 header span), never to a physical time +#: difference. It absorbs float64 representation, not physics (CONTRA-007). +TIME_REPRESENTATION_ALLOWANCE_S = 1e-3 + +#: §2.6 `OBSERVED`: `CTR` declares `cts/sec`. Read and checked, never assumed (F-07). +DECLARED_RATE_UNIT = "cts/sec" + + +def mjd_to_timestamp(mjd: float) -> Timestamp: + """MJD → a domain Timestamp, UTC. + + `OBSERVED` (§2.5): `TSTART = 61017.0000988685` is 2025-12-08. The conversion is exact + arithmetic on the declared epoch, not a library guess: MJD 40587 is 1970-01-01. + """ + from datetime import datetime, timezone + + seconds = (mjd - MJD_UNIX_EPOCH) * SECONDS_PER_DAY + moment = datetime.fromtimestamp(seconds, tz=timezone.utc) + return Timestamp(moment.isoformat().replace("+00:00", "Z")) + + +@dataclass(frozen=True) +class Band: + """One energy band of one detector, with edges read from its `EXTNAME`.""" + + detector: str + low_kev: float + high_kev: float + extname: str + mjd: tuple[float, ...] + rate: tuple[float, ...] + stat_err: tuple[float, ...] + + @property + def is_total(self) -> bool: + return (self.low_kev, self.high_kev) == TOTAL_BAND[family_of(self.detector)] + + @property + def quantity(self) -> str: + """The instrument's own term, carrying the band it belongs to. + + The energy range is part of what was measured: a rate in 20–40 keV and a rate in + 40–60 keV are different quantities, and a name that omitted the band would make them + indistinguishable once serialised. + """ + return f"count_rate_{self.low_kev:g}_{self.high_kev:g}_kev" + + +@dataclass(frozen=True) +class BandLightCurves: + """A parsed `lightcurve_*.fits`: five bands of one detector.""" + + detector: str + source_path: Path + source_digest: Digest + bands: tuple[Band, ...] + header_tstart: float + header_tstop: float + + @property + def instrument_id(self) -> Identifier: + return Identifier(f"hel1os-{self.detector.lower()}") + + @property + def family(self) -> str: + return family_of(self.detector) + + def band(self, low_kev: float, high_kev: float) -> Band: + for candidate in self.bands: + if (candidate.low_kev, candidate.high_kev) == (low_kev, high_kev): + return candidate + _fits.fail( + "F-10", f"/{self.source_path.name}", + f"no band {low_kev}-{high_kev} keV in {self.detector}; present: " + f"{[(b.low_kev, b.high_kev) for b in self.bands]}", + ) + + def observations(self, *, source_id: Identifier, ingest_time: Timestamp | None, + ) -> Iterator[Observation]: + """One Observation per (sample, band), streamed. + + The unit is `cts/s` — declared by the archive, not inferred — and is deliberately + different from the SoLEXS lightcurve's `counts`. Emitting both as one unit is the + cross-instrument mismatch F-07 exists to prevent. + """ + for band in self.bands: + for when, value in zip(band.mjd, band.rate): + yield Observation( + source_id=source_id, + instrument_id=self.instrument_id, + quantity=band.quantity, + unit="cts/s", + valid_time=mjd_to_timestamp(when), + ingest_time=ingest_time, + value=value, + source_digest=self.source_digest, + ) + + +def parse(path: Path, digest: Digest, *, detector: str) -> BandLightCurves: + """Parse one band-lightcurve product. Every §2.6 validation, no repair.""" + from astropy.io import fits + + source = path.name + family = family_of(detector) + allowed = BANDS_KEV[family] + + try: + opened = fits.open(path) + except OSError as exc: + _fits.fail("F-01", f"/{source}", f"not readable as FITS: {exc}") + + bands: list[Band] = [] + with opened as hdul: + extensions = [hdu for hdu in hdul[1:]] + if len(extensions) != EXPECTED_BAND_COUNT: + _fits.fail( + "F-10", f"/{source}", + f"{len(extensions)} band HDUs, expected exactly {EXPECTED_BAND_COUNT} " + f"(§2.6). Names present: {[h.name for h in extensions]}", + ) + + header_tstart = float( + _fits.keyword(extensions[0].header, "TSTART", source=source, + hdu_name=extensions[0].name) + ) + header_tstop = float( + _fits.keyword(extensions[0].header, "TSTOP", source=source, + hdu_name=extensions[0].name) + ) + + for hdu in extensions: + match = BAND_EXTNAME.match(hdu.name) + if match is None: + _fits.fail( + "F-10", f"/{source}#{hdu.name}", + f"EXTNAME {hdu.name!r} does not encode a band; §2.6 records the form " + f"_LC_BAND_KEV_TO_KEV", + ) + low, high = float(match.group("lo")), float(match.group("hi")) + if (low, high) not in allowed: + _fits.fail( + "F-10", f"/{source}#{hdu.name}", + f"band {low}-{high} keV is not in the {family} allowlist " + f"{list(allowed)}. §2.6: never silently accept an unknown band — an " + f"unlisted band would attribute a rate to the wrong energy range.", + ) + if match.group("det").lower() != detector.lower(): + _fits.fail( + "F-10", f"/{source}#{hdu.name}", + f"EXTNAME names detector {match.group('det')!r}, the file is " + f"{detector!r}", + ) + + ctr_unit = _fits.declared_unit(hdu, "CTR", source=source, hdu_name=hdu.name) + if ctr_unit != DECLARED_RATE_UNIT: + _fits.fail( + "F-07", f"/{source}#{hdu.name}/CTR", + f"CTR declares unit {ctr_unit!r}; §2.6 records a declared rate " + f"{DECLARED_RATE_UNIT!r}. Emitting a rate over any other declaration " + f"would assume a convention the product does not state.", + ) + + mjd = tuple( + float(v) for v in _fits.column(hdu, "MJD", source=source, hdu_name=hdu.name) + ) + rate = tuple( + float(v) for v in _fits.column(hdu, "CTR", source=source, hdu_name=hdu.name) + ) + stat_err = tuple( + float(v) + for v in _fits.column(hdu, "STAT_ERR", source=source, hdu_name=hdu.name) + ) + _validate_band(mjd, rate, source, hdu.name) + + bands.append( + Band( + detector=detector, + low_kev=low, + high_kev=high, + extname=hdu.name, + mjd=mjd, + rate=rate, + stat_err=stat_err, + ) + ) + + found = {(b.low_kev, b.high_kev) for b in bands} + if found != set(allowed): + _fits.fail( + "F-10", f"/{source}", + f"bands present {sorted(found)} do not match the {family} allowlist " + f"{sorted(allowed)}", + ) + + return BandLightCurves( + detector=detector, + source_path=path, + source_digest=digest, + bands=tuple(bands), + header_tstart=header_tstart, + header_tstop=header_tstop, + ) + + +def _validate_band(mjd, rate, source: str, hdu_name: str) -> None: + """§2.6: `MJD` strictly increasing, `CTR ≥ 0`. + + Strictly increasing here, unlike the housekeeping product where §2.8 r4 **removed** that + requirement after the archive falsified it. The two are different measurements: a light + curve is a binned series, housekeeping is telemetry in arrival order. Applying one + product's rule to the other is how a real defect gets tolerated or a real value rejected. + """ + for index, value in enumerate(mjd): + if not math.isfinite(value): + _fits.fail("F-16", f"/{source}#{hdu_name}/MJD[{index}]", + f"MJD[{index}] is not finite ({value!r})") + for index in range(1, len(mjd)): + if mjd[index] <= mjd[index - 1]: + _fits.fail( + "F-16", f"/{source}#{hdu_name}/MJD[{index}]", + f"MJD is not strictly increasing at row {index}: " + f"{mjd[index - 1]!r} then {mjd[index]!r}", + ) + for index, value in enumerate(rate): + if math.isnan(value): + # §2.6 declares no missing-value sentinel for CTR, unlike SoLEXS §2.1 which + # declares NaN. Treating a NaN here as "absent" would import another + # instrument's convention; the spec does not authorise it, so it terminates. + _fits.fail( + "F-07", f"/{source}#{hdu_name}/CTR[{index}]", + f"CTR[{index}] is NaN. §2.6 declares no missing-value sentinel for this " + f"product; SoLEXS §2.1 declares one, and importing that convention across " + f"instruments is exactly what F-07 forbids.", + ) + if value < 0: + _fits.fail("F-19", f"/{source}#{hdu_name}/CTR[{index}]", + f"CTR[{index}] is {value!r}; a negative rate is physically impossible") + + +__all__ = [ + "BANDS_KEV", + "BAND_EXTNAME", + "Band", + "BandLightCurves", + "DECLARED_RATE_UNIT", + "EXPECTED_BAND_COUNT", + "MJD_UNIX_EPOCH", + "SECONDS_PER_DAY", + "TIME_REPRESENTATION_ALLOWANCE_S", + "TOTAL_BAND", + "mjd_to_timestamp", + "parse", +] diff --git a/contexts/ingest/parsers/hel1os/orbit.py b/contexts/ingest/parsers/hel1os/orbit.py new file mode 100644 index 0000000..bbf437c --- /dev/null +++ b/contexts/ingest/parsers/hel1os/orbit.py @@ -0,0 +1,233 @@ +"""Orbit identity and version precedence — `SPEC-parsers@r7` §4. + +§4 is titled *"Version-Selection Policy (HEL1OS) — **naive ingestion made impossible**"*, and +the reason is measured rather than hypothetical: 46 overlapping orbit pairs exist, and their +`evt.fits` SHA-256 values **differ**, so they are genuine reprocessings rather than byte +copies. Deduplicating by content hash would silently keep both. + +IDENTIFICATION +-------------- +`HLS_(?P\\d{8})_(?P\\d{6})_(?P\\d+)sec_lev1_V(?P\\d{3})` + +`OBSERVED` distribution: V111 ×371, V211 ×16, V112 ×3, V311 ×1 — re-measured against the +archive by `test_the_version_distribution_matches_the_specification`, which found exactly +those counts across 391 orbits. + +**The three digits are undocumented in the archive and are treated as opaque** (§8 A-1). +Nothing here decodes them into major/minor/patch, because inventing that decomposition would +make V112 > V211 or the reverse depending on a reading nobody has authority for. They compare +as one integer, which is what §4 specifies. + +PRECEDENCE, DETERMINISTIC, IN ORDER +----------------------------------- +1. Higher `ver` integer wins. +2. Tie → longer `dur` wins (more coverage). +3. Tie → later header `DATE` (processing date) wins. +4. Still tied → **F-14 terminate. Never coin-flip.** + +WHAT THIS MODULE DELIBERATELY DOES NOT DO +------------------------------------------ +It does not merge orbits, and it offers no function that concatenates them. §4's mandatory +mechanism is a **minute-level coverage map** built *"before emitting T3/T4/T5"* — the +canonical minute tables — and §4.4 requires that *"the merge function MUST accept the +coverage map as a required argument. There is no API that concatenates orbit files +directly."* + +Those tables are the write path's, M3/E5/#19 (`grid`, `write`; row 19 depends on 17, 18 and +14). Building a coverage map here would require the minute grid this issue does not own, and +shipping a merge without it would create exactly the API §4 forbids. So #18 supplies the +*rules* — identity and precedence — and `test_no_orbit_merge_api_exists` asserts that the +concatenation §4 forbids has not appeared. +""" + +from __future__ import annotations + +import re +from dataclasses import dataclass +from datetime import datetime, timezone + +from contexts.ingest.parsers import _fits +from domain.values import Identifier, Timestamp + +#: §4's identification regex, verbatim. +ORBIT_STEM = re.compile( + r"^HLS_(?P\d{8})_(?P\d{6})_(?P\d+)sec_lev1_V(?P\d{3})$" +) + +#: The four detectors, and the two families they fall into (§2.6, §2.7). +CZT_DETECTORS = ("czt1", "czt2") +CDTE_DETECTORS = ("cdte1", "cdte2") +DETECTORS = CZT_DETECTORS + CDTE_DETECTORS + + +def family_of(detector: str) -> str: + """`czt` or `cdte`. Fails on anything else rather than guessing a family. + + The family decides the PHA channel space — 341 for CZT, 511 for CdTe (§2.7 r5) — so a + wrong family is a wrong channel space, which F-11 exists to prevent. + """ + lowered = detector.lower() + if lowered in CZT_DETECTORS: + return "czt" + if lowered in CDTE_DETECTORS: + return "cdte" + _fits.fail( + "F-07", f"/{detector}", + f"detector {detector!r} is in neither HEL1OS family; §2.6 and §2.7 name exactly " + f"{list(DETECTORS)}", + ) + + +@dataclass(frozen=True) +class OrbitId: + """One HEL1OS orbit archive, identified from its directory name (§4). + + `version` is an opaque integer. `processing_date` is the header `DATE`, supplied by a + caller that has read the product — this module reads no file, so precedence rule 3 takes + it as an argument rather than discovering it. + """ + + date: str + start: str + duration_s: int + version: int + processing_date: Timestamp | None = None + + @property + def start_epoch(self) -> float: + """Seconds since the Unix epoch for the declared start, from the stem alone. + + Derived from the directory name rather than from a header, because `overlaps` must + be answerable before any file is opened — §4 resolves coverage across orbits, and + opening 391 archives to compare intervals would defeat the point. + """ + moment = datetime.strptime(f"{self.date}{self.start}", "%Y%m%d%H%M%S").replace( + tzinfo=timezone.utc + ) + return moment.timestamp() + + @property + def stem(self) -> str: + return ( + f"HLS_{self.date}_{self.start}_{self.duration_s}sec_lev1_V{self.version:03d}" + ) + + @property + def instrument_id(self) -> Identifier: + return Identifier("hel1os") + + def detector_id(self, detector: str) -> Identifier: + """`hel1os-czt1`. The detector is part of the instrument identity (ADR-0003). + + HEL1OS carries four physically distinct detectors in two families with different + channel spaces; collapsing them to `hel1os` would make a CZT row and a CdTe row + indistinguishable, which is the confusion F-11 forbids. + """ + family_of(detector) + return Identifier(f"hel1os-{detector.lower()}") + + +def parse_stem(name: str) -> OrbitId: + """Read an orbit identity from a directory name. F-18 on anything that is not one.""" + match = ORBIT_STEM.match(name) if isinstance(name, str) else None + if match is None: + _fits.fail( + "F-18", f"/{name}", + f"{name!r} is not a HEL1OS orbit stem; §4 identifies orbits as " + f"HLS___sec_lev1_V", + ) + return OrbitId( + date=match.group("date"), + start=match.group("start"), + duration_s=int(match.group("dur")), + version=int(match.group("ver")), + ) + + +def precedence(first: OrbitId, second: OrbitId) -> OrbitId: + """Which of two orbits wins, by §4's four rules in order. + + Returns the winner. Raises F-14 when all four are exhausted — §4 is explicit that the + fourth outcome is termination: **"Never coin-flip."** A tie that reached rule 4 means two + products claim the same coverage with identical version, duration and processing date, + and nothing in the archive distinguishes them; choosing either would be a decision the + data does not support. + """ + if first.version != second.version: + return first if first.version > second.version else second + + if first.duration_s != second.duration_s: + return first if first.duration_s > second.duration_s else second + + if first.processing_date is not None and second.processing_date is not None: + if first.processing_date.instant != second.processing_date.instant: + return ( + first + if first.processing_date.instant > second.processing_date.instant + else second + ) + + _fits.fail( + "F-14", f"/{first.stem}", + f"version precedence is unresolved between {first.stem!r} and {second.stem!r}: " + f"equal version V{first.version:03d}, equal duration {first.duration_s}s, and " + f"processing dates " + f"{'absent' if first.processing_date is None or second.processing_date is None else 'equal'}" + f". §4 rule 4 terminates rather than choosing — never coin-flip.", + ) + + +def rule_applied(first: OrbitId, second: OrbitId) -> str: + """Which of §4's rules decided, for the resolution log §4.2 requires. + + §4.2: *"log every resolution … with winner, losers, and rule invoked."* The log itself is + written by the merge step (#19); naming the rule is part of the precedence decision and + so belongs with it. + """ + if first.version != second.version: + return "rule-1-higher-version" + if first.duration_s != second.duration_s: + return "rule-2-longer-duration" + if ( + first.processing_date is not None + and second.processing_date is not None + and first.processing_date.instant != second.processing_date.instant + ): + return "rule-3-later-processing-date" + return "rule-4-terminate" + + +def overlaps(first: OrbitId, second: OrbitId) -> bool: + """Whether two orbits claim any second in common. + + §4 names two classes, both `OBSERVED` in the archive and both present locally: + Class A — identical interval, different version. + Class B — partial overlap, different start and duration. + + Class B is why §4 says file-level selection is *insufficient*: each file covers seconds + the other lacks. This predicate reports the overlap; it does not resolve it, because + resolving Class B needs the minute-level coverage map that belongs to #19. + + §4 defines no overlap measure. Over the stem-declared intervals of the 391-orbit archive + this predicate reports 50 overlapping pairs, where §4 records 46 and the Milestone VI engine + recorded 49; the difference is recorded, not reconciled (CONTRA-007 Observation F). + """ + first_start, second_start = first.start_epoch, second.start_epoch + return ( + first_start < second_start + second.duration_s + and second_start < first_start + first.duration_s + ) + + +__all__ = [ + "CDTE_DETECTORS", + "CZT_DETECTORS", + "DETECTORS", + "ORBIT_STEM", + "OrbitId", + "family_of", + "overlaps", + "parse_stem", + "precedence", + "rule_applied", +] diff --git a/contexts/ingest/parsers/hel1os/spectra.py b/contexts/ingest/parsers/hel1os/spectra.py new file mode 100644 index 0000000..ed92769 --- /dev/null +++ b/contexts/ingest/parsers/hel1os/spectra.py @@ -0,0 +1,326 @@ +"""HEL1OS `hel1os_{czt,cdte}_spectra_{det}.fits` — per-detector PHA spectra. §2.7. + +TWO CHANNEL SPACES, VALIDATED AGAINST A FAMILY-SPECIFIC ALLOWLIST +------------------------------------------------------------------ +§2.7 as amended at r5: HEL1OS has **two detector families with different PHA channel +spaces**, and `DETCHANS` is validated against a family allowlist, *"never a single scalar"*. + + CZT (czt1, czt2) DETCHANS 341 CHANTYPE PHA + CdTe (cdte1, cdte2) DETCHANS 511 CHANTYPE PHA + +An unlisted `(family, DETCHANS)` pair terminates via **F-07**, exactly as an unlisted band +terminates via F-10. + +And a third space exists one module away: SoLEXS is **340 PI** channels at 1 s (§2.2). +§2.7 is explicit — *"`PI ≠ PHA` (gain-corrected vs raw pulse height) and 340 ≠ 341. No v2 +code may treat these as a common channel space (F-11)."* This module cannot merge them +because it knows nothing of SoLEXS, and a test asserts the three counts stay distinct. + +EPOCH RESOLUTION R-1 (§2.7, AMENDED r4 AND r7) +----------------------------------------------- +The r0 framing read the conflict between column `unit='s'` and header `TSTART`-as-MJD as an +*epoch* ambiguity. It is not. The `unit='s'` declaration is correct, and both metadata +statements are true and **compose**: + + absolute_time = header TSTART (MJD → UTC) + column TSTART (offset seconds) + +Three hypotheses are evaluated in order, and **only if all fail does F-06 terminate**: + + H3 relative seconds from header TSTART col[0] == 0 exactly, and the column span + agrees with the header span within one + EXPOSURE bin + ε_t (r7) + H1 MJD days |col[0] − header TSTART| × 86400 ≤ 1 s + H2 Unix seconds |col[0] − unix(header TSTART)| ≤ 1 s + +H3 is tested first because `unit='s'` literally declares seconds. H1 and H2 are retained +because a reprocessed product could legitimately switch to an absolute epoch, and silently +misreading one would be worse than an extra branch. + +**The resolved hypothesis and its residual are recorded** on the parse result, as §2.7 +requires them to be recorded in T7 provenance — the recording itself is #19's, the +measurement is this parser's. + +`col_span` AND THE COMPARISON PRECISION — §2.7 AS AMENDED AT r7 +---------------------------------------------------------------- +r4 left `col_span` undefined and stated the H3 bound exactly. r7 defines both (CONTRA-007 +Defect A): + + col_span = column TSTOP[last] − column TSTART[0] the covered interval + header_span = (header TSTOP − header TSTART) × 86400 s + accept H3 iff col[0] == 0 exactly and |col_span − header_span| ≤ EXPOSURE + ε_t + +`ε_t` is §5.1's time-representation allowance, 1 ms. It absorbs float64 representation and +nothing physical: on the reference orbit `(61017.49963590554 − 61017.0000988685) × 86400` +evaluates to `43160.00000007916`, 79 ns over an exact one-bin bound, and on other valid products +the excess reaches about a microsecond. A genuine disagreement is at least a further 20 s bin, +and a two-bin disagreement fails by a wide margin. The value is the contract's, not this +module's — an allowance chosen here would be the silent tolerance r7 exists to rule out. +""" + +from __future__ import annotations + +import math +from collections.abc import Iterator +from dataclasses import dataclass +from pathlib import Path + +from contexts.ingest.parsers import _fits +from contexts.ingest.parsers.hel1os.lc import ( + MJD_UNIX_EPOCH, + SECONDS_PER_DAY, + TIME_REPRESENTATION_ALLOWANCE_S, + mjd_to_timestamp, +) +from contexts.ingest.parsers.hel1os.orbit import family_of +from domain.values import Digest, Identifier, Timestamp + +#: §2.7 r5's family allowlist. Never a single scalar. +DETCHANS_BY_FAMILY: dict[str, int] = {"czt": 341, "cdte": 511} + +#: The three incommensurable channel spaces F-11 names, kept here so a test can assert they +#: stay distinct. SoLEXS's 340 is stated as a number, not imported: importing it would create +#: the very coupling between channel spaces the rule forbids. +INCOMMENSURABLE_CHANNEL_SPACES = { + ("solexs", "PI"): 340, + ("czt", "PHA"): 341, + ("cdte", "PHA"): 511, +} + +CHANTYPE = "PHA" +#: §2.7 `OBSERVED`: `HDUCLAS3='COUNT'` — singular, unlike SoLEXS's `'COUNTS'`. +HDUCLAS3 = "COUNT" + + +@dataclass(frozen=True) +class EpochResolution: + """The outcome of R-1, recorded rather than assumed (§2.7).""" + + hypothesis: str + residual_s: float + header_tstart_mjd: float + exposure_s: float + + def absolute_time(self, column_tstart: float) -> Timestamp: + """Compose the resolved offset onto the header epoch.""" + if self.hypothesis == "H3": + return mjd_to_timestamp( + self.header_tstart_mjd + column_tstart / SECONDS_PER_DAY + ) + if self.hypothesis == "H1": + return mjd_to_timestamp(column_tstart) + return mjd_to_timestamp(MJD_UNIX_EPOCH + column_tstart / SECONDS_PER_DAY) + + +@dataclass(frozen=True) +class Spectrum: + """One ~20-second PHA spectrum.""" + + spec_num: int + tstart: float + tstop: float + exposure: float + counts: tuple[float, ...] + stat_err: tuple[float, ...] + valid_time: Timestamp + + @property + def total_counts(self) -> float: + return math.fsum(self.counts) + + +@dataclass(frozen=True) +class SpectraHeader: + detector: str + family: str + detchans: int + chantype: str + hduclas3: str + tstart_mjd: float + tstop_mjd: float + rows: int + + +@dataclass(frozen=True) +class Spectra: + """A validated spectra product. Rows are streamed, never held.""" + + header: SpectraHeader + source_path: Path + source_digest: Digest + channel_map: tuple[int, ...] + epoch: EpochResolution + + @property + def instrument_id(self) -> Identifier: + return Identifier(f"hel1os-{self.header.detector.lower()}") + + def spectra(self) -> Iterator[Spectrum]: + """Yield each spectrum lazily from a freshly opened file (E5 §12).""" + from astropy.io import fits + + source = self.source_path.name + with fits.open(self.source_path) as hdul: + table = _fits.hdu(hdul, "SPECTRUM", source=source) + spec_nums = _fits.column(table, "SPEC_NUM", source=source, hdu_name="SPECTRUM") + tstarts = _fits.column(table, "TSTART", source=source, hdu_name="SPECTRUM") + tstops = _fits.column(table, "TSTOP", source=source, hdu_name="SPECTRUM") + exposures = _fits.column(table, "EXPOSURE", source=source, hdu_name="SPECTRUM") + counts = _fits.column(table, "COUNTS", source=source, hdu_name="SPECTRUM") + errors = _fits.column(table, "STAT_ERR", source=source, hdu_name="SPECTRUM") + + for index in range(len(tstarts)): + row = tuple(float(v) for v in counts[index]) + for channel, value in enumerate(row): + if value < 0: + _fits.fail( + "F-19", f"/{source}#SPECTRUM/COUNTS[{index}][{channel}]", + f"count {value!r} is negative, which is physically impossible", + ) + exposure = float(exposures[index]) + if exposure < 0: + _fits.fail("F-19", f"/{source}#SPECTRUM/EXPOSURE[{index}]", + f"EXPOSURE[{index}] is {exposure!r}") + yield Spectrum( + spec_num=int(spec_nums[index]), + tstart=float(tstarts[index]), + tstop=float(tstops[index]), + exposure=exposure, + counts=row, + stat_err=tuple(float(v) for v in errors[index]), + valid_time=self.epoch.absolute_time(float(tstarts[index])), + ) + + +def resolve_epoch(column_tstart, column_tstop, exposures, header_tstart: float, + header_tstop: float, source: str) -> EpochResolution: + """R-1, in the order §2.7 mandates. F-06 only if all three hypotheses fail.""" + first = float(column_tstart[0]) + exposure = float(exposures[0]) + header_span = (header_tstop - header_tstart) * SECONDS_PER_DAY + + # ── H3: relative seconds from the header epoch (the observed convention) ───── + # §2.7 r7: col_span runs from the first bin's start to the last bin's end, and the one-bin + # bound carries §5.1's ε_t — representation, never physics (CONTRA-007 Defect A). + column_span = float(column_tstop[-1]) - first + residual = abs(column_span - header_span) + if first == 0.0 and residual <= exposure + TIME_REPRESENTATION_ALLOWANCE_S: + return EpochResolution("H3", residual, header_tstart, exposure) + + # ── H1: MJD days ──────────────────────────────────────────────────────────── + h1_residual = abs(first - header_tstart) * SECONDS_PER_DAY + if h1_residual <= 1.0: + return EpochResolution("H1", h1_residual, header_tstart, exposure) + + # ── H2: Unix seconds ──────────────────────────────────────────────────────── + header_unix = (header_tstart - MJD_UNIX_EPOCH) * SECONDS_PER_DAY + h2_residual = abs(first - header_unix) + if h2_residual <= 1.0: + return EpochResolution("H2", h2_residual, header_tstart, exposure) + + _fits.fail( + "F-06", f"/{source}#SPECTRUM/TSTART[0]", + f"epoch resolution R-1 failed for all three hypotheses. column TSTART[0]={first!r}, " + f"header TSTART={header_tstart!r} MJD, column span={column_span!r}s, header " + f"span={header_span!r}s, EXPOSURE={exposure!r}s. H3 residual={residual!r}s, " + f"H1={h1_residual!r}s, H2={h2_residual!r}s. §5 admits no ambiguous time.", + ) + + +def parse(path: Path, digest: Digest, *, detector: str) -> Spectra: + """Validate one spectra product and resolve its epoch. Row data is not read here.""" + from astropy.io import fits + + source = path.name + family = family_of(detector) + expected_detchans = DETCHANS_BY_FAMILY[family] + + try: + opened = fits.open(path) + except OSError as exc: + _fits.fail("F-01", f"/{source}", f"not readable as FITS: {exc}") + + with opened as hdul: + table = _fits.hdu(hdul, "SPECTRUM", source=source) + header = table.header + + _fits.expect(header, "CHANTYPE", CHANTYPE, source=source, + hdu_name="SPECTRUM", rule="F-07") + _fits.expect(header, "HDUCLAS3", HDUCLAS3, source=source, + hdu_name="SPECTRUM", rule="F-07") + + detchans = int(_fits.keyword(header, "DETCHANS", source=source, + hdu_name="SPECTRUM")) + if detchans != expected_detchans: + _fits.fail( + "F-07", f"/{source}#SPECTRUM/DETCHANS", + f"DETCHANS is {detchans} for family {family!r}, which the §2.7 r5 allowlist " + f"gives as {expected_detchans}. An unlisted (family, DETCHANS) pair " + f"terminates: {DETCHANS_BY_FAMILY}", + ) + + channel_column = _fits.column(table, "CHANNEL", source=source, hdu_name="SPECTRUM") + channel_map = _constant_channel_map(channel_column, detchans, source) + + tstart_mjd = float(_fits.keyword(header, "TSTART", source=source, + hdu_name="SPECTRUM")) + tstop_mjd = float(_fits.keyword(header, "TSTOP", source=source, + hdu_name="SPECTRUM")) + rows = int(_fits.keyword(header, "NAXIS2", source=source, hdu_name="SPECTRUM")) + + epoch = resolve_epoch( + _fits.column(table, "TSTART", source=source, hdu_name="SPECTRUM"), + _fits.column(table, "TSTOP", source=source, hdu_name="SPECTRUM"), + _fits.column(table, "EXPOSURE", source=source, hdu_name="SPECTRUM"), + tstart_mjd, tstop_mjd, source, + ) + + built = SpectraHeader( + detector=detector, + family=family, + detchans=detchans, + chantype=CHANTYPE, + hduclas3=HDUCLAS3, + tstart_mjd=tstart_mjd, + tstop_mjd=tstop_mjd, + rows=rows, + ) + + return Spectra( + header=built, + source_path=path, + source_digest=digest, + channel_map=channel_map, + epoch=epoch, + ) + + +def _constant_channel_map(column, detchans: int, source: str) -> tuple[int, ...]: + """F-08: the CHANNEL vector must be identical across rows. + + Same rule as SoLEXS §2.2, and the reason is the same: a varying map means a channel index + does not mean the same thing in every row. The two are checked separately in separate + modules, because sharing the check would require sharing a channel space. + """ + first = tuple(int(v) for v in column[0]) + if len(first) != detchans: + _fits.fail("F-08", f"/{source}#SPECTRUM/CHANNEL[0]", + f"CHANNEL vector has {len(first)} entries, DETCHANS declares {detchans}") + for index in range(1, len(column)): + if tuple(int(v) for v in column[index]) != first: + _fits.fail("F-08", f"/{source}#SPECTRUM/CHANNEL[{index}]", + f"CHANNEL vector at row {index} differs from row 0") + return first + + +__all__ = [ + "CHANTYPE", + "DETCHANS_BY_FAMILY", + "EpochResolution", + "HDUCLAS3", + "INCOMMENSURABLE_CHANNEL_SPACES", + "Spectra", + "SpectraHeader", + "Spectrum", + "parse", + "resolve_epoch", +] diff --git a/contexts/ingest/parsers/solexs/gti.py b/contexts/ingest/parsers/solexs/gti.py index 38c093b..4845e82 100644 --- a/contexts/ingest/parsers/solexs/gti.py +++ b/contexts/ingest/parsers/solexs/gti.py @@ -43,7 +43,7 @@ from dataclasses import dataclass from pathlib import Path -from contexts.ingest.parsers.solexs import _fits +from contexts.ingest.parsers import _fits from contexts.ingest.parsers.solexs.lc import unix_to_timestamp from domain.values import Digest, Identifier, Timestamp diff --git a/contexts/ingest/parsers/solexs/lc.py b/contexts/ingest/parsers/solexs/lc.py index e64cbb6..f528957 100644 --- a/contexts/ingest/parsers/solexs/lc.py +++ b/contexts/ingest/parsers/solexs/lc.py @@ -41,7 +41,7 @@ from datetime import datetime, timezone from pathlib import Path -from contexts.ingest.parsers.solexs import _fits +from contexts.ingest.parsers import _fits from domain.entities import Observation from domain.values import Digest, Identifier, Timestamp diff --git a/contexts/ingest/parsers/solexs/pi.py b/contexts/ingest/parsers/solexs/pi.py index c819bc5..a8c6b81 100644 --- a/contexts/ingest/parsers/solexs/pi.py +++ b/contexts/ingest/parsers/solexs/pi.py @@ -48,7 +48,7 @@ from dataclasses import dataclass from pathlib import Path -from contexts.ingest.parsers.solexs import _fits +from contexts.ingest.parsers import _fits from contexts.ingest.parsers.solexs.lc import LightCurve, unix_to_timestamp from domain.values import Digest, Identifier, Timestamp diff --git a/contexts/ingest/tests/hel1os_fixtures.py b/contexts/ingest/tests/hel1os_fixtures.py new file mode 100644 index 0000000..951562b --- /dev/null +++ b/contexts/ingest/tests/hel1os_fixtures.py @@ -0,0 +1,264 @@ +"""FITS fixtures with the layouts `SPEC-parsers@r7` §2.5–§2.9 record. + +Real FITS written by astropy, carrying the headers, EXTNAMEs, column names and declared units +the archive carries — so the parsers can be exercised without the 132 GB corpus, and so a +*violating* product can be built deliberately. A fail-loud rule nobody has watched fire is a +comment. + +Every default is a value `OBSERVED` in the spec or in the archive: `TSTART = 61017.0000988685` +(2025-12-08), `EXPOSURE = 20.0` s, `DETCHANS` 341 for CZT and 511 for CdTe, `CHANTYPE='PHA'`, +`HDUCLAS3='COUNT'` (singular), band EXTNAMEs of the `_LC_BAND_KEV_TO_KEV` form, +`CTR` declared `cts/sec`, lowercase `tstart`/`tstop` in the GTI products, the housekeeping +columns and units of §2.8, and `utc-isot` strings rounded to the millisecond. + +Row counts are small by default. None of the real counts is validated against a fixed +expectation by any HEL1OS rule, unlike SoLEXS's F-17, so a short fixture exercises the same +code paths. +""" + +from __future__ import annotations + +from datetime import datetime, timezone +from pathlib import Path + +import numpy as np +from astropy.io import fits + +#: 2025-12-08T00:00:08.5Z — the reference orbit's start (§2.5, §2.7 `OBSERVED`). +TSTART_MJD = 61017.0000988685 +EXPOSURE_S = 20.0 +SECONDS_PER_DAY = 86_400.0 + +CZT_BANDS = ((20.0, 40.0), (40.0, 60.0), (60.0, 80.0), (80.0, 150.0), (18.0, 160.0)) +CDTE_BANDS = ((5.0, 20.0), (20.0, 30.0), (30.0, 40.0), (40.0, 60.0), (1.8, 90.0)) +BANDS = {"czt": CZT_BANDS, "cdte": CDTE_BANDS} +DETCHANS = {"czt": 341, "cdte": 511} + +#: §2.8's decisive housekeeping columns: (FITS format, declared unit) as the archive carries them. +HK_COLUMNS: dict[str, tuple[str, str | None]] = { + "cdte1pilectr": ("K", None), "cdte2pilectr": ("K", None), + "czt1satctr1": ("K", None), "czt2satctr1": ("K", None), + "czthvmon": ("D", "V"), "cdtehvmon": ("D", "V"), + "czt1enth": ("D", "keV"), "czt2enth": ("D", None), + "cdte1enerthr": ("D", None), "cdte2enerthr": ("D", None), + "czt1temp": ("D", "degC"), "czt2temp": ("D", "degC"), + "cdte1temp": ("D", "degC"), "cdte2temp": ("D", "degC"), + "czt1hotpix": ("K", None), "czt2hotpix": ("K", None), + "czt1hotpixcnt": ("K", None), "czt2hotpixcnt": ("K", None), + "czt1hotpixthr": ("K", None), "czt2hotpixthr": ("K", None), + "czt1hotpixlgcstat": ("K", None), "czt2hotpixlgcstat": ("K", None), + "czt1bunpxctr": ("K", None), "czt2bunpxctr": ("K", None), + "fehkstat": ("K", None), + "czt1ctr": ("D", "c/s"), "czt2ctr": ("D", "c/s"), + "cdte1ctr": ("D", "c/s"), "cdte2ctr": ("D", "c/s"), + "sunradeg": ("D", None), "sundecdeg": ("D", None), + "sun2yawdeg": ("D", None), "sun2rolldeg": ("D", None), "sun2pitchdeg": ("D", None), + "l0dhobt": ("D", None), + "l0utcyr": ("J", None), "l0utcmon": ("J", None), "l0utcdy": ("J", None), + "l0utchr": ("J", None), "l0utcmin": ("J", None), "l0utcsc": ("J", None), + "l0utcmsc": ("J", None), +} + + +def band_extname(detector: str, low: float, high: float) -> str: + return f"{detector.upper()}_LC_BAND_{low:.2f}KEV_TO_{high:.2f}KEV" + + +def isot_of(mjd: float) -> str: + """The `utc-isot` string the archive writes for an mjd: the instant to the millisecond.""" + seconds = round((mjd - 40587.0) * SECONDS_PER_DAY, 3) + moment = datetime.fromtimestamp(seconds, tz=timezone.utc) + return moment.strftime("%Y-%m-%dT%H:%M:%S.%f")[:23] + + +def write_lightcurve(path: Path, *, detector: str = "czt1", rows: int = 120, + bands: tuple[tuple[float, float], ...] | None = None, + rates: dict[tuple[float, float], list[float]] | None = None, + mjd: list[float] | None = None, + extnames: list[str] | None = None, + ctr_unit: str | None = "cts/sec") -> Path: + """A `lightcurve_.fits` per §2.6: five band HDUs, `MJD`/`ISOT`/`CTR`/`STAT_ERR`.""" + family = "czt" if detector.lower().startswith("czt") else "cdte" + bands = bands if bands is not None else BANDS[family] + if mjd is None: + mjd = [TSTART_MJD + i / SECONDS_PER_DAY for i in range(rows)] + + hdus = [fits.PrimaryHDU()] + for index, (low, high) in enumerate(bands): + values = (rates or {}).get((low, high), [10.0 + i for i in range(len(mjd))]) + name = extnames[index] if extnames else band_extname(detector, low, high) + table = fits.BinTableHDU.from_columns([ + fits.Column(name="MJD", format="D", unit="MJD", array=np.array(mjd, dtype=float)), + fits.Column(name="ISOT", format="30A", unit="UT", + array=np.array(["2025-12-08T00:00:00.0"] * len(mjd))), + fits.Column(name="CTR", format="D", unit=ctr_unit, + array=np.array(values, dtype=float)), + fits.Column(name="STAT_ERR", format="D", unit="cts/sec", + array=np.ones(len(mjd), dtype=float)), + ], name=name) + table.header["TSTART"] = TSTART_MJD + table.header["TSTOP"] = mjd[-1] + hdus.append(table) + + path.parent.mkdir(parents=True, exist_ok=True) + fits.HDUList(hdus).writeto(path, overwrite=True) + return path + + +def write_spectra(path: Path, *, detector: str = "czt1", rows: int = 20, + detchans: int | None = None, chantype: str = "PHA", + hduclas3: str = "COUNT", col_tstart: list[float] | None = None, + channel_map: list[list[int]] | None = None, + counts: np.ndarray | None = None, + header_tstart: float = TSTART_MJD, + header_tstop: float | None = None, + overrides: dict | None = None) -> Path: + """A `hel1os__spectra_.fits` per §2.7 — Type II PHA, H3 by construction.""" + family = "czt" if detector.lower().startswith("czt") else "cdte" + detchans = detchans if detchans is not None else DETCHANS[family] + if col_tstart is None: + col_tstart = [i * EXPOSURE_S for i in range(rows)] + col_tstop = [t + EXPOSURE_S for t in col_tstart] + if header_tstop is None: + # One EXPOSURE bin beyond the last bin's end — the relationship §2.7's H3 admits. + header_tstop = header_tstart + (col_tstop[-1] + EXPOSURE_S) / SECONDS_PER_DAY + if channel_map is None: + channel_map = [list(range(detchans))] * rows + if counts is None: + counts = np.tile(np.arange(detchans, dtype=float), (rows, 1)) + + table = fits.BinTableHDU.from_columns([ + fits.Column(name="SPEC_NUM", format="I", array=np.arange(rows, dtype=np.int16)), + fits.Column(name="CHANNEL", format=f"{detchans}J", + array=np.array(channel_map, dtype=np.int32)), + fits.Column(name="COUNTS", format=f"{detchans}D", unit="cts", array=counts), + fits.Column(name="STAT_ERR", format=f"{detchans}D", array=np.ones_like(counts)), + fits.Column(name="ROWID", format="12A", array=np.array(["r"] * rows)), + fits.Column(name="TSTART", format="D", unit="s", + array=np.array(col_tstart, dtype=float)), + fits.Column(name="TSTOP", format="D", unit="s", + array=np.array(col_tstop, dtype=float)), + fits.Column(name="EXPOSURE", format="D", unit="s", + array=np.full(rows, EXPOSURE_S)), + ], name="SPECTRUM") + + for key, value in {"CHANTYPE": chantype, "HDUCLAS3": hduclas3, "HDUCLAS4": "TYPE:II", + "DETCHANS": detchans, "TSTART": header_tstart, + "TSTOP": header_tstop}.items(): + table.header[key] = value + for key, value in (overrides or {}).items(): + if value is None: + del table.header[key] + else: + table.header[key] = value + + path.parent.mkdir(parents=True, exist_ok=True) + fits.HDUList([fits.PrimaryHDU(), table]).writeto(path, overwrite=True) + return path + + +def write_gti(path: Path, *, detector: str = "czt1", + intervals: list[tuple[float, float]] | None = None) -> Path: + """A `gti.fits` per §2.9 — LOWERCASE column names, unlike SoLEXS.""" + if intervals is None: + intervals = [(TSTART_MJD, TSTART_MJD + 0.5)] + table = fits.BinTableHDU.from_columns([ + fits.Column(name="tstart", format="D", + array=np.array([a for a, _ in intervals], dtype=float)), + fits.Column(name="tstop", format="D", + array=np.array([b for _, b in intervals], dtype=float)), + ], name=f"GTI_{detector.upper()}") + path.parent.mkdir(parents=True, exist_ok=True) + fits.HDUList([fits.PrimaryHDU(), table]).writeto(path, overwrite=True) + return path + + +def write_hk(path: Path, *, mjd: list[float] | None = None, + suninfov: list[int] | None = None, rows: int = 50, + header_tstart: float | None = None, + header_tstop: float | None = None, + units: dict[str, str | None] | None = None, + values: dict[str, list] | None = None, + omit: tuple[str, ...] = ()) -> Path: + """An `hk.fits` per §2.8 — HDU `HLSHK`, every decisive column, telemetry order preserved.""" + if mjd is None: + mjd = [TSTART_MJD + i / SECONDS_PER_DAY for i in range(rows)] + if suninfov is None: + suninfov = [1] * len(mjd) + + columns = [ + fits.Column(name="mjd", format="D", array=np.array(mjd, dtype=float)), + fits.Column(name="suninfov", format="I", + array=np.array(suninfov, dtype=np.int16)), + ] + for name, (form, unit) in HK_COLUMNS.items(): + if name in omit: + continue + unit = (units or {}).get(name, unit) + dtype = np.int64 if form == "K" else (np.int32 if form == "J" else float) + data = (values or {}).get(name, list(range(len(mjd)))) + columns.append(fits.Column(name=name, format=form, unit=unit, + array=np.array(data, dtype=dtype))) + + table = fits.BinTableHDU.from_columns(columns, name="HLSHK") + table.header["TSTART"] = header_tstart if header_tstart is not None else min(mjd) + table.header["TSTOP"] = header_tstop if header_tstop is not None else max(mjd) + + path.parent.mkdir(parents=True, exist_ok=True) + fits.HDUList([fits.PrimaryHDU(), table]).writeto(path, overwrite=True) + return path + + +def write_events(path: Path, *, rows: int = 10, + detectors: tuple[str, ...] = ("cdte1", "cdte2", "czt1", "czt2"), + detnam: dict[str, str] | None = None, + energies: list[float] | None = None, + mjd: list[float] | None = None, + isot: list[str] | None = None, + ener_unit: str | None = "keV", + omit: tuple[str, ...] = (), + header_tstart: float = TSTART_MJD, + header_tstop: float = TSTART_MJD + 0.5) -> Path: + """An `evt.fits` per §2.5 — four detector HDUs, `ener` in keV, CZT `pix`/`offsetchn`.""" + names = {"cdte1": "CDTE1-EVENTS", "cdte2": "CDTE2-EVENTS", + "czt1": "CZT1-EVENTS", "czt2": "CZT2-EVENTS"} + default_detnam = {"cdte1": "CdTe1", "cdte2": "CdTe2", "czt1": "CZT1", "czt2": "CZT2"} + detnam = {**default_detnam, **(detnam or {})} + if mjd is None: + mjd = [TSTART_MJD + i / SECONDS_PER_DAY for i in range(rows)] + rows = len(mjd) + if energies is None: + energies = [20.0 + i for i in range(rows)] + if isot is None: + isot = [isot_of(value) for value in mjd] + + hdus = [fits.PrimaryHDU()] + for detector in detectors: + columns = [ + fits.Column(name="mjd", format="D", array=np.array(mjd, dtype=float)), + fits.Column(name="hlsobt", format="D", unit="s", + array=np.arange(rows, dtype=float)), + fits.Column(name="currtemp", format="D", unit="degC", + array=np.full(rows, 12.0)), + fits.Column(name="chn", format="I", array=np.arange(rows, dtype=np.int16)), + fits.Column(name="ener", format="D", unit=ener_unit, + array=np.array(energies, dtype=float)), + fits.Column(name="recnum", format="J", array=np.arange(rows, dtype=np.int32)), + fits.Column(name="utc-isot", format="23A", array=np.array(isot)), + ] + if detector.startswith("czt"): + columns += [ + fits.Column(name="pix", format="B", array=np.arange(rows, dtype=np.uint8)), + fits.Column(name="offsetchn", format="I", + array=np.arange(rows, dtype=np.int16)), + ] + columns = [c for c in columns if c.name not in omit] + table = fits.BinTableHDU.from_columns(columns, name=names[detector]) + table.header["DETNAM"] = detnam[detector] + table.header["TSTART"] = header_tstart + table.header["TSTOP"] = header_tstop + hdus.append(table) + + path.parent.mkdir(parents=True, exist_ok=True) + fits.HDUList(hdus).writeto(path, overwrite=True) + return path diff --git a/contexts/ingest/tests/test_hel1os_parsers.py b/contexts/ingest/tests/test_hel1os_parsers.py new file mode 100644 index 0000000..75cf599 --- /dev/null +++ b/contexts/ingest/tests/test_hel1os_parsers.py @@ -0,0 +1,903 @@ +"""HEL1OS parsers — per-field vs spec, no imputation, and §4 version resolution. + +`| 18 | 18 | E5 | 900 | L | 17 | per-field vs spec | no imputation | parse fixtures |` + +Every assertion cites the clause of `SPEC-parsers@r7` it enforces. Each fail-loud rule is +fired against a product built to violate exactly it — the archive contains no malformed +HEL1OS product, so a violating fixture is the only way to watch a rule reject anything. + +r7 adds §5.1's time-representation allowance `ε_t`, so the comparisons that carry it are tested +from both sides: representation noise inside `ε_t` is admitted, and a disagreement beyond it +still terminates. §2.5's non-decreasing event rule is tested as written; the archive falsifies +it and the falsification is recorded OPEN in CONTRA-008 rather than absorbed here. +""" + +from __future__ import annotations + +import math +from pathlib import Path + +import numpy as np +import pytest + +from contexts.ingest.parsers.hel1os import events, gti, hk, lc, orbit, spectra +from contexts.ingest.tests import hel1os_fixtures as fx +from domain.errors import ContractViolation +from domain.values import Digest, Identifier, Timestamp + +DIGEST = Digest("a" * 64) + + +def violation(caught) -> str: + return caught.value.message + + +# ═══════════════════════════════════════════ §4 — orbit identity and precedence + + +def test_the_orbit_stem_is_parsed_by_the_specified_regex(): + parsed = orbit.parse_stem("HLS_20251208_000008_43178sec_lev1_V111") + assert (parsed.date, parsed.start, parsed.duration_s, parsed.version) == ( + "20251208", "000008", 43178, 111) + assert parsed.stem == "HLS_20251208_000008_43178sec_lev1_V111" + + +@pytest.mark.parametrize("bad", [ + "HLS_20251208_000008_43178sec_lev1", "AL1_SLX_L1_20240514_v1.0", + "HLS_20251208_000008_43178sec_lev1_V11", "", None, +]) +def test_a_name_that_is_not_an_orbit_stem_is_refused(bad): + with pytest.raises(ContractViolation) as caught: + orbit.parse_stem(bad) + assert violation(caught).startswith("F-18") + + +def test_rule_1_higher_version_wins(): + """§4 precedence 1. Class A: identical interval, different version.""" + a = orbit.parse_stem("HLS_20251208_000008_43178sec_lev1_V111") + b = orbit.parse_stem("HLS_20251208_000008_43178sec_lev1_V211") + assert orbit.precedence(a, b) is b + assert orbit.precedence(b, a) is b + assert orbit.rule_applied(a, b) == "rule-1-higher-version" + + +def test_the_version_digits_are_opaque_not_decomposed(): + """§8 A-1: the three digits are undocumented and compared as one integer. + + V112 vs V211: a major/minor reading would make 1.1.2 lose to 2.1.1 for a reason nobody + has authority for. As integers 112 < 211, which is what §4 specifies — and the point is + that the rule is stated, not inferred. + """ + low = orbit.parse_stem("HLS_20251208_000008_43178sec_lev1_V112") + high = orbit.parse_stem("HLS_20251208_000008_43178sec_lev1_V211") + assert orbit.precedence(low, high) is high + + +def test_rule_2_longer_duration_wins_on_a_version_tie(): + a = orbit.parse_stem("HLS_20251207_120003_43195sec_lev1_V111") + b = orbit.parse_stem("HLS_20251207_120003_42570sec_lev1_V111") + assert orbit.precedence(a, b) is a + assert orbit.rule_applied(a, b) == "rule-2-longer-duration" + + +def test_rule_3_later_processing_date_wins(): + base = orbit.parse_stem("HLS_20251208_000008_43178sec_lev1_V111") + early = orbit.OrbitId(base.date, base.start, base.duration_s, base.version, + Timestamp("2025-12-10T00:00:00Z")) + late = orbit.OrbitId(base.date, base.start, base.duration_s, base.version, + Timestamp("2025-12-20T00:00:00Z")) + assert orbit.precedence(early, late) is late + assert orbit.rule_applied(early, late) == "rule-3-later-processing-date" + + +def test_rule_4_terminates_and_never_coin_flips(): + """§4: "Still tied → F-14 terminate. Never coin-flip." + + The most important rule in §4, because the alternative is a silent, arbitrary choice + between two products that claim the same coverage. + """ + a = orbit.parse_stem("HLS_20251208_000008_43178sec_lev1_V111") + b = orbit.parse_stem("HLS_20251208_000008_43178sec_lev1_V111") + with pytest.raises(ContractViolation) as caught: + orbit.precedence(a, b) + assert violation(caught).startswith("F-14") + assert "never coin-flip" in violation(caught) + assert orbit.rule_applied(a, b) == "rule-4-terminate" + + +def test_an_equal_processing_date_still_terminates(): + base = orbit.parse_stem("HLS_20251208_000008_43178sec_lev1_V111") + same = Timestamp("2025-12-20T00:00:00Z") + with pytest.raises(ContractViolation): + orbit.precedence( + orbit.OrbitId(base.date, base.start, base.duration_s, base.version, same), + orbit.OrbitId(base.date, base.start, base.duration_s, base.version, same), + ) + + +def test_class_b_partial_overlap_is_detected(): + """§4: Class B — different start and duration, each covering seconds the other lacks. + + Detected, not resolved: resolving it needs the minute-level coverage map, which §4 ties + to emitting T3/T4/T5 and which belongs to M3/E5/#19. + """ + a = orbit.parse_stem("HLS_20251207_120003_43195sec_lev1_V211") + b = orbit.parse_stem("HLS_20251207_121028_42570sec_lev1_V111") + assert orbit.overlaps(a, b) + assert orbit.overlaps(b, a) + + +def test_non_overlapping_orbits_are_not_reported_as_overlapping(): + a = orbit.parse_stem("HLS_20251207_000006_43190sec_lev1_V111") + b = orbit.parse_stem("HLS_20251208_000008_43178sec_lev1_V111") + assert not orbit.overlaps(a, b) + + +def test_no_orbit_merge_api_exists(): + """§4.4: "There is no API that concatenates orbit files directly." + + The structural guarantee that naive ingestion cannot occur. #18 supplies the precedence + rules; the merge that consumes them must take the coverage map as a required argument, + and it belongs to #19. Shipping a merge here without that map would create precisely the + API §4 forbids. + """ + for forbidden in ("merge", "concat", "concatenate", "combine", "coverage_map", + "resolve_all", "ingest"): + assert not hasattr(orbit, forbidden), f"orbit.py grew {forbidden!r}" + + +@pytest.mark.parametrize("detector, family", [ + ("czt1", "czt"), ("czt2", "czt"), ("cdte1", "cdte"), ("cdte2", "cdte"), +]) +def test_each_detector_maps_to_its_family(detector, family): + assert orbit.family_of(detector) == family + + +@pytest.mark.parametrize("bad", ["czt3", "sdd2", "solexs", "cdte", ""]) +def test_an_unknown_detector_has_no_family(bad): + with pytest.raises(ContractViolation) as caught: + orbit.family_of(bad) + assert violation(caught).startswith("F-07") + + +def test_the_detector_is_part_of_the_instrument_identity(): + parsed = orbit.parse_stem("HLS_20251208_000008_43178sec_lev1_V111") + assert parsed.detector_id("CZT1") == Identifier("hel1os-czt1") + assert parsed.detector_id("cdte2") == Identifier("hel1os-cdte2") + + +# ═══════════════════════════════════════════ §2.6 — band light curves + + +def test_the_bands_are_parsed_from_extname(tmp_path): + product = lc.parse(fx.write_lightcurve(tmp_path / "lc.fits"), DIGEST, detector="czt1") + assert [(b.low_kev, b.high_kev) for b in product.bands] == list(fx.CZT_BANDS) + assert product.bands[0].extname == "CZT1_LC_BAND_20.00KEV_TO_40.00KEV" + + +def test_the_cdte_family_has_its_own_bands(tmp_path): + product = lc.parse(fx.write_lightcurve(tmp_path / "lc.fits", detector="cdte1"), + DIGEST, detector="cdte1") + assert [(b.low_kev, b.high_kev) for b in product.bands] == list(fx.CDTE_BANDS) + + +def test_an_unknown_band_is_never_silently_accepted(tmp_path): + """F-10. §2.6: an unlisted band would attribute a rate to the wrong energy range.""" + bands = ((20.0, 40.0), (40.0, 60.0), (60.0, 80.0), (80.0, 150.0), (200.0, 400.0)) + with pytest.raises(ContractViolation) as caught: + lc.parse(fx.write_lightcurve(tmp_path / "lc.fits", bands=bands), DIGEST, + detector="czt1") + assert violation(caught).startswith("F-10") + + +def test_a_cdte_band_in_a_czt_file_is_refused(tmp_path): + """The families' allowlists are separate; 5-20 keV is CdTe's, not CZT's.""" + with pytest.raises(ContractViolation) as caught: + lc.parse(fx.write_lightcurve(tmp_path / "lc.fits", bands=fx.CDTE_BANDS), + DIGEST, detector="czt1") + assert violation(caught).startswith("F-10") + + +def test_an_extname_that_encodes_no_band_is_refused(tmp_path): + names = ["CZT1_LIGHTCURVE"] + [fx.band_extname("czt1", lo, hi) + for lo, hi in fx.CZT_BANDS[1:]] + with pytest.raises(ContractViolation) as caught: + lc.parse(fx.write_lightcurve(tmp_path / "lc.fits", extnames=names), DIGEST, + detector="czt1") + assert violation(caught).startswith("F-10") + + +def test_the_wrong_number_of_band_hdus_is_refused(tmp_path): + with pytest.raises(ContractViolation) as caught: + lc.parse(fx.write_lightcurve(tmp_path / "lc.fits", bands=fx.CZT_BANDS[:4]), + DIGEST, detector="czt1") + assert violation(caught).startswith("F-10") + + +def test_the_total_band_is_identified_not_assumed_by_position(tmp_path): + """HDU order is not a contract, so the total is recognised by its edges.""" + product = lc.parse(fx.write_lightcurve(tmp_path / "lc.fits"), DIGEST, detector="czt1") + totals = [b for b in product.bands if b.is_total] + assert len(totals) == 1 + assert (totals[0].low_kev, totals[0].high_kev) == (18.0, 160.0) + + +def test_the_rate_unit_is_declared_and_differs_from_solexs(tmp_path): + """F-07: `CTR` is a rate in cts/sec; SoLEXS `.lc` is undeclared counts. Never shared.""" + product = lc.parse(fx.write_lightcurve(tmp_path / "lc.fits"), DIGEST, detector="czt1") + observed = list(product.observations(source_id=Identifier("issdc-pradan"), + ingest_time=None)) + assert {o.unit for o in observed} == {"cts/s"} + + from contexts.ingest.parsers.solexs import lc as solexs_lc + assert solexs_lc.UNIT == "counts" + assert solexs_lc.UNIT != "cts/s" + + +def test_the_ctr_unit_must_be_declared_as_a_rate(tmp_path): + """§2.6: `CTR` declares `cts/sec`; the parser reads that rather than assuming it (F-07).""" + with pytest.raises(ContractViolation) as caught: + lc.parse(fx.write_lightcurve(tmp_path / "lc.fits", ctr_unit="counts"), DIGEST, + detector="czt1") + assert violation(caught).startswith("F-07") + + +def test_an_undeclared_ctr_unit_is_refused(tmp_path): + with pytest.raises(ContractViolation) as caught: + lc.parse(fx.write_lightcurve(tmp_path / "lc.fits", ctr_unit=None), DIGEST, + detector="czt1") + assert violation(caught).startswith("F-07") + + +def test_a_negative_rate_is_physically_impossible(tmp_path): + rates = {(20.0, 40.0): [1.0, -2.0] + [3.0] * 118} + with pytest.raises(ContractViolation) as caught: + lc.parse(fx.write_lightcurve(tmp_path / "lc.fits", rates=rates), DIGEST, + detector="czt1") + assert violation(caught).startswith("F-19") + + +def test_a_nan_rate_is_refused_rather_than_read_as_absent(tmp_path): + """NO IMPUTATION, and no imported convention either. + + SoLEXS §2.1 declares NaN as its missing-data sentinel. §2.6 declares none for `CTR`. + Treating a NaN here as "observed to be absent" would import one instrument's convention + into another, which is exactly what F-07 forbids — so it terminates instead. + """ + rates = {(20.0, 40.0): [1.0, math.nan] + [3.0] * 118} + with pytest.raises(ContractViolation) as caught: + lc.parse(fx.write_lightcurve(tmp_path / "lc.fits", rates=rates), DIGEST, + detector="czt1") + assert violation(caught).startswith("F-07") + assert "across instruments" in violation(caught) + + +def test_a_non_increasing_time_axis_is_refused(tmp_path): + """§2.6: MJD strictly increasing — unlike §2.8 housekeeping, where r4 removed it.""" + mjd = [fx.TSTART_MJD + i / 86400.0 for i in range(120)] + mjd[10] = mjd[9] + with pytest.raises(ContractViolation) as caught: + lc.parse(fx.write_lightcurve(tmp_path / "lc.fits", mjd=mjd), DIGEST, + detector="czt1") + assert violation(caught).startswith("F-16") + + +def test_the_mjd_epoch_is_converted_exactly(tmp_path): + """MJD 40587 = the Unix epoch. §2.5 anchors TSTART 61017.0000988685 to 2025-12-08.""" + assert lc.MJD_UNIX_EPOCH == 40587.0 + assert str(lc.mjd_to_timestamp(fx.TSTART_MJD)).startswith("2025-12-08T00:00:08") + + +# ═══════════════════════════════════════════ §2.7 — spectra and R-1 + + +@pytest.mark.parametrize("detector, detchans", [ + ("czt1", 341), ("czt2", 341), ("cdte1", 511), ("cdte2", 511), +]) +def test_detchans_is_validated_against_the_family_allowlist(tmp_path, detector, detchans): + """§2.7 r5: never a single scalar.""" + product = spectra.parse( + fx.write_spectra(tmp_path / "s.fits", detector=detector), DIGEST, detector=detector) + assert product.header.detchans == detchans + assert spectra.DETCHANS_BY_FAMILY[product.header.family] == detchans + + +def test_a_czt_file_declaring_the_cdte_channel_count_is_refused(tmp_path): + """An unlisted (family, DETCHANS) pair terminates via F-07.""" + with pytest.raises(ContractViolation) as caught: + spectra.parse(fx.write_spectra(tmp_path / "s.fits", detector="czt1", detchans=511), + DIGEST, detector="czt1") + assert violation(caught).startswith("F-07") + + +@pytest.mark.parametrize("detchans", [340, 342, 512]) +def test_any_other_channel_count_is_refused(tmp_path, detchans): + """340 is SoLEXS's PI space — the nearest miss, and the one F-11 names.""" + with pytest.raises(ContractViolation) as caught: + spectra.parse(fx.write_spectra(tmp_path / "s.fits", detchans=detchans), + DIGEST, detector="czt1") + assert violation(caught).startswith("F-07") + + +def test_the_three_channel_spaces_stay_distinct(): + """F-11: SoLEXS PI(340), CZT PHA(341), CdTe PHA(511) are incommensurable.""" + spaces = spectra.INCOMMENSURABLE_CHANNEL_SPACES + assert spaces[("solexs", "PI")] == 340 + assert spaces[("czt", "PHA")] == 341 + assert spaces[("cdte", "PHA")] == 511 + assert len(set(spaces.values())) == 3 + + from contexts.ingest.parsers.solexs import pi as solexs_pi + assert solexs_pi.DETCHANS == 340 + assert solexs_pi.DETCHANS not in spectra.DETCHANS_BY_FAMILY.values() + + +def test_the_chantype_is_pha_not_pi(tmp_path): + """PI is gain-corrected, PHA is raw pulse height. Not the same space.""" + with pytest.raises(ContractViolation) as caught: + spectra.parse(fx.write_spectra(tmp_path / "s.fits", chantype="PI"), + DIGEST, detector="czt1") + assert violation(caught).startswith("F-07") + + +def test_hduclas3_is_count_singular(tmp_path): + """§2.7 `OBSERVED`: 'COUNT', where SoLEXS §2.2 declares 'COUNTS'. Read, not assumed.""" + assert spectra.HDUCLAS3 == "COUNT" + with pytest.raises(ContractViolation) as caught: + spectra.parse(fx.write_spectra(tmp_path / "s.fits", hduclas3="COUNTS"), + DIGEST, detector="czt1") + assert violation(caught).startswith("F-07") + + +def test_r1_resolves_h3_on_the_observed_convention(tmp_path): + """§2.7 r4: relative seconds from the header epoch, tested first because unit='s'.""" + product = spectra.parse(fx.write_spectra(tmp_path / "s.fits"), DIGEST, detector="czt1") + assert product.epoch.hypothesis == "H3" + assert product.epoch.exposure_s == 20.0 + + +def test_r1_records_its_residual(tmp_path): + """§2.7: "The resolved hypothesis and its residual are recorded in T7 provenance." + + The recording is #19's; the measurement is this parser's, and it must be available. + """ + product = spectra.parse(fx.write_spectra(tmp_path / "s.fits"), DIGEST, detector="czt1") + assert product.epoch.residual_s >= 0 + assert math.isfinite(product.epoch.residual_s) + + +def test_r1_falls_through_to_h1_for_mjd_columns(tmp_path): + """H1 retained because a reprocessed product could switch to an absolute epoch.""" + columns = [fx.TSTART_MJD + i / 86400.0 for i in range(20)] + product = spectra.parse( + fx.write_spectra(tmp_path / "s.fits", col_tstart=columns), DIGEST, detector="czt1") + assert product.epoch.hypothesis == "H1" + + +def test_r1_falls_through_to_h2_for_unix_seconds(tmp_path): + """H2 retained for the same reason as H1.""" + unix0 = (fx.TSTART_MJD - 40587.0) * 86400.0 + columns = [unix0 + i * 20.0 for i in range(20)] + product = spectra.parse( + fx.write_spectra(tmp_path / "s.fits", col_tstart=columns), DIGEST, detector="czt1") + assert product.epoch.hypothesis == "H2" + + +def test_r1_terminates_when_every_hypothesis_fails(tmp_path): + """F-06: §5 admits no ambiguous time.""" + columns = [500_000.0 + i * 20.0 for i in range(20)] + with pytest.raises(ContractViolation) as caught: + spectra.parse(fx.write_spectra(tmp_path / "s.fits", col_tstart=columns), + DIGEST, detector="czt1") + assert violation(caught).startswith("F-06") + assert "R-1" in violation(caught) + + +def test_h3_requires_the_first_offset_to_be_exactly_zero(tmp_path): + """§2.7: "col[0] == 0 exactly". A near-zero offset is not H3, and falls through.""" + columns = [0.5 + i * 20.0 for i in range(20)] + with pytest.raises(ContractViolation) as caught: + spectra.parse(fx.write_spectra(tmp_path / "s.fits", col_tstart=columns), + DIGEST, detector="czt1") + assert violation(caught).startswith("F-06") + + +def test_a_two_bin_span_disagreement_still_fails_h3(tmp_path): + """The MJD-precision allowance covers 79 ns of float error, not a real disagreement. + + A header span two bins longer than the columns cover is a genuine mismatch and must not + be absorbed. It falls through H3, then H1 and H2, and terminates. + """ + header_tstop = fx.TSTART_MJD + (20 * 20.0 + 3 * 20.0) / 86400.0 + with pytest.raises(ContractViolation) as caught: + spectra.parse(fx.write_spectra(tmp_path / "s.fits", header_tstop=header_tstop), + DIGEST, detector="czt1") + assert violation(caught).startswith("F-06") + + +def test_a_varying_channel_map_is_refused(tmp_path): + """F-08, same rule as SoLEXS §2.2 but checked separately — sharing it would share a space.""" + varying = [list(range(341))] * 19 + [list(range(1, 342))] + with pytest.raises(ContractViolation) as caught: + spectra.parse(fx.write_spectra(tmp_path / "s.fits", channel_map=varying), + DIGEST, detector="czt1") + assert violation(caught).startswith("F-08") + + +def test_the_spectra_stream_rather_than_materialise(tmp_path): + import inspect + assert inspect.isgeneratorfunction(spectra.Spectra.spectra) + product = spectra.parse(fx.write_spectra(tmp_path / "s.fits"), DIGEST, detector="czt1") + first = next(iter(product.spectra())) + assert len(first.counts) == 341 + assert first.exposure == 20.0 + + +def test_a_negative_spectral_count_is_refused(tmp_path): + counts = np.tile(np.arange(341, dtype=float), (20, 1)) + counts[3][7] = -1.0 + with pytest.raises(ContractViolation) as caught: + list(spectra.parse(fx.write_spectra(tmp_path / "s.fits", counts=counts), + DIGEST, detector="czt1").spectra()) + assert violation(caught).startswith("F-19") + + +# ═══════════════════════════════════════════ §2.8 — housekeeping + + +def test_archive_order_is_preserved_exactly(tmp_path): + """§2.8 r4, binding: the parser performs NO sorting. + + "The parser is a lossless representation of the archive — reading and transforming are + separate acts, and a parser that silently reorders is no longer a faithful reader." + """ + jittered = [fx.TSTART_MJD + i / 86400.0 for i in range(20)] + jittered[5], jittered[6] = jittered[6], jittered[5] + + product = hk.parse(fx.write_hk(tmp_path / "hk.fits", mjd=jittered), DIGEST) + assert list(product.mjd) == jittered + assert list(product.mjd) != sorted(jittered) + + +def test_no_sorting_function_exists(): + for forbidden in ("sort", "sorted", "reorder", "chronological"): + assert not any(forbidden in name.lower() for name in dir(hk)), ( + f"hk.py exposes {forbidden!r}; §2.8 places chronological_sort outside the " + f"parser layer" + ) + + +def test_an_inversion_is_recorded_not_rejected(tmp_path): + """§2.8 r4: the non-decreasing requirement was falsified by the archive and removed.""" + jittered = [fx.TSTART_MJD + i / 86400.0 for i in range(20)] + jittered[5], jittered[6] = jittered[6], jittered[5] + + product = hk.parse(fx.write_hk(tmp_path / "hk.fits", mjd=jittered), DIGEST) + assert product.inversions.n_out_of_order == 1 + assert product.inversions.max_backward_step_s > 0 + + +def test_no_jitter_threshold_is_defined(): + """§2.8: "The magnitude is reported, never compared against an invented tolerance." + + So `InversionStats` has no limit, no `is_acceptable`, and nothing to compare against. + """ + stats = hk.InversionStats(rows=10, steps=9, n_out_of_order=3, max_backward_step_s=0.9) + public = [name for name in dir(stats) if not name.startswith("_")] + for forbidden in ("threshold", "limit", "tolerance", "acceptable", "passes", "within"): + assert not any(forbidden in name.lower() for name in public), ( + f"InversionStats exposes {forbidden!r}; §2.8 reports the magnitude and defines " + f"no threshold to compare it against" + ) + assert set(public) == { + "max_backward_step_ms", "max_backward_step_s", "n_out_of_order", "rows", "steps" + } + + +def test_a_duplicate_timestamp_is_a_defect(tmp_path): + """§2.8 r4 kept this one: an inversion is jitter, a duplicate is a genuine defect.""" + duplicated = [fx.TSTART_MJD + i / 86400.0 for i in range(20)] + duplicated[7] = duplicated[6] + with pytest.raises(ContractViolation) as caught: + hk.parse(fx.write_hk(tmp_path / "hk.fits", mjd=duplicated), DIGEST) + assert violation(caught).startswith("F-16") + assert "duplicate" in violation(caught) + + +def test_the_header_span_must_contain_the_measurements(tmp_path): + """§2.8 r7 assigns this check F-06; the bound is widened by ε_t and by nothing else.""" + with pytest.raises(ContractViolation) as caught: + hk.parse(fx.write_hk(tmp_path / "hk.fits", header_tstop=fx.TSTART_MJD), DIGEST) + assert violation(caught).startswith("F-06") + + +def test_suninfov_is_a_first_class_flag(tmp_path): + """§2.8 binding: data outside Sun-in-FOV is not solar signal.""" + flags = [1] * 10 + [0] * 10 + mjd = [fx.TSTART_MJD + i / 86400.0 for i in range(len(flags))] + product = hk.parse( + fx.write_hk(tmp_path / "hk.fits", mjd=mjd, suninfov=flags), DIGEST) + assert product.sun_in_fov(0) is True + assert product.sun_in_fov(15) is False + assert list(product.suninfov) == flags + + +@pytest.mark.parametrize("bad", [2, -1, 7]) +def test_a_suninfov_outside_the_declared_domain_is_refused(tmp_path, bad): + """Coercing it would decide whether data is solar signal on the parser's authority.""" + flags = [1] * 19 + [bad] + mjd = [fx.TSTART_MJD + i / 86400.0 for i in range(len(flags))] + with pytest.raises(ContractViolation) as caught: + hk.parse(fx.write_hk(tmp_path / "hk.fits", mjd=mjd, suninfov=flags), DIGEST) + assert violation(caught).startswith("F-07") + + +def test_the_czt2enth_unit_assumption_is_recorded(tmp_path): + """§8 A-4: the assumption travels with the data rather than living in a docstring.""" + product = hk.parse(fx.write_hk(tmp_path / "hk.fits"), DIGEST) + assert product.assumptions + assert any("czt2enth" in a for a in product.assumptions) + + +def test_the_phase_1a_columns_are_captured(tmp_path): + """§2.8's decisive columns: pile-up, saturation, thermal, pointing.""" + product = hk.parse(fx.write_hk(tmp_path / "hk.fits"), DIGEST) + for name in hk.PILEUP + hk.SATURATION + hk.THERMAL: + assert name in product.columns + + +# ═══════════════════════════════════════════ §2.9 — GTI + + +def test_columns_are_read_case_insensitively(tmp_path): + """§2.9 / §8 A-2: lowercase `tstart`/`tstop`, where SoLEXS uses uppercase START/STOP. + + A case-sensitive lookup would fail as F-04 and look like a missing column rather than a + spelling difference between instruments. + """ + product = gti.parse(fx.write_gti(tmp_path / "g.fits"), DIGEST, detector="czt1") + assert len(product.intervals) == 1 + assert product.detector_active + + from contexts.ingest.parsers.solexs import gti as solexs_gti + assert solexs_gti is not gti + + +def test_an_empty_gti_is_legal(tmp_path): + """F-12, the single deliberate non-terminating rule.""" + product = gti.parse(fx.write_gti(tmp_path / "g.fits", intervals=[]), DIGEST, + detector="czt1") + assert product.detector_active is False + assert product.intervals == () + + +def test_no_exposure_equality_rule_is_imported_from_solexs(tmp_path): + """§2.3's Σ(STOP−START+1) == EXPOSURE rests on 1 s sampling and a declared EXPOSURE. + + §2.9 declares neither, so importing the rule would apply one instrument's convention to + another — F-07's whole subject. + """ + product = gti.parse(fx.write_gti(tmp_path / "g.fits"), DIGEST, detector="czt1") + assert product.live_time_s > 0 + assert not hasattr(product, "declared_exposure") + + +def test_a_bound_that_is_not_an_mjd_is_refused(tmp_path): + """§2.9 leaves the unit undeclared; a bound outside the plausible range is not + reinterpreted as another epoch.""" + with pytest.raises(ContractViolation) as caught: + gti.parse(fx.write_gti(tmp_path / "g.fits", intervals=[(1.7e9, 1.7e9 + 100)]), + DIGEST, detector="czt1") + assert violation(caught).startswith("F-05") + + +def test_an_inverted_interval_is_refused(tmp_path): + with pytest.raises(ContractViolation) as caught: + gti.parse(fx.write_gti(tmp_path / "g.fits", + intervals=[(fx.TSTART_MJD + 0.5, fx.TSTART_MJD)]), + DIGEST, detector="czt1") + assert violation(caught).startswith("F-09") + + +# ═══════════════════════════════════════════ §2.5 — events + + +def test_all_four_detector_hdus_are_required(tmp_path): + """F-03: three detectors' events are a different measurement, not a smaller one.""" + with pytest.raises(ContractViolation) as caught: + events.parse(fx.write_events(tmp_path / "e.fits", + detectors=("czt1", "czt2", "cdte1")), DIGEST) + assert violation(caught).startswith("F-03") + + +def test_detnam_must_match_the_hdu_it_labels(tmp_path): + with pytest.raises(ContractViolation) as caught: + events.parse(fx.write_events(tmp_path / "e.fits", detnam={"czt1": "CZT2"}), DIGEST) + assert violation(caught).startswith("F-03") + + +def test_the_event_reader_emits_no_observations(tmp_path): + """§2.5, binding: the parser MUST expose an event reader but MUST NOT ingest events + into the canonical minute tables. + + Enforced by absence — there is no `observations()` to call, so events cannot flow into + #19's tables by default. + """ + product = events.parse(fx.write_events(tmp_path / "e.fits"), DIGEST) + assert not hasattr(product, "observations") + assert not hasattr(events, "observations") + for detector in product.detectors: + assert not hasattr(detector, "observations") + + +def test_event_energy_is_already_calibrated_in_kev(tmp_path): + """§2.5: unlike SoLEXS, HEL1OS ships keV. Both facts are true and must not be merged.""" + product = events.parse(fx.write_events(tmp_path / "e.fits"), DIGEST) + first = next(iter(product.events("czt1"))) + assert first.energy_kev == 20.0 + assert events.ENERGY_UNIT == "keV" + + from contexts.ingest.parsers.solexs import pi as solexs_pi + assert not hasattr(solexs_pi, "ENERGY_UNIT") + + +def test_a_non_positive_event_energy_is_refused(tmp_path): + """§2.5 validates `ener > 0`.""" + energies = [20.0] * 9 + [0.0] + with pytest.raises(ContractViolation) as caught: + list(events.parse(fx.write_events(tmp_path / "e.fits", energies=energies), + DIGEST).events("czt1")) + assert violation(caught).startswith("F-19") + + +def test_events_stream_rather_than_materialise(tmp_path): + import inspect + assert inspect.isgeneratorfunction(events.EventList.events) + + +# ═══════════════════════════════════════════ §5.1 — the time-representation allowance + + +def test_the_allowance_is_the_contract_value(): + """§5.1: ε_t = 1 ms, and it is the contract's value rather than each module's. + + hk, spectra and events all compare an MJD-derived time against a bound; each imports the + one constant, so no module can quietly hold a tolerance of its own. + """ + assert lc.TIME_REPRESENTATION_ALLOWANCE_S == 1e-3 + for module in (spectra, hk, events): + source = Path(module.__file__).read_text() + assert "TIME_REPRESENTATION_ALLOWANCE_S" in source + assert "math.ulp" not in source + + +def test_h3_admits_representation_error_within_the_allowance(tmp_path): + """§2.7 r7: the one-bin bound carries ε_t, so float64 noise does not reject a valid product.""" + header_tstop = fx.TSTART_MJD + (20 * 20.0 + 20.0 + 0.4e-3) / 86400.0 + product = spectra.parse(fx.write_spectra(tmp_path / "s.fits", header_tstop=header_tstop), + DIGEST, detector="czt1") + assert product.epoch.hypothesis == "H3" + + +def test_h3_refuses_a_disagreement_beyond_the_allowance(tmp_path): + """ε_t is representation, not physics: 3 ms past the bin bound is a disagreement.""" + header_tstop = fx.TSTART_MJD + (20 * 20.0 + 20.0 + 3e-3) / 86400.0 + with pytest.raises(ContractViolation) as caught: + spectra.parse(fx.write_spectra(tmp_path / "s.fits", header_tstop=header_tstop), + DIGEST, detector="czt1") + assert violation(caught).startswith("F-06") + + +def test_the_hk_header_span_admits_representation_error(tmp_path): + """§2.8 r7: a bound a fraction of a millisecond inside the data is representation noise.""" + mjd = [fx.TSTART_MJD + i / 86400.0 for i in range(20)] + product = hk.parse( + fx.write_hk(tmp_path / "hk.fits", mjd=mjd, + header_tstop=max(mjd) - 0.4e-3 / 86400.0), DIGEST) + assert product.inversions.rows == 20 + + +def test_the_hk_header_span_refuses_a_real_excursion(tmp_path): + mjd = [fx.TSTART_MJD + i / 86400.0 for i in range(20)] + with pytest.raises(ContractViolation) as caught: + hk.parse(fx.write_hk(tmp_path / "hk.fits", mjd=mjd, + header_tstop=max(mjd) - 3e-3 / 86400.0), DIGEST) + assert violation(caught).startswith("F-06") + + +# ═══════════════════════════════════════════ §2.8 — every decisive column, with its unit + + +def test_every_decisive_column_is_captured_with_its_archive_dtype(tmp_path): + """§2.8's decisive columns, by the archive's spelling (CONTRA-007 Observation E).""" + product = hk.parse(fx.write_hk(tmp_path / "hk.fits"), DIGEST) + for name in hk.CAPTURED: + assert name in product.columns, name + for name in ("czt1hotpixcnt", "czt2bunpxctr", "fehkstat", "cdte1pilectr"): + assert isinstance(product.columns[name][0], int) + assert isinstance(product.columns["czthvmon"][0], float) + + +def test_the_stated_units_are_checked_not_assumed(tmp_path): + product = hk.parse(fx.write_hk(tmp_path / "hk.fits"), DIGEST) + for name, stated in hk.STATED_UNITS.items(): + assert product.units[name] == stated + assert product.units["cdte1enerthr"] is None + + +def test_a_declared_unit_contradicting_section_2_8_is_refused(tmp_path): + with pytest.raises(ContractViolation) as caught: + hk.parse(fx.write_hk(tmp_path / "hk.fits", units={"czthvmon": "mV"}), DIGEST) + assert violation(caught).startswith("F-07") + + +def test_czt2enth_takes_the_czt1enth_unit_under_a4(tmp_path): + """§8 A-4: undeclared, so czt1enth's keV is applied and the assumption travels.""" + product = hk.parse(fx.write_hk(tmp_path / "hk.fits"), DIGEST) + assert product.units["czt2enth"] == "keV" + assert any("czt2enth" in assumption for assumption in product.assumptions) + + +def test_a_contradictory_czt2enth_unit_is_refused(tmp_path): + """A-4 covers an undeclared unit, never two disagreeing declarations of one quantity.""" + with pytest.raises(ContractViolation) as caught: + hk.parse(fx.write_hk(tmp_path / "hk.fits", units={"czt2enth": "meV"}), DIGEST) + assert violation(caught).startswith("F-07") + + +def test_a_missing_decisive_column_is_refused(tmp_path): + with pytest.raises(ContractViolation) as caught: + hk.parse(fx.write_hk(tmp_path / "hk.fits", omit=("fehkstat",)), DIGEST) + assert violation(caught).startswith("F-04") + + +def test_only_the_two_czt_temperatures_must_be_finite(tmp_path): + """§2.8 names `czt1temp`/`czt2temp`. The CdTe temperatures are carried, NaN included.""" + rows = 20 + product = hk.parse( + fx.write_hk(tmp_path / "hk.fits", rows=rows, + values={"cdte1temp": [math.nan] * rows}), DIGEST) + assert math.isnan(product.columns["cdte1temp"][0]) + + with pytest.raises(ContractViolation) as caught: + hk.parse(fx.write_hk(tmp_path / "hk.fits", rows=rows, + values={"czt1temp": [math.nan] * rows}), DIGEST) + assert violation(caught).startswith("F-16") + + +# ═══════════════════════════════════════════ §2.5 r7 — the event row rules + + +def test_event_energy_must_be_declared_in_kev(tmp_path): + """§2.5: HEL1OS ships calibrated keV; the declaration is read, not assumed (F-07).""" + with pytest.raises(ContractViolation) as caught: + events.parse(fx.write_events(tmp_path / "e.fits", ener_unit="MeV"), DIGEST) + assert violation(caught).startswith("F-07") + + +def test_czt_event_hdus_require_pix_and_offsetchn(tmp_path): + with pytest.raises(ContractViolation) as caught: + events.parse(fx.write_events(tmp_path / "e.fits", omit=("pix",)), DIGEST) + assert violation(caught).startswith("F-04") + + +def test_czt_events_expose_the_pixel_and_cdte_events_do_not(tmp_path): + """§2.5 lists `pix`/`offsetchn` for CZT only; a reader that dropped them would be lossy.""" + product = events.parse(fx.write_events(tmp_path / "e.fits"), DIGEST) + czt = next(iter(product.events("czt1"))) + cdte = next(iter(product.events("cdte1"))) + assert czt.pixel is not None and czt.offset_channel is not None + assert cdte.pixel is None and cdte.offset_channel is None + + +def test_an_event_outside_the_header_span_is_refused(tmp_path): + """§2.5 r7: the span check, widened by ε_t and by nothing else.""" + beyond = [fx.TSTART_MJD + 0.6] + with pytest.raises(ContractViolation) as caught: + list(events.parse(fx.write_events(tmp_path / "e.fits", mjd=beyond), + DIGEST).events("czt1")) + assert violation(caught).startswith("F-06") + + +def test_v_evt_2_accepts_the_millisecond_rounding_the_archive_writes(tmp_path): + """`utc-isot` is the mjd instant to the millisecond, so it disagrees by ≤ 0.5 ms.""" + product = events.parse(fx.write_events(tmp_path / "e.fits"), DIGEST) + assert sum(1 for _ in product.events("czt1")) == 10 + + +def test_v_evt_2_refuses_a_disagreeing_isot(tmp_path): + mjd = [fx.TSTART_MJD + i / 86400.0 for i in range(5)] + drifted = [fx.isot_of(value + 0.05 / 86400.0) for value in mjd] + with pytest.raises(ContractViolation) as caught: + list(events.parse(fx.write_events(tmp_path / "e.fits", mjd=mjd, isot=drifted), + DIGEST).events("czt1")) + assert violation(caught).startswith("F-06") + assert "V-EVT-2" in violation(caught) + + +def test_an_unparseable_isot_is_refused(tmp_path): + mjd = [fx.TSTART_MJD + i / 86400.0 for i in range(3)] + with pytest.raises(ContractViolation) as caught: + list(events.parse(fx.write_events(tmp_path / "e.fits", mjd=mjd, + isot=["not-an-instant"] * 3), + DIGEST).events("czt1")) + assert violation(caught).startswith("F-06") + + +def test_a_decreasing_event_timestamp_is_refused(tmp_path): + """§2.5 as written, enforced although the archive falsifies it (CONTRA-008, OPEN).""" + mjd = [fx.TSTART_MJD + i / 86400.0 for i in range(6)] + mjd[4] = mjd[2] + with pytest.raises(ContractViolation) as caught: + list(events.parse(fx.write_events(tmp_path / "e.fits", mjd=mjd), DIGEST).events("czt1")) + assert violation(caught).startswith("F-16") + assert "CONTRA-008" in violation(caught) + + +def test_equal_event_timestamps_are_permitted(tmp_path): + """§2.5 says non-decreasing, not strictly increasing: events share a packet timestamp.""" + mjd = [fx.TSTART_MJD, fx.TSTART_MJD, fx.TSTART_MJD + 1 / 86400.0] + product = events.parse(fx.write_events(tmp_path / "e.fits", mjd=mjd), DIGEST) + assert sum(1 for _ in product.events("czt1")) == 3 + + +def test_rows_before_a_violation_are_yielded_and_valid(tmp_path): + """Streaming means the consumer sees the valid prefix, then the rule terminates the run.""" + mjd = [fx.TSTART_MJD + i / 86400.0 for i in range(6)] + mjd[4] = mjd[2] + stream = events.parse(fx.write_events(tmp_path / "e.fits", mjd=mjd), DIGEST).events("czt1") + seen = [] + with pytest.raises(ContractViolation): + for event in stream: + seen.append(event) + assert len(seen) == 4 + + +# ═══════════════════════════════════════════ traceability and purity + + +@pytest.mark.parametrize("product", ["lc", "spectra", "gti", "hk", "events"]) +def test_every_product_carries_the_digest_of_its_bytes(tmp_path, product): + """ADR-0005: every parsed entity traceable to its originating archive and digest.""" + built = { + "lc": lambda: lc.parse(fx.write_lightcurve(tmp_path / "a.fits"), DIGEST, + detector="czt1"), + "spectra": lambda: spectra.parse(fx.write_spectra(tmp_path / "b.fits"), DIGEST, + detector="czt1"), + "gti": lambda: gti.parse(fx.write_gti(tmp_path / "c.fits"), DIGEST, + detector="czt1"), + "hk": lambda: hk.parse(fx.write_hk(tmp_path / "d.fits"), DIGEST), + "events": lambda: events.parse(fx.write_events(tmp_path / "e.fits"), DIGEST), + }[product]() + assert built.source_digest == DIGEST + + +def test_no_hel1os_parser_reads_a_clock(): + """TIS §0.4: ingest_time is stamped at the acquisition boundary and nowhere else.""" + for module in (lc, spectra, gti, hk, events, orbit): + source = Path(module.__file__).read_text() + for forbidden in ("datetime.now", "utcnow", "time.time", "boundary.stamp"): + assert forbidden not in source, f"{module.__name__} reads a clock" + + +def test_no_hel1os_parser_interpolates_or_repairs(): + for module in (lc, spectra, gti, hk, events): + names = {name.lower() for name in dir(module)} + for forbidden in ("interpolate", "interp", "smooth", "fillna", "impute", + "resample", "repair", "ffill", "bfill"): + assert not any(forbidden in name for name in names), ( + f"{module.__name__} exposes {forbidden!r}") + + +def test_neither_instrument_imports_the_other(): + """F-07 and F-11 made structural: the two parser packages are independent. + + The only shared code is `_fits`, which knows nothing of either instrument's conventions. + """ + import contexts.ingest.parsers.solexs as solexs_pkg + + for module in (lc, spectra, gti, hk, events, orbit): + assert "parsers.solexs" not in Path(module.__file__).read_text() + for name in ("lc", "pi", "gti"): + source = Path(getattr(solexs_pkg, name).__file__).read_text() + assert "parsers.hel1os" not in source diff --git a/specs/contradictions/CONTRA-007.md b/specs/contradictions/CONTRA-007.md new file mode 100644 index 0000000..edae985 --- /dev/null +++ b/specs/contradictions/CONTRA-007.md @@ -0,0 +1,191 @@ +--- +id: CONTRA-007 +title: HEL1OS time comparisons are undecidable at float64 resolution +status: active +state: CLOSED +resolving_revision: r7 +supersedes: [] +superseded_by: null +source: null +source_date: 2026-09-12 +origin: authored +--- + +> **Authored, not migrated.** CONTRA-001 … CONTRA-006 are carried verbatim from +> `artifacts/v2/phase05/`. This record has no such source: it was raised by the governed +> re-implementation of the HEL1OS parsers (M3/E5 Issue #18) and its archive-wide verification. +> Cite it by ID; do not restate it ([ADR-0013](../../adr/ADR-0013.md)). + +# CONTRA-007 — HEL1OS time comparisons are undecidable at float64 resolution + +**Status: CLOSED 2026-09-12 — resolved by [SPEC-parsers](../parsers/SPEC-parsers.md) r7 §5.1.** + +**Defects A, B and C are one defect in three places:** the contract compares times derived from +MJD values without saying at what precision, and the values it compares are IEEE-754 float64. +Implemented exactly as written, three checks reject valid archive products by margins of one or +two representable steps. r7 adds `ε_t = 1 ms` (§5.1) and defines the comparisons; it changes no +physical rule and no statistic enters the contract. + +**Observations D, E and F are not falsifications** and carry no amendment. They are recorded +because #18's implementation had to resolve each of them, and an unrecorded resolution is an +invented convention. + +## How the defect reached a governed implementation + +**CONTRADICTION-006 Defect A ruled on the same numerical problem for §2.8 and fixed it in code +only** — *"The specification text is unchanged — the epsilon exists solely to prevent IEEE754 +boundary artifacts"* — and the Milestone V parsers applied the same `_FLOAT_EPS_S = 1e-3` to +§2.7 R-1 without the contract ever recording it. #18 was written from the contract text, as the +contract requires, and therefore reproduced the rejections that ruling had already resolved. + +**An implementation-only fix to a contract defect does not survive a re-implementation.** That +is the reason r7 moves the value into the text rather than repeating the code fix, and it is the +one generalisation this record makes. + +## DEFECT A — §2.7 R-1 H3 rejects 647 of 1,564 spectra products + +**The rule (r4):** accept H3 iff `col[0] == 0` exactly **and** `abs(col_span − header_span) ≤ one +EXPOSURE bin`. `col_span` is not defined anywhere in the contract. + +**`OBSERVED`, archive-wide (391 orbits, 1,564 spectra products):** + +| | Products | +|---|---| +| Resolve to H3 | 917 | +| Terminate at F-06 | **647** | +| Terminate for any other reason | 0 | + +Every one of the 647 is H3 in shape: `col TSTART[0] == 0.0` exactly, a uniform 20 s `EXPOSURE`, +and a span that agrees with the header to **between 6.4×10⁻⁷ s and 1.36×10⁻⁶ s past the one-bin +bound** — at most about **two representable float64 steps** of an MJD value at this epoch +(≈ 6.3×10⁻⁷ s each). The specification's own reference orbit is one of the products that fails: +`(61017.49963590554 − 61017.0000988685) × 86400` evaluates to `43160.00000007916`, 79 ns over an +exact bound, while §2.7 records that orbit as resolving to H3. + +**Two readings of `col_span` were available**, and they disagree on that reference orbit: +`TSTART[last] − TSTART[0]` is two bins short of the header span, `TSTOP[last] − TSTART[0]` is one. +Only the second reproduces the contract's own recorded verdict, and it is the one the Milestone V +parser used. + +### Amendment (APPLIED as r7) +1. **§2.7** — define `col_span = column TSTOP[last] − column TSTART[0]` and `header_span`, and + state the H3 bound as `≤ one EXPOSURE bin + ε_t`. +2. **§5.1** — define `ε_t`, its scope and what it is not. + +## DEFECT B — §2.8 header-span consistency rejects 94 of 391 housekeeping products + +**The rule (r4):** *"the global `mjd` range lies within the header `TSTART`/`TSTOP`"* — no +precision, no rule id. + +**`OBSERVED`, archive-wide:** 94 of 391 products terminate. **Every violation is exactly one +representable float64 step** (6.286×10⁻⁷ s): a boundary timestamp equals its header bound to +physical precision and differs from it by one step. None exceeds 1 µs; none is within four orders +of magnitude of a physically meaningful excursion, which would be seconds. + +This is the identical measurement CONTRADICTION-006 Defect A recorded (60 orbits above `TSTOP`, +40 below `TSTART`, all ~6×10⁻⁷ s). The count differs because that ruling measured the legacy +build's orbit set; the property is the same. + +### Amendment (APPLIED as r7) +**§2.8** — state header-span consistency as `TSTART − ε_t ≤ min(mjd)` and `max(mjd) ≤ TSTOP + ε_t`, +and assign it **F-06**, the id the Milestone V parser already raised. + +## DEFECT C — §2.5's span check has the same defect, and V-EVT-2 has no rule at all + +**The rules (r0):** *"`mjd` within `[TSTART,TSTOP]`"*, and *"`utc-isot` is a cross-check +(V-EVT-2)"* — named, with no comparison and no tolerance. + +**`OBSERVED`, archive-wide (1,564 detector HDUs, 1,352,158,522 event rows):** + +| Measurement | Result | +|---|---| +| HDUs whose `min(mjd)` falls below header `TSTART` | 182 — **all exactly one float64 step** (6.29×10⁻⁷ s) | +| HDUs whose `max(mjd)` rises above header `TSTOP` | 270 — **all exactly one float64 step** | +| HDUs exceeding either bound by ≥ 1 ms | **0** | +| `utc-isot` strings that fail to parse | **0** | +| `instant(mjd) − instant(utc-isot)` | **within ±0.5000 ms, every row** | +| Rows with `ener ≤ 0` | **0** | + +`utc-isot` is the `mjd` instant **rounded to the millisecond it carries**: the agreement is +exactly half the string's own resolution, in both directions, across every row of the corpus. +That is a measurement of the product, so the rule follows the value rather than a chosen number. + +### Amendment (APPLIED as r7) +**§2.5** — the span check gains `ε_t` and **F-06**; V-EVT-2 becomes +`|instant(mjd) − instant(utc-isot)| ≤ ½·r_isot + ε_t`, where `r_isot` is the resolution the string +itself carries, with **F-06** for a violation or an unparseable string. + +## OBSERVATION D — the five bands of one detector need not share a time axis + +§2.6 records `NAXIS2 ≈ 43171` and does not say whether a detector's five band HDUs are sampled +together. **They frequently are not:** of 1,564 light-curve products, **794 carry bands of unequal +length** — every CdTe product (782 of 782) and 12 CZT products. On the reference orbit CdTe1's +bands hold 43,154 / 43,171 / 43,163 / 43,133 / 43,171 samples. + +A consumer assuming one axis per detector would misalign up to 38 samples. **No amendment:** §2.6 +neither states nor denies a shared axis, and the parser keeps each band's own `MJD` column, which +is the only reading the measurement supports. Recorded so the property is pinned by a test rather +than rediscovered. + +## OBSERVATION E — §2.8's detector-health names are abbreviations + +§2.8 lists `czt{1,2}hotpix`, `hotpixcnt`, `hotpixthr`, `hotpixlgcstat`, `bunpxctr`. **The archive +carries no unprefixed column of those names.** All 391 orbits carry an identical 62-column +`HLSHK` table in which they appear as `czt1hotpixcnt`, `czt2hotpixthr`, `czt1bunpxctr` and so on — +which is also how §3 T4 names them (`czt1hotpixcnt`, `czt2hotpixcnt`). + +**No amendment:** the prefixed spelling is what both the archive and §3 carry, so §2.8's +abbreviation resolves without ambiguity. The reading is recorded rather than left implicit. +`czt1enth` declares `keV` and `czt2enth` declares no unit in **391 of 391** orbits, exactly as +§2.8 and §8 A-4 describe. + +## OBSERVATION F — three different counts of overlapping orbit pairs + +§4 records **46** time-overlapping orbit pairs; the Milestone VI version-resolution engine +recorded **49**; #18's `overlaps` predicate, applied to the stem-declared interval +(`start`, `start + dur`) of all 391 orbits, reports **50**. + +§4 defines no overlap measure, and the three numbers are consistent with three different ones. +**No amendment, and nothing is reconciled here:** #18 detects overlap and resolves precedence; +the authority on coverage is the minute-level map §4 assigns to the write path (M3/E5/#19), which +will produce the only count that governs ingestion. + +## OBSERVATION G — housekeeping inversions are far larger archive-wide than on the reference orbit + +§2.8 records the reference orbit's telemetry jitter — 424 of 9,513 steps decreasing, **max backward +step 892.4 ms** — and A-12 obliges Milestone VIII to report the distribution across all 391 orbits. +Parsing every product for #18 produced that measurement as a by-product: + +| `max_backward_step_s` across the 389 housekeeping products that parse | Value | +|---|---| +| median | **0.0 s** (most orbits have no inversion at all) | +| maximum | **1,153.4 s** | + +**No amendment, and no threshold.** §2.8 is explicit that the magnitude is reported and never +compared against an invented tolerance, so the parser records `n_out_of_order` and +`max_backward_step_s` and compares them to nothing. The figure is recorded here because it is +evidence for A-12, and because "sub-second packet-arrival jitter" — the reading the reference +orbit supports — does not describe a 1,153 s step. Whether that is jitter, a telemetry gap, or +something else is a Milestone VIII scientific question, and this record asserts no mechanism. + +## What is NOT affected + +Nothing physical. `ε_t` never touches a physical time difference: §2.8's inversion statistics +remain unthresholded, §2.6's `MJD` stays strictly increasing, GTI durations are unchanged, and +F-09's exact SoLEXS identity is untouched. The 20 fail-loud rule ids are unchanged. SoLEXS +parsing is unchanged in behaviour. + +**Duplicate HK timestamps remain F-16.** Two orbits (`HLS_20260201_120005_43198sec_lev1_V111`, +`HLS_20260202_000005_43183sec_lev1_V111`) terminate on duplicate `mjd`, exactly as +CONTRADICTION-006 Defect B was ruled, and exactly the two the V&V plan lists as known archive +defects. No amendment; they are archive-quality findings. + +**The §2.5 non-decreasing rule is NOT amended here.** Its archive-wide falsification is a separate +question with a separate answer, recorded OPEN in [CONTRA-008](CONTRA-008.md). + +## Measurements, and where they live + +Every number above is a measurement of the archive, recorded here and **not** in the contract — +the discipline CONTRADICTION-005 Defect B established: *measurements belong in the profile and the +contradiction record; invariants belong in the contract.* r7 encodes one value, `ε_t`, and states +its derivation without importing a statistic. diff --git a/specs/contradictions/CONTRA-008.md b/specs/contradictions/CONTRA-008.md new file mode 100644 index 0000000..34c6be0 --- /dev/null +++ b/specs/contradictions/CONTRA-008.md @@ -0,0 +1,103 @@ +--- +id: CONTRA-008 +title: §2.5's non-decreasing event rule is falsified by every event HDU in the archive +status: active +state: OPEN +resolving_revision: null +supersedes: [] +superseded_by: null +source: null +source_date: 2026-09-12 +origin: authored +--- + +> **Authored, not migrated**, and **OPEN**: it records a falsification and a proposed amendment +> that has **not** been applied. The parser enforces §2.5 as written until the owner rules. +> Cite it by ID; do not restate it ([ADR-0013](../../adr/ADR-0013.md)). + +# CONTRA-008 — §2.5's non-decreasing event rule is falsified by every event HDU + +**Status: OPEN 2026-09-12.** Raised by M3/E5 Issue #18, which implemented §2.5's row rules and +ran them against the whole corpus. + +**The rule (§2.5, unamended since r0):** *"Validation: … `mjd` non-decreasing …"* + +**`OBSERVED`, archive-wide (391 orbits, 1,564 detector HDUs, 1,352,158,522 rows):** + +| Measurement | Result | +|---|---| +| HDUs containing at least one backward step | **1,564 of 1,564 (100%)** | +| Orbits affected | **391 of 391** | +| Backward steps, total | **37,912,843** | +| Backward steps per HDU | min 4 · median 23,216 · max 229,417 | +| Largest backward step | min 0.2 ms · median 6.13 s · **max 1,156.36 s** | +| Consecutive rows sharing a timestamp | 1,265,162,302 of 1,352,156,958 steps (93.6%) | + +**Enforced as written, no real event stream completes.** The reader terminates at F-16 on the +first decrease of every HDU in the archive, having yielded only rows that satisfied every other +§2.5 rule. That is the behaviour #18 ships, because weakening a frozen rule in code is the one +thing the contract forbids: *"a deviation requires a logged amendment, not a code change."* + +## Why this was not seen before + +The Milestone V parser implemented the same check, but only under `load_columns=True`, and the +Milestone VII build never set it: `evt.fits` is 85.7 GB and events are deliberately outside the +canonical tables, so the rows were never read. **The rule had never been executed against the +archive.** This is the same pattern as CONTRADICTION-001, -003, -004 A and -005 A — a property +asserted from a single reading and never checked against the population — with one difference: +here the property was never checked at all. + +## What the measurement says about the data + +The three facts that bear on the ruling, and nothing beyond them: + +1. **The rows are not a sorted index.** 93.6% of consecutive rows share a timestamp, and the + backward steps are frequent (median 23,216 per HDU) rather than isolated. +2. **The backward steps are large.** A median largest-step of 6.13 s and a maximum of 1,156 s are + not sub-second packet jitter; they are far beyond the scale §2.8 r4 recorded for housekeeping + telemetry (max 892.4 ms on the reference orbit). +3. **Every other §2.5 rule holds on every row**: all four detector HDUs present, `DETNAM` correct, + `ener > 0` on all 1.35 billion rows, `mjd` inside the header span to within one float64 step + (CONTRA-007 Defect C), and `utc-isot` agreeing with `mjd` to within half the millisecond it + carries. + +**No mechanism is asserted here.** Whether event rows are written per detector packet, per +telemetry frame, or in some other archive order is not established by this measurement, and this +record does not guess. + +## Proposed amendment (NOT applied — the owner's ruling is required) + +The shape that already exists in the contract for exactly this situation is §2.8 r4, where the +owner removed HK's non-decreasing requirement, **declined** a proposal to sort inside the parser, +and required the parser to remain a lossless reader while recording inversion statistics. + +1. **§2.5** — remove `mjd` non-decreasing for event lists. Replace it with: `mjd` finite; `mjd` + within the header span (already r7); **inversion statistics recorded, never thresholded**, as + §2.8 requires for housekeeping. +2. **The parser MUST NOT sort.** §2.8 r4's general v2 principle holds: reading and transforming + are separate acts, and a consumer needing chronological order invokes the documented + out-of-parser utility. +3. **§8** — a new assumption scoping what the inversion statistics do and do not establish, + matching A-12's treatment of housekeeping jitter. +4. **No statistic enters the contract.** The measurements stay in this record. + +**My recommendation is amendment**, on the reasoning that the alternative disposition — treating +this as an archive defect, as CONTRADICTION-005 Defect C treated 12 SoLEXS days — would discard +the event streams of **all 391 orbits**, which is not an archive defect but a convention the +contract never established. §2.5 already calls `mjd` the canonical representation of three +redundant ones and offers `utc-isot` as its cross-check; ordering was assumed, not measured. + +**But this is the owner's call, not mine,** and for the same reason CONTRADICTION-006 Defect B +was: it turns on whether a backward step in an event list is telemetry or corruption. The r4 +housekeeping ruling answers that question one way for a different product, and a 1,156 s step is +large enough that the answer should not be inherited without a decision. + +## Until it is ruled + +- §2.5 stands as written; the parser enforces it and fails closed (F-16). +- The event reader is exposed and every other row rule is exercised on real data. +- `tests/integration/test_hel1os_parse_with_platform.py` pins the falsification: it asserts that a + real stream yields valid rows and then terminates at F-16, so the day the rule changes, the test + changes with it. +- No canonical table is affected: §2.5 keeps events out of T1–T7, and the write path (M3/E5/#19) + does not read them. diff --git a/specs/contradictions/index.md b/specs/contradictions/index.md index 8c1b4ff..56437cd 100644 --- a/specs/contradictions/index.md +++ b/specs/contradictions/index.md @@ -12,5 +12,7 @@ a record showing only accepted outcomes hides its most informative part. | [CONTRA-004](CONTRA-004.md) | CLOSED | `r4` | two HEL1OS parser-level rules are falsified | | [CONTRA-005](CONTRA-005.md) | CLOSED | `r5` | the archive-wide build falsifies three frozen rules | | [CONTRA-006](CONTRA-006.md) | CLOSED | `—` | two §2.8 HK checks falsified by the archive-wide rebuild | +| [CONTRA-007](CONTRA-007.md) | CLOSED | `r7` | HEL1OS time comparisons are undecidable at float64 resolution | +| [CONTRA-008](CONTRA-008.md) | OPEN | `—` | §2.5's non-decreasing event rule is falsified by every event HDU in the archive | Governing specification: [SPEC-parsers](../parsers/SPEC-parsers.md). diff --git a/specs/parsers/SPEC-parsers.md b/specs/parsers/SPEC-parsers.md index 49d8181..bda2f02 100644 --- a/specs/parsers/SPEC-parsers.md +++ b/specs/parsers/SPEC-parsers.md @@ -2,8 +2,8 @@ id: SPEC-parsers title: Real FITS parser specification (contract) status: active -revision: r6 -revisions: [r0, r1, r2, r3, r4, r5, r6] +revision: r7 +revisions: [r0, r1, r2, r3, r4, r5, r6, r7] supersedes: [] superseded_by: null source: artifacts/v2/phase05/PARSER_SPECIFICATION.md @@ -13,7 +13,9 @@ source_date: 2026-07-17 > Carried verbatim from `artifacts/v2/phase05/PARSER_SPECIFICATION.md`. This is a **contract**: the parsers > implement it, and a deviation requires a logged amendment, not a code change. > Amendments are adjudicated in [contradiction records](../contradictions/index.md). -> Current revision **r6**; history in the Revision History section below. +> Current revision **r7**; history in the Revision History section below. +> r0–r6 are carried from the source artifact. **r7 is authored in-tree** and has no counterpart there; +> every r7 change is marked *AMENDED r7* at the point of use, and the text it replaces is quoted, not deleted. # Phase 0.5.2 — Real FITS Parser Specification (CONTRACT) @@ -132,7 +134,8 @@ The implication is the strong half: it forbids the dangerous direction while ass **Columns:** `mjd` (D), `hlsobt` (D, `s`), `currtemp` (D, `degC`), `chn` (I), `ener` (D, **`keV`**), `recnum` (J), `utc-isot` (23A); **CZT additionally** `pix` (B), `offsetchn` (I). **Timestamps:** header `TSTART/TSTOP` in **MJD** (61017.0000988685 = 2025-12-08); columns provide `mjd`, spacecraft `hlsobt`, and an ISO string — **three redundant representations**. Canonical = `mjd`; `utc-isot` is a cross-check (V-EVT-2). **Units:** `ener` is **already energy-calibrated in keV** (unlike SoLEXS). `currtemp` is a per-event detector temperature. -**Validation:** all 4 HDUs present (F-03); `DETNAM` matches EXTNAME; `mjd` non-decreasing; `ener>0`; `mjd` within `[TSTART,TSTOP]`. +**Validation (AMENDED r7 — see §10 / CONTRA-007 Defect C):** all 4 HDUs present (F-03); `DETNAM` matches EXTNAME; `mjd` non-decreasing *(retained exactly as written — its archive-wide falsification is recorded OPEN in CONTRA-008 and is not amended here)*; `ener>0`; **`mjd` within `[TSTART − ε_t, TSTOP + ε_t]`** (§5.1; a violation → F-06); **V-EVT-2:** `|instant(mjd) − instant(utc-isot)| ≤ ½·r_isot + ε_t`, where `r_isot` is the resolution the `utc-isot` string itself carries (10⁻³ s for the observed 23-character form `YYYY-MM-DDTHH:MM:SS.sss`); an unparseable string or a violation → F-06. +*r6 text, superseded by r7:* "`mjd` within `[TSTART,TSTOP]`" — no comparison precision and no rule id; and V-EVT-2 named as a cross-check without a comparison rule. **Volume:** 85.7 GB total. **Not required for the canonical tables** — retained for Phase 1a pile-up/gain work. The 0.5.2 parser MUST expose an event reader but MUST NOT ingest events into the canonical minute tables. ### 2.6 HEL1OS `lightcurve_{czt,cdte}{1,2}.fits` — band light curves @@ -167,12 +170,14 @@ An unlisted `(family, DETCHANS)` pair **terminates via F-07** — exactly as an > > | # | Hypothesis | Accept iff | > |---|---|---| -> | **H3** | **relative seconds from header `TSTART`** *(the observed convention)* | `col[0] == 0` exactly **and** `abs(col_span − header_span)` ≤ one `EXPOSURE` bin | +> | **H3** | **relative seconds from header `TSTART`** *(the observed convention)* | `col[0] == 0` exactly **and** `abs(col_span − header_span)` ≤ one `EXPOSURE` bin **+ ε_t** *(amended r7; r4 text: "≤ one `EXPOSURE` bin")* | > | H1 | MJD days | `abs(col[0] − header TSTART) × 86400` ≤ 1 s | > | H2 | Unix seconds | `abs(col[0] − unix(header TSTART))` ≤ 1 s | > > **H3 is tested first** because `unit='s'` literally declares seconds. **H1 and H2 are retained for future compatibility** — a reprocessed product could legitimately switch to an absolute epoch, and silently mis-reading it would be worse than an extra branch. The resolved hypothesis and its residual are recorded in T7 provenance. > +> **Definitions (AMENDED r7 — see §10 / CONTRA-007 Defect A).** `col_span = column TSTOP[last] − column TSTART[0]` — the covered interval, from the first bin's start to the last bin's end. `header_span = (header TSTOP − header TSTART) × 86400 s`. `EXPOSURE` is the declared bin width. `ε_t` is the §5.1 time-representation allowance: it is added to the one-bin bound and never widens it to a second bin. The r4 text left `col_span` undefined and stated the bound exactly; on valid archive products the float64 MJD subtraction alone exceeds an exact one-bin bound (CONTRA-007). H1's and H2's 1 s bounds are unchanged. +> > `OBSERVED` (`hel1os_czt_spectra_czt1.fits`, orbit `HLS_20251208_000008`): col `TSTART` = `[0.0, 20.0, 40.0, …, 43120.0]`, uniform `EXPOSURE` = 20.0 s, header span 43,160.0 s → **H3**. **Note the cross-instrument asymmetry (do not conflate):** SoLEXS = 340 **PI** channels @ 1 s; HEL1OS = 341 **PHA** channels @ ~20 s. `PI ≠ PHA` (gain-corrected vs raw pulse height) and 340 ≠ 341. **No v2 code may treat these as a common channel space** (F-11). @@ -193,7 +198,7 @@ An unlisted `(family, DETCHANS)` pair **terminates via F-07** — exactly as an **Validation (AMENDED r4 — see §10 / CONTRADICTION-004 Defect B).** The strict **non-decreasing** requirement is **REMOVED**: it is falsified by the archive. `mjd` is a *measurement* written in telemetry-arrival order, not a sorted index. Validation is now: - `mjd` **finite**; - `mjd` **unique** (duplicates remain F-16 — a repeated timestamp is a genuine defect); -- **header-span consistency**: the global `mjd` range lies within the header `TSTART`/`TSTOP`; +- **header-span consistency** *(AMENDED r7 — see §10 / CONTRA-007 Defect B)*: `TSTART − ε_t ≤ min(mjd)` and `max(mjd) ≤ TSTOP + ε_t` (§5.1); a violation → F-06. *r4 text, superseded by r7:* "the global `mjd` range lies within the header `TSTART`/`TSTOP`" — no comparison precision and no rule id; - **inversion statistics recorded** (not thresholded): `n_out_of_order` and `max_backward_step_s` in T7 provenance; - `czt1temp`/`czt2temp` finite; `suninfov ∈ {0,1}`. @@ -318,6 +323,17 @@ One row per parsed source file: `src_file`, `src_sha256` (must equal 0.5.1 manif **F-12 is the single deliberate non-terminating rule** and is enumerated here so the exception is explicit rather than discovered. +### 5.1 Time-representation allowance `ε_t` *(added r7 — see §10 / CONTRA-007)* + +**`ε_t = 1 ms`** (1.157407407×10⁻⁸ MJD day). + +Where this contract compares a HEL1OS time derived from an MJD value — or a column offset composed onto one — against a bound or against another representation of the same instant, the comparison is made with `ε_t`: **§2.5** span and V-EVT-2, **§2.7** R-1 H3, **§2.8** header-span consistency. **Nowhere else.** + +- **What it absorbs:** numerical representation only. Header MJD keywords and time columns are IEEE-754 float64. One representable step of an MJD value is ≈ 6.3×10⁻⁷ s between MJD 32,768 and 65,535 (≈ 1.3×10⁻⁶ s from MJD 65,536), and time-offset columns are written at 10⁻⁶ s. Such values cannot be compared below the microsecond scale; a strict inequality there tests float representation, not data validity (CONTRADICTION-006 Defect A). +- **What it is not:** a physical tolerance, a jitter allowance, or a timing-accuracy claim. It never applies to physical time differences — §2.8 inversion statistics, GTI durations, §2.6's strictly increasing `MJD`, or F-09's exact SoLEXS identity — and it does not alter H1's or H2's 1 s bounds. +- **Why 1 ms:** it is the value the owner approved as the representation slack in CONTRADICTION-006 Defect A — applied there in code only, with this text left unchanged — and the value the Milestone V–VII parsers applied to R-1. r7 moves it from code into the contract. It lies about three orders of magnitude above the representation scale and at least three below the physical scales the checks protect: the 1 s light-curve cadence, the 20 s spectral bin, and the seconds-or-more of a wrong epoch, day or origin. A disagreement of 1 ms or more is reported, never absorbed. +- **No archive statistic is encoded here.** The measurements that motivated r7 are recorded in CONTRA-007. + --- ## 6. Validation Protocol (≥3 manually inspected days) @@ -484,3 +500,15 @@ Original contract, grounded in structure-only schema discovery of the real archi **Unchanged.** All rules, all other schemas, every measurement. r6 changes contract prose to match an owner-accepted table shape; it alters no value and adds no rule id. **Disposition.** T3 deviation **RATIFIED**. Milestone VII **CLOSED**; dataset version **FROZEN**. **Milestone VIII is now the final validation milestone: it shall discharge A-8, A-11, A-12, A-13, A-14 and resolve CONTRADICTION-003 through archive-wide scientific validation.** + +### r7 — 2026-09-12 (maintainer-approved direction D1/D1b for M3/E5 Issue #18; authored in-tree; raised by CONTRA-007) + +**Trigger.** The governed re-implementation of the HEL1OS parsers (M3/E5 Issue #18) was verified against all 391 orbits. Implemented from the r6 text alone, R-1 H3 rejected **647 of 1,564** spectra products and §2.8 header-span consistency rejected **94 of 391** housekeeping products — every rejection within about two representable float64 steps of its bound, and none by as much as 2×10⁻⁶ s. CONTRADICTION-006 Defect A had approved a 1 ms representation slack for §2.8 as an implementation-only fix and left this text unchanged, and the Milestone V parsers carried the same slack in R-1; a faithful implementation of the contract text therefore reproduced the rejections. The defect is the contract's silence about comparison precision, not the data. + +**Changes.** §5.1 added: the time-representation allowance `ε_t = 1 ms`, its scope (§2.5, §2.7 R-1 H3, §2.8 header span only) and what it is not. §2.7 R-1 H3's bound becomes one `EXPOSURE` bin + `ε_t`, with `col_span` and `header_span` defined. §2.8 header-span consistency is stated as inequalities with `ε_t` and assigned F-06. §2.5 `mjd` within `[TSTART, TSTOP]` gains `ε_t` and F-06, and V-EVT-2 receives a comparison rule (`≤ ½·r_isot + ε_t`, F-06). Every replaced clause is quoted beside its amendment. + +**Unchanged.** Every other rule, schema and policy. §2.5 `mjd` non-decreasing is **not** amended: its archive-wide falsification is recorded OPEN in CONTRA-008 and awaits a ruling. No archive statistic enters the contract. The 20 fail-loud rule ids are untouched; r7 assigns F-06 to §2.5 and §2.8 time checks that carried no id. + +**Bound value.** The direction — make these comparisons numerically well defined through a recorded amendment rather than a silent tolerance — was approved by the maintainer. The value `ε_t = 1 ms` is carried from CONTRADICTION-006 Defect A and is submitted for maintainer review with the Issue #18 pull request. + +**Disposition.** **CONTRA-007: CLOSED** by this revision. **CONTRA-008: OPEN.** diff --git a/specs/parsers/index.md b/specs/parsers/index.md index aefafb5..6539815 100644 --- a/specs/parsers/index.md +++ b/specs/parsers/index.md @@ -4,9 +4,9 @@ The contract the instrument parsers implement. | ID | Revision | Title | | --- | --- | --- | -| [SPEC-parsers](SPEC-parsers.md) | r6 | Real FITS parser specification (contract) | +| [SPEC-parsers](SPEC-parsers.md) | r7 | Real FITS parser specification (contract) | -Revisions present: `r0`, `r1`, `r2`, `r3`, `r4`, `r5`, `r6`. +Revisions present: `r0`, `r1`, `r2`, `r3`, `r4`, `r5`, `r6`, `r7`. Each revision beyond `r0` was forced by an adjudicated contradiction; see [contradiction records](../contradictions/index.md). diff --git a/tests/integration/test_hel1os_parse_with_platform.py b/tests/integration/test_hel1os_parse_with_platform.py new file mode 100644 index 0000000..a777efc --- /dev/null +++ b/tests/integration/test_hel1os_parse_with_platform.py @@ -0,0 +1,506 @@ +"""Real HEL1OS products through the whole platform, coexisting with SoLEXS. + +`| 18 | 18 | E5 | 900 | L | 17 | per-field vs spec | no imputation | parse fixtures |` + + #16 ISSDC adapter the acquisition boundary the archive arrives through + #10 provenance kernel mints every digest, records the derivation + #15 ingest contract supplies the one sanctioned clock read + #18 HEL1OS parsers canonicalise five product types + #12 domain model holds the rows and both times + #11 contract schemas validate the serialised Observations + #14 manifest records the Tier 0 orbit as referenced, never deposited + #13 import rules say Ingest may do this and reach no other context + +COEXISTENCE WITHOUT LEAKAGE IS THE POINT OF THIS FILE +------------------------------------------------------ +SoLEXS and HEL1OS disagree on nearly every convention — counts vs rates, undeclared vs +declared units, Unix seconds vs MJD, PI(340) vs PHA(341)/PHA(511), uppercase vs lowercase +column names. F-07 and F-11 exist because assuming a shared convention across them is the +specific error that has already been made. + +So the tests below put both instruments' Observations in one collection and assert that +every distinguishing fact survives: instrument identity, unit, quantity name, channel space +and time encoding. Nothing merges them, and nothing can. + +Real-corpus tests skip where the archive is absent (STD-12, E5 §17); the guard tests for +`.fits` **products**, never a directory — `aux/cztdis/*.txt` files are tracked in git while +the 132 GB of FITS are not, which is the trap E5 §16 names by name. +""" + +from __future__ import annotations + +import json +from pathlib import Path + +import jsonschema +import pytest +from referencing import Registry, Resource + +from contexts.ingest.parsers.hel1os import events, gti, hk, lc, orbit, spectra +from contexts.ingest.parsers.solexs import lc as solexs_lc +from contexts.ingest.parsers.solexs import pi as solexs_pi +from contexts.ingest.tests import hel1os_fixtures as fx +from contexts.ingest.tests import solexs_fixtures as sfx +from domain.entities import Observation +from domain.errors import ContractViolation +from domain.invariants import observation_is_wellformed +from domain.values import Digest, Identifier, Timestamp +from kernel.provenance import ( + Digest as KernelDigest, + ProvenanceStore, + begin_run, + digest_file, +) + +REPO_ROOT = Path(__file__).resolve().parents[2] +CONTRACTS = REPO_ROOT / "contracts" +HEL1OS_ROOT = REPO_ROOT / "research" / "data" / "aditya_l1" / "real_l1_v1" / "hel1os" +SOLEXS_ROOT = REPO_ROOT / "research" / "data" / "aditya_l1" / "real_l1_v1" / "solexs" + +#: The reference orbit §2.5–2.9 record their OBSERVED values against. +REFERENCE_ORBIT = "HLS_20251208_000008_43178sec_lev1_V111" + + +def real_orbit_present() -> bool: + """Products, not directories — E5 §16. + + `aux/cztdis/*.txt` pixel maps are tracked in git while the FITS are not, so the orbit + directory exists in a clean checkout and an `isdir` guard would pass where the data does + not exist. That exact guard silently disabled 188 tests once. + """ + if not HEL1OS_ROOT.is_dir(): + return False + return any(p for p in HEL1OS_ROOT.rglob("*.fits") if not p.name.startswith("._")) + + +real_only = pytest.mark.skipif( + not real_orbit_present(), reason="real HEL1OS orbit archive not extracted" +) +solexs_only = pytest.mark.skipif( + not (SOLEXS_ROOT.is_dir() + and any(p for p in SOLEXS_ROOT.rglob("*.gz") if not p.name.startswith("._"))), + reason="real SoLEXS archive not extracted", +) + + +def registry() -> Registry: + built = Registry() + for path in sorted(CONTRACTS.glob("*.schema.json")): + if path.name.startswith("._"): + continue + built = Resource.from_contents(json.loads(path.read_text())) @ built + return built + + +def validator_for(name: str): + schema = json.loads((CONTRACTS / f"{name}.schema.json").read_text()) + return jsonschema.Draft202012Validator(schema, registry=registry()) + + +def orbit_dir() -> Path: + return HEL1OS_ROOT / REFERENCE_ORBIT / "2025" / "12" / "08" / REFERENCE_ORBIT + + +# ═══════════════════════════════════════════ the whole path, on fixtures + + +def test_hel1os_rates_become_contract_valid_observations(tmp_path): + """Parse → domain → contract → provenance, with the digest carried throughout.""" + store = ProvenanceStore(tmp_path / "store") + path = fx.write_lightcurve(tmp_path / "lightcurve_czt1.fits", rows=60) + minted = digest_file(path) + registered = store.put_file(path) + + curve = lc.parse(path, Digest(minted.hex), detector="czt1") + observations = list(curve.observations( + source_id=Identifier("issdc-pradan"), + ingest_time=Timestamp("2025-12-20T09:00:00Z"), + )) + + # Five bands × 60 samples. + assert len(observations) == 5 * 60 + + validator = validator_for("observation") + for observation in observations[:40]: + validator.validate(observation.to_dict()) + assert observation_is_wellformed(observation) + + assert {str(o.source_digest) for o in observations} == {minted.hex} + assert store.has_artifact(KernelDigest(minted.hex)) + assert registered.digest.hex == minted.hex + + run = begin_run(context="ingest", event="parse") + canonical = store.put_bytes( + json.dumps([o.to_dict() for o in observations[:5]], sort_keys=True).encode()) + store.record(run, inputs=[minted], outputs=[canonical.digest]) + assert minted in store.ancestors(canonical.digest) + + +def test_the_orbit_becomes_a_valid_tier_0_manifest(tmp_path): + """ADR-0023: an orbit archive is referenced, never deposited.""" + path = fx.write_lightcurve(tmp_path / "lightcurve_czt1.fits", rows=10) + minted = digest_file(path) + + manifest = { + "kind": "dataset", "digest": minted.hex, "tier": 0, + "recorded_at": "2025-12-20T09:00:00Z", + "retention": {"class": "permanent"}, + "retrieval": {"provider": "ISSDC PRADAN", + "locator": f"hel1os/{REFERENCE_ORBIT}", + "requires_credentials": True}, + } + validator_for("manifest").validate(manifest) + + redistributing = dict(manifest) + redistributing["deposition"] = { + "provider": "Zenodo", "url": "https://zenodo.org/records/1", "doi": None} + assert list(validator_for("manifest").iter_errors(redistributing)) + + +def test_the_parsers_stay_inside_the_ingest_import_rule(): + """ADR-0026, against the shipped policy.""" + from tools.gates.imports import POLICIES, run + + report, code = run(POLICIES) + assert code == 0, report.violations + + ingest = next(p for p in POLICIES if p.package == "contexts.ingest") + internal = ingest.allow & {"contracts", "domain", "kernel", "contexts", "apps", + "tools", "registry", "tests"} + assert internal == {"contracts", "domain", "kernel"} + assert "contexts" not in ingest.allow + + +# ═══════════════════════════════════════════ coexistence without leakage + + +def test_both_instruments_coexist_in_one_collection(tmp_path): + """SoLEXS and HEL1OS Observations side by side, each keeping its own conventions. + + Every distinguishing fact survives: instrument identity, unit, and quantity name. If any + collapsed, a rate would be summed with a count and nothing downstream would notice. + """ + solexs_path = sfx.write_lc(sfx.sdd2(tmp_path) / "s.lc", rows=86_400) + solexs = solexs_lc.parse(solexs_path, Digest(digest_file(solexs_path).hex)) + solexs_rows = [ + o for i, o in enumerate(solexs.observations( + source_id=Identifier("issdc-pradan"), ingest_time=None)) if i < 50 + ] + + hel1os_path = fx.write_lightcurve(tmp_path / "lightcurve_czt1.fits", rows=10) + hel1os = lc.parse(hel1os_path, Digest(digest_file(hel1os_path).hex), detector="czt1") + hel1os_rows = list(hel1os.observations( + source_id=Identifier("issdc-pradan"), ingest_time=None)) + + both: list[Observation] = solexs_rows + hel1os_rows + validator = validator_for("observation") + for observation in both: + validator.validate(observation.to_dict()) + + instruments = {str(o.instrument_id) for o in both} + assert instruments == {"solexs-sdd2", "hel1os-czt1"} + + # The unit distinguishes them, and the two are never the same string. + assert {o.unit for o in solexs_rows} == {"counts"} + assert {o.unit for o in hel1os_rows} == {"cts/s"} + assert not ({o.unit for o in solexs_rows} & {o.unit for o in hel1os_rows}) + + # And so does the quantity: HEL1OS carries its band, SoLEXS carries none. + assert {o.quantity for o in solexs_rows} == {"counts"} + assert all(o.quantity.startswith("count_rate_") for o in hel1os_rows) + + +def test_the_three_channel_spaces_never_merge(tmp_path): + """F-11: SoLEXS PI(340), CZT PHA(341), CdTe PHA(511). + + Three spaces, three modules, three separately-validated declarations. Stacking any two + would fabricate a channel correspondence that does not exist. + """ + czt = spectra.parse(fx.write_spectra(tmp_path / "czt.fits", detector="czt1"), + Digest("a" * 64), detector="czt1") + cdte = spectra.parse(fx.write_spectra(tmp_path / "cdte.fits", detector="cdte1"), + Digest("b" * 64), detector="cdte1") + solexs = solexs_pi.parse( + sfx.write_pi(sfx.sdd2(tmp_path) / "s.pi", rows=86_400), Digest("c" * 64)) + + assert len(czt.channel_map) == 341 + assert len(cdte.channel_map) == 511 + assert len(solexs.channel_map) == 340 + assert len({len(czt.channel_map), len(cdte.channel_map), len(solexs.channel_map)}) == 3 + + # PI is not PHA, and the two modules say so independently. + assert solexs.header.chantype == "PI" + assert czt.header.chantype == "PHA" == cdte.header.chantype + + +def test_the_two_time_encodings_do_not_cross(tmp_path): + """SoLEXS states Unix seconds; HEL1OS states MJD. Both resolve to the same UTC. + + Reading one with the other's epoch is a ~49-year error of the kind F-05 exists for, so + each parser converts with its own declared encoding and the results are compared as + instants rather than as numbers. + """ + solexs_instant = solexs_lc.unix_to_timestamp(1_715_644_800.0) + hel1os_instant = lc.mjd_to_timestamp(fx.TSTART_MJD) + + assert str(solexs_instant) == "2024-05-14T00:00:00Z" + assert str(hel1os_instant).startswith("2025-12-08T00:00:08") + assert solexs_lc.MJDREFI_UNIX == 40587 + assert lc.MJD_UNIX_EPOCH == 40587.0 + + # The same instant expressed both ways agrees. + unix_of_mjd = (fx.TSTART_MJD - 40587.0) * 86_400.0 + assert lc.mjd_to_timestamp(fx.TSTART_MJD) == solexs_lc.unix_to_timestamp(unix_of_mjd) + + +# ═══════════════════════════════════════════ the real corpus + + +@real_only +def test_the_version_distribution_matches_the_specification(): + """§4 `OBSERVED`: V111 ×371, V211 ×16, V112 ×3, V311 ×1 across 391 orbits.""" + counts: dict[int, int] = {} + stems = [] + for path in HEL1OS_ROOT.iterdir(): + if not path.is_dir() or path.name.startswith("._"): + continue + parsed = orbit.parse_stem(path.name) + stems.append(parsed) + counts[parsed.version] = counts.get(parsed.version, 0) + 1 + + assert len(stems) == 391 + assert counts == {111: 371, 211: 16, 112: 3, 311: 1} + + +@real_only +def test_both_overlap_classes_exist_and_resolve(): + """§4's Class A and Class B, both present in the archive. + + Class A resolves by precedence at file level. Class B is *detected* here and resolved by + the minute-level coverage map, which §4 ties to emitting T3/T4/T5 — #19's work. + """ + a1 = orbit.parse_stem("HLS_20251208_000008_43178sec_lev1_V111") + a2 = orbit.parse_stem("HLS_20251208_000008_43178sec_lev1_V211") + assert (HEL1OS_ROOT / a1.stem).is_dir() and (HEL1OS_ROOT / a2.stem).is_dir() + assert orbit.overlaps(a1, a2) + assert orbit.precedence(a1, a2) is a2 + + b1 = orbit.parse_stem("HLS_20251207_120003_43195sec_lev1_V211") + b2 = orbit.parse_stem("HLS_20251207_121028_42570sec_lev1_V111") + assert (HEL1OS_ROOT / b1.stem).is_dir() and (HEL1OS_ROOT / b2.stem).is_dir() + assert orbit.overlaps(b1, b2) + assert b1.start_epoch != b2.start_epoch and b1.duration_s != b2.duration_s + + +@real_only +def test_the_real_band_lightcurves_carry_the_specified_bands(): + """§2.6 `OBSERVED`, re-measured: five bands per detector, edges from EXTNAME.""" + for detector, expected in (("czt1", fx.CZT_BANDS), ("cdte1", fx.CDTE_BANDS)): + family = "czt" if detector.startswith("czt") else "cdte" + path = orbit_dir() / family / f"lightcurve_{detector}.fits" + product = lc.parse(path, Digest(digest_file(path).hex), detector=detector) + + assert [(b.low_kev, b.high_kev) for b in product.bands] == list(expected) + assert sum(1 for b in product.bands if b.is_total) == 1 + + # §2.6 states `NAXIS2 ≈ 43171` — approximate, and the archive shows why the word + # matters. Measured on this orbit: + # CZT1 43171, 43171, 43171, 43171, 43171 — one shared axis + # CdTe1 43154, 43171, 43163, 43133, 43171 — FIVE DIFFERENT LENGTHS + # So the bands of one CdTe detector do NOT share a time axis. §2.6 neither states + # nor denies this; the parser keeps each band's own MJD column, which is the only + # reading that survives the measurement. Recorded as CONTRA-007 Observation D. + for band in product.bands: + assert 43_000 <= len(band.mjd) <= 43_200, ( + f"{detector} {band.extname}: {len(band.mjd)} samples") + assert len(band.mjd) == len(band.rate) == len(band.stat_err), ( + f"{detector} {band.extname}: columns of unequal length") + + +@real_only +def test_cdte_bands_do_not_share_a_time_axis(): + """CONTRA-007 Observation D, asserted rather than assumed. + + A consumer building a per-detector band × time array would silently misalign up to 38 + samples if it assumed one axis. The parser gives each band its own, and this pins the + property so a future change to either the archive or the parser is visible. + """ + path = orbit_dir() / "cdte" / "lightcurve_cdte1.fits" + product = lc.parse(path, Digest(digest_file(path).hex), detector="cdte1") + lengths = {band.extname: len(band.mjd) for band in product.bands} + assert len(set(lengths.values())) > 1, lengths + + czt_path = orbit_dir() / "czt" / "lightcurve_czt1.fits" + czt = lc.parse(czt_path, Digest(digest_file(czt_path).hex), detector="czt1") + assert len({len(band.mjd) for band in czt.bands}) == 1 + + +@real_only +def test_the_real_spectra_resolve_to_h3_with_both_channel_spaces(): + """§2.7 `OBSERVED`: 341 for CZT, 511 for CdTe, 2157 rows, epoch H3.""" + for detector, detchans in (("czt1", 341), ("cdte1", 511)): + family = "czt" if detector.startswith("czt") else "cdte" + path = orbit_dir() / family / f"hel1os_{family}_spectra_{detector}.fits" + product = spectra.parse(path, Digest(digest_file(path).hex), detector=detector) + + assert product.header.detchans == detchans + assert product.header.chantype == "PHA" + assert product.header.hduclas3 == "COUNT" + assert product.header.rows == 2157 + assert product.epoch.hypothesis == "H3" + assert product.epoch.exposure_s == 20.0 + + first = next(iter(product.spectra())) + assert len(first.counts) == detchans + + +@real_only +def test_the_real_housekeeping_reproduces_the_recorded_inversion_statistics(): + """§2.8 `OBSERVED`: 9,514 rows; 424 of 9,513 steps decrease; max backward 892.4 ms; + 0 duplicates; range within the header span. + + Every figure re-measured. §2.8 reports these and defines no threshold, so the assertions + are on the measurements themselves rather than on any tolerance. + """ + path = orbit_dir() / "aux" / "hk.fits" + product = hk.parse(path, Digest(digest_file(path).hex)) + stats = product.inversions + + assert stats.rows == 9_514 + assert stats.steps == 9_513 + assert stats.n_out_of_order == 424 + assert round(stats.max_backward_step_ms, 1) == 892.4 + assert len(set(product.mjd)) == len(product.mjd) + assert product.header_tstart <= min(product.mjd) + assert max(product.mjd) <= product.header_tstop + + +@real_only +def test_the_real_housekeeping_is_not_sorted(): + """§2.8 r4, binding: the parser preserves archive order exactly and performs no sorting. + + 424 inversions survive, which is the observable proof that nothing reordered them. + """ + path = orbit_dir() / "aux" / "hk.fits" + product = hk.parse(path, Digest(digest_file(path).hex)) + assert list(product.mjd) != sorted(product.mjd) + assert product.inversions.n_out_of_order == 424 + + +@real_only +def test_the_real_gti_products_use_lowercase_columns(): + """§2.9 / §8 A-2: lowercase `tstart`/`tstop`, unlike SoLEXS's uppercase.""" + for detector in ("czt1", "czt2", "cdte1", "cdte2"): + path = orbit_dir() / "aux" / f"gti{detector}.fits" + product = gti.parse(path, Digest(digest_file(path).hex), detector=detector) + assert product.detector_active + assert len(product.intervals) == 1 + assert product.live_time_s > 43_000 + + +@real_only +def test_the_real_event_list_has_all_four_detectors(): + """§2.5 `OBSERVED`: four detector HDUs, 1.3–1.6 M rows each, `ener` in keV.""" + path = orbit_dir() / "events" / "evt.fits" + product = events.parse(path, Digest(digest_file(path).hex)) + + assert {d.detector for d in product.detectors} == {"czt1", "czt2", "cdte1", "cdte2"} + for detector in product.detectors: + assert 1_300_000 <= detector.rows <= 1_600_000 + + first = next(iter(product.events("czt1"))) + assert first.energy_kev > 0 + assert str(first.valid_time).startswith("2025-12-08T") + + # §2.5: exposed, and deliberately not ingestible into the canonical tables. + assert not hasattr(product, "observations") + + +@real_only +@solexs_only +def test_both_real_instruments_coexist_without_leakage(): + """The coexistence claim, on real data from both instruments. + + SoLEXS 2024-05-14 (Unix seconds, PI 340, counts) and HEL1OS 2025-12-08 (MJD, PHA 341, + cts/s) parsed by separate modules into one collection, with every distinguishing fact + intact and no shared convention anywhere. + """ + solexs_path = (SOLEXS_ROOT / "AL1_SLX_L1_20240514_v1.0" / "AL1_SLX_L1_20240514_v1.0" + / "SDD2" / "AL1_SOLEXS_20240514_SDD2_L1.lc.gz") + solexs = solexs_lc.parse(solexs_path, Digest(digest_file(solexs_path).hex)) + + hel1os_path = orbit_dir() / "czt" / "lightcurve_czt1.fits" + hel1os = lc.parse(hel1os_path, Digest(digest_file(hel1os_path).hex), detector="czt1") + + solexs_rows = [ + o for i, o in enumerate(solexs.observations( + source_id=Identifier("issdc-pradan"), ingest_time=None)) if i < 20 + ] + hel1os_rows = [] + for observation in hel1os.observations(source_id=Identifier("issdc-pradan"), + ingest_time=None): + hel1os_rows.append(observation) + if len(hel1os_rows) == 20: + break + + validator = validator_for("observation") + for observation in solexs_rows + hel1os_rows: + validator.validate(observation.to_dict()) + + assert {str(o.instrument_id) for o in solexs_rows} == {"solexs-sdd2"} + assert {str(o.instrument_id) for o in hel1os_rows} == {"hel1os-czt1"} + assert {o.unit for o in solexs_rows} == {"counts"} + assert {o.unit for o in hel1os_rows} == {"cts/s"} + + # Different missions, different years, different digests — nothing shared but the schema. + assert str(solexs_rows[0].valid_time).startswith("2024-05-14") + assert str(hel1os_rows[0].valid_time).startswith("2025-12-08") + assert solexs_rows[0].source_digest != hel1os_rows[0].source_digest + +@real_only +def test_the_real_housekeeping_carries_every_decisive_column_with_its_unit(): + """§2.8's decisive columns on the reference orbit, with the units the archive declares. + + The detector-health columns exist only with a detector prefix (`czt1hotpixcnt`), which is + how §3 T4 names them and how §2.8's abbreviation resolves — recorded as CONTRA-007 + Observation E rather than inferred silently. + """ + path = orbit_dir() / "aux" / "hk.fits" + product = hk.parse(path, Digest(digest_file(path).hex)) + + for name in hk.CAPTURED: + assert name in product.columns, name + assert product.units["czthvmon"] == "V" + assert product.units["cdtehvmon"] == "V" + assert product.units["czt1enth"] == "keV" + assert {product.units[name] for name in hk.THERMAL} == {"degC"} + assert {product.units[name] for name in hk.RATES} == {"c/s"} + + # §8 A-4: czt2enth declares no unit in the archive, so czt1enth's is applied and recorded. + assert product.units["czt2enth"] == "keV" + assert any("czt2enth" in assumption for assumption in product.assumptions) + assert isinstance(product.columns["czt1hotpixcnt"][0], int) + + +@real_only +def test_real_event_rows_pass_the_r7_rules_until_the_archive_falsifies_section_2_5(): + """The §2.5 row rules on real events, and the falsification they run into. + + Every row the archive yields here satisfies r7's span check and V-EVT-2 — `utc-isot` agrees + with `mjd` to the millisecond it carries. The stream then terminates at F-16 because `mjd` + steps backward, which §2.5 forbids and the archive does anyway: CONTRA-008, OPEN. The rule + is enforced as written, so this test pins the falsification rather than hiding it. + """ + path = orbit_dir() / "events" / "evt.fits" + product = events.parse(path, Digest(digest_file(path).hex)) + + seen = 0 + with pytest.raises(ContractViolation) as caught: + for event in product.events("czt1"): + seen += 1 + assert event.energy_kev > 0 + assert event.pixel is not None + assert seen > 0 + assert caught.value.message.startswith("F-16") + assert "CONTRA-008" in caught.value.message +