diff --git a/autolens/weak/dataset.py b/autolens/weak/dataset.py index a349cadd3..737c14fbd 100644 --- a/autolens/weak/dataset.py +++ b/autolens/weak/dataset.py @@ -15,8 +15,12 @@ the same one pinned by ``PyAutoGalaxy`` PR #366: column 0 is :math:`\\gamma_2`, column 1 is :math:`\\gamma_1`, and the (y, x) galaxy positions are accessible via ``shear_yx.grid``. """ +import csv +from pathlib import Path from typing import List, Optional, Union +import numpy as np + import autoarray as aa from autogalaxy.util.shear_field import ShearYX2DIrregular @@ -28,6 +32,8 @@ def __init__( shear_yx: ShearYX2DIrregular, noise_map: Union[float, aa.ArrayIrregular, List[float]], name: str = "", + redshifts: Optional[Union[aa.ArrayIrregular, List[float]]] = None, + is_reduced: bool = False, ): """ A weak-lensing shear catalogue: a ``ShearYX2DIrregular`` shear field plus a per-galaxy noise map. @@ -45,6 +51,15 @@ def __init__( name Optional label, mirroring ``PointDataset.name``. Used by downstream fitting code to pair this dataset with a corresponding model component when multiple datasets are fitted simultaneously. + redshifts + Optional per-galaxy source redshifts (one value per galaxy). Stored for provenance and for + future per-galaxy lensing-efficiency (sigma_crit) scaling; the current fit assumes a single + effective source plane (the tracer's source redshift). + is_reduced + Whether the stored values are *reduced* shears ``g = gamma / (1 - kappa)`` — what real surveys + measure from galaxy ellipticities — rather than the shear ``gamma`` itself. ``FitWeak`` computes + the matching model quantity, so set this to match how the catalogue was measured. Simulated + datasets from ``SimulatorShearYX`` default to plain shear. """ self.name = name @@ -72,6 +87,18 @@ def __init__( self.noise_map = noise_map + if redshifts is not None and not isinstance(redshifts, aa.ArrayIrregular): + redshifts = aa.ArrayIrregular(values=list(redshifts)) + + if redshifts is not None and len(redshifts) != n_galaxies: + raise ValueError( + f"WeakDataset.redshifts has length {len(redshifts)} but shear_yx has " + f"{n_galaxies} entries; the two must match." + ) + + self.redshifts = redshifts + self.is_reduced = is_reduced + @property def positions(self) -> aa.Grid2DIrregular: """The (y, x) sky positions of the source galaxies the shear is measured at.""" @@ -90,8 +117,246 @@ def info(self) -> str: f"n_galaxies : {self.n_galaxies}\n" f"shear_yx : {self.shear_yx}\n" f"noise_map : {self.noise_map}\n" + f"redshifts : {self.redshifts}\n" + f"is_reduced : {self.is_reduced}\n" + ) + + @classmethod + def from_arrays( + cls, + positions, + gamma_1, + gamma_2, + noise_map: Optional[Union[float, List[float]]] = None, + weights: Optional[List[float]] = None, + redshifts: Optional[List[float]] = None, + is_reduced: bool = True, + name: str = "", + ) -> "WeakDataset": + """ + Build a ``WeakDataset`` from plain per-galaxy arrays — the shared entry point of the + catalogue loaders (:meth:`from_csv`, :meth:`from_fits`). + + Exactly one of ``noise_map`` / ``weights`` must be given: real catalogues quote either a + per-galaxy shear standard deviation or an inverse-variance weight, converted here via + ``sigma = weights**-0.5``. + + Parameters + ---------- + positions + The ``(N, 2)`` ``(y, x)`` galaxy positions in arc-seconds. + gamma_1 + The ``(N,)`` gamma_1 (axis-aligned) shear components. + gamma_2 + The ``(N,)`` gamma_2 (diagonal) shear components. + noise_map + Per-galaxy shear standard deviation (scalar broadcasts). + weights + Per-galaxy inverse-variance weights, the common alternative in lensing catalogues. + redshifts + Optional per-galaxy source redshifts. + is_reduced + Whether the components are reduced shears ``g`` (the default here — real catalogues + measure ellipticities, i.e. reduced shear) or plain shear ``gamma``. + """ + if (noise_map is None) == (weights is None): + raise ValueError( + "WeakDataset.from_arrays requires exactly one of noise_map or weights." + ) + + if weights is not None: + weights = np.asarray(weights, dtype=float) + if np.any(weights <= 0.0): + raise ValueError( + "WeakDataset.from_arrays weights must all be positive to convert " + "to noise via sigma = weights**-0.5." + ) + noise_map = list(weights**-0.5) + + grid = aa.Grid2DIrregular(values=np.asarray(positions, dtype=float)) + values = np.stack( + [np.asarray(gamma_2, dtype=float), np.asarray(gamma_1, dtype=float)], + axis=1, + ) + shear_yx = ShearYX2DIrregular(values=values, grid=grid) + + return cls( + shear_yx=shear_yx, + noise_map=noise_map, + name=name, + redshifts=redshifts, + is_reduced=is_reduced, + ) + + def to_csv(self, file_path: Union[str, Path]): + """ + Write this dataset to ``file_path`` as a CSV with one row per galaxy. + + Columns: ``name, y, x, gamma_2, gamma_1, noise`` plus ``redshift`` when this dataset + carries per-galaxy redshifts. Whether the components are reduced shear is not a column — + pass ``is_reduced`` when loading (:meth:`from_csv` defaults to this dataset's convention + only if you tell it to). + """ + fieldnames = ["name", "y", "x", "gamma_2", "gamma_1", "noise"] + if self.redshifts is not None: + fieldnames.append("redshift") + + values = np.asarray(self.shear_yx) + positions = np.asarray(self.positions) + noise = np.asarray(self.noise_map) + + with open(file_path, "w", newline="") as f: + writer = csv.DictWriter(f, fieldnames=fieldnames) + writer.writeheader() + for i in range(self.n_galaxies): + row = { + "name": self.name, + "y": positions[i, 0], + "x": positions[i, 1], + "gamma_2": values[i, 0], + "gamma_1": values[i, 1], + "noise": noise[i], + } + if self.redshifts is not None: + row["redshift"] = np.asarray(self.redshifts)[i] + writer.writerow(row) + + @classmethod + def from_csv( + cls, + file_path: Union[str, Path], + is_reduced: bool = True, + name: Optional[str] = None, + ) -> "WeakDataset": + """ + Load a ``WeakDataset`` from a CSV with columns ``name, y, x, gamma_2, gamma_1, noise`` + (optionally ``redshift``), as written by :meth:`to_csv`. + + Parameters + ---------- + file_path + Path to the CSV file. + is_reduced + Whether the catalogue's components are reduced shears ``g`` (the default — real + catalogues measure ellipticities) or plain shear ``gamma``. + name + When given, only rows whose ``name`` column matches are loaded. + """ + positions, gamma_1, gamma_2, noise, redshifts = [], [], [], [], [] + row_name = name or "" + + with open(file_path, newline="") as f: + reader = csv.DictReader(f) + has_redshift = "redshift" in (reader.fieldnames or []) + for row in reader: + if name is not None and row["name"] != name: + continue + row_name = row["name"] + positions.append((float(row["y"]), float(row["x"]))) + gamma_1.append(float(row["gamma_1"])) + gamma_2.append(float(row["gamma_2"])) + noise.append(float(row["noise"])) + if has_redshift: + redshifts.append(float(row["redshift"])) + + if len(positions) == 0: + raise ValueError( + f"CSV file {str(file_path)!r} contained no WeakDataset rows" + + (f" with name {name!r}." if name is not None else ".") + ) + + return cls.from_arrays( + positions=positions, + gamma_1=gamma_1, + gamma_2=gamma_2, + noise_map=noise, + redshifts=redshifts if redshifts else None, + is_reduced=is_reduced, + name=row_name, ) + @classmethod + def from_fits( + cls, + file_path: Union[str, Path], + y_col: str, + x_col: str, + gamma_1_col: str, + gamma_2_col: str, + noise_col: Optional[str] = None, + weight_col: Optional[str] = None, + redshift_col: Optional[str] = None, + hdu: int = 1, + is_reduced: bool = True, + name: str = "", + ) -> "WeakDataset": + """ + Load a ``WeakDataset`` from a FITS binary table — the standard distribution format of + real shear catalogues — with explicit column-name mapping. + + Exactly one of ``noise_col`` / ``weight_col`` must be given (inverse-variance weights are + converted via ``sigma = weight**-0.5``). Positions are read as-is: convert RA/Dec to + projected arc-second offsets about the cluster centre *before* writing the table, or use + a catalogue whose columns are already tangent-plane coordinates. + + Parameters + ---------- + file_path + Path to the FITS file. + y_col, x_col + Names of the position columns (arc-seconds, (y, x) convention). + gamma_1_col, gamma_2_col + Names of the shear-component columns. + noise_col + Name of the per-galaxy shear standard-deviation column. + weight_col + Name of the per-galaxy inverse-variance weight column. + redshift_col + Optional name of a per-galaxy source-redshift column. + hdu + The FITS extension holding the binary table. + is_reduced + Whether the components are reduced shears ``g`` (the default — real catalogues + measure ellipticities) or plain shear ``gamma``. + """ + from astropy.io import fits as astropy_fits + + if (noise_col is None) == (weight_col is None): + raise ValueError( + "WeakDataset.from_fits requires exactly one of noise_col or weight_col." + ) + + with astropy_fits.open(file_path) as hdul: + table = hdul[hdu].data + + positions = np.stack( + [np.asarray(table[y_col], dtype=float), np.asarray(table[x_col], dtype=float)], + axis=1, + ) + + return cls.from_arrays( + positions=positions, + gamma_1=np.asarray(table[gamma_1_col], dtype=float), + gamma_2=np.asarray(table[gamma_2_col], dtype=float), + noise_map=( + list(np.asarray(table[noise_col], dtype=float)) + if noise_col is not None + else None + ), + weights=( + list(np.asarray(table[weight_col], dtype=float)) + if weight_col is not None + else None + ), + redshifts=( + list(np.asarray(table[redshift_col], dtype=float)) + if redshift_col is not None + else None + ), + is_reduced=is_reduced, + name=name, + ) + def extent_from(self, buffer: float = 0.1) -> List[float]: """The axis-aligned bounding box of the source-galaxy positions, padded by ``buffer`` on each side.""" positions = self.positions diff --git a/autolens/weak/fit.py b/autolens/weak/fit.py index 2389278e1..2fe111e72 100644 --- a/autolens/weak/fit.py +++ b/autolens/weak/fit.py @@ -42,10 +42,28 @@ def __init__(self, dataset: WeakDataset, tracer): @cached_property def model_shear(self) -> ShearYX2DIrregular: - """The model shear field evaluated at the galaxy positions, via ``LensCalc``.""" - return LensCalc.from_tracer(self.tracer).shear_yx_2d_via_hessian_from( + """ + The model shear field evaluated at the galaxy positions, via ``LensCalc``. + + When the dataset is marked ``is_reduced`` (real catalogues measure galaxy ellipticities, + i.e. the reduced shear) the model quantity is ``g = gamma / (1 - kappa)``, with the + convergence from the same Hessian primitive — so data and model always live in the same + space. + """ + lens_calc = LensCalc.from_tracer(self.tracer) + + shear = lens_calc.shear_yx_2d_via_hessian_from(grid=self.dataset.positions) + + if not getattr(self.dataset, "is_reduced", False): + return shear + + convergence = lens_calc.convergence_2d_via_hessian_from( grid=self.dataset.positions ) + return ShearYX2DIrregular( + values=np.asarray(shear) / (1.0 - np.asarray(convergence))[:, None], + grid=self.dataset.positions, + ) @property def residual_map(self) -> np.ndarray: diff --git a/autolens/weak/simulator.py b/autolens/weak/simulator.py index 292fd62ae..8e848699f 100644 --- a/autolens/weak/simulator.py +++ b/autolens/weak/simulator.py @@ -27,7 +27,12 @@ class SimulatorShearYX: - def __init__(self, noise_sigma: float = 0.3, seed: Optional[int] = None): + def __init__( + self, + noise_sigma: float = 0.3, + seed: Optional[int] = None, + reduced: bool = False, + ): """ Simulator for weak-lensing shear catalogues. @@ -40,9 +45,14 @@ def __init__(self, noise_sigma: float = 0.3, seed: Optional[int] = None): seed Optional seed for the random number generator. When set, both the noise and the random-position draws (for ``via_tracer_random_positions_from``) are reproducible. + reduced + When ``True`` the simulator produces *reduced* shears ``g = gamma / (1 - kappa)`` — the quantity + real surveys measure from galaxy ellipticities — and marks the output dataset ``is_reduced``. + The default plain-shear mode matches the earlier tutorials in the weak-lensing series. """ self.noise_sigma = float(noise_sigma) self.seed = seed + self.reduced = reduced self._rng = np.random.default_rng(seed) def via_tracer_from( @@ -72,6 +82,13 @@ def via_tracer_from( true_shear = self._true_shear_yx_from(tracer=tracer, grid=grid) + if self.reduced: + convergence = self._convergence_from(tracer=tracer, grid=grid) + true_shear = ShearYX2DIrregular( + values=np.asarray(true_shear) / (1.0 - np.asarray(convergence))[:, None], + grid=grid, + ) + if self.noise_sigma > 0.0: noise = self._rng.normal( loc=0.0, scale=self.noise_sigma, size=true_shear.shape @@ -86,7 +103,12 @@ def via_tracer_from( values=[self.noise_sigma] * len(grid) ) - return WeakDataset(shear_yx=shear_yx, noise_map=noise_map, name=name) + return WeakDataset( + shear_yx=shear_yx, + noise_map=noise_map, + name=name, + is_reduced=self.reduced, + ) def via_tracer_random_positions_from( self, @@ -125,6 +147,19 @@ def via_tracer_random_positions_from( grid = aa.Grid2DIrregular(values=positions) return self.via_tracer_from(tracer=tracer, grid=grid, name=name) + @staticmethod + def _convergence_from(tracer, grid: aa.Grid2DIrregular): + """ + Evaluate the convergence of ``tracer`` on ``grid`` via the same Hessian primitive and + input-dispatch rules as ``_true_shear_yx_from``. + """ + method = getattr(tracer, "convergence_2d_via_hessian_from", None) + if method is not None: + return method(grid=grid) + if hasattr(tracer, "deflections_between_planes_from"): + return LensCalc.from_tracer(tracer).convergence_2d_via_hessian_from(grid=grid) + return LensCalc.from_mass_obj(tracer).convergence_2d_via_hessian_from(grid=grid) + @staticmethod def _true_shear_yx_from(tracer, grid: aa.Grid2DIrregular): """ diff --git a/test_autolens/weak/test_dataset_io.py b/test_autolens/weak/test_dataset_io.py new file mode 100644 index 000000000..bf17ba91f --- /dev/null +++ b/test_autolens/weak/test_dataset_io.py @@ -0,0 +1,186 @@ +import numpy as np +import pytest + +import autolens as al + +from autolens.weak.fit import FitWeak + + +def _sis_tracer(einstein_radius=1.6, centre=(0.0, 0.0)): + lens = al.Galaxy( + redshift=0.5, + mass=al.mp.IsothermalSph(centre=centre, einstein_radius=einstein_radius), + ) + return al.Tracer(galaxies=[lens, al.Galaxy(redshift=1.0)]) + + +def _dataset(**kwargs): + return al.SimulatorShearYX( + noise_sigma=kwargs.pop("noise_sigma", 0.3), seed=1, **kwargs + ).via_tracer_random_positions_from( + tracer=_sis_tracer(), n_galaxies=20, grid_extent=3.0, name="io_test" + ) + + +def test__from_arrays__noise_or_weights_exclusively(): + positions = [(1.0, 0.0), (0.0, 1.0)] + + dataset = al.WeakDataset.from_arrays( + positions=positions, + gamma_1=[0.1, -0.1], + gamma_2=[0.05, 0.02], + weights=[25.0, 100.0], + ) + + # sigma = weights**-0.5 + np.testing.assert_allclose(np.asarray(dataset.noise_map), [0.2, 0.1]) + assert dataset.is_reduced is True # loader default: real catalogues are reduced + + # [gamma_2, gamma_1] storage convention round-trips. + np.testing.assert_allclose(np.asarray(dataset.shear_yx)[:, 1], [0.1, -0.1]) + np.testing.assert_allclose(np.asarray(dataset.shear_yx)[:, 0], [0.05, 0.02]) + + with pytest.raises(ValueError): + al.WeakDataset.from_arrays( + positions=positions, gamma_1=[0.1, -0.1], gamma_2=[0.05, 0.02] + ) + with pytest.raises(ValueError): + al.WeakDataset.from_arrays( + positions=positions, + gamma_1=[0.1, -0.1], + gamma_2=[0.05, 0.02], + noise_map=[0.3, 0.3], + weights=[25.0, 100.0], + ) + + +def test__csv_round_trip(tmp_path): + dataset = _dataset() + dataset.redshifts = None + + file_path = tmp_path / "catalogue.csv" + dataset.to_csv(file_path) + + loaded = al.WeakDataset.from_csv(file_path, is_reduced=False) + + assert loaded.name == "io_test" + assert loaded.is_reduced is False + np.testing.assert_allclose(np.asarray(loaded.shear_yx), np.asarray(dataset.shear_yx)) + np.testing.assert_allclose(np.asarray(loaded.positions), np.asarray(dataset.positions)) + np.testing.assert_allclose(np.asarray(loaded.noise_map), np.asarray(dataset.noise_map)) + assert loaded.redshifts is None + + +def test__csv_round_trip__with_redshifts(tmp_path): + dataset = al.WeakDataset.from_arrays( + positions=[(1.0, 0.0), (0.0, 1.0)], + gamma_1=[0.1, -0.1], + gamma_2=[0.05, 0.02], + noise_map=[0.3, 0.25], + redshifts=[1.2, 0.9], + name="with_z", + ) + + file_path = tmp_path / "catalogue.csv" + dataset.to_csv(file_path) + loaded = al.WeakDataset.from_csv(file_path) + + np.testing.assert_allclose(np.asarray(loaded.redshifts), [1.2, 0.9]) + assert loaded.is_reduced is True + + +def test__from_fits__column_mapping_and_weights(tmp_path): + from astropy.io import fits as astropy_fits + from astropy.table import Table + + table = Table( + { + "y_arcsec": [1.0, -1.0, 0.5], + "x_arcsec": [0.0, 1.0, -0.5], + "e1": [0.1, -0.05, 0.02], + "e2": [0.03, 0.07, -0.01], + "weight": [25.0, 100.0, 4.0], + "z_source": [1.1, 0.8, 1.4], + } + ) + file_path = tmp_path / "catalogue.fits" + astropy_fits.BinTableHDU(table).writeto(file_path) + + dataset = al.WeakDataset.from_fits( + file_path=file_path, + y_col="y_arcsec", + x_col="x_arcsec", + gamma_1_col="e1", + gamma_2_col="e2", + weight_col="weight", + redshift_col="z_source", + name="fits_test", + ) + + assert dataset.n_galaxies == 3 + assert dataset.is_reduced is True + np.testing.assert_allclose(np.asarray(dataset.noise_map), [0.2, 0.1, 0.5]) + np.testing.assert_allclose(np.asarray(dataset.shear_yx)[:, 1], [0.1, -0.05, 0.02]) + np.testing.assert_allclose(np.asarray(dataset.redshifts), [1.1, 0.8, 1.4]) + + +def test__json_round_trip__preserves_new_fields(tmp_path): + dataset = al.WeakDataset.from_arrays( + positions=[(1.0, 0.0), (0.0, 1.0)], + gamma_1=[0.1, -0.1], + gamma_2=[0.05, 0.02], + noise_map=[0.3, 0.25], + redshifts=[1.2, 0.9], + is_reduced=True, + name="json_test", + ) + + file_path = tmp_path / "dataset.json" + al.output_to_json(obj=dataset, file_path=file_path) + loaded = al.from_json(file_path=file_path) + + assert loaded.is_reduced is True + np.testing.assert_allclose(np.asarray(loaded.redshifts), [1.2, 0.9]) + + +def test__reduced_shear__round_trip_and_differs_from_shear(): + """ + A reduced-shear simulation fitted by its truth tracer must round-trip to zero residuals + (FitWeak computes the model reduced shear when the dataset declares is_reduced), and the + reduced/plain quantities must differ where the convergence is non-negligible. + """ + tracer = _sis_tracer() + + dataset_reduced = al.SimulatorShearYX( + noise_sigma=0.0, seed=2, reduced=True + ).via_tracer_random_positions_from(tracer=tracer, n_galaxies=50, grid_extent=3.0) + + assert dataset_reduced.is_reduced is True + + fit = FitWeak(dataset=dataset_reduced, tracer=tracer) + np.testing.assert_allclose(fit.residual_map, 0.0, atol=1e-10) + + dataset_plain = al.SimulatorShearYX( + noise_sigma=0.0, seed=2, reduced=False + ).via_tracer_random_positions_from(tracer=tracer, n_galaxies=50, grid_extent=3.0) + + # For an SIS, |g| = |gamma| / (1 - kappa) with kappa = theta_E / 2r > 0, so the + # reduced values must be strictly larger in magnitude at every galaxy. + reduced_mag = np.linalg.norm(np.asarray(dataset_reduced.shear_yx), axis=1) + plain_mag = np.linalg.norm(np.asarray(dataset_plain.shear_yx), axis=1) + assert (reduced_mag > plain_mag).all() + + +def test__reduced_mismatch__truth_tracer_no_longer_perfect(): + """A plain-shear dataset fitted as if reduced (or vice versa) must NOT round-trip — guards + against the flag being ignored somewhere in the chain.""" + tracer = _sis_tracer() + + dataset_plain = al.SimulatorShearYX( + noise_sigma=0.0, seed=3, reduced=False + ).via_tracer_random_positions_from(tracer=tracer, n_galaxies=30, grid_extent=3.0) + + dataset_plain.is_reduced = True # deliberately mislabel + + fit = FitWeak(dataset=dataset_plain, tracer=tracer) + assert np.abs(fit.residual_map).max() > 1e-3