#!/usr/bin/env python3 """Controlled FITS integrity observations. Only generated temporary fixtures are read.""" import argparse import hashlib import importlib.metadata import io import json import platform import struct import tempfile import warnings from datetime import datetime, timezone from pathlib import Path import astropy import numpy as np from astropy.io import fits VALUES = [10, 20, 30, 40] COMMENT = "Hesela deterministic fixture" STATUS = {0: "mismatch", 1: "valid", 2: "absent"} def sha(data): return hashlib.sha256(data).hexdigest() def sum32(data): """Independent accumulator for these 4-byte-aligned, uncompressed fixtures.""" if len(data) % 4: raise ValueError("Expected complete 32-bit words") total = sum(word[0] for word in struct.iter_unpack(">I", data)) while total >> 32: total = (total & 0xffffffff) + (total >> 32) return total def seal(raw, extension=False): cls = fits.ImageHDU if extension else fits.PrimaryHDU hdu = cls.fromstring(raw) hdu.add_datasum(when=COMMENT) hdu.add_checksum(when=COMMENT, override_datasum=True) # Fixture-only serialization: one 2880-byte header and one data record. header = hdu.header.tostring().encode("ascii") assert len(header) == 2880 and len(raw) == 5760 result = header + raw[2880:] assert sum32(result) == 0xffffffff return result def fixture(extension=False): hdu = (fits.ImageHDU if extension else fits.PrimaryHDU)( np.array(VALUES, dtype=np.int32)) hdu.header["BSCALE"] = 1 hdu.header["BZERO"] = 0 if extension: hdu.header["EXTNAME"] = "SCI" stream = io.BytesIO() if extension: # ImageHDU.writeto inserts a primary HDU; select the generated extension. fits.HDUList([fits.PrimaryHDU(), hdu]).writeto(stream) raw = stream.getvalue()[2880:] else: hdu.writeto(stream) raw = stream.getvalue() return seal(raw, extension) def replace_card(raw, keyword, replacement): header = bytearray(raw[:2880]) offsets = [i for i in range(0, 2880, 80) if header[i:i + 8].decode("ascii").strip() == keyword] assert len(offsets) == 1 encoded = (fits.Card(keyword, replacement).image.encode("ascii") if replacement is not None else b" " * 80) header[offsets[0]:offsets[0] + 80] = encoded return bytes(header) + raw[2880:] def fixtures(): baseline = fixture() flipped = bytearray(baseline) flipped[2883] ^= 1 padding = bytearray(baseline) # Byte 19 in the data record is padding, beyond the four 4-byte pixels. padding[2880 + 19] ^= 1 swapped = baseline[:2880] + baseline[2884:2888] + baseline[2880:2884] + baseline[2888:] scaled = replace_card(baseline, "BSCALE", 2) absent = replace_card(replace_card(baseline, "CHECKSUM", None), "DATASUM", None) extension = fixture(extension=True) corrupt_extension = bytearray(extension) corrupt_extension[2883] ^= 1 return [ ("baseline", "Unchanged image", baseline, baseline, 1), ("bit_flip", "One data bit flipped", bytes(flipped), baseline, 1), ("padding_bit_flip", "One padding bit flipped", bytes(padding), baseline, 1), ("word_swap", "Two 32-bit pixels exchanged", swapped, baseline, 1), ("scale_changed", "BSCALE changed to 2", scaled, baseline, 1), ("scale_resealed", "Changed scale, sums regenerated", seal(scaled), baseline, 1), ("keywords_absent", "Both checksum cards removed", absent, baseline, 1), ("late_corruption", "Only second HDU corrupted", baseline + bytes(corrupt_extension), baseline + extension, 2), ("extension_removed", "Whole second HDU removed", baseline, baseline + extension, 2), ] def observe(case_id, label, raw, reference, expected_hdus, root): path = root / (case_id + ".fits") path.write_bytes(raw) with warnings.catch_warnings(record=True) as default_warnings: warnings.simplefilter("always") with fits.open(path, memmap=False, lazy_load_hdus=False) as hdus: decoded = [hdu.data.tolist() if hdu.data is not None else None for hdu in hdus] with warnings.catch_warnings(record=True) as checked_warnings: warnings.simplefilter("always") with fits.open(path, checksum=True, memmap=False, lazy_load_hdus=False, do_not_scale_image_data=True) as hdus: checks = [] for index, hdu in enumerate(hdus): info = hdu.fileinfo() start, data = info["hdrLoc"], info["datLoc"] end = data + info["datSpan"] checks.append({ "index": index, "datasum_code": hdu.verify_datasum(), "checksum_code": hdu.verify_checksum(), "datasum_card": hdu.header.get("DATASUM"), "checksum_card": hdu.header.get("CHECKSUM"), "data_sha256": sha(raw[data:end]), "raw_data_sum32": sum32(raw[data:end]), "raw_hdu_sum32": sum32(raw[start:end]), }) hdus.verify("exception") # A separate read-only open isolates data loading from scaling and checksum=True. with fits.open(path, mode="readonly", memmap=False, lazy_load_hdus=False, do_not_scale_image_data=True, checksum=False) as hdus: read_states = [] for index, hdu in enumerate(hdus): before = [hdu.verify_datasum(), hdu.verify_checksum()] stored_pixels = hdu.data.tolist() if hdu.data is not None else None after = [hdu.verify_datasum(), hdu.verify_checksum()] read_states.append({ "index": index, "before_data_access": before, "after_data_access": after, "stored_pixels": stored_pixels, }) after_read = path.read_bytes() assert after_read == raw, "Read-only observation changed its input" return { "id": case_id, "label": label, "bytes": len(raw), "sha256": sha(raw), "reference_sha256": sha(reference), "byte_identical_to_reference": raw == reference, "expected_hdus": expected_hdus, "observed_hdus": len(checks), "structure_verify": "returned_without_exception", "default_open": "returned", "default_warning_count": len(default_warnings), "checksum_open": "returned", "checksum_warnings": [str(w.message).strip() for w in checked_warnings], "decoded_pixels": decoded, "hdu_checks": checks, "read_state_checks": read_states, "after_read_sha256": sha(after_read), "all_present_hdus_valid": all( c["datasum_code"] == 1 and c["checksum_code"] == 1 for c in checks), } def run(): with tempfile.TemporaryDirectory(prefix="hesela-fits-") as name: rows = [observe(*item, Path(name)) for item in fixtures()] return { "schemaVersion": "1.1", "kind": "controlled-software-experiment", "script_sha256": sha(Path(__file__).read_bytes()), "generatedAt": datetime.now(timezone.utc).isoformat(), "environment": { "python": platform.python_version(), "platform": platform.platform(), "astropy": astropy.__version__, "numpy": np.__version__, "dependencies": {name: importlib.metadata.version(name) for name in ["astropy", "numpy", "pyerfa", "astropy-iers-data", "packaging", "PyYAML"]}, }, "fixture": {"pixels": VALUES, "bitpix": 32, "bscale": 1, "bzero": 0, "header_bytes": 2880, "data_record_bytes": 2880, "compressed": False, "astronomical_observation": False}, "statusCodes": {str(k): v for k, v in STATUS.items()}, "cases": rows, "limits": [ "Constructed images, not sampled corruption or astronomical observations.", "Astropy version-scoped behavior; CFITSIO is not executed.", "No compression, tables, remote storage, concurrent writes or power loss.", "Hashes compared with a separate generated reference, not proof of origin.", "The independent accumulator only accepts complete 32-bit words.", ], } def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--output", type=Path, required=True) args = parser.parse_args() # Exclusive creation prevents overwriting an earlier result. with args.output.open("x", encoding="utf-8") as target: json.dump(run(), target, indent=2) target.write("\n") if __name__ == "__main__": main()