diff --git a/.github/workflows/constitution.yml b/.github/workflows/constitution.yml index ae34fdd..13c1672 100644 --- a/.github/workflows/constitution.yml +++ b/.github/workflows/constitution.yml @@ -95,7 +95,10 @@ jobs: with: python-version: "3.12" - - run: pip install pytest + # numpy + astropy added by M3/E5/#17: the SoLEXS products are FITS, and the parsers' + # unit tests build real FITS fixtures to fire each fail-loud rule against. Without + # these the whole parser suite would error at import rather than run. + - run: pip install pytest numpy astropy # `domain` added by M2/E4/#12 and `contexts` by M3/E5/#15, each placing its unit tests # beside the code they test. Without the path here those tests would exist and never @@ -134,7 +137,9 @@ jobs: # domain types against the normative contracts (ADR-0019), which cannot be done # without a validator. The test skips rather than fails if the library is absent, so # omitting this would have turned a cross-component check into a silent skip. - - run: pip install pytest jsonschema + # astropy + numpy added by M3/E5/#17 for the same reason as the unit job; jsonschema + # was added by #12 to validate the domain against the normative contracts. + - run: pip install pytest jsonschema numpy astropy - name: Integration tests run: pytest tests/integration -q -rs diff --git a/contexts/ingest/parsers/__init__.py b/contexts/ingest/parsers/__init__.py new file mode 100644 index 0000000..e0904af --- /dev/null +++ b/contexts/ingest/parsers/__init__.py @@ -0,0 +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. + +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 +no parser registry and no dispatch: a caller selects the module for the product it holds. +""" diff --git a/contexts/ingest/parsers/solexs/__init__.py b/contexts/ingest/parsers/solexs/__init__.py new file mode 100644 index 0000000..676f4e6 --- /dev/null +++ b/contexts/ingest/parsers/solexs/__init__.py @@ -0,0 +1,17 @@ +"""SoLEXS parsers — `.lc`, `.pi` and `.gti` (M3/E5/#17). + +Each implements its section of `SPEC-parsers@r6`, which is a **contract**: a deviation +requires a logged amendment, not a code change. + + lc §2.1 total-band counts per second. NaN is data, not an error. + pi §2.2 340-channel PI spectra per second. Ordinal channels; keV is prohibited. + gti §2.3 inclusive second-marks; Σ(STOP−START+1) == EXPOSURE, exactly. + +`.hk` (§2.4, v1.1 only, 5 of 436 days) is NOT implemented here: TIS E5 §5 lists this epic's +SoLEXS modules as `parsers/solexs/{lc,pi,gti}` and assigns no `.hk` module to any issue. The +gap is reported in the Issue #17 completion record rather than filled by guessing at scope. +""" + +from contexts.ingest.parsers.solexs import gti, lc, pi + +__all__ = ["gti", "lc", "pi"] diff --git a/contexts/ingest/parsers/solexs/_fits.py b/contexts/ingest/parsers/solexs/_fits.py new file mode 100644 index 0000000..7d0d5be --- /dev/null +++ b/contexts/ingest/parsers/solexs/_fits.py @@ -0,0 +1,96 @@ +"""Strict FITS access. Every lookup fails loud; none has a default. + +`SPEC-parsers@r6` §5 opens with the rule this module exists to make unbreakable: + +> **No default-on-missing-key** (`header.get(K, default)` is banned for any physically +> meaningful key — that idiom is the direct cause of v1's thirty-sprint failure). + +The concrete failure was F-05: v1 defaulted `MJDREF` to 58484 where the archive declares +40587, producing a **~49-year** timestamp error that nothing detected. A default is not a +convenience — it is a fabricated measurement wearing the shape of a real one. + +So there is no `get` here. `keyword()` raises. `column()` raises. A caller who wants a +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. +""" + +from __future__ import annotations + +from typing import Any + +from domain.errors import ContractViolation + + +def fail(rule: str, pointer: str, message: str) -> None: + """Terminate with a diagnostic naming file, HDU and expectation (§5). + + `ContractViolation` is the TIS §0.2 class for *"input fails its schema or an invariant"*. + The F-rule identifier is carried in the message rather than in a new exception type, + because §0.2 permits exactly five failure classes and no others. + """ + raise ContractViolation(pointer, f"{rule}: {message}") + + +def hdu(hdul: Any, name: str, *, source: str) -> Any: + """Fetch an HDU **by name**. F-02: HDU order is not a contract. + + v1 read HDUs positionally. The archive's ordering happens to be stable, which is exactly + what makes a positional read dangerous: it works until it does not, and then it reads a + different HDU rather than failing. + """ + try: + return hdul[name] + except KeyError: + present = [x.name for x in hdul] + fail("F-02", f"/{source}", f"HDU {name!r} absent; present by name: {present}") + + +def keyword(header: Any, key: str, *, source: str, hdu_name: str) -> Any: + """Fetch a header keyword. Absent is F-04/F-05 territory, never a default.""" + if key not in header: + fail( + "F-05", + f"/{source}#{hdu_name}/{key}", + f"keyword {key!r} absent from HDU {hdu_name!r}. A default here is how v1 " + f"produced a ~49-year timestamp error (§5).", + ) + return header[key] + + +def expect( + header: Any, key: str, expected: Any, *, source: str, hdu_name: str, rule: str +) -> Any: + """Fetch a keyword and require an exact value.""" + actual = keyword(header, key, source=source, hdu_name=hdu_name) + if actual != expected: + fail( + rule, + f"/{source}#{hdu_name}/{key}", + f"{key} is {actual!r}, expected {expected!r}", + ) + return actual + + +def column(table: Any, name: str, *, source: str, hdu_name: str) -> Any: + """Fetch a column by name, case-insensitively. F-04. + + Case-insensitive because FITS column names are not case-sensitive by the standard, and a + case difference is a spelling of the same column rather than a different one. Absence is + still fatal: *"v1 sought `RATE`; real is `COUNTS`"* is F-04's stated origin, and a parser + that silently produced no column would have reproduced it. + """ + wanted = name.upper() + for candidate in table.columns.names: + if candidate.upper() == wanted: + return table.data[candidate] + fail( + "F-04", + f"/{source}#{hdu_name}/{name}", + f"column {name!r} absent; present: {list(table.columns.names)}", + ) + + +__all__ = ["column", "expect", "fail", "hdu", "keyword"] diff --git a/contexts/ingest/parsers/solexs/gti.py b/contexts/ingest/parsers/solexs/gti.py new file mode 100644 index 0000000..38c093b --- /dev/null +++ b/contexts/ingest/parsers/solexs/gti.py @@ -0,0 +1,266 @@ +"""SoLEXS `.gti.gz` — good time intervals. + +Implements `SPEC-parsers@r6` §2.3. + +THE INCLUSIVE CONVENTION, AND WHY THE TOLERANCE IS ZERO +-------------------------------------------------------- +§2.3 as amended at r1: `START`/`STOP` are **inclusive second-marks** at 1-second sampling, so + + live_time(interval) = STOP − START + 1 + +`OBSERVED` on 2024-05-14 SDD2: 5 intervals, `Σ(STOP−START+1) = 86395.0 s`, equal to the +declared `EXPOSURE` **with exactly zero error**. The pre-amendment exclusive reading produced +"4 gaps of ~2 s", which was an artifact of the assumption rather than a property of the data. + +F-09 requires that equality to be **exact, tolerance 0 s**. §5 states why: *"the relation is +definitional, and a tolerance would re-admit the ambiguity that produced CONTRADICTION-001."* +A tolerance here would not absorb noise — there is no noise in a sum of integers — it would +absorb a convention error, which is precisely the thing that must not be absorbed. + +**The convention is verified for the implementation target only** (§2.3, §8 A-8). It is not +promoted to archive truth. Milestone VIII must verify it across all 436 archives, and §2.3 is +explicit that *"any deviation is a scientific finding and MUST terminate validation (never +tolerate, never widen)"*. So this parser fails on deviation and reports the arithmetic; it +does not widen. + +EXPOSURE IS A STRING IN THE ARCHIVE +----------------------------------- +§2.3 records it: HDU1 `EXPOSURE` is `'86395.0'`, a string, while the primary `TSTART`/`TSTOP` +are ISO-8601 strings. The parser converts explicitly and fails loud if the conversion does +not succeed — it does not coerce, and it does not fall back to a computed value, because a +computed fallback would make F-09 compare a number against itself. + +AN EMPTY GTI IS LEGAL +--------------------- +F-12 is *"the single deliberate non-terminating rule"*. `NAXIS2 == 0` means the detector was +inactive — §1.1 records SDD1 at zero rows on 2024-05-14 with empty-string `TSTART`/`TSTOP` — +and the correct response is `detector_active = False` and zero intervals, not an error. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass +from pathlib import Path + +from contexts.ingest.parsers.solexs import _fits +from contexts.ingest.parsers.solexs.lc import unix_to_timestamp +from domain.values import Digest, Identifier, Timestamp + + +@dataclass(frozen=True) +class Interval: + """One good time interval, with inclusive second-marks (§2.3 r1).""" + + start: float + stop: float + + def __post_init__(self) -> None: + if not (math.isfinite(self.start) and math.isfinite(self.stop)): + _fits.fail("F-16", "/GTI", f"interval bound not finite: {self}") + if self.start >= self.stop: + _fits.fail( + "F-09", "/GTI", + f"START {self.start!r} is not less than STOP {self.stop!r}; §2.3 requires " + f"START float: + """`STOP − START + 1` — the inclusive convention, amended at r1.""" + return self.stop - self.start + 1 + + @property + def start_utc(self) -> Timestamp: + return unix_to_timestamp(self.start) + + @property + def stop_utc(self) -> Timestamp: + return unix_to_timestamp(self.stop) + + def covers(self, second: float) -> bool: + """Inclusive at both ends, which is the whole content of the r1 amendment.""" + return self.start <= second <= self.stop + + +@dataclass(frozen=True) +class GoodTimeIntervals: + """A parsed `.gti` product. + + `detector_active` is False exactly when the archive declared zero rows — the F-12 path. + It is a recorded state, not an inference: nothing here decides a detector was inactive + from the absence of something else. + """ + + detector: str + source_path: Path + source_digest: Digest + intervals: tuple[Interval, ...] + declared_exposure: float | None + detector_active: bool + + @property + def instrument_id(self) -> Identifier: + return Identifier(f"solexs-{self.detector.lower()}") + + @property + def live_time(self) -> float: + """Σ(STOP − START + 1) over all intervals.""" + return sum(interval.live_time for interval in self.intervals) + + def excluded_seconds(self, tstart: float, rows: int) -> tuple[int, ...]: + """Day-offsets not covered by any interval. + + §2.3 `OBSERVED` on 2024-05-14 SDD2: [0, 5, 30072, 30078, 83951] — second 0 because + data begins at 00:00:01, plus four isolated 1-second dropouts. + + Offered here so the day-assembly layer (#19) can evaluate the `NaN ⇒ GTI-excluded` + implication that §2.1 explicitly places outside the single-file parsers. + """ + covered = set() + for interval in self.intervals: + first = int(interval.start - tstart) + last = int(interval.stop - tstart) + covered.update(range(first, last + 1)) + return tuple(offset for offset in range(rows) if offset not in covered) + + +def parse(path: Path, digest: Digest) -> GoodTimeIntervals: + """Parse one `.gti.gz` product. Every §2.3 validation, exact where the spec says exact.""" + 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: + primary = _fits.hdu(hdul, "PRIMARY", source=source).header + gti_hdu = _fits.hdu(hdul, "GTI", source=source) + gti = gti_hdu.header + + rows = int(_fits.keyword(gti, "NAXIS2", source=source, hdu_name="GTI")) + detector = _detector_from_path(path, source) + + # ── F-12: the single deliberate non-terminating rule ──────────────────── + if rows == 0: + return GoodTimeIntervals( + detector=detector, + source_path=path, + source_digest=digest, + intervals=(), + declared_exposure=None, + detector_active=False, + ) + + starts = _fits.column(gti_hdu, "START", source=source, hdu_name="GTI") + stops = _fits.column(gti_hdu, "STOP", source=source, hdu_name="GTI") + intervals = tuple( + Interval(float(a), float(b)) for a, b in zip(starts, stops) + ) + + _validate_ordering(intervals, source) + + declared = _declared_exposure(gti, source) + summed = sum(interval.live_time for interval in intervals) + if summed != declared: + # F-09, exact. §5: a tolerance would re-admit the ambiguity that produced + # CONTRADICTION-001, and the sum of integer second-counts has no noise to absorb. + _fits.fail( + "F-09", f"/{source}#GTI/EXPOSURE", + f"Σ(STOP−START+1) = {summed!r} s but EXPOSURE declares {declared!r} s. " + f"§2.3 requires exact equality with tolerance 0 s. The inclusive convention " + f"is verified for 2024-05-14 SDD2 only (§8 A-8); a deviation is a scientific " + f"finding and must terminate validation, never widen the tolerance.", + ) + + _validate_within_day(primary, intervals, source) + + return GoodTimeIntervals( + detector=detector, + source_path=path, + source_digest=digest, + intervals=intervals, + declared_exposure=declared, + detector_active=True, + ) + + +def _detector_from_path(path: Path, source: str) -> str: + """SDD1 or SDD2, from the directory the archive places the product in (§1.1). + + Read from the path because the `.gti` HDU carries no `FILTER` keyword — unlike `.lc` and + `.pi`, which do. Verified against the archive rather than assumed, and a path that does + not name a detector fails rather than defaulting to one. + """ + for part in path.parts[::-1]: + if part.upper() in ("SDD1", "SDD2"): + return part.upper() + _fits.fail( + "F-18", f"/{source}", + f"no SDD1/SDD2 directory in {path}; §1.1 places every SoLEXS product under one, and " + f"the .gti HDU carries no FILTER keyword to fall back on", + ) + + +def _declared_exposure(header, source: str) -> float: + """§2.3: HDU1 `EXPOSURE` is a **string** in the archive. Converted explicitly.""" + raw = _fits.keyword(header, "EXPOSURE", source=source, hdu_name="GTI") + try: + return float(raw) + except (TypeError, ValueError): + _fits.fail( + "F-09", f"/{source}#GTI/EXPOSURE", + f"EXPOSURE is {raw!r}, which is not a number. §2.3 records it as the string " + f"'86395.0'; a value that will not convert cannot be compared against F-09.", + ) + + +def _validate_ordering(intervals: tuple[Interval, ...], source: str) -> None: + """§2.3: rows sorted and non-overlapping.""" + for index in range(1, len(intervals)): + previous, current = intervals[index - 1], intervals[index] + if current.start <= previous.stop: + _fits.fail( + "F-09", f"/{source}#GTI/START[{index}]", + f"interval {index} starts at {current.start!r}, at or before the previous " + f"stop {previous.stop!r}. Overlapping intervals would double-count live time.", + ) + + +def _validate_within_day(primary, intervals: tuple[Interval, ...], source: str) -> None: + """§2.3: all intervals within [OBS_DATE 00:00, 23:59:59]. + + The primary `TSTART`/`TSTOP` are ISO-8601 strings (§2.3), so they are compared as + instants rather than as numbers. A non-parsing value is F-05 rather than a skipped check: + silently not checking is the failure mode this whole specification exists to prevent. + """ + raw_start = _fits.keyword(primary, "TSTART", source=source, hdu_name="PRIMARY") + raw_stop = _fits.keyword(primary, "TSTOP", source=source, hdu_name="PRIMARY") + try: + declared_start = Timestamp(str(raw_start).replace("+00:00", "Z")) + declared_stop = Timestamp(str(raw_stop).replace("+00:00", "Z")) + except Exception: + _fits.fail( + "F-05", f"/{source}#PRIMARY/TSTART", + f"primary TSTART/TSTOP are {raw_start!r}/{raw_stop!r}; §2.3 records them as " + f"ISO-8601 strings such as '2024-05-14T00:00:01+00:00'", + ) + + for index, interval in enumerate(intervals): + if interval.start_utc.instant < declared_start.instant: + _fits.fail( + "F-09", f"/{source}#GTI/START[{index}]", + f"interval {index} starts at {interval.start_utc} — before the declared " + f"TSTART {declared_start}", + ) + if interval.stop_utc.instant > declared_stop.instant: + _fits.fail( + "F-09", f"/{source}#GTI/STOP[{index}]", + f"interval {index} stops at {interval.stop_utc} — after the declared " + f"TSTOP {declared_stop}", + ) + + +__all__ = ["GoodTimeIntervals", "Interval", "parse"] diff --git a/contexts/ingest/parsers/solexs/lc.py b/contexts/ingest/parsers/solexs/lc.py new file mode 100644 index 0000000..e64cbb6 --- /dev/null +++ b/contexts/ingest/parsers/solexs/lc.py @@ -0,0 +1,298 @@ +"""SoLEXS `.lc.gz` — the total-band count time series. + +Implements `SPEC-parsers@r6` §2.1. That section is a contract, not a description: a deviation +requires a logged amendment, not a code change. + +THE TWO THINGS THIS PARSER MUST GET RIGHT +------------------------------------------ +**The values are counts, not a rate.** The HDU is named `RATE` and the values are counts per +1-second bin. `HDUCLAS3='COUNTS'` is what declares this, and F-07 requires the parser to use +`HDUCLAS3` and **never** `EXTNAME` to decide semantics. Reading the name would divide by a +second that was never applied and publish a number that is right by coincidence of magnitude. + +**NaN is data.** §2.1, as amended at r2: *"SoLEXS `COUNTS` uses NaN as the missing-data +sentinel… NaN values represent absent measurements. Zero remains a valid physical count and +MUST NOT be treated as missing."* The binding parser behaviour is that NaN values pass +through **unchanged**, are **never** imputed, **never** converted to zero, **never** removed. + +Here that maps exactly onto the domain: `Observation.value = None` means *observed to be +absent* (ADR-0017, L-07), and `0.0` means a measured zero. The mapping is +`NaN → None`, `0.0 → 0.0`, and the two are tested against each other in both directions +because collapsing them is the specific error that would be undetectable downstream. + +**Finiteness of `COUNTS` is deliberately NOT validated** (§2.1, r2): NaN is data, not an +error. What *is* validated is `TIME` finiteness — F-16 — because a NaN timestamp defeats a +monotonicity test silently, all NaN comparisons being False. + +WHAT THIS PARSER DOES NOT DO +---------------------------- +It does not check `NaN ⇒ GTI-excluded`. §2.1 places that at the day-assembly layer, which +needs both `.lc` and `.gti`: *"enforced at the day-assembly layer (Milestone VII), not inside +the single-file `.lc` parser."* It does not aggregate to a minute grid — §3's T1 contract is +the write path's, M3/E5/#19. It reads a clock never; `ingest_time` is supplied by the +acquisition that produced the file. +""" + +from __future__ import annotations + +import math +from collections.abc import Iterator +from dataclasses import dataclass +from datetime import datetime, timezone +from pathlib import Path + +from contexts.ingest.parsers.solexs import _fits +from domain.entities import Observation +from domain.values import Digest, Identifier, Timestamp + +#: Unix epoch as Modified Julian Date. `MJDREFI=40587` is what the archive declares, and +#: F-05 fires on anything else — v1 defaulted 58484 and produced a ~49-year error. +MJDREFI_UNIX = 40587 + +#: §2.1: one row per second of a UTC day. +EXPECTED_ROWS = 86_400 + +#: The quantity name is the instrument's own term. `HDUCLAS3` declares COUNTS, and the unit is +#: counts accumulated in one `TIMEDEL` second — not a rate, despite `EXTNAME='RATE'` (F-07). +QUANTITY = "counts" +UNIT = "counts" + + +@dataclass(frozen=True) +class LightCurveHeader: + """Every keyword §2.1 lists under "Metadata to capture", read strictly. + + Split across PRIMARY and the `RATE` extension exactly as the archive stores them — + `OBSERVED` on the 2024-05-14 SDD2 product. Reading a keyword from the wrong header would + require a fallback, and a fallback is the banned idiom (§5). + """ + + mission: str + telescope: str + instrument: str + origin: str + creator: str + filename: str + obs_date: str + obs_id: str + processing_date: str + detector: str + tstart: float + tstop: float + timedel: float + numband: str + rows: int + + +@dataclass(frozen=True) +class LightCurve: + """A parsed `.lc` product: its header, and the second-by-second series. + + `counts` holds `None` where the archive holds NaN. The conversion happens once, here, so + that no downstream consumer can encounter a raw NaN and decide for itself what it meant. + """ + + header: LightCurveHeader + source_path: Path + source_digest: Digest + times: tuple[float, ...] + counts: tuple[float | None, ...] + + @property + def n_finite(self) -> int: + """Seconds carrying a measurement. §2.1 `OBSERVED`: 86,395 on 2024-05-14 SDD2.""" + return sum(1 for value in self.counts if value is not None) + + @property + def n_absent(self) -> int: + """Seconds observed to be absent. §2.1 `OBSERVED`: 5 on 2024-05-14 SDD2.""" + return len(self.counts) - self.n_finite + + @property + def absent_offsets(self) -> tuple[int, ...]: + """Day-offsets of the absent seconds. §2.1 records [0, 5, 30072, 30078, 83951].""" + return tuple(i for i, value in enumerate(self.counts) if value is None) + + @property + def instrument_id(self) -> Identifier: + """`solexs-sdd2`. + + The detector is part of the instrument identity, not a label on it: SDD1 and SDD2 are + two physical SDD detectors, and §1.1 records that they carry different products — + SDD1 supplies GTI only. Collapsing them to `solexs` would make an SDD1 row and an + SDD2 row indistinguishable (ADR-0003: an Instrument is a sensor). + """ + return Identifier(f"solexs-{self.header.detector.lower()}") + + def observations(self, *, source_id: Identifier, ingest_time: Timestamp | None + ) -> Iterator[Observation]: + """One Observation per second, streamed (E5 §12). + + A generator, not a list: §2.2 records ~472 MB/day decompressed for the sibling `.pi` + product and the same discipline applies here — a day must not be materialised twice. + + `ingest_time` is supplied by the acquisition, never read from a clock here. `None` is + legitimate and means exactly what ADR-0022 says: unknown, predates bitemporal capture. + """ + for offset, (when, value) in enumerate(zip(self.times, self.counts)): + yield Observation( + source_id=source_id, + instrument_id=self.instrument_id, + quantity=QUANTITY, + unit=UNIT, + valid_time=unix_to_timestamp(when), + ingest_time=ingest_time, + value=value, + source_digest=self.source_digest, + quality_flags=() if value is not None else ("no_data",), + ) + del offset + + +def unix_to_timestamp(seconds: float) -> Timestamp: + """Unix seconds UTC → a domain Timestamp. + + §2.1 establishes the epoch from the file rather than by assumption: `MJDREFI=40587`, + `MJDREFF=0` → MJD 40587 = 1970-01-01 = the Unix epoch, with `TIMESYS='UTC'`, + `TIMEUNIT='s'`, `TIMZERO=0`. Verified there against `TSTART=1715644800.0` = + 2024-05-14T00:00:00Z, and re-verified by this module's tests. + """ + moment = datetime.fromtimestamp(seconds, tz=timezone.utc) + return Timestamp(moment.isoformat().replace("+00:00", "Z")) + + +def _validate_time_axis(times, source: str) -> None: + """F-16 and §2.1's validation list, in the order a failure is most diagnostic. + + Finiteness first: a NaN timestamp would defeat the monotonicity check silently, because + every comparison against NaN is False. §2.1 calls this out explicitly. + """ + for index, value in enumerate(times): + if not math.isfinite(value): + _fits.fail( + "F-16", f"/{source}#RATE/TIME[{index}]", + f"TIME[{index}] is not finite ({value!r}). A NaN timestamp defeats the " + f"monotonicity test silently — every comparison against NaN is False.", + ) + for index in range(1, len(times)): + step = times[index] - times[index - 1] + if step != 1: + _fits.fail( + "F-16", f"/{source}#RATE/TIME[{index}]", + f"TIME step at row {index} is {step!r}s, expected exactly 1s " + f"(TIMEDEL=1, §2.1)", + ) + + +def parse(path: Path, digest: Digest) -> LightCurve: + """Parse one `.lc.gz` product. Every §2.1 validation, no repair. + + `digest` is the content address the kernel minted for these exact bytes (ADR-0005). It is + carried onto every Observation, so a published number resolves to the file it came from. + """ + 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: + primary = _fits.hdu(hdul, "PRIMARY", source=source).header + rate_hdu = _fits.hdu(hdul, "RATE", source=source) + rate = rate_hdu.header + + # ── F-07: semantics come from HDUCLAS3, never from EXTNAME ────────────── + _fits.expect(rate, "HDUCLAS1", "LIGHTCURVE", source=source, hdu_name="RATE", rule="F-07") + _fits.expect(rate, "HDUCLAS2", "TOTAL", source=source, hdu_name="RATE", rule="F-07") + _fits.expect(rate, "HDUCLAS3", "COUNTS", source=source, hdu_name="RATE", rule="F-07") + + # ── F-05: the epoch is read, never defaulted ──────────────────────────── + _fits.expect(rate, "MJDREFI", MJDREFI_UNIX, source=source, hdu_name="RATE", rule="F-05") + _fits.expect(rate, "MJDREFF", 0, source=source, hdu_name="RATE", rule="F-05") + _fits.expect(rate, "TIMESYS", "UTC", source=source, hdu_name="RATE", rule="F-05") + _fits.expect(rate, "TIMEUNIT", "s", source=source, hdu_name="RATE", rule="F-05") + _fits.expect(rate, "TIMZERO", 0, source=source, hdu_name="RATE", rule="F-05") + + timedel = _fits.keyword(rate, "TIMEDEL", source=source, hdu_name="RATE") + if timedel != 1: + _fits.fail("F-05", f"/{source}#RATE/TIMEDEL", + f"TIMEDEL is {timedel!r}s; §2.1 declares 1s bins") + + rows = _fits.keyword(rate, "NAXIS2", source=source, hdu_name="RATE") + if rows != EXPECTED_ROWS: + _fits.fail("F-17", f"/{source}#RATE/NAXIS2", + f"NAXIS2 is {rows}, expected {EXPECTED_ROWS} for a SoLEXS daily " + f"product — a truncated product must not be parsed as a whole day") + + times = [float(v) for v in _fits.column(rate_hdu, "TIME", source=source, hdu_name="RATE")] + raw_counts = _fits.column(rate_hdu, "COUNTS", source=source, hdu_name="RATE") + + _validate_time_axis(times, source) + + tstart = float(_fits.keyword(rate, "TSTART", source=source, hdu_name="RATE")) + if times and tstart != times[0]: + _fits.fail("F-06", f"/{source}#RATE/TSTART", + f"TSTART {tstart!r} does not equal TIME[0] {times[0]!r}") + + # ── NaN passes through as "observed to be absent". Never zero, never dropped ── + counts: list[float | None] = [] + for index, value in enumerate(raw_counts): + number = float(value) + if math.isnan(number): + counts.append(None) + continue + if number < 0: + # F-19 covers negative counts only, and is inherently NaN-safe because + # `NaN < 0` is False — which is why this check follows the NaN branch. + _fits.fail("F-19", f"/{source}#RATE/COUNTS[{index}]", + f"COUNTS[{index}] is {number!r}; a negative count is physically " + f"impossible") + counts.append(number) + + detector = _fits.keyword(rate, "FILTER", source=source, hdu_name="RATE") + instrument = _fits.keyword(primary, "INSTRUME", source=source, hdu_name="PRIMARY") + if instrument != "SoLEXS": + _fits.fail("F-07", f"/{source}#PRIMARY/INSTRUME", + f"INSTRUME is {instrument!r}, expected 'SoLEXS'") + + header = LightCurveHeader( + mission=_fits.keyword(primary, "MISSION", source=source, hdu_name="PRIMARY"), + telescope=_fits.keyword(primary, "TELESCOP", source=source, hdu_name="PRIMARY"), + instrument=instrument, + origin=_fits.keyword(primary, "ORIGIN", source=source, hdu_name="PRIMARY"), + creator=_fits.keyword(primary, "CREATOR", source=source, hdu_name="PRIMARY"), + filename=_fits.keyword(primary, "FILENAME", source=source, hdu_name="PRIMARY"), + obs_date=str(_fits.keyword(primary, "OBS_DATE", source=source, hdu_name="PRIMARY")), + obs_id=str(_fits.keyword(primary, "OBS_ID", source=source, hdu_name="PRIMARY")), + processing_date=str( + _fits.keyword(primary, "DATE", source=source, hdu_name="PRIMARY") + ), + detector=detector, + tstart=tstart, + tstop=float(_fits.keyword(rate, "TSTOP", source=source, hdu_name="RATE")), + timedel=float(timedel), + numband=str(_fits.keyword(rate, "NUMBAND", source=source, hdu_name="RATE")), + rows=int(rows), + ) + + return LightCurve( + header=header, + source_path=path, + source_digest=digest, + times=tuple(times), + counts=tuple(counts), + ) + + +__all__ = [ + "EXPECTED_ROWS", + "LightCurve", + "LightCurveHeader", + "MJDREFI_UNIX", + "QUANTITY", + "UNIT", + "parse", + "unix_to_timestamp", +] diff --git a/contexts/ingest/parsers/solexs/pi.py b/contexts/ingest/parsers/solexs/pi.py new file mode 100644 index 0000000..c819bc5 --- /dev/null +++ b/contexts/ingest/parsers/solexs/pi.py @@ -0,0 +1,306 @@ +"""SoLEXS `.pi.gz` — per-second spectra. The scientific core. + +Implements `SPEC-parsers@r6` §2.2: OGIP Type II PHA, `DETCHANS=340`, `CHANTYPE='PI'`, +one 340-channel spectrum per second, `NAXIS2 = 86400`. + +THE BINDING PROHIBITION +----------------------- +§2.2 states it as a `CRITICAL GAP`, `OBSERVED`: + +> **no RMF/ARF response file exists anywhere in the archive.** SoLEXS PI-channel → keV +> conversion is therefore **impossible from archive contents alone**. **Binding rule: no v2 +> artifact may state a SoLEXS energy in keV until a response file is acquired.** + +`CHANNEL` is an **ordinal PI index, not energy**. This module therefore exposes channel +indices and never an energy, has no keV field, no conversion function and no calibration +constant — and `test_no_energy_is_stated_anywhere` asserts that none appears. The prohibition +is enforced by absence rather than by comment, because a comment is not a constraint. + +STREAMING IS A DESIGN REQUIREMENT, NOT A PREFERENCE +---------------------------------------------------- +§2.2, marked design-binding: `NAXIS1 = 5468 B × 86400` ≈ **472 MB/day decompressed**, and 436 +days ≈ **206 GB**. *"The parser MUST stream day-by-day and MUST NOT hold multiple days of +`COUNTS` in memory (16 GB RAM)."* + +So `parse` returns a header-and-validation result, and `spectra()` yields rows lazily from a +freshly opened file. Nothing here accumulates a day's `COUNTS`. + +`CHANNEL` IS READ ONCE AND DISCARDED +------------------------------------ +§2.2: *"`CHANNEL` is 235 MB/day of pure redundancy — read once, validate constant, discard."* +F-08 fires if the vector is not identical across rows, because a varying channel map voids +the assumption that a channel index means the same thing in every row. The constant map is +kept once, on the parse result. + +WHAT THIS PARSER DOES NOT DO +---------------------------- +No aggregation — §3's T2 contract belongs to the write path (#19). No merging with any other +channel space: F-11 forbids combining SoLEXS PI(340) with HEL1OS CZT PHA(341) or CdTe +PHA(511), and this module cannot do so because it knows nothing of HEL1OS. The cross-product +check `.pi TSTART[0] == .lc TSTART` (V-PI-3, F-06) is offered as a function taking both, since +a single-file parser cannot perform it alone. +""" + +from __future__ import annotations + +import math +from collections.abc import Iterator +from dataclasses import dataclass +from pathlib import Path + +from contexts.ingest.parsers.solexs import _fits +from contexts.ingest.parsers.solexs.lc import LightCurve, unix_to_timestamp +from domain.values import Digest, Identifier, Timestamp + +#: §2.2, `OBSERVED`. Not a configurable expectation: a product declaring a different channel +#: count is a different product, and F-11 forbids treating two channel spaces as one. +DETCHANS = 340 + +#: One spectrum per second of a UTC day. +EXPECTED_ROWS = 86_400 + + +@dataclass(frozen=True) +class Spectrum: + """One second's 340-channel spectrum. + + `counts` is per-channel counts accumulated over `exposure` seconds. Channels are ordinal + PI indices; the class carries no energy and there is nowhere to put one. + """ + + tstart: float + telapse: float + spec_num: int + exposure: float + counts: tuple[float, ...] + + @property + def valid_time(self) -> Timestamp: + return unix_to_timestamp(self.tstart) + + @property + def total_counts(self) -> float: + return math.fsum(self.counts) + + +@dataclass(frozen=True) +class SpectraHeader: + """The OGIP Type II PHA declarations §2.2 requires, read strictly.""" + + mission: str + telescope: str + instrument: str + origin: str + creator: str + filename: str + obs_date: str + obs_id: str + processing_date: str + detector: str + detchans: int + chantype: str + poisserr: bool + areascal: float + corrscal: float + rows: int + + +@dataclass(frozen=True) +class Spectra: + """A validated `.pi` product. Rows are streamed, never held. + + `channel_map` is the constant ordinal channel vector, validated once (F-08) and kept once + rather than 86,400 times. + """ + + header: SpectraHeader + source_path: Path + source_digest: Digest + channel_map: tuple[int, ...] + first_tstart: float + + @property + def instrument_id(self) -> Identifier: + return Identifier(f"solexs-{self.header.detector.lower()}") + + def spectra(self) -> Iterator[Spectrum]: + """Yield each second's spectrum, one at a time, from a freshly opened file. + + A generator over a re-opened file rather than a retained handle: §2.2's memory + constraint is about a day never being materialised, and a consumer that abandons the + iteration should not leave a 472 MB mapping alive. + """ + 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) + tstarts = _fits.column(table, "TSTART", source=source, hdu_name="SPECTRUM") + telapses = _fits.column(table, "TELAPSE", source=source, hdu_name="SPECTRUM") + spec_nums = _fits.column(table, "SPEC_NUM", 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") + + 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", + ) + yield Spectrum( + tstart=float(tstarts[index]), + telapse=float(telapses[index]), + spec_num=int(spec_nums[index]), + exposure=float(exposures[index]), + counts=row, + ) + + +def parse(path: Path, digest: Digest) -> Spectra: + """Validate one `.pi.gz` product and capture its constant channel map. + + The row data is **not** read here. §2.2 makes streaming design-binding, and a `parse` that + returned 86,400 spectra would have already violated it before any caller could choose. + """ + 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: + primary = _fits.hdu(hdul, "PRIMARY", source=source).header + table = _fits.hdu(hdul, "SPECTRUM", source=source) + spectrum = table.header + + # ── OGIP Type II PHA, declared not assumed (§2.2) ─────────────────────── + _fits.expect(spectrum, "HDUCLAS1", "SPECTRUM", source=source, + hdu_name="SPECTRUM", rule="F-07") + _fits.expect(spectrum, "HDUCLAS3", "COUNTS", source=source, + hdu_name="SPECTRUM", rule="F-07") + _fits.expect(spectrum, "HDUCLAS4", "TYPE:II", source=source, + hdu_name="SPECTRUM", rule="F-07") + _fits.expect(spectrum, "CHANTYPE", "PI", source=source, + hdu_name="SPECTRUM", rule="F-07") + + detchans = int(_fits.keyword(spectrum, "DETCHANS", source=source, + hdu_name="SPECTRUM")) + if detchans != DETCHANS: + _fits.fail( + "F-11", f"/{source}#SPECTRUM/DETCHANS", + f"DETCHANS is {detchans}, expected {DETCHANS}. A different channel count is " + f"a different channel space, and F-11 forbids treating two as one.", + ) + + # Channel space before row count, deliberately. A product declaring 341 channels is + # a different channel space (F-11) whether or not it also has the right number of + # rows, and reporting the completeness defect first would name the less specific + # problem. Both checks run; only the order of reporting is chosen. + channel_column = _fits.column(table, "CHANNEL", source=source, hdu_name="SPECTRUM") + channel_map = _constant_channel_map(channel_column, detchans, source) + + rows = int(_fits.keyword(spectrum, "NAXIS2", source=source, hdu_name="SPECTRUM")) + if rows != EXPECTED_ROWS: + _fits.fail("F-17", f"/{source}#SPECTRUM/NAXIS2", + f"NAXIS2 is {rows}, expected {EXPECTED_ROWS}") + + tstarts = _fits.column(table, "TSTART", source=source, hdu_name="SPECTRUM") + first_tstart = float(tstarts[0]) + + exposures = _fits.column(table, "EXPOSURE", source=source, hdu_name="SPECTRUM") + for index, value in enumerate(exposures): + if float(value) < 0: + _fits.fail("F-19", f"/{source}#SPECTRUM/EXPOSURE[{index}]", + f"EXPOSURE[{index}] is {float(value)!r}; a negative exposure is " + f"physically impossible") + + instrument = _fits.keyword(primary, "INSTRUME", source=source, hdu_name="PRIMARY") + if instrument != "SoLEXS": + _fits.fail("F-07", f"/{source}#PRIMARY/INSTRUME", + f"INSTRUME is {instrument!r}, expected 'SoLEXS'") + + header = SpectraHeader( + mission=_fits.keyword(primary, "MISSION", source=source, hdu_name="PRIMARY"), + telescope=_fits.keyword(primary, "TELESCOP", source=source, hdu_name="PRIMARY"), + instrument=instrument, + origin=_fits.keyword(primary, "ORIGIN", source=source, hdu_name="PRIMARY"), + creator=_fits.keyword(primary, "CREATOR", source=source, hdu_name="PRIMARY"), + filename=_fits.keyword(primary, "FILENAME", source=source, hdu_name="PRIMARY"), + obs_date=str(_fits.keyword(primary, "OBS_DATE", source=source, + hdu_name="PRIMARY")), + obs_id=str(_fits.keyword(primary, "OBS_ID", source=source, hdu_name="PRIMARY")), + processing_date=str(_fits.keyword(primary, "DATE", source=source, + hdu_name="PRIMARY")), + detector=_fits.keyword(spectrum, "FILTER", source=source, hdu_name="SPECTRUM"), + detchans=detchans, + chantype="PI", + poisserr=bool(_fits.keyword(spectrum, "POISSERR", source=source, + hdu_name="SPECTRUM")), + areascal=float(_fits.keyword(spectrum, "AREASCAL", source=source, + hdu_name="SPECTRUM")), + corrscal=float(_fits.keyword(spectrum, "CORRSCAL", source=source, + hdu_name="SPECTRUM")), + rows=rows, + ) + + return Spectra( + header=header, + source_path=path, + source_digest=digest, + channel_map=channel_map, + first_tstart=first_tstart, + ) + + +def _constant_channel_map(column, detchans: int, source: str) -> tuple[int, ...]: + """F-08: every row's CHANNEL vector must be identical. + + Read once and discarded (§2.2 — 235 MB/day of pure redundancy). A varying map would void + the assumption that channel *n* means the same thing in every row, which is the + assumption every downstream spectral operation rests on. + """ + 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. The channel map must be " + f"constant, or a channel index does not mean the same thing in every row.", + ) + return first + + +def check_epoch_agrees_with_lightcurve(spectra: Spectra, lightcurve: LightCurve) -> None: + """V-PI-3 / F-06: `.pi TSTART[0]` must equal the `.lc` `TSTART` for the same day. + + A cross-product check, so it takes both products rather than living inside either parser + — the same placement §2.1 uses for the `NaN ⇒ GTI-excluded` implication. Two products of + one day disagreeing about when the day began is an ambiguous time, and §5 admits no + ambiguous time. + """ + if spectra.first_tstart != lightcurve.header.tstart: + _fits.fail( + "F-06", f"/{spectra.source_path.name}#SPECTRUM/TSTART[0]", + f".pi TSTART[0] is {spectra.first_tstart!r} but .lc TSTART is " + f"{lightcurve.header.tstart!r} for the same day. V-PI-3 requires equality.", + ) + + +__all__ = [ + "DETCHANS", + "EXPECTED_ROWS", + "Spectra", + "SpectraHeader", + "Spectrum", + "check_epoch_agrees_with_lightcurve", + "parse", +] diff --git a/contexts/ingest/tests/solexs_fixtures.py b/contexts/ingest/tests/solexs_fixtures.py new file mode 100644 index 0000000..89d57e7 --- /dev/null +++ b/contexts/ingest/tests/solexs_fixtures.py @@ -0,0 +1,196 @@ +"""FITS fixtures with the layouts `SPEC-parsers@r6` §2.1–§2.3 record. + +Real FITS files, written by astropy, carrying the headers and columns the archive carries — +built so the parsers can be exercised without the 3.5 GB corpus, and so a *violating* product +can be constructed deliberately. The archive contains no malformed product to test against, +and a fail-loud rule nobody has watched fire is a comment. + +Every default here is a value `OBSERVED` in the spec, not a plausible-looking invention: +`MJDREFI=40587`, `TIMEDEL=1`, `NAXIS2=86400`, `DETCHANS=340`, `CHANTYPE='PI'`, +`HDUCLAS3='COUNTS'`, `FILTER='SDD2'`, and the 2024-05-14 epoch `TSTART=1715644800.0`. + +The fixtures are deliberately **short** by default — 60 seconds rather than 86,400 — because +a test that writes a full day writes 472 MB for `.pi`. Where a test needs the real row count +it asks for it, and where it needs to prove the row-count check fires it asks for a wrong one. +""" + +from __future__ import annotations + +import math +from pathlib import Path + +import numpy as np +from astropy.io import fits + +#: 2024-05-14T00:00:00Z — the epoch §2.1 verifies `TSTART` against. +FLARE_DAY_TSTART = 1_715_644_800.0 + +#: §2.1 `OBSERVED` on 2024-05-14 SDD2: NaN at these day-offsets, identical to the +#: GTI-excluded set. Reproduced here so a fixture can carry the same shape at small scale. +ABSENT_OFFSETS = (0, 5, 30_072, 30_078, 83_951) + + +def primary(**overrides): + """The PRIMARY header, with every keyword §2.1 lists under "Metadata to capture".""" + header = fits.Header() + values = { + "MISSION": "ADITYA-L1", "TELESCOP": "AL1", "INSTRUME": "SoLEXS", + "ORIGIN": "ISSDC", "CREATOR": "solexs_pipeline-1.4", + "FILENAME": "AL1_SOLEXS_20240514_SDD2_L1.lc", "OBS_DATE": "2024-05-14", + "OBS_ID": "AL1_SLX_L1_20240514", "DATE": "2024-06-16", + } + values.update(overrides) + for key, value in values.items(): + if value is not None: + header[key] = value + return header + + +def write_lc(path: Path, *, rows: int = 86_400, tstart: float = FLARE_DAY_TSTART, + counts: list[float] | None = None, times: list[float] | None = None, + rate_overrides: dict | None = None, + primary_overrides: dict | None = None, + extname: str = "RATE") -> Path: + """A `.lc` product per §2.1.""" + if times is None: + times = [tstart + i for i in range(rows)] + if counts is None: + counts = [100.0 + i for i in range(rows)] + + table = fits.BinTableHDU.from_columns([ + fits.Column(name="TIME", format="D", array=np.array(times, dtype=float)), + fits.Column(name="COUNTS", format="D", array=np.array(counts, dtype=float)), + ], name=extname) + + header = table.header + for key, value in { + "HDUCLAS1": "LIGHTCURVE", "HDUCLAS2": "TOTAL", "HDUCLAS3": "COUNTS", + "MJDREFI": 40587, "MJDREFF": 0, "TIMESYS": "UTC", "TIMEUNIT": "s", + "TIMEDEL": 1, "TIMZERO": 0, "TSTART": tstart, + "TSTOP": tstart + rows - 1, "FILTER": "SDD2", "NUMBAND": "4", + }.items(): + header[key] = value + for key, value in (rate_overrides or {}).items(): + if value is None: + del header[key] + else: + header[key] = value + + path.parent.mkdir(parents=True, exist_ok=True) + fits.HDUList([ + fits.PrimaryHDU(header=primary(**(primary_overrides or {}))), table, + ]).writeto(path, overwrite=True) + return path + + +def write_gti(path: Path, *, intervals: list[tuple[float, float]] | None = None, + exposure: str | float | None = None, + tstart: str = "2024-05-14T00:00:00+00:00", + tstop: str = "2024-05-14T23:59:59+00:00", + gti_overrides: dict | None = None) -> Path: + """A `.gti` product per §2.3, with `EXPOSURE` as the string the archive stores. + + `exposure` defaults to `Σ(STOP−START+1)`, the inclusive convention amended at r1, so a + fixture is F-09-consistent unless a test deliberately makes it otherwise. + """ + if intervals is None: + intervals = [(FLARE_DAY_TSTART + 1, FLARE_DAY_TSTART + 86_399)] + if exposure is None: + exposure = str(float(sum(stop - start + 1 for start, stop in intervals))) + + starts = np.array([a for a, _ in intervals], dtype=float) + stops = np.array([b for _, b in intervals], dtype=float) + table = fits.BinTableHDU.from_columns([ + fits.Column(name="START", format="D", array=starts), + fits.Column(name="STOP", format="D", array=stops), + ], name="GTI") + if exposure is not None: + table.header["EXPOSURE"] = exposure + for key, value in (gti_overrides or {}).items(): + if value is None: + del table.header[key] + else: + table.header[key] = value + + head = primary(FILENAME="AL1_SOLEXS_20240514_SDD2_L1.gti") + if tstart is not None: + head["TSTART"] = tstart + if tstop is not None: + head["TSTOP"] = tstop + + path.parent.mkdir(parents=True, exist_ok=True) + fits.HDUList([fits.PrimaryHDU(header=head), table]).writeto(path, overwrite=True) + return path + + +def write_pi(path: Path, *, rows: int = 60, detchans: int = 340, + tstart: float = FLARE_DAY_TSTART, + channel_map: list[list[int]] | None = None, + counts: np.ndarray | None = None, + exposures: list[float] | None = None, + spectrum_overrides: dict | None = None) -> Path: + """A `.pi` product per §2.2 — OGIP Type II PHA.""" + 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)) + if exposures is None: + exposures = [1.0] * rows + + table = fits.BinTableHDU.from_columns([ + fits.Column(name="TSTART", format="D", unit="s", + array=np.array([tstart + i for i in range(rows)], dtype=float)), + fits.Column(name="TELAPSE", format="D", unit="s", + array=np.ones(rows, dtype=float)), + fits.Column(name="SPEC_NUM", format="J", array=np.arange(rows, dtype=np.int32)), + fits.Column(name="CHANNEL", format=f"{detchans}K", + array=np.array(channel_map, dtype=np.int64)), + fits.Column(name="COUNTS", format=f"{detchans}D", array=counts), + fits.Column(name="EXPOSURE", format="D", unit="s", + array=np.array(exposures, dtype=float)), + ], name="SPECTRUM") + + header = table.header + for key, value in { + "HDUCLAS1": "SPECTRUM", "HDUCLAS2": "TOTAL", "HDUCLAS3": "COUNTS", + "HDUCLAS4": "TYPE:II", "CHANTYPE": "PI", "DETCHANS": detchans, + "POISSERR": False, "AREASCAL": 1.0, "CORRSCAL": 1.0, "FILTER": "SDD2", + }.items(): + header[key] = value + for key, value in (spectrum_overrides or {}).items(): + if value is None: + del header[key] + else: + header[key] = value + + path.parent.mkdir(parents=True, exist_ok=True) + fits.HDUList([ + fits.PrimaryHDU(header=primary(FILENAME="AL1_SOLEXS_20240514_SDD2_L1.pi")), table, + ]).writeto(path, overwrite=True) + return path + + +def sdd2(root: Path) -> Path: + directory = root / "AL1_SLX_L1_20240514_v1.0" / "AL1_SLX_L1_20240514_v1.0" / "SDD2" + directory.mkdir(parents=True, exist_ok=True) + return directory + + +def sdd1(root: Path) -> Path: + directory = root / "AL1_SLX_L1_20240514_v1.0" / "AL1_SLX_L1_20240514_v1.0" / "SDD1" + directory.mkdir(parents=True, exist_ok=True) + return directory + + +def nan_at(rows: int, offsets: tuple[int, ...]) -> list[float]: + """Counts with NaN at the given offsets — the archive's missing-data sentinel (§2.1 r2). + + Zero appears deliberately at an offset that is *not* absent, so any test asserting on + absence also has a measured zero present to be confused with. + """ + values = [float(i) for i in range(rows)] + values[1] = 0.0 + for offset in offsets: + if offset < rows: + values[offset] = math.nan + return values diff --git a/contexts/ingest/tests/test_solexs_parsers.py b/contexts/ingest/tests/test_solexs_parsers.py new file mode 100644 index 0000000..93421bc --- /dev/null +++ b/contexts/ingest/tests/test_solexs_parsers.py @@ -0,0 +1,589 @@ +"""SoLEXS parsers — per-field vs spec, and the no-imputation property. + +`| 17 | 17 | E5 | 900 | L | 15,12 | per-field vs spec | no imputation | parse fixtures |` + +Every assertion below cites the clause of `SPEC-parsers@r6` it enforces. Two halves: + +**Per-field vs spec** — each declared keyword, column, unit and count is checked, and each +fail-loud rule (§5) is fired deliberately against a product built to violate exactly it. The +archive contains no malformed product, so a violating fixture is the only way to watch a rule +reject anything. + +**No imputation** — the Property column. §2.1 as amended at r2 is binding: NaN passes through +unchanged, is never imputed, never converted to zero, never removed. The property is +quantified over every position and both directions, because collapsing NaN and 0.0 is the one +error that would be undetectable in every downstream artifact. +""" + +from __future__ import annotations + +import math +from pathlib import Path + +import numpy as np +import pytest + +from contexts.ingest.parsers.solexs import gti, lc, pi +from contexts.ingest.tests import solexs_fixtures as fx +from domain.errors import ContractViolation +from domain.values import Digest, Identifier, Timestamp + +DIGEST = Digest("a" * 64) +SOURCE = Identifier("issdc-pradan") + + +#: Fixtures are written UNCOMPRESSED. astropy reads both, the parsers never inspect the +#: extension, and gzipping every fixture costs 1.26 s against 0.01 s — 126× for no coverage. +#: The real `.gz` path is exercised by `test_a_gzipped_product_parses_identically` below and +#: by every real-corpus test, which read the archive's own compressed products. +def parsed_lc(tmp_path, **kwargs): + return lc.parse(fx.write_lc(fx.sdd2(tmp_path) / "x.lc", **kwargs), DIGEST) + + +def parsed_gti(tmp_path, **kwargs): + return gti.parse(fx.write_gti(fx.sdd2(tmp_path) / "x.gti", **kwargs), DIGEST) + + +def parsed_pi(tmp_path, **kwargs): + """A `.pi` product for the NEGATIVE cases, deliberately short. + + Short is safe here because §2.2's declarations and channel space are checked before the + row count, so a fixture violating F-07/F-08/F-11 fails on the rule under test. Positive + cases use `full_day_pi`, which is full size and written once per session. + """ + return pi.parse(fx.write_pi(fx.sdd2(tmp_path) / "x.pi", **kwargs), DIGEST) + + +@pytest.fixture(scope="session") +def full_day_pi(tmp_path_factory): + """One real-size `.pi`: 86,400 spectra × 340 channels. + + Session-scoped because writing it costs ~11 s. §2.2 records ~472 MB/day decompressed, so + this is the one place the suite pays that price, and it pays it once. + """ + directory = tmp_path_factory.mktemp("pi_full") + return pi.parse(fx.write_pi(directory / "full.pi", rows=86_400), DIGEST) + + +def violation(caught) -> str: + return caught.value.message + + +# ═══════════════════════════════════════════════════ §2.1 — .lc per-field + + +def test_the_lightcurve_parses_every_declared_field(tmp_path): + """§2.1's "Metadata to capture" list, read from the headers that actually carry it.""" + product = parsed_lc(tmp_path) + header = product.header + assert header.mission == "ADITYA-L1" + assert header.telescope == "AL1" + assert header.instrument == "SoLEXS" + assert header.origin == "ISSDC" + assert header.creator == "solexs_pipeline-1.4" + assert header.obs_date == "2024-05-14" + assert header.detector == "SDD2" + assert header.timedel == 1.0 + assert header.numband == "4" + assert header.rows == 86_400 + + +def test_the_epoch_is_unix_seconds_utc(tmp_path): + """§2.1: MJDREFI=40587 → MJD 40587 = 1970-01-01 = the Unix epoch. + + Verified against the archive's own anchor: TSTART=1715644800.0 = 2024-05-14T00:00:00Z. + """ + product = parsed_lc(tmp_path) + assert product.header.tstart == 1_715_644_800.0 + assert str(lc.unix_to_timestamp(product.header.tstart)) == "2024-05-14T00:00:00Z" + assert lc.MJDREFI_UNIX == 40587 + + +def test_the_values_are_counts_not_a_rate(tmp_path): + """F-07: semantics come from HDUCLAS3, never from EXTNAME. + + The HDU is called RATE and holds counts per 1-second bin. Reading the name would divide + by a second that was never applied. + """ + product = parsed_lc(tmp_path) + assert lc.QUANTITY == "counts" + assert lc.UNIT == "counts" + observations = list(product.observations(source_id=SOURCE, ingest_time=None)) + assert {o.unit for o in observations} == {"counts"} + assert {o.quantity for o in observations} == {"counts"} + + +def test_the_detector_is_part_of_the_instrument_identity(tmp_path): + """SDD1 and SDD2 are two physical detectors carrying different products (§1.1).""" + assert parsed_lc(tmp_path).instrument_id == Identifier("solexs-sdd2") + + +@pytest.mark.parametrize( + "key, value, rule", + [ + ("MJDREFI", 58484, "F-05"), # v1's default — a ~49-year error + ("MJDREFF", 1, "F-05"), + ("TIMESYS", "TT", "F-05"), + ("TIMEUNIT", "d", "F-05"), + ("TIMZERO", 100, "F-05"), + ("TIMEDEL", 2, "F-05"), + ("HDUCLAS1", "SPECTRUM", "F-07"), + ("HDUCLAS2", "NET", "F-07"), + ("HDUCLAS3", "RATE", "F-07"), + ], +) +def test_a_wrong_declaration_fails_loud(tmp_path, key, value, rule): + with pytest.raises(ContractViolation) as caught: + parsed_lc(tmp_path, rate_overrides={key: value}) + assert violation(caught).startswith(rule) + + +@pytest.mark.parametrize("key", ["MJDREFI", "TIMESYS", "TIMEUNIT", "TSTART", "NUMBAND"]) +def test_an_absent_keyword_is_never_defaulted(tmp_path, key): + """§5's opening rule: `header.get(K, default)` is banned for any meaningful key. + + v1 defaulted MJDREF and produced a ~49-year timestamp error that nothing detected. + """ + with pytest.raises(ContractViolation) as caught: + parsed_lc(tmp_path, rate_overrides={key: None}) + assert "F-05" in violation(caught) + + +def test_the_hdu_is_found_by_name_not_by_index(tmp_path): + """F-02: HDU order is not a contract.""" + with pytest.raises(ContractViolation) as caught: + parsed_lc(tmp_path, extname="LIGHTCURVE") + assert violation(caught).startswith("F-02") + assert "RATE" in violation(caught) + + +def test_a_truncated_product_is_not_parsed_as_a_whole_day(tmp_path): + """F-17: NAXIS2 must be 86400 for a SoLEXS daily product.""" + with pytest.raises(ContractViolation) as caught: + lc.parse(fx.write_lc(fx.sdd2(tmp_path) / "d.lc", rows=86_399, + rate_overrides={"NAXIS2": 86_399}), DIGEST) + assert violation(caught).startswith("F-17") + + +def test_a_full_day_of_rows_is_accepted(tmp_path): + """The positive case for F-17, so the check is not simply always-failing.""" + product = parsed_lc(tmp_path) + assert product.header.rows == lc.EXPECTED_ROWS == 86_400 + + +def test_a_nan_timestamp_is_rejected_before_monotonicity_is_tested(tmp_path): + """F-16, and §2.1's stated reason: every comparison against NaN is False, so a NaN + timestamp would defeat the monotonicity test silently.""" + times = [fx.FLARE_DAY_TSTART + i for i in range(86_400)] + times[7] = math.nan + with pytest.raises(ContractViolation) as caught: + parsed_lc(tmp_path, times=times) + assert violation(caught).startswith("F-16") + assert "not finite" in violation(caught) + + +@pytest.mark.parametrize("mutation", ["duplicate", "backwards", "gap"]) +def test_a_broken_time_axis_is_rejected(tmp_path, mutation): + """F-16: Δ must be exactly 1 s. A gap is not repaired — §2.1 declares TIMEDEL=1.""" + times = [fx.FLARE_DAY_TSTART + i for i in range(86_400)] + if mutation == "duplicate": + times[10] = times[9] + elif mutation == "backwards": + times[10] = times[9] - 5 + else: + times[10] = times[9] + 4 + with pytest.raises(ContractViolation) as caught: + parsed_lc(tmp_path, times=times) + assert violation(caught).startswith("F-16") + + +def test_tstart_must_equal_the_first_timestamp(tmp_path): + """F-06: a header and its data disagreeing about when the day began is an ambiguous time.""" + with pytest.raises(ContractViolation) as caught: + parsed_lc(tmp_path, rate_overrides={"TSTART": fx.FLARE_DAY_TSTART + 1}) + assert violation(caught).startswith("F-06") + + +def test_a_negative_count_is_physically_impossible(tmp_path): + """F-19.""" + counts = [10.0] * 86_400 + counts[3] = -1.0 + with pytest.raises(ContractViolation) as caught: + parsed_lc(tmp_path, counts=counts) + assert violation(caught).startswith("F-19") + + +def test_the_negative_count_check_is_nan_safe(tmp_path): + """§2.1: F-19 is inherently NaN-safe because `NaN < 0` is False. + + Asserted because the ordering of the two branches is what makes it so: if the sign check + ran before the NaN branch it would still pass, but for the wrong reason. + """ + counts = [10.0] * 86_400 + counts[3] = math.nan + product = parsed_lc(tmp_path, counts=counts) + assert product.counts[3] is None + + +def test_a_wrong_instrument_is_refused(tmp_path): + with pytest.raises(ContractViolation) as caught: + parsed_lc(tmp_path, primary_overrides={"INSTRUME": "HEL1OS"}) + assert violation(caught).startswith("F-07") + + +# ═══════════════════════════════════════════ NO IMPUTATION — the Property column + + +def test_nan_becomes_absent_and_never_zero(tmp_path): + """§2.1 r2, binding: NaN is the missing-data sentinel; zero remains a valid count. + + Both directions, in one product: the NaN positions become `None`, the measured zero stays + `0.0`, and neither becomes the other. Collapsing them is the error that would be + undetectable in every downstream artifact. + """ + counts = fx.nan_at(86_400, (0, 5, 30)) + product = parsed_lc(tmp_path, counts=counts) + + assert product.absent_offsets == (0, 5, 30) + for offset in (0, 5, 30): + assert product.counts[offset] is None + assert product.counts[1] == 0.0, "a measured zero must not become absent" + assert product.counts[1] is not None + + +def test_no_value_is_invented_at_any_position(tmp_path): + """Quantified over every position: parsed[i] is None iff the archive held NaN there.""" + counts = fx.nan_at(86_400, (0, 5, 30_072, 30_078, 83_951)) + product = parsed_lc(tmp_path, counts=counts) + + assert len(product.counts) == len(counts) + for index, (raw, parsed) in enumerate(zip(counts, product.counts)): + if math.isnan(raw): + assert parsed is None, f"position {index}: NaN was imputed to {parsed!r}" + else: + assert parsed == raw, f"position {index}: {raw!r} became {parsed!r}" + + +def test_absent_seconds_are_never_removed(tmp_path): + """§2.1: NaN values are never removed. The row count is preserved exactly.""" + counts = fx.nan_at(86_400, fx.ABSENT_OFFSETS) + product = parsed_lc(tmp_path, counts=counts) + + # §2.1 `OBSERVED` on 2024-05-14 SDD2: 86,395 finite, 5 absent, at these exact offsets. + assert len(product.counts) == 86_400 + assert product.absent_offsets == fx.ABSENT_OFFSETS + assert product.n_absent == 5 + assert product.n_finite == 86_395 + assert product.n_finite + product.n_absent == 86_400 + + +def test_absence_survives_into_the_observation(tmp_path): + """The domain's `None` means observed-to-be-absent (ADR-0017, L-07) — the same meaning. + + The mapping is checked at the boundary where it would be easiest to lose: a generator + that filtered or filled would produce a shorter or a different sequence. + """ + counts = fx.nan_at(86_400, (0, 5, 30)) + product = parsed_lc(tmp_path, counts=counts) + observations = list(product.observations(source_id=SOURCE, ingest_time=None)) + + assert len(observations) == 86_400 + assert [o.value for o in observations][:6] == [None, 0.0, 2.0, 3.0, 4.0, None] + assert observations[0].value_is_absent + assert not observations[1].value_is_absent + assert observations[0].quality_flags == ("no_data",) + assert observations[1].quality_flags == () + + +def test_no_interpolation_smoothing_or_repair_function_exists(): + """The prohibition enforced by absence rather than by comment. + + A module that *can* interpolate will interpolate. Nothing in these parsers can. + """ + for module in (lc, pi, gti): + names = {name.lower() for name in dir(module)} + for forbidden in ("interpolate", "interp", "smooth", "fillna", "fill", + "impute", "resample", "repair", "ffill", "bfill"): + assert not any(forbidden in name for name in names), ( + f"{module.__name__} exposes {forbidden!r}" + ) + + +# ═══════════════════════════════════════════════════ §2.3 — .gti per-field + + +def test_the_interval_convention_is_inclusive(tmp_path): + """§2.3 r1: live_time = STOP − START + 1. + + The amendment's whole content. Under the exclusive reading a 59-second interval would + measure 58 s and the F-09 equality could never hold. + """ + product = parsed_gti(tmp_path, intervals=[(fx.FLARE_DAY_TSTART + 1, fx.FLARE_DAY_TSTART + 59)], exposure="59.0") + assert product.intervals[0].live_time == 59.0 + assert product.live_time == 59.0 + + +def test_the_exposure_equality_is_exact(tmp_path): + """F-09, tolerance 0 s. §5: the relation is definitional and a tolerance would re-admit + the ambiguity that produced CONTRADICTION-001.""" + with pytest.raises(ContractViolation) as caught: + parsed_gti(tmp_path, intervals=[(fx.FLARE_DAY_TSTART + 1, fx.FLARE_DAY_TSTART + 59)], exposure="59.5") + assert violation(caught).startswith("F-09") + assert "tolerance 0" in violation(caught) + + +def test_even_a_one_second_discrepancy_is_fatal(tmp_path): + """The exclusive reading, off by exactly the number of intervals. Not tolerated.""" + with pytest.raises(ContractViolation) as caught: + parsed_gti(tmp_path, intervals=[(fx.FLARE_DAY_TSTART + 1, fx.FLARE_DAY_TSTART + 59)], exposure="58.0") + assert violation(caught).startswith("F-09") + + +def test_exposure_is_read_as_the_string_the_archive_stores(tmp_path): + """§2.3: HDU1 EXPOSURE is a string ('86395.0'). Converted explicitly, never coerced.""" + product = parsed_gti(tmp_path, intervals=[(fx.FLARE_DAY_TSTART + 1, fx.FLARE_DAY_TSTART + 59)], exposure="59.0") + assert product.declared_exposure == 59.0 + + +def test_an_unconvertible_exposure_fails_rather_than_falling_back(tmp_path): + """A computed fallback would make F-09 compare a number against itself.""" + with pytest.raises(ContractViolation) as caught: + parsed_gti(tmp_path, intervals=[(fx.FLARE_DAY_TSTART + 1, fx.FLARE_DAY_TSTART + 59)], exposure="unknown") + assert violation(caught).startswith("F-09") + + +def test_an_empty_gti_is_legal_and_marks_the_detector_inactive(tmp_path): + """F-12 — the single deliberate non-terminating rule. §1.1 records SDD1 at zero rows.""" + path = fx.write_gti(fx.sdd1(tmp_path) / "x.gti", intervals=[], exposure="0.0") + product = gti.parse(path, DIGEST) + assert product.detector_active is False + assert product.intervals == () + assert product.declared_exposure is None + assert product.detector == "SDD1" + + +def test_overlapping_intervals_are_refused(tmp_path): + """Overlap would double-count live time.""" + with pytest.raises(ContractViolation) as caught: + parsed_gti(tmp_path, intervals=[(fx.FLARE_DAY_TSTART + 1, fx.FLARE_DAY_TSTART + 59), + (fx.FLARE_DAY_TSTART + 50, fx.FLARE_DAY_TSTART + 99)]) + assert violation(caught).startswith("F-09") + + +def test_an_inverted_interval_is_refused(tmp_path): + with pytest.raises(ContractViolation) as caught: + parsed_gti(tmp_path, intervals=[(fx.FLARE_DAY_TSTART + 59, fx.FLARE_DAY_TSTART + 1)]) + assert violation(caught).startswith("F-09") + + +def test_an_interval_outside_the_declared_day_is_refused(tmp_path): + """§2.3: all intervals within [OBS_DATE 00:00, 23:59:59].""" + with pytest.raises(ContractViolation) as caught: + parsed_gti(tmp_path, intervals=[(fx.FLARE_DAY_TSTART - 100, fx.FLARE_DAY_TSTART - 41)]) + assert violation(caught).startswith("F-09") + + +def test_the_detector_comes_from_the_path_because_gti_carries_no_filter(tmp_path): + """Read from the directory §1.1 places the product in, and failed on rather than guessed.""" + assert parsed_gti(tmp_path).detector == "SDD2" + stray = tmp_path / "nowhere" / "x.gti" + fx.write_gti(stray) + with pytest.raises(ContractViolation) as caught: + gti.parse(stray, DIGEST) + assert violation(caught).startswith("F-18") + + +def test_excluded_seconds_are_reported_for_the_day_assembly_layer(tmp_path): + """§2.1 places `NaN ⇒ GTI-excluded` at day assembly (#19), which needs this.""" + product = parsed_gti( + tmp_path, + intervals=[(fx.FLARE_DAY_TSTART + 1, fx.FLARE_DAY_TSTART + 4), + (fx.FLARE_DAY_TSTART + 6, fx.FLARE_DAY_TSTART + 59)], + exposure="58.0", + ) + assert product.excluded_seconds(fx.FLARE_DAY_TSTART, 60) == (0, 5) + + +def test_interval_coverage_is_inclusive_at_both_ends(tmp_path): + interval = gti.Interval(1000.0, 1058.0) + assert interval.covers(1000.0) and interval.covers(1058.0) + assert not interval.covers(999.0) and not interval.covers(1059.0) + + +# ═══════════════════════════════════════════════════ §2.2 — .pi per-field + + +def test_the_spectra_product_declares_ogip_type_ii_pha(full_day_pi): + product = full_day_pi + assert product.header.detchans == pi.DETCHANS == 340 + assert product.header.chantype == "PI" + assert product.header.poisserr is False + assert product.header.areascal == 1.0 + assert product.header.corrscal == 1.0 + + +def test_the_channel_map_is_constant_and_kept_once(full_day_pi): + """§2.2: CHANNEL is 235 MB/day of pure redundancy — read once, validate, discard.""" + product = full_day_pi + assert product.channel_map == tuple(range(340)) + assert len(product.channel_map) == 340 + + +def test_a_varying_channel_map_voids_the_assumption_and_fails(tmp_path): + """F-08: a channel index must mean the same thing in every row.""" + varying = [list(range(340))] * 59 + [list(range(1, 341))] # short: F-08 precedes F-17 + with pytest.raises(ContractViolation) as caught: + parsed_pi(tmp_path, channel_map=varying) + assert violation(caught).startswith("F-08") + + +@pytest.mark.parametrize("detchans", [341, 511, 339]) +def test_a_different_channel_count_is_a_different_channel_space(tmp_path, detchans): + """F-11: SoLEXS PI(340), HEL1OS CZT PHA(341) and CdTe PHA(511) are incommensurable. + + 341 and 511 are chosen deliberately — they are the two spaces F-11 names. + """ + with pytest.raises(ContractViolation) as caught: + parsed_pi(tmp_path, detchans=detchans) + assert violation(caught).startswith("F-11") + + +def test_a_wrong_chantype_is_refused(tmp_path): + """PI is gain-corrected pulse-invariant; PHA is not the same space.""" + with pytest.raises(ContractViolation) as caught: + parsed_pi(tmp_path, spectrum_overrides={"CHANTYPE": "PHA"}) + assert violation(caught).startswith("F-07") + + +def test_the_spectra_stream_rather_than_materialise(full_day_pi): + """§2.2 is design-binding: ~472 MB/day decompressed; a day must not be held.""" + import inspect + + assert inspect.isgeneratorfunction(pi.Spectra.spectra) + first = next(iter(full_day_pi.spectra())) + assert len(first.counts) == 340 + assert first.valid_time == Timestamp("2024-05-14T00:00:00Z") + + +def test_every_spectrum_carries_its_second(full_day_pi): + seconds = [] + for index, spectrum in enumerate(full_day_pi.spectra()): + seconds.append(spectrum.tstart) + if index == 4: + break + assert seconds == [fx.FLARE_DAY_TSTART + i for i in range(5)] + + +def test_a_negative_spectral_count_is_refused(tmp_path): + counts = np.tile(np.arange(340, dtype=float), (86_400, 1)) + counts[4][17] = -1.0 + with pytest.raises(ContractViolation) as caught: + list(parsed_pi(tmp_path, rows=86_400, counts=counts).spectra()) + assert violation(caught).startswith("F-19") + + +def test_a_negative_exposure_is_refused(tmp_path): + exposures = [1.0] * 86_400 + exposures[9] = -1.0 + with pytest.raises(ContractViolation) as caught: + parsed_pi(tmp_path, rows=86_400, exposures=exposures) + assert violation(caught).startswith("F-19") + + +def test_the_pi_epoch_must_agree_with_the_lightcurve(tmp_path, full_day_pi): + """V-PI-3 / F-06: `.pi TSTART[0]` must equal the `.lc` TSTART for the same day.""" + curve = parsed_lc(tmp_path) + pi.check_epoch_agrees_with_lightcurve(full_day_pi, curve) + + shifted = parsed_lc(tmp_path, tstart=fx.FLARE_DAY_TSTART + 1) + with pytest.raises(ContractViolation) as caught: + pi.check_epoch_agrees_with_lightcurve(full_day_pi, shifted) + assert violation(caught).startswith("F-06") + + +# ═══════════════════════════════════════════ §2.2 — the keV prohibition + + +def test_no_energy_is_stated_anywhere(): + """§2.2, binding: no v2 artifact may state a SoLEXS energy in keV. + + No RMF/ARF exists anywhere in the archive, so PI-channel → keV is impossible from archive + contents alone. Channels are ordinal indices. Enforced by absence: there is no field, no + function and no constant to state an energy with. + """ + for holder in (pi, pi.Spectra, pi.Spectrum, pi.SpectraHeader): + names = {name.lower() for name in dir(holder)} + for forbidden in ("kev", "energy", "rmf", "arf", "response", "calibrat"): + assert not any(forbidden in name for name in names), ( + f"{getattr(holder, '__name__', holder)} exposes {forbidden!r}; §2.2 forbids " + f"stating a SoLEXS energy until a response file is acquired" + ) + + source = Path(pi.__file__).read_text().lower() + assert "def to_kev" not in source + assert "channel_energy" not in source + + +def test_the_channel_map_is_ordinal_not_physical(full_day_pi): + """The channel vector is 0..339 — indices, carrying no physical meaning by themselves.""" + product = full_day_pi + assert product.channel_map[0] == 0 + assert product.channel_map[-1] == 339 + assert all(isinstance(c, int) for c in product.channel_map) + + +# ═══════════════════════════════════════════ provenance and the clock + + +def test_every_product_carries_the_digest_of_the_bytes_it_came_from(tmp_path, full_day_pi): + """Every Observation must be traceable to the original acquired artifact (ADR-0005).""" + assert parsed_lc(tmp_path).source_digest == DIGEST + assert parsed_gti(tmp_path).source_digest == DIGEST + assert full_day_pi.source_digest == DIGEST + + +def test_the_observation_carries_the_source_digest(tmp_path): + product = parsed_lc(tmp_path) + for observation in product.observations(source_id=SOURCE, ingest_time=None): + assert observation.source_digest == DIGEST + break + + +def test_no_parser_reads_a_clock(): + """TIS §0.4: `ingest_time` is stamped at the acquisition boundary and nowhere else.""" + for module in (lc, pi, gti): + source = Path(module.__file__).read_text() + for forbidden in ("datetime.now", "utcnow", "time.time", "stamp()"): + assert forbidden not in source, f"{module.__name__} reads a clock via {forbidden}" + + +def test_the_ingest_time_is_supplied_not_invented(tmp_path): + """The parser propagates what the acquisition stamped, including a legitimate None.""" + product = parsed_lc(tmp_path) + stamped = Timestamp("2024-06-16T09:00:00Z") + + with_time = next(iter(product.observations(source_id=SOURCE, ingest_time=stamped))) + without = next(iter(product.observations(source_id=SOURCE, ingest_time=None))) + assert with_time.ingest_time == stamped + assert without.ingest_time is None + + +def test_a_gzipped_product_parses_identically(tmp_path): + """The archive ships `.gz`; the fixtures above are uncompressed for speed. + + This is the one place the two are compared, so the shortcut is proved harmless rather + than assumed to be: the same bytes, compressed and not, must parse to the same product. + """ + import gzip + import shutil + + plain = fx.write_lc(fx.sdd2(tmp_path) / "plain.lc", rows=86_400) + zipped = plain.with_suffix(".lc.gz") + with open(plain, "rb") as raw, gzip.open(zipped, "wb") as out: + shutil.copyfileobj(raw, out) + + from_plain = lc.parse(plain, DIGEST) + from_gzip = lc.parse(zipped, DIGEST) + assert from_plain.counts == from_gzip.counts + assert from_plain.times == from_gzip.times + assert from_plain.header == from_gzip.header diff --git a/tests/architecture/test_context_imports.py b/tests/architecture/test_context_imports.py index 6556898..f3e0d2f 100644 --- a/tests/architecture/test_context_imports.py +++ b/tests/architecture/test_context_imports.py @@ -64,6 +64,12 @@ CONTEXT_NAMES = [name for name, _, _ in CONTEXT_RULES] +#: This repository's own package roots. An `allow` entry outside this set is a third-party +#: library, which ADR-0026 does not govern and does not grant by default. +INTERNAL_ROOTS = frozenset( + {"kernel", "domain", "contexts", "contracts", "apps", "tools", "tests", "registry"} +) + def load_gate(root: Path): """Load the analyser with its repository root re-pointed at a fixture tree. @@ -124,11 +130,41 @@ def test_the_real_policy_set_covers_all_six_contexts(): for name, allow, citation in CONTEXT_RULES: policy = declared.get(f"contexts.{name}") assert policy is not None, f"contexts.{name} has no import policy ({citation})" - assert policy.allow == allow, ( - f"contexts.{name} permits {sorted(policy.allow)}, ADR-0026 grants {sorted(allow)}" + internal = policy.allow & INTERNAL_ROOTS + assert internal == allow, ( + f"contexts.{name} permits internal {sorted(internal)}, " + f"ADR-0026 grants {sorted(allow)}" ) +def test_every_third_party_grant_is_named_and_is_not_a_context(): + """A context may be granted a library. It may never be granted another context. + + Comparing only the internal roots above would let a third-party grant appear unnoticed, + so every non-internal grant is enumerated here as well. Adding one is then two deliberate + edits — POLICIES and this list — rather than one silent widening. M3/E5/#17 granted + `astropy` to Ingest because the SoLEXS products are FITS; ADR-0026 grants none by default + (STD-11), which is what #13 recorded when it declined to pre-grant any. + """ + from tools.gates.imports import GOVERNED_ROOTS, POLICIES + + granted: dict[str, set[str]] = {} + for policy in POLICIES: + if not policy.package.startswith("contexts."): + continue + third_party = policy.allow - INTERNAL_ROOTS + assert not (third_party & set(GOVERNED_ROOTS)), ( + f"{policy.package} grants {sorted(third_party & set(GOVERNED_ROOTS))}, which is " + f"a governed root of this repository, not a third-party library" + ) + if third_party: + granted[policy.package] = third_party + + assert granted == {"contexts.ingest": {"astropy"}}, ( + f"third-party grants changed: {granted}" + ) + + def test_no_context_policy_grants_the_contexts_root(): """ADR-0026 rule 5. Granting `contexts` would permit every context to import every other.""" from tools.gates.imports import POLICIES diff --git a/tests/integration/test_ingest_contract_with_platform.py b/tests/integration/test_ingest_contract_with_platform.py index 5ae730b..655d2d4 100644 --- a/tests/integration/test_ingest_contract_with_platform.py +++ b/tests/integration/test_ingest_contract_with_platform.py @@ -332,7 +332,13 @@ def test_the_ingest_context_imports_only_what_adr_0026_grants(): ingest = next(p for p in POLICIES if p.package == "contexts.ingest") assert ingest.populated is True, "the ingest policy still declares itself unpopulated" - assert ingest.allow == frozenset({"contracts", "domain", "kernel"}) + + # Internal grants are exactly ADR-0026's. `astropy` was added by M3/E5/#17 because the + # SoLEXS products are FITS; it is a library, not a context, which is the part that matters. + internal = ingest.allow & {"contracts", "domain", "kernel", "contexts", "apps", + "tools", "registry", "tests"} + assert internal == {"contracts", "domain", "kernel"} + assert "contexts" not in ingest.allow report, code = run(POLICIES) assert code == 0, report.violations diff --git a/tests/integration/test_issdc_pradan_acquire.py b/tests/integration/test_issdc_pradan_acquire.py index 1cc56ef..4f32609 100644 --- a/tests/integration/test_issdc_pradan_acquire.py +++ b/tests/integration/test_issdc_pradan_acquire.py @@ -338,7 +338,10 @@ def test_the_adapter_lives_within_the_ingest_import_rule(): assert code == 0, report.violations ingest = next(p for p in POLICIES if p.package == "contexts.ingest") - assert ingest.allow == frozenset({"contracts", "domain", "kernel"}) + internal = ingest.allow & {"contracts", "domain", "kernel", "contexts", "apps", + "tools", "registry", "tests"} + assert internal == {"contracts", "domain", "kernel"} + assert "contexts" not in ingest.allow, "Ingest must reach no other context" assert ingest.populated is True diff --git a/tests/integration/test_solexs_parse_with_platform.py b/tests/integration/test_solexs_parse_with_platform.py new file mode 100644 index 0000000..e2300a1 --- /dev/null +++ b/tests/integration/test_solexs_parse_with_platform.py @@ -0,0 +1,410 @@ +"""Real SoLEXS products become valid Observations, across the whole platform. + +`| 17 | 17 | E5 | 900 | L | 15,12 | per-field vs spec | no imputation | parse fixtures |` + +The path an Aditya-L1 X-ray measurement now takes, with each step performed by the component +that owns it: + + #16 ISSDC adapter resolves a date to the day's archive and its products + #10 provenance kernel mints the digest of every product, and records the derivation + #15 ingest contract admits the acquisition; supplies the one sanctioned clock read + #17 SoLEXS parsers canonicalise a product into second-by-second Observations + #12 domain model holds the row, and refuses a fabricated absence + #11 contract schemas validate the serialised Observation + #14 manifest records the Tier 0 archive as referenced, never deposited + #13 import rules say Ingest may do all of this and reach no other context + +TWO CLASSES OF TEST +------------------- +**Fixture tests** run everywhere. They use FITS products written by the test with the layout +`SPEC-parsers@r6` §2.1–2.3 records. + +**Real-corpus tests** are marked `real_only` and skip where the archive is absent (STD-12, +E5 §17). They re-measure the spec's own `OBSERVED` values against the 2024-05-14 SDD2 +products — the X8.7 flare day — so a claim in the specification and the behaviour of this +code are checked against each other rather than against nothing. + +The guard tests for `.gz` **products**, never for a directory: five v1.1 archive directories +exist in a clean checkout because their quicklook PNGs are tracked, while the 3.5 GB of +science products are not. +""" + +from __future__ import annotations + +import json +import math +from pathlib import Path + +import jsonschema +import pytest +from referencing import Registry, Resource + +from contexts.ingest import verify_conformance +from contexts.ingest.adapters.issdc_pradan import IssdcPradanAdapter +from contexts.ingest.parsers.solexs import gti, lc, pi +from contexts.ingest.tests import solexs_fixtures as fx +from domain.entities import Observation +from domain.invariants import ( + absence_survives_serialisation, + ingest_time_is_not_backfilled, + 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" +REAL_ARCHIVE = REPO_ROOT / "research" / "data" / "aditya_l1" / "real_l1_v1" / "solexs" + +FLARE_DAY = "2024-05-14" +#: §2.1 / §2.3 `OBSERVED` on 2024-05-14 SDD2. Re-measured by the real-corpus tests below. +SPEC_ABSENT_OFFSETS = (0, 5, 30_072, 30_078, 83_951) +SPEC_EXPOSURE = 86_395.0 + + +def real_products_present() -> bool: + if not REAL_ARCHIVE.is_dir(): + return False + return any(p for p in REAL_ARCHIVE.rglob("*.gz") if not p.name.startswith("._")) + + +real_only = pytest.mark.skipif( + not real_products_present(), 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 real_product(kind: str) -> Path: + stem = REAL_ARCHIVE / "AL1_SLX_L1_20240514_v1.0" / "AL1_SLX_L1_20240514_v1.0" + return stem / "SDD2" / f"AL1_SOLEXS_20240514_SDD2_L1.{kind}.gz" + + +# ═══════════════════════════════════════════ the full path, on fixtures + + +@pytest.fixture +def acquired(tmp_path): + """A conforming acquisition of a fixture archive carrying real-shaped products.""" + root = tmp_path / "archive" + directory = fx.sdd2(root) + fx.write_lc(directory / "AL1_SOLEXS_20240514_SDD2_L1.lc.gz", + counts=fx.nan_at(86_400, SPEC_ABSENT_OFFSETS)) + fx.write_gti( + directory / "AL1_SOLEXS_20240514_SDD2_L1.gti.gz", + intervals=_intervals_excluding(SPEC_ABSENT_OFFSETS), + ) + store = ProvenanceStore(tmp_path / "store") + adapter = IssdcPradanAdapter(root, store) + return store, adapter, verify_conformance(adapter, selector=FLARE_DAY) + + +def _intervals_excluding(offsets: tuple[int, ...]) -> list[tuple[float, float]]: + """GTI intervals covering the day except the given second-offsets. + + Built to the inclusive convention (§2.3 r1) so `Σ(STOP−START+1)` equals the covered + second count exactly — the fixture is F-09-consistent by construction rather than by a + tolerance. + """ + excluded = set(offsets) + intervals: list[tuple[float, float]] = [] + start: int | None = None + for offset in range(86_400): + if offset in excluded: + if start is not None: + intervals.append((fx.FLARE_DAY_TSTART + start, + fx.FLARE_DAY_TSTART + offset - 1)) + start = None + elif start is None: + start = offset + if start is not None: + intervals.append((fx.FLARE_DAY_TSTART + start, fx.FLARE_DAY_TSTART + 86_399)) + return intervals + + +def test_an_acquired_product_becomes_contract_valid_observations(acquired, tmp_path): + """The whole path: acquire → digest → parse → Observation → contract → provenance.""" + store, adapter, acquisition = acquired + product = adapter.archive_for(FLARE_DAY) / "AL1_SLX_L1_20240514_v1.0" / "SDD2" / \ + "AL1_SOLEXS_20240514_SDD2_L1.lc.gz" + + # The kernel mints the product's digest — the parser never could (ADR-0005). + minted = digest_file(product) + curve = lc.parse(product, Digest(minted.hex)) + + observations = product_observations(curve, acquisition) + assert len(observations) == 86_400 + + validator = validator_for("observation") + for observation in observations[:200]: + validator.validate(observation.to_dict()) + assert observation_is_wellformed(observation) + + # Every row traces to the exact bytes it was parsed from. + assert {str(o.source_digest) for o in observations} == {minted.hex} + + # And the derivation is recorded, so a published number is answerable to the archive. + run = begin_run(context="ingest", event="parse") + canonical = store.put_bytes( + json.dumps([o.to_dict() for o in observations[:10]], sort_keys=True).encode() + ) + store.record(run, inputs=[KernelDigest(str(acquisition.artifact.digest))], + outputs=[canonical.digest]) + assert KernelDigest(str(acquisition.artifact.digest)) in store.ancestors(canonical.digest) + + +def product_observations(curve, acquisition) -> list[Observation]: + return list(curve.observations( + source_id=acquisition.provenance.source_id, + ingest_time=acquisition.provenance.ingest_time, + )) + + +def test_absence_survives_the_whole_path(acquired, tmp_path): + """The no-imputation property, checked at the far end rather than at the parser. + + A NaN second becomes `value = None`, validates against the contract as `null`, survives + serialisation, and is still distinguishable from a measured zero. Any layer that filled + it would break exactly one of those four. + """ + _, adapter, acquisition = acquired + product = adapter.archive_for(FLARE_DAY) / "AL1_SLX_L1_20240514_v1.0" / "SDD2" / \ + "AL1_SOLEXS_20240514_SDD2_L1.lc.gz" + curve = lc.parse(product, Digest(digest_file(product).hex)) + observations = product_observations(curve, acquisition) + + validator = validator_for("observation") + for offset in SPEC_ABSENT_OFFSETS: + absent = observations[offset] + assert absent.value is None + assert absent.value_is_absent + document = absent.to_dict() + validator.validate(document) + assert document["value"] is None + assert absence_survives_serialisation(absent) + + measured_zero = observations[1] + assert measured_zero.value == 0.0 + assert not measured_zero.value_is_absent + validator.validate(measured_zero.to_dict()) + assert measured_zero.to_dict()["value"] == 0.0 + + +def test_the_ingest_time_comes_from_the_acquisition_not_the_parser(acquired, tmp_path): + """ADR-0004/0022: the parser propagates a stamp; it never reads a clock or backfills.""" + _, adapter, acquisition = acquired + product = adapter.archive_for(FLARE_DAY) / "AL1_SLX_L1_20240514_v1.0" / "SDD2" / \ + "AL1_SOLEXS_20240514_SDD2_L1.lc.gz" + curve = lc.parse(product, Digest(digest_file(product).hex)) + + stamped = product_observations(curve, acquisition)[100] + assert stamped.ingest_time == acquisition.provenance.ingest_time + assert ingest_time_is_not_backfilled(stamped, Timestamp("2026-05-01T00:00:00Z")) + + historical = next(iter(curve.observations(source_id=Identifier("issdc-pradan"), + ingest_time=None))) + assert historical.ingest_time is None + validator_for("observation").validate(historical.to_dict()) + + +def test_the_valid_time_is_the_second_the_sun_was_observed(acquired, tmp_path): + """The epoch resolved from the file (§2.1), not assumed.""" + _, adapter, acquisition = acquired + product = adapter.archive_for(FLARE_DAY) / "AL1_SLX_L1_20240514_v1.0" / "SDD2" / \ + "AL1_SOLEXS_20240514_SDD2_L1.lc.gz" + curve = lc.parse(product, Digest(digest_file(product).hex)) + observations = product_observations(curve, acquisition) + + assert str(observations[0].valid_time) == "2024-05-14T00:00:00Z" + assert str(observations[60_660].valid_time) == "2024-05-14T16:51:00Z" # the X8.7 peak + assert str(observations[-1].valid_time) == "2024-05-14T23:59:59Z" + + +def test_the_manifest_still_refuses_to_redistribute_the_parsed_archive(acquired): + """Parsing does not change the tier: the bytes stay referenced (ADR-0023, STD-23).""" + _, _, acquisition = acquired + manifest = { + "kind": "dataset", "digest": acquisition.artifact.digest.hex, "tier": 0, + "recorded_at": str(acquisition.provenance.ingest_time), + "retention": {"class": "permanent"}, + "retrieval": acquisition.artifact.retrieval.to_dict(), + } + 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, and the one grant #17 added. + + `astropy` is now permitted to Ingest — reading FITS needs a FITS reader. Nothing else + changed, and no context root was granted. + """ + 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") + assert ingest.allow == frozenset({"contracts", "domain", "kernel", "astropy"}) + assert "contexts" not in ingest.allow + + +# ═══════════════════════════════════════════ the real corpus + + +@real_only +def test_the_real_flare_day_lightcurve_reproduces_the_specification(): + """§2.1 `OBSERVED`, re-measured: 86,400 rows, 5 NaN at exactly the recorded offsets. + + The specification states these as observations of the archive. This asserts the parser + reproduces them, so a drift in either the code or the archive is a red test rather than a + quiet disagreement between a document and reality. + """ + path = real_product("lc") + curve = lc.parse(path, Digest(digest_file(path).hex)) + + assert curve.header.rows == 86_400 + assert curve.header.tstart == 1_715_644_800.0 + assert str(lc.unix_to_timestamp(curve.header.tstart)) == "2024-05-14T00:00:00Z" + assert curve.absent_offsets == SPEC_ABSENT_OFFSETS + assert curve.n_absent == 5 + assert curve.n_finite == 86_395 + assert curve.instrument_id == Identifier("solexs-sdd2") + + +@real_only +def test_the_real_flare_day_gti_satisfies_the_exact_exposure_equality(): + """§2.3 / F-09: Σ(STOP−START+1) == EXPOSURE, tolerance 0 s. + + The inclusive convention's whole justification, measured on the day it was verified for. + """ + path = real_product("gti") + intervals = gti.parse(path, Digest(digest_file(path).hex)) + + assert intervals.detector_active + assert len(intervals.intervals) == 5 + assert intervals.live_time == SPEC_EXPOSURE == 86_395.0 + assert intervals.declared_exposure == SPEC_EXPOSURE + assert intervals.live_time == intervals.declared_exposure + + +@real_only +def test_the_nan_set_equals_the_gti_excluded_set_on_the_reference_day(): + """§2.1 `OBSERVED`: on 2024-05-14 SDD2 the NaN set is exactly the GTI-excluded set. + + The r5 amendment withdrew *equality* as an archive-wide claim and replaced it with the + implication `NaN ⇒ GTI-excluded`. This asserts only what §2.1 states for **this day**, + and is deliberately not generalised: §2.3 §8 A-8 requires Milestone VIII to check all 436 + archives, and *"any deviation is a scientific finding and MUST terminate validation"*. + + Evaluating the implication across a day is the day-assembly layer's job (#19); this is + the single-day observation the specification records. + """ + lc_path, gti_path = real_product("lc"), real_product("gti") + curve = lc.parse(lc_path, Digest(digest_file(lc_path).hex)) + intervals = gti.parse(gti_path, Digest(digest_file(gti_path).hex)) + + excluded = intervals.excluded_seconds(curve.header.tstart, curve.header.rows) + assert curve.absent_offsets == excluded == SPEC_ABSENT_OFFSETS + assert curve.n_finite == intervals.declared_exposure + + +@real_only +def test_the_real_sdd1_gti_is_empty_and_that_is_legal(): + """F-12, the single deliberate non-terminating rule. §1.1 records SDD1 at zero rows.""" + path = REAL_ARCHIVE / "AL1_SLX_L1_20240514_v1.0" / "AL1_SLX_L1_20240514_v1.0" / \ + "SDD1" / "AL1_SOLEXS_20240514_SDD1_L1.gti.gz" + intervals = gti.parse(path, Digest(digest_file(path).hex)) + + assert intervals.detector_active is False + assert intervals.intervals == () + assert intervals.detector == "SDD1" + + +@real_only +def test_the_real_spectra_declare_340_ordinal_channels(): + """§2.2: DETCHANS=340, CHANTYPE='PI', constant channel map, and no energy anywhere.""" + path = real_product("pi") + spectra = pi.parse(path, Digest(digest_file(path).hex)) + + assert spectra.header.detchans == 340 + assert spectra.header.chantype == "PI" + assert spectra.channel_map == tuple(range(340)) + assert spectra.header.poisserr is False + + first = next(iter(spectra.spectra())) + assert len(first.counts) == 340 + assert str(first.valid_time) == "2024-05-14T00:00:00Z" + + +@real_only +def test_the_real_spectra_and_lightcurve_agree_on_the_epoch(): + """V-PI-3 / F-06, on the real products.""" + lc_path, pi_path = real_product("lc"), real_product("pi") + curve = lc.parse(lc_path, Digest(digest_file(lc_path).hex)) + spectra = pi.parse(pi_path, Digest(digest_file(pi_path).hex)) + pi.check_epoch_agrees_with_lightcurve(spectra, curve) + assert spectra.first_tstart == curve.header.tstart + + +@real_only +def test_real_observations_validate_against_the_contract(): + """Real Aditya-L1 X-ray measurements, as contract-valid Observations. + + A sample rather than all 86,400: validating every row against a JSON Schema would take + minutes and prove the same thing. The sample spans the day and includes every absent + second the archive declares, which are the rows most likely to be mishandled. + """ + path = real_product("lc") + curve = lc.parse(path, Digest(digest_file(path).hex)) + observations = list(curve.observations( + source_id=Identifier("issdc-pradan"), + ingest_time=Timestamp("2024-06-16T09:00:00Z"), + )) + + assert len(observations) == 86_400 + validator = validator_for("observation") + sample = [observations[i] for i in (0, 1, 5, 30_072, 30_078, 60_660, 83_951, 86_399)] + for observation in sample: + validator.validate(observation.to_dict()) + assert observation_is_wellformed(observation) + + # The absent seconds are absent, and nothing else is. + for offset in SPEC_ABSENT_OFFSETS: + assert observations[offset].value is None + assert observations[60_660].value is not None, "the flare peak minute is measured" + + +@real_only +def test_the_flare_peak_carries_a_real_measurement(): + """2024-05-14T16:51 UTC — the GOES X8.7 peak, and SALVAGE-001's reference event. + + The assertion is deliberately weak on physics: that the second is measured, finite and + non-negative. §2.2 forbids stating a SoLEXS energy, and no response file exists, so any + claim beyond "a count was recorded here" would exceed what the archive supports. + """ + path = real_product("lc") + curve = lc.parse(path, Digest(digest_file(path).hex)) + + peak_offset = 16 * 3600 + 51 * 60 + value = curve.counts[peak_offset] + assert value is not None + assert math.isfinite(value) and value >= 0 + assert str(lc.unix_to_timestamp(curve.times[peak_offset])) == "2024-05-14T16:51:00Z" diff --git a/tools/gates/imports.py b/tools/gates/imports.py index 39595a0..492c8dd 100644 --- a/tools/gates/imports.py +++ b/tools/gates/imports.py @@ -323,8 +323,15 @@ def run(policies: list[Policy]) -> tuple[Report, int]: # staleness check did its job on the way: with the modules added and the flag still False # the gate went red — "declared unpopulated but contains 7 module(s)" — which is exactly # what M2/E4/#13 built it for. The rule now binds against real code rather than a fixture. + # + # `astropy` granted by M3/E5/#17. The SoLEXS products are FITS, and reading FITS needs a + # FITS reader — hand-rolling one to avoid a dependency would be a far larger risk to + # scientific fidelity than the dependency is to the architecture. This is the path #13 + # described when it declined to pre-grant third-party roots: "a context that later needs + # `astropy` should add one reviewable line rather than inherit a blanket permission + # nobody voted for (STD-11)." This is that line, and it grants astropy to Ingest only. Policy(package="contexts.ingest", populated=True, - allow=frozenset({"contracts", "domain", "kernel"})), + allow=frozenset({"contracts", "domain", "kernel", "astropy"})), # R2 — Curation. Freezes observations into digest-addressed releases; the kernel is what # mints those digests (ADR-0005, ADR-0006).