"""Run the threshold aperture on the committed real sector 56 products. This script reads two local files, ``research/tess-timing/data/real-run/tpf.fits`` and ``research/tess-timing/data/real-run/lightcurve.fits``, with the existing :mod:`tess_timing.aperture` functions :func:`~tess_timing.aperture.read_target_pixel_file`, :func:`~tess_timing.aperture.reference_image`, :func:`~tess_timing.aperture.select_aperture`, and :func:`~tess_timing.aperture.extract_light_curve`, and the existing :func:`tess_timing.aperture_compare.compare_to_mission`. The default aperture configuration and the default 175 quality mask are used unchanged. The committed outputs are: research/tess-timing/data/real-aperture-run/real_aperture.json research/tess-timing/data/real-aperture-run/real_aperture.csv docs/TESS-REAL-APERTURE.md Run from the repository root to regenerate them: python research/tess-timing/scripts/real_aperture.py Recompute into a temporary directory under ``.scratch`` and compare the JSON, CSV, and Markdown against the committed files: python research/tess-timing/scripts/real_aperture.py --check No network request is made. No detrending, transit fit, folding, or timing measurement is performed. The comparison is a diagnostic on one common quality mask, not a validation of the SPOC pipeline, and it does not claim that a transit is present. """ from __future__ import annotations import argparse import csv import json import math import shutil import sys import tempfile from pathlib import Path from typing import Any PACKAGE_ROOT = Path(__file__).resolve().parents[1] REPO_ROOT = PACKAGE_ROOT.parents[1] SRC = PACKAGE_ROOT / "src" if str(SRC) not in sys.path: sys.path.insert(0, str(SRC)) # isort: off from astropy.io import fits from tess_timing.aperture import ( DEFAULT_APERTURE_CONFIG, ApertureError, extract_light_curve, read_target_pixel_file, reference_image, select_aperture, ) from tess_timing.aperture_compare import compare_to_mission from tess_timing.lightcurve import ( DEFAULT_QUALITY_MASK, LightCurveError, load_lightcurve, portable_path, sha256_file, software_versions, ) # isort: on TASK = "jdp-real-aperture" SCHEMA_VERSION = 1 LABEL = "real SPOC sector 56 product comparison" TPF_PATH = PACKAGE_ROOT / "data" / "real-run" / "tpf.fits" LIGHTCURVE_PATH = PACKAGE_ROOT / "data" / "real-run" / "lightcurve.fits" DEFAULT_JSON_PATH = PACKAGE_ROOT / "data" / "real-aperture-run" / "real_aperture.json" DEFAULT_CSV_PATH = PACKAGE_ROOT / "data" / "real-aperture-run" / "real_aperture.csv" DEFAULT_MARKDOWN_PATH = REPO_ROOT / "docs" / "TESS-REAL-APERTURE.md" CHECK_ROOT = REPO_ROOT / ".scratch" / "real-aperture-check" EXPECTED_TPF_SHA256 = ( "2b5bd6ee4747cfe0c5405c869b871064fd74e4b5f0521c4e6be9fdc2b5c96030" ) EXPECTED_LIGHTCURVE_SHA256 = ( "1b76b4a4b73e24954fa6e29a7e28a174685113ecdbad1bf10ef0e7b766ce99e9" ) EXPECTED_TARGET_NAME = "HD 209458 b" EXPECTED_TIC_ID = 420814525 EXPECTED_SECTOR = 56 EXPECTED_MASK = DEFAULT_QUALITY_MASK.value EXPECTED_CAMERA = 1 EXPECTED_CCD = 2 CSV_COLUMNS = ( "time_btjd", "aperture_flux_normalized", "aperture_flux_err_normalized", "sap_flux_normalized", "pdcsap_flux_normalized", ) LIMITS = [ ( "This is a real SPOC product comparison with no detrending, no transit " "fit, and no timing measurement." ), ( "The comparison is a diagnostic on one common quality mask, not a " "validation of the SPOC pipeline." ), ( "No transit, depth, duration, period, ephemeris, or mid-transit time " "is measured or claimed." ), ( "The aperture series and each mission series are divided by their " "median on the common mask before the difference metrics and the CSV." ), ( "The CSV has one row per common-mask cadence, in time order. The " "aperture uncertainty is the extracted uncertainty divided by the same " "aperture median." ), ] class RealApertureError(RuntimeError): """The input products or the recorded outputs are not acceptable.""" def _display(path: Path) -> str: return portable_path(path) def _header_int(header: fits.Header, key: str) -> int: value = header.get(key) if value is None: raise RealApertureError(f"the primary FITS header has no {key} value") try: return int(value) except (TypeError, ValueError) as error: raise RealApertureError( f"the primary FITS header {key} value {value!r} is not an integer" ) from error def _preflight(path: Path, expected_sha256: str, label: str) -> dict[str, Any]: """Verify the digest and primary header before photometry.""" path = Path(path) if not path.exists(): raise RealApertureError(f"{label} input file does not exist: {path}") digest = sha256_file(path) if digest != expected_sha256: raise RealApertureError( f"{label} SHA-256 is {digest}; expected {expected_sha256}. The file " "is not the committed sector 56 product." ) try: with fits.open(path, memmap=False) as hdul: header = hdul[0].header tic_id = _header_int(header, "TICID") sector = _header_int(header, "SECTOR") camera = _header_int(header, "CAMERA") ccd = _header_int(header, "CCD") object_name = header.get("OBJECT") synthetic = bool( header.get("SYNTH", False) or str(header.get("ORIGIN", "")).lower() == "synthetic" ) except OSError as error: raise RealApertureError( f"could not open {label} FITS file {path}: {error}" ) from error if tic_id != EXPECTED_TIC_ID: raise RealApertureError( f"{label} primary header TICID is {tic_id}; expected {EXPECTED_TIC_ID}" ) if sector != EXPECTED_SECTOR: raise RealApertureError( f"{label} primary header SECTOR is {sector}; expected {EXPECTED_SECTOR}" ) if synthetic: raise RealApertureError(f"the {label} FITS header marks the input as synthetic") return { "tic_id": tic_id, "sector": sector, "camera": camera, "ccd": ccd, "object": str(object_name) if object_name is not None else None, "synthetic": synthetic, "sha256": digest, } def _comparison_block(raw: dict[str, Any]) -> dict[str, Any]: return { "n_points": raw["n_points"], "pearson_r": raw["pearson_r"], "r_squared": raw["r_squared"], "mean_difference": raw["mean_difference"], "median_difference": raw["median_difference"], "rms_difference": raw["rms_difference"], "rms_difference_ppm": raw["rms_difference_ppm"], "max_absolute_difference": raw["max_abs_difference"], "correlation_defined": raw["correlation_defined"], } def build_record( tpf_path: Path = TPF_PATH, lightcurve_path: Path = LIGHTCURVE_PATH, ) -> tuple[dict[str, Any], Any, Any]: """Build the JSON record and return it with the series for the CSV.""" tpf_path = Path(tpf_path) lightcurve_path = Path(lightcurve_path) tpf_header = _preflight(tpf_path, EXPECTED_TPF_SHA256, "target pixel file") _preflight(lightcurve_path, EXPECTED_LIGHTCURVE_SHA256, "light curve") config = DEFAULT_APERTURE_CONFIG try: tpf = read_target_pixel_file(tpf_path) selection = select_aperture(reference_image(tpf), config) extracted = extract_light_curve(tpf, selection, mask=EXPECTED_MASK) light_curve = load_lightcurve(lightcurve_path, mask=EXPECTED_MASK) comparison = compare_to_mission(extracted, light_curve, config) except (ApertureError, LightCurveError, OSError) as error: raise RealApertureError(str(error)) from error if tpf.sha256 != EXPECTED_TPF_SHA256: raise RealApertureError("the loader read a different target pixel file digest") if light_curve.input_sha256 != EXPECTED_LIGHTCURVE_SHA256: raise RealApertureError("the loader read a different light-curve digest") if tpf.synthetic or light_curve.synthetic: raise RealApertureError("one of the loaded products is marked synthetic") if tpf.metadata.get("TICID") not in (None, EXPECTED_TIC_ID): raise RealApertureError("the target pixel file TICID metadata disagrees") if extracted.mask.value != EXPECTED_MASK or light_curve.mask.value != EXPECTED_MASK: raise RealApertureError("the quality mask is not the expected 175") raw_comparison = comparison.metrics for key in ("sap", "pdcsap"): block = raw_comparison[key] if block["pearson_r"] is None: raise RealApertureError( f"{key} correlation is undefined; the comparison cannot be recorded" ) record: dict[str, Any] = { "schema_version": SCHEMA_VERSION, "task": TASK, "status": "OK", "label": LABEL, "synthetic": False, "target": { "name": EXPECTED_TARGET_NAME, "tic_id": EXPECTED_TIC_ID, "sector": EXPECTED_SECTOR, "camera": tpf_header["camera"], "ccd": tpf_header["ccd"], }, "inputs": { "tpf": { "path": _display(tpf.path), "filename": tpf.path.name, "sha256": tpf.sha256, "size_bytes": tpf.size_bytes, "synthetic": tpf.synthetic, "shape": { "cadences": tpf.cadence_count, "rows": tpf.shape[0], "cols": tpf.shape[1], }, "metadata": dict(tpf.metadata), }, "lightcurve": { "path": _display(light_curve.input_path), "filename": light_curve.input_path.name, "sha256": light_curve.input_sha256, "size_bytes": light_curve.input_size_bytes, "synthetic": light_curve.synthetic, }, }, "quality_mask": { "mask_value": extracted.mask.value, "mask_bits": extracted.mask.value, "bits": list(extracted.mask.bits), "bit_table": extracted.mask.bit_records(), "policy": "A cadence is excluded when (QUALITY & mask_value) != 0.", }, "config": { "threshold_fraction": config.threshold_fraction, "inner_radius_pixels": config.inner_radius, "outer_radius_pixels": config.outer_radius, "time_tolerance_seconds": config.time_tolerance_seconds, }, "aperture": selection.summary(), "cadence_counts": dict(extracted.counts), "retained_cadence_count": int(extracted.counts["retained"]), "comparison": { "common_mask": raw_comparison["common_mask"], "common_cadence_count": comparison.count, "matched_pairs": raw_comparison["matched_pairs"], "unmatched_aperture_cadences": raw_comparison[ "unmatched_aperture_cadences" ], "unmatched_mission_cadences": raw_comparison["unmatched_mission_cadences"], "max_time_offset_seconds": raw_comparison["max_time_offset_seconds"], "time_tolerance_seconds": raw_comparison["time_tolerance_seconds"], "normalization": raw_comparison["normalization"], "medians": { "aperture_raw_flux": raw_comparison["aperture_median_raw_flux"], "sap_raw_flux": raw_comparison["sap_median_raw_flux"], "pdcsap_raw_flux": raw_comparison["pdcsap_median_raw_flux"], }, "sap": _comparison_block(raw_comparison["sap"]), "pdcsap": _comparison_block(raw_comparison["pdcsap"]), }, "software": software_versions(), "limits": LIMITS, } return record, extracted, comparison def write_json(record: dict[str, Any], path: Path) -> Path: path = Path(path) path.parent.mkdir(parents=True, exist_ok=True) path.write_text( json.dumps(record, indent=2, ensure_ascii=False) + "\n", encoding="utf-8", ) return path def write_csv(extracted: Any, comparison: Any, path: Path) -> Path: """Write one row per common-mask cadence with the documented schema.""" path = Path(path) path.parent.mkdir(parents=True, exist_ok=True) aperture_median = comparison.metrics.get("aperture_median_raw_flux") if aperture_median is None or not math.isfinite(float(aperture_median)): raise RealApertureError("the aperture median is not finite") if float(aperture_median) == 0.0: raise RealApertureError("the aperture median is zero; normalization is undefined") indices = comparison.extracted_index with path.open("w", encoding="utf-8", newline="") as handle: writer = csv.writer(handle, lineterminator="\n") writer.writerow(CSV_COLUMNS) for position in range(comparison.count): index = int(indices[position]) writer.writerow( [ f"{comparison.time[position]:.12f}", f"{comparison.aperture_normalized[position]:.12g}", f"{extracted.flux_err[index] / float(aperture_median):.12g}", f"{comparison.sap_normalized[position]:.12g}", f"{comparison.pdcsap_normalized[position]:.12g}", ] ) return path def _metric_text(value: object) -> str: if value is None: return "undefined" return f"{float(value):.6f}" def _science_text(value: object) -> str: if value is None: return "undefined" return f"{float(value):.6g}" def markdown(record: dict[str, Any]) -> str: target = record["target"] aperture = record["aperture"] counts = record["cadence_counts"] comparison = record["comparison"] sap = comparison["sap"] pdcsap = comparison["pdcsap"] bits_text = ", ".join(str(bit) for bit in record["quality_mask"]["bits"]) csv_schema = ", ".join(f"`{name}`" for name in CSV_COLUMNS) lines = [ "# Real TESS aperture product comparison", "", ( "Status: real SPOC sector 56 product comparison. A threshold " "aperture was run on the committed target pixel file and compared " "with the mission SAP and PDCSAP columns. No detrending, no " "transit fit, and no timing measurement was performed." ), "", f"- Target: {target['name']}", f"- TIC identifier: {target['tic_id']}", f"- Sector: {target['sector']}", f"- Camera / CCD: {target['camera']} / {target['ccd']}", f"- Target pixel file: `{record['inputs']['tpf']['path']}`", f"- Target pixel file SHA-256: `{record['inputs']['tpf']['sha256']}`", f"- Light curve: `{record['inputs']['lightcurve']['path']}`", f"- Light curve SHA-256: `{record['inputs']['lightcurve']['sha256']}`", f"- Quality mask: {record['quality_mask']['mask_value']} (bits {bits_text})", "- Synthetic: no", ( "- Machine-readable record: " "`research/tess-timing/data/real-aperture-run/real_aperture.json`" ), ( "- Per-cadence table: " "`research/tess-timing/data/real-aperture-run/real_aperture.csv`" ), "", "## Aperture", "", "| Quantity | Value |", "| --- | ---: |", f"| Aperture pixels | {aperture['aperture_pixel_count']} |", f"| Annulus pixels | {aperture['annulus_pixel_count']} |", ( "| Threshold fraction of peak excess | " f"{aperture['threshold_fraction_of_peak_excess']} |" ), f"| Threshold value (electrons/s) | {aperture['threshold_flux']:.6f} |", f"| Background (electrons/s) | {aperture['background_flux']:.6f} |", ( "| Centroid (row, column) | " f"{aperture['centroid_row']:.6f}, {aperture['centroid_col']:.6f} |" ), ( "| Aperture touches image edge | " f"{'yes' if aperture['aperture_touches_image_edge'] else 'no'} |" ), "", ( "The aperture configuration is the module default: threshold " f"fraction {record['config']['threshold_fraction']}, annulus radii " f"{record['config']['inner_radius_pixels']} to " f"{record['config']['outer_radius_pixels']} pixels, and pairing " f"tolerance {record['config']['time_tolerance_seconds']} s." ), "", "## Cadence counts", "", "| Quantity | Cadences |", "| --- | ---: |", f"| Input rows | {counts['input_cadences']:,} |", f"| Passed the quality mask | {counts['quality_pass']:,} |", f"| Retained by the aperture extraction | {counts['retained']:,} |", f"| Rejected by the quality mask | {counts['rejected_by_quality_mask']:,} |", f"| Rejected for a non-finite aperture pixel | {counts['rejected_nonfinite_pixel']:,} |", f"| Rejected for a zero-weight aperture pixel | {counts['rejected_zero_weight_pixel']:,} |", f"| Rejected for the background annulus | {counts['rejected_background']:,} |", f"| Common-mask cadences | {comparison['common_cadence_count']:,} |", "", ( "The common mask requires the same quality mask on both inputs and " "a valid aperture extraction. The comparison paired " f"{comparison['matched_pairs']:,} cadences with a maximum time " f"offset of {comparison['max_time_offset_seconds']} s." ), "", "## Comparison", "", ( "Both the aperture series and each mission series are divided by " "their median on the common mask. The difference is aperture minus " "mission. The RMS difference is reported in ppm of the normalized " "series." ), "", "| Series | Points | Pearson r | RMS difference (ppm) | Max absolute difference |", "| --- | ---: | ---: | ---: | ---: |", ( f"| SAP | {sap['n_points']:,} | " f"{_metric_text(sap['pearson_r'])} | " f"{_metric_text(sap['rms_difference_ppm'])} | " f"{_science_text(sap['max_absolute_difference'])} |" ), ( f"| PDCSAP | {pdcsap['n_points']:,} | " f"{_metric_text(pdcsap['pearson_r'])} | " f"{_metric_text(pdcsap['rms_difference_ppm'])} | " f"{_science_text(pdcsap['max_absolute_difference'])} |" ), "", "## CSV schema", "", f"Columns, in order: {csv_schema}.", "", ( "The CSV has one header row and one row per common-mask cadence, in " "time order. Flux and uncertainty columns are normalized by the " "aperture median; SAP and PDCSAP are normalized by their own " "medians on the common mask." ), "", "## Reproduction", "", "```", "python research/tess-timing/scripts/real_aperture.py", "python research/tess-timing/scripts/real_aperture.py --check", "```", "", "## Limits", "", ] for limit in record["limits"]: lines.append(f"- {limit}") lines.append("") return "\n".join(lines) def write_markdown(record: dict[str, Any], path: Path) -> Path: path = Path(path) path.parent.mkdir(parents=True, exist_ok=True) path.write_text(markdown(record), encoding="utf-8") return path def _compare(recorded: object, recomputed: object, prefix: str = "") -> list[str]: issues: list[str] = [] if isinstance(recorded, dict) and isinstance(recomputed, dict): for key in sorted(set(recorded) | set(recomputed)): location = f"{prefix}.{key}" if prefix else key if key not in recorded: issues.append(f"{location}: missing from recorded value") elif key not in recomputed: issues.append(f"{location}: missing from recomputed value") else: issues.extend(_compare(recorded[key], recomputed[key], location)) return issues if isinstance(recorded, list) and isinstance(recomputed, list): if len(recorded) != len(recomputed): issues.append( f"{prefix}: recorded {len(recorded)} item(s); " f"recomputed {len(recomputed)}" ) return issues for index, (left, right) in enumerate(zip(recorded, recomputed, strict=True)): issues.extend(_compare(left, right, f"{prefix}[{index}]")) return issues if isinstance(recorded, bool) or isinstance(recomputed, bool): if recorded != recomputed: issues.append(f"{prefix}: recorded {recorded!r}; recomputed {recomputed!r}") return issues if isinstance(recorded, (int, float)) and isinstance(recomputed, (int, float)): if not math.isclose( float(recorded), float(recomputed), rel_tol=1e-12, abs_tol=1e-12 ): issues.append(f"{prefix}: recorded {recorded!r}; recomputed {recomputed!r}") return issues if recorded != recomputed: issues.append(f"{prefix}: recorded {recorded!r}; recomputed {recomputed!r}") return issues def check_record( json_path: Path = DEFAULT_JSON_PATH, csv_path: Path = DEFAULT_CSV_PATH, markdown_path: Path = DEFAULT_MARKDOWN_PATH, ) -> int: """Recompute into ``.scratch`` and compare with the committed outputs.""" json_path = Path(json_path) csv_path = Path(csv_path) markdown_path = Path(markdown_path) if not json_path.exists(): raise RealApertureError(f"recorded JSON does not exist: {json_path}") if not csv_path.exists(): raise RealApertureError(f"recorded CSV does not exist: {csv_path}") try: recorded = json.loads(json_path.read_text(encoding="utf-8")) except (OSError, json.JSONDecodeError) as error: raise RealApertureError(f"could not read recorded JSON {json_path}: {error}") from error CHECK_ROOT.mkdir(parents=True, exist_ok=True) temporary = Path(tempfile.mkdtemp(prefix="check-", dir=CHECK_ROOT)) try: recomputed, extracted, comparison = build_record(TPF_PATH, LIGHTCURVE_PATH) temporary_json = write_json(recomputed, temporary / json_path.name) temporary_csv = write_csv(extracted, comparison, temporary / csv_path.name) temporary_markdown = write_markdown( recomputed, temporary / (markdown_path.name), ) issues = _compare(recorded, recomputed) if csv_path.read_bytes() != temporary_csv.read_bytes(): issues.append("CSV: recorded bytes differ from the recomputed bytes") if markdown_path.exists() and ( markdown_path.read_bytes() != temporary_markdown.read_bytes() ): issues.append("Markdown: recorded bytes differ from the recomputed bytes") if not temporary_json.exists(): issues.append("JSON: recomputed file was not written") finally: shutil.rmtree(temporary, ignore_errors=True) if issues: raise RealApertureError( "recorded outputs do not match the recomputed values: " + "; ".join(issues[:10]) + (f"; {len(issues) - 10} more difference(s)" if len(issues) > 10 else "") ) return 0 def _parser() -> argparse.ArgumentParser: parser = argparse.ArgumentParser( description=( "Run the threshold aperture on the committed real sector 56 target " "pixel file and compare it with the mission SAP and PDCSAP " "columns. No network request is made. No detrending, transit fit, " "or timing measurement is performed." ) ) parser.add_argument( "--check", action="store_true", help=( "recompute into a temporary directory under .scratch and compare " "the JSON, CSV, and Markdown with the committed outputs" ), ) parser.add_argument( "--json-out", type=Path, default=None, help=f"output JSON path (default: {_display(DEFAULT_JSON_PATH)})", ) parser.add_argument( "--csv-out", type=Path, default=None, help=f"output CSV path (default: {_display(DEFAULT_CSV_PATH)})", ) parser.add_argument( "--markdown-out", type=Path, default=None, help=f"output Markdown path (default: {_display(DEFAULT_MARKDOWN_PATH)})", ) return parser def main(argv: list[str] | None = None) -> int: arguments = _parser().parse_args(argv) if arguments.check: if ( arguments.json_out is not None or arguments.csv_out is not None or arguments.markdown_out is not None ): print( "error: --check cannot be combined with an output-path option", file=sys.stderr, ) return 2 try: check_record() except (RealApertureError, OSError) as error: print(f"error: {error}", file=sys.stderr) return 1 print( f"status=OK check=record tpf_sha256={EXPECTED_TPF_SHA256} " f"lightcurve_sha256={EXPECTED_LIGHTCURVE_SHA256} " f"tic_id={EXPECTED_TIC_ID} sector={EXPECTED_SECTOR} " f"mask_bits={EXPECTED_MASK}" ) return 0 json_path = arguments.json_out or DEFAULT_JSON_PATH csv_path = arguments.csv_out or DEFAULT_CSV_PATH markdown_path = arguments.markdown_out or DEFAULT_MARKDOWN_PATH try: record, extracted, comparison = build_record(TPF_PATH, LIGHTCURVE_PATH) write_json(record, json_path) write_csv(extracted, comparison, csv_path) write_markdown(record, markdown_path) except (RealApertureError, OSError) as error: print(f"error: {error}", file=sys.stderr) return 1 sap = record["comparison"]["sap"] pdcsap = record["comparison"]["pdcsap"] print( f"status=OK label=real-spoc-sector-56 " f"target={EXPECTED_TARGET_NAME} tic_id={EXPECTED_TIC_ID} " f"sector={EXPECTED_SECTOR} mask_bits={EXPECTED_MASK} " f"aperture_pixels={record['aperture']['aperture_pixel_count']} " f"annulus_pixels={record['aperture']['annulus_pixel_count']} " f"threshold_flux={record['aperture']['threshold_flux']} " f"retained={record['retained_cadence_count']} " f"common={record['comparison']['common_cadence_count']} " f"r_sap={_metric_text(sap['pearson_r'])} " f"r_pdcsap={_metric_text(pdcsap['pearson_r'])} " f"wrote={_display(json_path)},{_display(csv_path)} " f"markdown={_display(markdown_path)}" ) return 0 if __name__ == "__main__": raise SystemExit(main())