from __future__ import annotations import json from pathlib import Path import struct from astropy.coordinates import SkyCoord import astropy.units as u import numpy as np import pandas as pd import pytest import generate_explorer_products as explorer HERE = Path(__file__).resolve().parent PRODUCTS = HERE / "products" FIGURES = HERE / "figures" @pytest.fixture(scope="session") def regenerated(tmp_path_factory: pytest.TempPathFactory) -> Path: output_root = tmp_path_factory.mktemp("regenerated-explorer") outputs = explorer.generate(output_root) assert len(outputs) == 21 return output_root def test_current_generator_reproduces_committed_products(regenerated: Path) -> None: for filename in explorer.DATA_FILENAMES: assert (regenerated / "products" / filename).read_bytes() == ( PRODUCTS / filename ).read_bytes() figure_paths = sorted((regenerated / "figures").glob("*")) assert len(figure_paths) == 16 pngs = [path for path in figure_paths if path.suffix == ".png"] assert len(pngs) == 8 for path in pngs: width, height = _png_dimensions(path) expected_width = 1340 if "_1col." in path.name else 2800 assert width == expected_width assert (width, height) == _png_dimensions(FIGURES / path.name) @pytest.fixture(scope="module") def latitude() -> pd.DataFrame: return pd.read_csv(PRODUCTS / "mw_latitude_profile.csv") @pytest.fixture(scope="module") def axis() -> pd.DataFrame: return pd.read_csv(PRODUCTS / "mw_m31_gc_signed_axis_profile.csv") @pytest.fixture(scope="module") def fields() -> pd.DataFrame: return pd.read_csv( PRODUCTS / "m31_field_profile_decomposition.csv", dtype={"obsid": str} ) def test_components_sum_without_cross_term( latitude: pd.DataFrame, axis: pd.DataFrame ) -> None: for table in (latitude, axis): assert np.allclose( table["em_total_kpc_cm-6"], table["em_halo_kpc_cm-6"] + table["em_disk_kpc_cm-6"], rtol=2.0e-12, atol=0.0, ) for quantity, suffix in ( ("o8", "lu"), ("direct", "absorbed_0p4_1p25_primary_fluxunit"), ("response_matched", "absorbed_0p4_1p25_primary_fluxunit"), ): assert np.allclose( table[f"{quantity}_total_{suffix}"], table[f"{quantity}_halo_{suffix}"] + table[f"{quantity}_disk_{suffix}"], rtol=2.0e-12, atol=0.0, ) def test_m31_center_coordinates_are_exact(axis: pd.DataFrame) -> None: center = axis.loc[np.isclose(axis["signed_axis_kpc"], 0.0)].iloc[0] assert center["ra_deg"] == pytest.approx(explorer.M31_RA_DEG, abs=1.0e-12) assert center["dec_deg"] == pytest.approx(explorer.M31_DEC_DEG, abs=1.0e-12) assert center["galactic_l_deg"] == pytest.approx(121.17432906153702, abs=1.0e-10) assert center["galactic_b_deg"] == pytest.approx(-21.573308798437054, abs=1.0e-10) coordinates = SkyCoord( axis["ra_deg"].to_numpy() * u.deg, axis["dec_deg"].to_numpy() * u.deg, frame="icrs", ).galactic assert np.allclose(coordinates.l.deg, axis["galactic_l_deg"], atol=2.0e-10) assert np.allclose(coordinates.b.deg, axis["galactic_b_deg"], atol=2.0e-10) def test_axis_endpoints_sign_and_tangent_radius(axis: pd.DataFrame) -> None: assert axis["signed_axis_kpc"].iloc[0] == -50.0 assert axis["signed_axis_kpc"].iloc[-1] == 50.0 expected_theta = np.rad2deg( np.arctan(axis["signed_axis_kpc"] / explorer.M31_DISTANCE_KPC) ) assert np.allclose(axis["angular_offset_deg"], expected_theta, atol=2.0e-14) assert np.allclose( explorer.M31_DISTANCE_KPC * np.tan(np.deg2rad(axis["angular_offset_deg"])), axis["signed_axis_kpc"], atol=2.0e-12, ) endpoint_coordinates = SkyCoord( axis.loc[[0, len(axis) - 1], "ra_deg"].to_numpy() * u.deg, axis.loc[[0, len(axis) - 1], "dec_deg"].to_numpy() * u.deg, ) signed, perpendicular = explorer.gnomonic_projection_kpc(endpoint_coordinates) assert np.allclose(signed, [-50.0, 50.0], atol=3.0e-11) assert np.allclose(perpendicular, 0.0, atol=3.0e-11) assert axis["galactic_l_deg"].iloc[-1] < axis["galactic_l_deg"].iloc[0] assert set(axis["axis_positive_direction"]) == {"toward Galactic center"} def test_mw_foreground_is_nearly_flat_from_zero_to_30_kpc(axis: pd.DataFrame) -> None: values = axis.set_index("signed_axis_kpc") center = values.loc[ 0.0, "response_matched_total_absorbed_0p4_1p25_primary_fluxunit" ] plus_30 = values.loc[ 30.0, "response_matched_total_absorbed_0p4_1p25_primary_fluxunit" ] assert abs(plus_30 / center - 1.0) < 0.002 def test_actual_field_predictions_match_frozen_csv(fields: pd.DataFrame) -> None: frozen = pd.read_csv(explorer.FROZEN_REFERENCE_CSV, dtype={"obsid": str}) joined = fields.merge( frozen[ [ "obsid", "direct_em_target_apec_absorbed_0p4_1p25_primary_fluxunit", "reference_response_matched_absorbed_0p4_1p25_primary_fluxunit", ] ], on="obsid", suffixes=("_new", "_frozen"), validate="one_to_one", ) for stem in ( "direct_em_target_apec_absorbed_0p4_1p25_primary_fluxunit", "reference_response_matched_absorbed_0p4_1p25_primary_fluxunit", ): assert np.allclose( joined[f"{stem}_new"], joined[f"{stem}_frozen"], rtol=5.0e-15 ) assert np.allclose( fields["actual_nh_recomputed_response_matched_primary_fluxunit"], fields["reference_response_matched_absorbed_0p4_1p25_primary_fluxunit"], rtol=5.0e-15, ) def test_observed_columns_and_aliases_are_not_relabeled(fields: pd.DataFrame) -> None: source = pd.read_csv(explorer.PUBLIC_CSV, dtype={"obsid": str}) assert list(fields.columns[: len(source.columns)]) == list(source.columns) assert list(fields["obsid"]) == list(source["obsid"]) for column in source.columns: if pd.api.types.is_numeric_dtype(source[column]): assert np.allclose( fields[column], source[column], rtol=2.0e-15, atol=0.0, equal_nan=True ) else: assert ( fields[column].fillna("").tolist() == source[column].fillna("").tolist() ) aliases = { "cgmsum_kt_keV": "mwhalo_kt_keV", "cgmsum_kt_err_keV": "mwhalo_kt_err_keV", "cgmsum_norm": "mwhalo_norm", "cgmsum_norm_err": "mwhalo_norm_err", "cgmsum_cov_kt_norm": "mwhalo_cov_kt_norm", } for alias, historical in aliases.items(): assert np.allclose(fields[alias], fields[historical], equal_nan=True) def test_nominal_aperture_radial_extent_is_explicit(fields: pd.DataFrame) -> None: expected = explorer.M31_DISTANCE_KPC * np.tan(np.deg2rad(15.0 / 60.0)) assert expected == pytest.approx(3.403413640147222) assert np.allclose(fields["nominal_aperture_radius_arcmin"], 15.0) assert np.allclose(fields["nominal_aperture_radial_extent_kpc"], expected) assert np.allclose( fields["rproj_aperture_min_kpc"], np.maximum(0.0, fields["rproj_kpc"] - expected), ) assert np.allclose(fields["rproj_aperture_max_kpc"], fields["rproj_kpc"] + expected) def test_zhang_profile_and_grayson_bin_only_semantics() -> None: templates = pd.read_csv(PRODUCTS / "conditional_m31_templates.csv") zhang = templates.loc[ templates["template_id"].eq("zhang2024_m31_mass_beta_profile") ] at_20 = zhang.loc[np.isclose(zhang["radius_kpc"], 20.0)].iloc[0] assert at_20["primary_central"] == pytest.approx(explorer.ZHANG_S20_PRIMARY) assert zhang["radius_kpc"].min() == 10.0 assert zhang["radius_kpc"].max() == 30.0 grayson = templates.loc[templates["template_id"].str.startswith("grayson2025_")] assert len(grayson) == 5 assert grayson["radius_kpc"].isna().all() assert (grayson["bin_low_kpc"] == 10.0).all() assert (grayson["bin_high_kpc"] == 30.0).all() assert set(grayson["representation"]) == { "single_10_to_30_kpc_bin_band_no_interpolation" } def _png_dimensions(path: Path) -> tuple[int, int]: with path.open("rb") as handle: signature = handle.read(24) assert signature[:8] == b"\x89PNG\r\n\x1a\n" return struct.unpack(">II", signature[16:24]) def test_all_outputs_exist_and_dual_sizes_are_exact() -> None: outputs = explorer.expected_output_paths(HERE) assert len(outputs) == 21 assert all(path.is_file() and path.stat().st_size > 100 for path in outputs) for stem in explorer.FIGURE_STEMS: one_col = FIGURES / f"{stem}_1col.png" two_col = FIGURES / f"{stem}_2col.png" assert _png_dimensions(one_col)[0] == 1340 assert _png_dimensions(two_col)[0] == 2800 assert (FIGURES / f"{stem}_1col.pdf").read_bytes().startswith(b"%PDF") assert (FIGURES / f"{stem}_2col.pdf").read_bytes().startswith(b"%PDF") def test_summary_contains_closure_formula_constants_and_estimators() -> None: summary = json.loads( (PRODUCTS / "explorer_summary.json").read_text(encoding="utf-8") ) assert summary["constants"]["o8_response_matched_closure"] == pytest.approx( 0.9498760469603238 ) assert summary["constants"]["nominal_aperture_radial_extent_kpc"] == pytest.approx( 3.403413640147222 ) assert ( "no cross term" in summary["formula"]["density_components"]["emission_measure"] ) assert set(summary["estimators_primary_fluxunit"]) == { "observed_cgmsum", "actual_nh_direct_em", "actual_nh_response_matched", "fixed_nh_geometry_only_direct", "fixed_nh_geometry_only_response_matched", } assert summary["estimators_primary_fluxunit"]["actual_nh_response_matched"][ "all_inverse_variance" ] == pytest.approx(1.6017123021968118) assert summary["expected_products"] == [ path.relative_to(HERE).as_posix() for path in explorer.expected_output_paths(HERE) ] def test_summary_estimators_recompute_from_authoritative_field_rows( fields: pd.DataFrame, ) -> None: summary = json.loads( (PRODUCTS / "explorer_summary.json").read_text(encoding="utf-8") ) assert fields["side"].value_counts().to_dict() == {"North/NW": 10, "South/SE": 4} error = fields[ "absorbed_flux_soft_0p4_1p25_staterr_erg_cm-2_s-1_arcmin-2" ].to_numpy(dtype=float) columns = { "observed_cgmsum": fields[ "absorbed_flux_soft_0p4_1p25_erg_cm-2_s-1_arcmin-2" ].to_numpy(dtype=float) / explorer.PRIMARY_FLUX_UNIT, "actual_nh_response_matched": fields[ "reference_response_matched_absorbed_0p4_1p25_primary_fluxunit" ].to_numpy(dtype=float), } for label, values in columns.items(): expected = summary["estimators_primary_fluxunit"][label] weights = 1.0 / np.square(error) all_mean = float(np.sum(weights * values) / np.sum(weights)) side_means = {} for side in ("North/NW", "South/SE"): selected = fields["side"].eq(side).to_numpy() side_means[side] = float( np.sum(weights[selected] * values[selected]) / np.sum(weights[selected]) ) assert expected["all_inverse_variance"] == pytest.approx(all_mean) assert expected["north_inverse_variance"] == pytest.approx( side_means["North/NW"] ) assert expected["south_inverse_variance"] == pytest.approx( side_means["South/SE"] ) assert expected["side_balanced"] == pytest.approx( 0.5 * (side_means["North/NW"] + side_means["South/SE"]) ) model = summary["estimators_primary_fluxunit"]["actual_nh_response_matched"] assert model["side_balanced"] != pytest.approx(model["all_inverse_variance"])