From 8e154998ff4af5d44dcd309f156ae8f69529ff16 Mon Sep 17 00:00:00 2001 From: John Franklin Crenshaw <41785729+jfcrenshaw@users.noreply.github.com> Date: Thu, 29 Feb 2024 00:57:43 -0800 Subject: [PATCH] v1.2 (#13) * Updates for v1.2 * Oops. Didn't add lsstV1/V2. * No longer forcing random_state. * Linting tests. * Linting --- .github/workflows/main.yml | 2 +- README.md | 10 +-- photerr/__init__.py | 12 ++- photerr/euclid.py | 21 ++--- photerr/lsstV1.py | 117 +++++++++++++++++++++++++++ photerr/lsstV2.py | 138 ++++++++++++++++++++++++++++++++ photerr/model.py | 91 +++++++++++++++------ photerr/params.py | 159 +++++++++++++++++++++++++------------ photerr/roman.py | 23 ++---- pyproject.toml | 50 ++++++------ tests/test_models.py | 21 ++++- tests/test_params.py | 18 +++++ 12 files changed, 519 insertions(+), 143 deletions(-) create mode 100644 photerr/lsstV1.py create mode 100644 photerr/lsstV2.py diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index aaca91d..cef693c 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -43,7 +43,7 @@ jobs: fail-fast: true matrix: os: ["ubuntu-latest", "macos-latest"] - python-version: ["3.8", "3.9", "3.10"] + python-version: ["3.10", "3.11", "3.12"] runs-on: ${{ matrix.os }} steps: diff --git a/README.md b/README.md index 57570cf..50361a6 100644 --- a/README.md +++ b/README.md @@ -88,8 +88,8 @@ If you are changing other dictionary-parameters at the same time (e.g. `nVisYr`, By default, PhotErr tries to use the provided information to calculate limiting magnitudes for you. If you would like to directly supply your own $5\sigma$ limits, you can do so using the `m5` parameter. -Note that PhotErr assumes these are single-visit limiting magnitudes. -If you want to supply coadded depths, you should also set `nYrObs=1` and `nVisYr={u: 1, g: 1, ...}`, so that the calculated coadded depths are equal to those you provided. +Note PhotErr assumes these are single-visit point-source limiting magnitudes. +If you want to supply coadded depths, you should also set `nYrObs=1` and `nVisYr=1`, so that the calculated coadded depths are equal to those you provided. ### *Handling non-detections* @@ -120,12 +120,12 @@ You can also calculate limiting magnitudes for apertures of a given size by pass ### Scaling the errors -If you want to scale up or scale down the errors in any band(s), you can use the keyword `scale`. +If you want to scale up or scale down the errors in any band(s), you can use the keyword `scale`. For example, `LsstErrorModel(scale={"u": 2, "y": 2})` will have all the same properties of the default error model, except the errors in the `u` and `y` bands will be doubled. This allows you to answer questions like "what happens to my science if the `u` band errors are doubled." Note this only scales the band-specific error. -The band-independent systematic error floor, `sigmaSys` is still the same, and so at high-SNR, near the systematic floor, the errors won't scale as you expect. +The band-independent systematic error floor, `sigmaSys` is still the same, and so at high-SNR near the systematic floor, the errors won't scale as you expect. ### *Other error models* @@ -169,7 +169,7 @@ In PhotErr, we do not make this approximation so that the error model generalize In addition, while the high-SNR model assumes photometric errors are Gaussian in magnitude space, we model errors as Gaussian in flux space. However, both of these high-SNR approximations can be restored with the keyword `highSNR=True`. -The LSST error model uses the parameters from [Ivezic (2019)](https://arxiv.org/abs/0805.2366). +The LSST error model uses parameters from [Ivezic (2019)](https://arxiv.org/abs/0805.2366), [Bianco 2022](https://pstn-054.lsst.io), and from [this Rubin systems engineering notebook](https://github.com/lsst-pst/syseng_throughputs/blob/main/notebooks/EvalReadnoise.ipynb). The Euclid and Roman error models follow [Graham (2020)](https://arxiv.org/abs/2004.07885) in setting $\gamma = 0.04$, which is nearly identical to the values for Rubin (which are all $\sim 0.039$). In addition to the random photometric error above, we include a system error of $\sigma_\text{sys} = 0.005$ which is added in quadrature to random error. Note that the system error can be changed using the keyword `sigmaSys`. diff --git a/photerr/__init__.py b/photerr/__init__.py index 202803d..7dd2d16 100644 --- a/photerr/__init__.py +++ b/photerr/__init__.py @@ -2,10 +2,20 @@ from importlib import metadata from .euclid import EuclidErrorModel, EuclidErrorParams -from .lsst import LsstErrorModel, LsstErrorParams +from .lsstV1 import LsstErrorModelV1, LsstErrorParamsV1 +from .lsstV2 import LsstErrorModelV2, LsstErrorParamsV2 from .model import ErrorModel from .params import ErrorParams from .roman import RomanErrorModel, RomanErrorParams +# set the version number __version__ = metadata.version(__package__) del metadata + +# alias the latest LSST error model +class LsstErrorParams(LsstErrorParamsV2): + ... + + +class LsstErrorModel(LsstErrorModelV2): + ... diff --git a/photerr/euclid.py b/photerr/euclid.py index 56995b3..4564805 100644 --- a/photerr/euclid.py +++ b/photerr/euclid.py @@ -1,6 +1,7 @@ """Photometric error model for Euclid.""" + from dataclasses import dataclass, field -from typing import Any, Dict +from typing import Any from photerr.model import ErrorModel from photerr.params import ErrorParams, param_docstring @@ -18,22 +19,10 @@ class EuclidErrorParams(ErrorParams): __doc__ += " Graham 2020 - https://arxiv.org/abs/2004.07885" nYrObs: float = 1.0 - nVisYr: Dict[str, float] = field( - default_factory=lambda: { - "Y": 1.0, - "J": 1.0, - "H": 1.0, - } - ) - gamma: Dict[str, float] = field( - default_factory=lambda: { - "Y": 0.04, - "J": 0.04, - "H": 0.04, - } - ) + nVisYr: dict[str, float] | float = 1.0 + gamma: dict[str, float] | float = 0.04 - m5: Dict[str, float] = field( + m5: dict[str, float] | float = field( default_factory=lambda: { "Y": 24.0, "J": 24.2, diff --git a/photerr/lsstV1.py b/photerr/lsstV1.py new file mode 100644 index 0000000..a7355c5 --- /dev/null +++ b/photerr/lsstV1.py @@ -0,0 +1,117 @@ +"""Photometric error model for LSST using defaults from Ivezic 2019.""" + +from dataclasses import dataclass, field +from typing import Any + +from photerr.model import ErrorModel +from photerr.params import ErrorParams, param_docstring + + +@dataclass +class LsstErrorParamsV1(ErrorParams): + """Parameters for the LSST photometric error model. + + Default values taken from pages 11, 12, 26 of Ivezic 2019. + """ + + __doc__ += param_docstring + + nYrObs: float = 10.0 + nVisYr: dict[str, float] | float = field( + default_factory=lambda: { + "u": 5.6, + "g": 8.0, + "r": 18.4, + "i": 18.4, + "z": 16.0, + "y": 16.0, + } + ) + gamma: dict[str, float] | float = field( + default_factory=lambda: { + "u": 0.038, + "g": 0.039, + "r": 0.039, + "i": 0.039, + "z": 0.039, + "y": 0.039, + } + ) + + m5: dict[str, float] | float = field(default_factory=lambda: {}) + + tvis: dict[str, float] | float = 30 + airmass: dict[str, float] | float = 1.2 + Cm: dict[str, float] | float = field( + default_factory=lambda: { + "u": 23.09, + "g": 24.42, + "r": 24.44, + "i": 24.32, + "z": 24.16, + "y": 23.73, + } + ) + dCmInf: dict[str, float] | float = field( + default_factory=lambda: { + "u": 0.62, + "g": 0.18, + "r": 0.10, + "i": 0.07, + "z": 0.05, + "y": 0.04, + } + ) + msky: dict[str, float] | float = field( + default_factory=lambda: { + "u": 22.99, + "g": 22.26, + "r": 21.20, + "i": 20.48, + "z": 19.60, + "y": 18.61, + } + ) + mskyDark: dict[str, float] | float = field( + default_factory=lambda: { + "u": 22.99, + "g": 22.26, + "r": 21.20, + "i": 20.48, + "z": 19.60, + "y": 18.61, + } + ) + theta: dict[str, float] | float = field( + default_factory=lambda: { + "u": 0.92, + "g": 0.87, + "r": 0.83, + "i": 0.80, + "z": 0.78, + "y": 0.76, + } + ) + km: dict[str, float] | float = field( + default_factory=lambda: { + "u": 0.491, + "g": 0.213, + "r": 0.126, + "i": 0.096, + "z": 0.069, + "y": 0.170, + } + ) + + tvisRef: float | None = 30 + + +class LsstErrorModelV1(ErrorModel): + """Photometric error model for Euclid.""" + + def __init__(self, **kwargs: Any) -> None: + """Create an LSST error model. + + Keyword arguments override default values in LsstErrorParams. + """ + super().__init__(LsstErrorParamsV1(**kwargs)) diff --git a/photerr/lsstV2.py b/photerr/lsstV2.py new file mode 100644 index 0000000..d8ae469 --- /dev/null +++ b/photerr/lsstV2.py @@ -0,0 +1,138 @@ +"""Photometric error model for LSST as of Feb 28, 2024.""" + +from dataclasses import dataclass, field +from typing import Any + +from photerr.model import ErrorModel +from photerr.params import ErrorParams, param_docstring + + +@dataclass +class LsstErrorParamsV2(ErrorParams): + """Parameters for the LSST photometric error model. + + Default values match v1.9 of the official LSST throughputs, including the + silver coatings on all three mirrors. The values for tvis and nVisYr were + taken from one of the OpSim runs that Lynne Jones said is likely to resemble + the final exposure times and visit distributions. + """ + + __doc__ += param_docstring + + nYrObs: float = 10.0 + nVisYr: dict[str, float] | float = field( + default_factory=lambda: { + "u": 6.3, + "g": 6.8, + "r": 18.1, + "i": 18.8, + "z": 16.6, + "y": 16.3, + } + ) + gamma: dict[str, float] | float = field( + default_factory=lambda: { + "u": 0.038, + "g": 0.039, + "r": 0.039, + "i": 0.039, + "z": 0.039, + "y": 0.039, + } + ) + + m5: dict[str, float] | float = field(default_factory=lambda: {}) + + tvis: dict[str, float] | float = field( + default_factory=lambda: { + "u": 38, + "g": 30, + "r": 30, + "i": 30, + "z": 30, + "y": 30, + } + ) + airmass: dict[str, float] | float = field( + default_factory=lambda: { + "u": 1.15, + "g": 1.15, + "r": 1.14, + "i": 1.15, + "z": 1.16, + "y": 1.16, + } + ) + Cm: dict[str, float] | float = field( + default_factory=lambda: { + "u": 22.968, + "g": 24.582, + "r": 24.602, + "i": 24.541, + "z": 24.371, + "y": 23.840, + } + ) + dCmInf: dict[str, float] | float = field( + default_factory=lambda: { + "u": 0.543, + "g": 0.088, + "r": 0.043, + "i": 0.028, + "z": 0.019, + "y": 0.016, + } + ) + msky: dict[str, float] | float = field( + default_factory=lambda: { + "u": 22.65, + "g": 22.07, + "r": 21.10, + "i": 19.99, + "z": 19.04, + "y": 18.28, + } + ) + mskyDark: dict[str, float] | float = field( + default_factory=lambda: { + "u": 23.052, + "g": 22.254, + "r": 21.198, + "i": 20.463, + "z": 19.606, + "y": 18.602, + } + ) + theta: dict[str, float] | float = field( + default_factory=lambda: { + "u": 1.22, + "g": 1.09, + "r": 1.02, + "i": 0.99, + "z": 1.01, + "y": 0.99, + } + ) + km: dict[str, float] | float = field( + default_factory=lambda: { + "u": 0.470, + "g": 0.213, + "r": 0.126, + "i": 0.096, + "z": 0.068, + "y": 0.171, + } + ) + + tvisRef: float | None = 30 + + +class LsstErrorModelV2(ErrorModel): + """Photometric error model for Euclid.""" + + def __init__(self, **kwargs: Any) -> None: + """Create an LSST error model. + + Keyword arguments override default values in LsstErrorParams. + """ + super().__init__(LsstErrorParamsV2(**kwargs)) diff --git a/photerr/model.py b/photerr/model.py index 97be8b9..adaebb6 100644 --- a/photerr/model.py +++ b/photerr/model.py @@ -1,5 +1,6 @@ """The photometric error model.""" -from typing import Any, Tuple, Union + +from typing import Any import numpy as np import pandas as pd @@ -51,7 +52,7 @@ def __init__(self, *args: ErrorParams, **kwargs: Any) -> None: self._params.update(**kwargs) # make a list of the bands - self._bands = [band for band in self._params.nVisYr.keys()] + self._bands = list(self._params.nVisYr.keys()) # calculate all of the 5-sigma limiting magnitudes self._calculate_m5() @@ -64,25 +65,46 @@ def params(self) -> ErrorParams: def _calculate_m5(self) -> None: """Calculate the single-visit 5-sigma limiting magnitudes. - Uses Eq. 6 from Ivezic 2019. - However, if m5 has been explicitly passed for any bands, those values are - used instead. + Uses Eq. 6 from Ivezic 2019 and Eq. 2-3 of Bianco 2022. + However, if m5 has been explicitly passed for any bands, + those values are used instead. """ - # calculate the m5 limiting magnitudes using Eq. 6 - m5 = { - band: self.params.Cm[band] - + 0.50 * (self.params.msky[band] - 21) - + 2.5 * np.log10(0.7 / self.params.theta[band]) - + 1.25 * np.log10(self.params.tvis / 30) - - self.params.km[band] * (self.params.airmass - 1) - for band in self.params.Cm - } + # calculate m5 for the bands with the requisite info + m5 = {} + for band in self.params.Cm: + # calculate the scaled visit time using Eq. 3 of Bianco 2022 + tau = ( + self.params.tvis[band] + / self.params.tvisRef + * 10 ** ((self.params.msky[band] - self.params.mskyDark[band]) / -2.5) + ) + + # calculate Cm exposure time adjustment using Eq. 2 from Bianco 2022 + dCmTau = self.params.dCmInf[band] - 1.25 * np.log10( + 1 + (10 ** (0.8 * self.params.dCmInf[band]) - 1) / tau + ) + + # calculate m5 using Eq. 6 from Ivezic 2019 and Eq. 2 from Bianco 2022 + m5[band] = ( + self.params.Cm[band] + + 0.50 * (self.params.msky[band] - 21) + + 2.5 * np.log10(0.7 / self.params.theta[band]) + + 1.25 * np.log10(self.params.tvis[band] / 30) + - self.params.km[band] * (self.params.airmass[band] - 1) + + dCmTau + ) + + # add the explicitly passed m5's to our dictionary m5.update(self.params.m5) + # save all the single-visit m5's together self._all_m5 = {band: m5[band] for band in self._bands} def _get_area_ratio_auto( - self, majors: np.ndarray, minors: np.ndarray, bands: list + self, + majors: np.ndarray, + minors: np.ndarray, + bands: list, ) -> np.ndarray: """Get the ratio between PSF area and galaxy aperture area for "auto" model. @@ -102,6 +124,8 @@ def _get_area_ratio_auto( """ # get the psf size for each band psf_size = np.array([self.params.theta[band] for band in bands]) + airmass = np.array([self.params.airmass[band] for band in bands]) + psf_size *= airmass**0.6 # convert PSF FWHM to Gaussian sigma psf_sig = psf_size / 2.355 @@ -118,7 +142,10 @@ def _get_area_ratio_auto( return A_ap / A_psf def _get_area_ratio_gaap( - self, majors: np.ndarray, minors: np.ndarray, bands: list + self, + majors: np.ndarray, + minors: np.ndarray, + bands: list, ) -> np.ndarray: """Get the ratio between PSF area and galaxy aperture area for "gaap" model. @@ -138,6 +165,8 @@ def _get_area_ratio_gaap( """ # get the psf size for each band psf_size = np.array([self.params.theta[band] for band in bands]) + airmass = np.array([self.params.airmass[band] for band in bands]) + psf_size *= airmass**0.6 # convert PSF FWHM to Gaussian sigma psf_sig = psf_size / 2.355 @@ -165,7 +194,11 @@ def _get_area_ratio_gaap( return A_ap / A_psf def _get_nsr_from_mags( - self, mags: np.ndarray, majors: np.ndarray, minors: np.ndarray, bands: list + self, + mags: np.ndarray, + majors: np.ndarray, + minors: np.ndarray, + bands: list, ) -> np.ndarray: """Calculate the noise-to-signal ratio. @@ -320,7 +353,7 @@ def _get_obs_and_errs( minors: np.ndarray, bands: list, rng: np.random.Generator, - ) -> Tuple[np.ndarray, np.ndarray]: + ) -> tuple[np.ndarray, np.ndarray]: """Calculate the noise-to-signal ratio. Uses Eqs 4 and 5 from Ivezic 2019. @@ -391,7 +424,9 @@ def _get_obs_and_errs( return obsMags, obsMagErrs def __call__( - self, catalog: pd.DataFrame, random_state: Union[np.random.Generator, int] + self, + catalog: pd.DataFrame, + random_state: np.random.Generator | int | None = None, ) -> pd.DataFrame: """Calculate photometric errors for the catalog and return in DataFrame. @@ -399,10 +434,10 @@ def __call__( ---------- catalog : pd.DataFrame The input catalog of galaxies in a pandas DataFrame. - random_sate : np.random.Generator or int + random_sate : np.random.Generator, int, or None The random state. Can either be a numpy random generator - (e.g. np.random.default_rng(42)), or an integer, which is then used - to seed np.random.default_rng. + (e.g. np.random.default_rng(42)), an integer (which is used + to seed np.random.default_rng), or None. Returns ------- @@ -411,10 +446,12 @@ def __call__( # set the rng if isinstance(random_state, np.random.Generator): rng = random_state - elif isinstance(random_state, int): + elif isinstance(random_state, int) or random_state is None: rng = np.random.default_rng(random_state) else: - raise TypeError("random_state must be a numpy random generator or an int.") + raise TypeError( + "random_state must be a numpy random generator, an int, or None." + ) # get the bands we will calculate errors for bands = [band for band in catalog.columns if band in self._bands] @@ -468,7 +505,10 @@ def __call__( return obsCatalog def getLimitingMags( - self, nSigma: float = 5, coadded: bool = True, aperture: float = 0 + self, + nSigma: float = 5, + coadded: bool = True, + aperture: float = 0, ) -> dict: """Return the limiting magnitudes in all bands. @@ -507,4 +547,5 @@ def getLimitingMags( ) def __repr__(self) -> str: # pragma: no cover + """Print the error model parameters.""" return "Photometric error model with parameters:\n\n" + str(self.params) diff --git a/photerr/params.py b/photerr/params.py index 5e6f9ec..d3d8386 100644 --- a/photerr/params.py +++ b/photerr/params.py @@ -1,39 +1,55 @@ """Parameter objects for the error model.""" + from __future__ import annotations from copy import deepcopy from dataclasses import InitVar, dataclass, field -from typing import Any, Dict, Optional, Union +from typing import Any import numpy as np param_docstring = """ + Note for all the dictionary parameters below, you can pass a float + instead and that value will be used for all bands. However at least + one of these parameters must be a dictionary so that the model can + infer the band names. + Parameters ---------- nYrObs : float Number of years of observations - nVisYr : dict + nVisYr : dict or float Mean number of visits per year in each band - gamma : dict + gamma : dict or float A band dependent parameter defined in Ivezic 2019 - m5 : dict + m5 : dict or float A dictionary of single visit 5-sigma limiting magnitudes. For any bands for which you pass a value in m5, this will be the 5-sigma limiting magnitude used, and any values for that band in Cm, msky, theta, and km will be ignored. - tvis : float - Exposure time in seconds for a single visit - airmass : float - The fiducial airmass - Cm : dict - A band dependent parameter defined in Ivezic 2019 - msky : dict - Median zenith sky brightness in each band, in AB mag / arcsec^2. - theta : dict - Median zenith seeing FWHM in arcseconds for each band - km : dict - Atmospheric extinction in each band + tvis : dict or float + Exposure time in seconds for a single visit in each band. + airmass : dict or float + The airmass in each band. + Cm : dict or float + System contribution to the limiting magnitude in each band that is + independent of observing conditions. + dCmInf : dict or float + The change in Cm in the limit of infinite visit time. This parameter + only matters if tvis differs from tvisRef. + msky : dict or float + Zenith sky brightness in each band, in AB mag / arcsec^2. + mskyDark : dict or float + Zenith dark sky brightness in each band, in AB mag / arcsec^2. + theta : dict or float + Median zenith seeing FWHM in arcseconds for each band. This corresponds + to seeingFwhmEff from OpSim runs. + km : dict or float + Atmospheric extinction at zenith in each band. + tvisRef : float + Reference exposure time used with tvis to scale read noise. This + must correspond to the visit time used for the calculation of Cm. sigmaSys : float; default=0.005 The irreducible error of the system in AB magnitudes. Sets the minimum photometric error. @@ -122,12 +138,14 @@ nVisYr and gamma for every band you wish to calculate errors for. In addition, you must provide either: - the single-visit 5-sigma limiting magnitude in the m5 dictionary - - tvis and airmass, plus per-band parameters in the Cm, msky, theta, and km - dictionaries, which are used to calculate the limiting magnitudes using - Eq. 6 from Ivezic 2019. + - tvis, airmass, Cm, dCmInf, msky, mskyDark, theta, and km, which are used + to calculate the limiting magnitudes using Eq. 6 of Ivezic 2019 and + Eq. 2-3 of Bianco 2022. Note if for any bands you pass a value in the m5 dictionary, the model will use - that value for that band, regardless of the values in Cm, msky, theta, and km. + that value for that band, regardless of the values in tvis, airmass, Cm, dCmInf, + msky, mskyDark, theta, and km. However you must still provide airmass and theta + if you wish to calculate errors for extended sources. When instantiating the ErrorParams object, it will determine which bands it has enough information to calculate errors for, and throw away all of the extraneous @@ -138,6 +156,7 @@ References ---------- Ivezic 2019 - https://arxiv.org/abs/0805.2366 + Bianco 2022 - https://pstn-054.lsst.io van den Busch 2020 - http://arxiv.org/abs/2007.01846 Kuijken 2019 - https://arxiv.org/abs/1902.11265 """ @@ -150,12 +169,15 @@ "nVisYr": [True, [int, float], [], False], "gamma": [True, [int, float], [], False], "m5": [True, [int, float], [], True], - "tvis": [False, [int, float, type(None)], [], False], - "airmass": [False, [int, float, type(None)], [], False], - "Cm": [True, [int, float], [], False], + "tvis": [True, [int, float], [], False], + "airmass": [True, [int, float], [], False], + "Cm": [True, [int, float], [], True], + "dCmInf": [True, [int, float], [], True], "msky": [True, [int, float], [], True], + "mskyDark": [True, [int, float], [], True], "theta": [True, [int, float], [], False], "km": [True, [int, float], [], False], + "tvisRef": [False, [int, float, type(None)], [], False], "sigmaSys": [False, [int, float], [], False], "sigLim": [False, [int, float], [], False], "ndMode": [False, [str], ["flag", "sigLim"], None], @@ -180,18 +202,21 @@ class ErrorParams: __doc__ += param_docstring nYrObs: float - nVisYr: Dict[str, float] - gamma: Dict[str, float] + nVisYr: dict[str, float] + gamma: dict[str, float] - m5: Dict[str, float] = field(default_factory=lambda: {}) + m5: dict[str, float] | float = field(default_factory=lambda: {}) - tvis: Optional[float] = None - airmass: Optional[float] = None - Cm: Dict[str, float] = field(default_factory=lambda: {}) - msky: Dict[str, float] = field(default_factory=lambda: {}) - theta: Dict[str, float] = field(default_factory=lambda: {}) - km: Dict[str, float] = field(default_factory=lambda: {}) + tvis: dict[str, float] | float = field(default_factory=lambda: {}) + airmass: dict[str, float] | float = field(default_factory=lambda: {}) + Cm: dict[str, float] | float = field(default_factory=lambda: {}) + dCmInf: dict[str, float] | float = field(default_factory=lambda: {}) + msky: dict[str, float] | float = field(default_factory=lambda: {}) + mskyDark: dict[str, float] | float = field(default_factory=lambda: {}) + theta: dict[str, float] | float = field(default_factory=lambda: {}) + km: dict[str, float] | float = field(default_factory=lambda: {}) + tvisRef: float | None = None sigmaSys: float = 0.005 sigLim: float = 0 @@ -207,20 +232,20 @@ class ErrorParams: decorrelate: bool = True highSNR: bool = False - scale: Dict[str, float] = field(default_factory=lambda: {}) + scale: dict[str, float] | float = field(default_factory=lambda: {}) errLoc: str = "after" - renameDict: InitVar[Optional[Dict[str, str]]] = None + renameDict: InitVar[dict[str, str] | None] = None validate: InitVar[bool] = True - def __post_init__( - self, renameDict: Union[Dict[str, str], None], validate: bool - ) -> None: + def __post_init__(self, renameDict: dict[str, str] | None, validate: bool) -> None: """Rename bands and remove duplicate parameters.""" - # if renameDict was provided, rename the bands - if renameDict is not None: - self.rename_bands(renameDict) + # make sure all dictionaries are dictionaries + self._convert_to_dict() + + # rename the bands + self.rename_bands(renameDict) # validate the parameters if validate: @@ -232,16 +257,43 @@ def __post_init__( # if using extended error types, if self.extendedType == "auto" or self.extendedType == "gaap": # make sure theta is provided for every band - if set(self.theta.keys()) != set(self.nVisYr.keys()): + if set(self.theta.keys()) != set(self.nVisYr.keys()) or set( + self.airmass.keys() + ) != set(self.airmass.keys()): raise ValueError( "If using one of the extended error types " "(i.e. extendedType == 'auto' or 'gaap'), " - "then theta must contain an entry for every band." + "then theta and airmass must be provided " + "for every band." ) # make sure that aMin < aMax elif self.aMin > self.aMax: raise ValueError("aMin must be less than aMax.") + def _convert_to_dict(self) -> None: + """For dict parameters that aren't dicts, convert to dicts.""" + # first loop over all the parameters and get a list of every band + bands = [] + for param in self.__dict__.values(): + if isinstance(param, dict): + bands.extend(list(param.keys())) + bands = list(set(bands)) + + if len(bands) == 0: + raise ValueError( + "At least one of the dictionary parameters " + "must actually be a dictionary." + ) + + # now loop over all the params and convert floats to dictionaries + for key, (is_dict, *_) in _val_dict.items(): + # get the parameter saved under the key + param = self.__dict__[key] + + # if it should be a dictionary but it's not + if is_dict and not isinstance(param, dict): + self.__dict__[key] = {band: param for band in bands} + def _clean_dictionaries(self) -> None: """Remove unnecessary info from all of the dictionaries. @@ -263,7 +315,7 @@ def _clean_dictionaries(self) -> None: for key, param in self.__dict__.items(): # get the parameters that are dictionaries if isinstance(param, dict) and key not in ignore_dicts: - if key != "theta": + if key != "theta" and key != "airmass": # remove bands that we have explicit m5's for self.__dict__[key] = { band: val for band, val in param.items() if band not in self.m5 @@ -285,11 +337,11 @@ def _clean_dictionaries(self) -> None: if len(all_bands) == 0: raise ValueError( "There are no bands left! You probably set one of the dictionary " - "parameters (i.e. nVisYr, gamma, m5, Cm, msky, theta, or km) such " - "that there was not enough info to calculate errors for any of the " - "photometric bands. Remember that for each band, you must have an " + "parameters (i.e. nVisYr, gamma, m5, Cm, dCmInf, msky, theta, or km) " + "such that there was not enough info to calculate errors for any of " + "the photometric bands. Remember that for each band, you must have an " "entry in nVisYr and gamma, plus either an entry in m5 or an entry " - "in all of Cm, msky, theta, and km." + "in all of Cm, dCmInf, msky, theta, and km." ) # finally, fill out the rest of the scale dictionary @@ -298,7 +350,7 @@ def _clean_dictionaries(self) -> None: for band in self.nVisYr } - def rename_bands(self, renameDict: Dict[str, str]) -> None: + def rename_bands(self, renameDict: dict[str, str]) -> None: """Rename the bands used in the error model. This method might be useful if you want to use the default parameters for an @@ -310,15 +362,20 @@ def rename_bands(self, renameDict: Dict[str, str]) -> None: A dictionary containing key: value pairs {old_name: new_name}. For example map={"u": "lsst_u"} will rename the u band to lsst_u. """ + if renameDict is None: + return + if not isinstance(renameDict, dict): + raise TypeError("renameDict must be a dict or None.") + # loop over every parameter for key, param in self.__dict__.items(): # get the parameters that are dictionaries if isinstance(param, dict): # rename bands in-place self.__dict__[key] = { - old_name - if old_name not in renameDict - else renameDict[old_name]: val + ( + old_name if old_name not in renameDict else renameDict[old_name] + ): val for old_name, val in param.items() } diff --git a/photerr/roman.py b/photerr/roman.py index c9c3c11..dedfc37 100644 --- a/photerr/roman.py +++ b/photerr/roman.py @@ -1,6 +1,7 @@ """Photometric error model for Roman.""" + from dataclasses import dataclass, field -from typing import Any, Dict +from typing import Any from photerr.model import ErrorModel from photerr.params import ErrorParams, param_docstring @@ -18,24 +19,10 @@ class RomanErrorParams(ErrorParams): __doc__ += " Graham 2020 - https://arxiv.org/abs/2004.07885" nYrObs: float = 1.0 - nVisYr: Dict[str, float] = field( - default_factory=lambda: { - "Y": 1.0, - "J": 1.0, - "H": 1.0, - "F": 1.0, - } - ) - gamma: Dict[str, float] = field( - default_factory=lambda: { - "Y": 0.04, - "J": 0.04, - "H": 0.04, - "F": 0.04, - } - ) + nVisYr: dict[str, float] | float = 1.0 + gamma: dict[str, float] | float = 0.04 - m5: Dict[str, float] = field( + m5: dict[str, float] | float = field( default_factory=lambda: { "Y": 26.9, "J": 26.95, diff --git a/pyproject.toml b/pyproject.toml index 4b9949d..0334f5a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,38 +1,38 @@ [tool.poetry] name = "photerr" -version = "1.1.1" +version = "1.2.0" description = "Photometric error model for astronomical imaging surveys" authors = ["John Franklin Crenshaw "] readme = "README.md" homepage = "https://github.com/jfcrenshaw/photerr" [tool.poetry.dependencies] -python = ">=3.8,<4.0" -numpy = "^1.23.2" -pandas = "^1.4.3" +python = ">=3.10,<4.0" +numpy = ">=1.23.2" +pandas = ">=1.4.3" [tool.poetry.dev-dependencies] -pytest = "^7.1" -pytest-cov = "^3.0.0" -black = "^22.6.0" -mypy = "^0.971" -isort = "^5.10.1" -flake8 = "^5.0.4" -flake8-black = "^0.3.3" -flake8-broken-line = "^0.5.0" -flake8-bugbear = "^22.8.23" -flake8-builtins = "^1.5.3" -flake8-comprehensions = "^3.10.0" -flake8-docstrings = "^1.6.0" -flake8-eradicate = "^1.3.0" -flake8-isort = "^4.2.0" -flake8-print = "^5.0.0" -flake8-pytest-style = "^1.6.0" -flake8-simplify = "^0.19.3" -flake8-tidy-imports = "^4.8.0" -jupyterlab = "^3.4.5" -matplotlib = "^3.5.3" -pre-commit = "^2.20.0" +pytest = ">=7.1" +pytest-cov = ">=3.0.0" +black = ">=22.6.0" +mypy = ">=0.971" +isort = ">=5.10.1" +flake8 = ">=5.0.4" +flake8-black = ">=0.3.3" +flake8-broken-line = ">=0.5.0" +flake8-bugbear = ">=22.8.23" +flake8-builtins = ">=1.5.3" +flake8-comprehensions = ">=3.10.0" +flake8-docstrings = ">=1.6.0" +flake8-eradicate = ">=1.3.0" +flake8-isort = ">=4.2.0" +flake8-print = ">=5.0.0" +flake8-pytest-style = ">=1.6.0" +flake8-simplify = ">=0.19.3" +flake8-tidy-imports = ">=4.8.0" +jupyterlab = ">=3.4.5" +matplotlib = ">=3.5.3" +pre-commit = ">=2.20.0" [build-system] requires = ["poetry-core>=1.0.0"] diff --git a/tests/test_models.py b/tests/test_models.py index d86454e..04933a7 100644 --- a/tests/test_models.py +++ b/tests/test_models.py @@ -10,6 +10,7 @@ ErrorModel, EuclidErrorModel, LsstErrorModel, + LsstErrorModelV1, LsstErrorParams, RomanErrorModel, ) @@ -341,6 +342,10 @@ def test_bad_instantiation() -> None: def test_other_models(data: pd.DataFrame) -> None: """Test instantiating other models and calculating errors.""" + + lsstData = LsstErrorModelV1()(data, 0) + assert lsstData.shape == (data.shape[0], data.shape[1] + 1) + euclidData = EuclidErrorModel()(data, 0) assert euclidData.shape == (data.shape[0], data.shape[1] + 1) @@ -374,7 +379,9 @@ def test_alternate_instantiation() -> None: @pytest.mark.parametrize("extendedType", ["point", "auto", "gaap"]) @pytest.mark.parametrize("highSNR", [True, False]) def test_mag_nsr_inversion( - extendedType: str, highSNR: bool, lsst_data: pd.DataFrame + extendedType: str, + highSNR: bool, + lsst_data: pd.DataFrame, ) -> None: """Test that the mag->nsr and nsr->mag methods are inverses.""" # get the arrays of data @@ -462,3 +469,15 @@ def test_scale(highSNR: bool) -> None: # but u band of low SNR is doubled np.allclose(ratio[1, 0], 2) + + +def test_limiting_mags() -> None: + """Compare V1 limiting mags to the values in Table 2 of Ivezic 2019.""" + # get the limiting mags from the error model + errM = LsstErrorModelV1(airmass=1) + m5 = errM.getLimitingMags(coadded=False) + + # compare to the Ivezic 2019 values + ivezic2019 = dict(u=23.78, g=24.81, r=24.35, i=23.92, z=23.34, y=22.45) + for band in m5: + assert np.isclose(m5[band], ivezic2019[band], rtol=1e-3) diff --git a/tests/test_params.py b/tests/test_params.py index d5f99b5..13d9032 100644 --- a/tests/test_params.py +++ b/tests/test_params.py @@ -134,3 +134,21 @@ def test_missing_theta() -> None: """Test fail if we have extended error but don't have theta for everyone.""" with pytest.raises(ValueError): LsstErrorParams(extendedType="auto", theta={"g": 0.1}, m5={"u": 23}) + + +def test_all_dicts_are_floats() -> None: + """Test that instantiation fails if all dictionaries are floats.""" + with pytest.raises(ValueError): + LsstErrorParams( + nVisYr=1, + gamma=1, + m5=1, + tvis=1, + airmass=1, + Cm=1, + dCmInf=1, + msky=1, + mskyDark=1, + theta=1, + km=1, + )