From 051bc6c5cd71da31eda557945d1847228ed38f3d Mon Sep 17 00:00:00 2001 From: John Franklin Crenshaw <41785729+jfcrenshaw@users.noreply.github.com> Date: Tue, 25 Apr 2023 22:35:41 -0700 Subject: [PATCH] v1.1.0 (#9) * Fixed most sigLim edge cases. * More edge cases; limiting mags with finite aperture. * Implemented error scaling. * Scale now only impacts band-specific error. Updated README. * Bumped minor version number. --- README.md | 25 +++++ photerr/model.py | 219 +++++++++++++++++++++++------------- photerr/params.py | 23 +++- pyproject.toml | 2 +- tests/test_models.py | 259 ++++++++++++++++++++++++++++++++++++++++--- tests/test_params.py | 4 +- 6 files changed, 431 insertions(+), 101 deletions(-) diff --git a/README.md b/README.md index a037c77..57570cf 100644 --- a/README.md +++ b/README.md @@ -36,6 +36,13 @@ Any extraneous columns in the DataFrame (e.g. a redshift column), remain in the *If compatibility with Astropy Tables, Ordered Dictionaries, etc., would be useful to you, let me know!* +You can also calculate limiting magnitudes: + +```python +errModel.getLimitingMags() # coadded point-source 5-sigma limits +errModel.getLimitingMags(nSigma=1, coadded=False) # single-image point-source 1-sigma limits +``` + # Tweaking the error model There are many parameters you can tweak to fine tune the error model. @@ -102,6 +109,24 @@ If `absFlux=True`, then the absolute value of all fluxes are taken before conver If combined with `sigLim=0`, this means every galaxy will have an observed flux in every band. This is useful if you do not want to worry about non-detections, but it results in a non-Gaussian error distribution for low-SNR sources. +### Errors for extended sources + +PhotErr can be used to calculate errors for extended sources as well. +You just have to pass `extendedType="auto"` or `extendedType="gaap"` to the constructor (see explanation below for the differences in these models). +PhotErr will then look for columns in the input DataFrame that correspond to the semi-major and -minor axes of the objects, corresponding to half-light radii in arcseconds. +By default, it looks for these in columns titled "major" and "minor", but you can change the names of these columns using the `majorCol` and `minorCol` keywords. + +You can also calculate limiting magnitudes for apertures of a given size by passing the `aperture` keyword to `errModel.getLimitingMags()` + +### Scaling the errors + +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. + ### *Other error models* In addition to `LsstErrorModel`, which comes with the LSST defaults, PhotErr includes `EuclidErrorModel` and `RomanErrorModel`, which come with the Euclid and Roman defaults, respectively. diff --git a/photerr/model.py b/photerr/model.py index f0a6349..97be8b9 100644 --- a/photerr/model.py +++ b/photerr/model.py @@ -56,9 +56,6 @@ def __init__(self, *args: ErrorParams, **kwargs: Any) -> None: # calculate all of the 5-sigma limiting magnitudes self._calculate_m5() - # calculate the limits for sigLim - self._sigLim = self.getLimitingMags(self._params.sigLim, coadded=True) - @property def params(self) -> ErrorParams: """The error model parameters in an ErrorParams objet.""" # noqa: D401 @@ -151,8 +148,8 @@ def _get_area_ratio_gaap( # convert galaxy half-light radii to Gaussian sigmas # this conversion factor comes from half-IQR -> sigma - majors /= 0.6745 - minors /= 0.6745 + majors = majors / 0.6745 + minors = minors / 0.6745 # calculate the area of the psf in each band A_psf = np.pi * psf_sig**2 @@ -167,7 +164,7 @@ def _get_area_ratio_gaap( # return their ratio return A_ap / A_psf - def _get_NSR( + def _get_nsr_from_mags( self, mags: np.ndarray, majors: np.ndarray, minors: np.ndarray, bands: list ) -> np.ndarray: """Calculate the noise-to-signal ratio. @@ -189,7 +186,7 @@ def _get_NSR( Returns ------- np.ndarray - The ratio of aperture size to PSF size for each band and galaxy. + The noise-to-signal ratio of each galaxy """ # get the 5-sigma limiting magnitudes for these bands m5 = np.array([self._all_m5[band] for band in bands]) @@ -197,6 +194,8 @@ def _get_NSR( gamma = np.array([self.params.gamma[band] for band in bands]) # and the number of visits per year nVisYr = np.array([self.params.nVisYr[band] for band in bands]) + # and the scales + scale = np.array([self.params.scale[band] for band in bands]) # calculate x as defined in the paper x = 10 ** (0.4 * (mags - m5)) @@ -218,6 +217,9 @@ def _get_NSR( A_ratio = 1 # type: ignore nsrRand *= np.sqrt(A_ratio) + # rescale the NSR + nsrRand *= scale + # get the irreducible system NSR if self.params.highSNR: nsrSys = self.params.sigmaSys @@ -229,13 +231,94 @@ def _get_NSR( return nsr + def _get_mags_from_nsr( + self, + nsr: np.ndarray, + majors: np.ndarray, + minors: np.ndarray, + bands: list, + coadded: bool = True, + ) -> np.ndarray: + """Calculate magnitudes that correspond to the given NSRs. + + Essentially inverts self._get_nsr_from_mags(). + + Parameters + ---------- + nsr : np.ndarray + The noise-to-signal ratios of the galaxies + majors : np.ndarray + The semi-major axes of the galaxies in arcseconds + minors : np.ndarray + The semi-minor axes of the galaxies in arcseconds + bands : list + The list of bands the galaxy is observed in + coadded : bool; default=True + If True, assumes NSR is after coaddition. + + Returns + ------- + np.ndarray + The magnitude corresponding to the NSR for each galaxy + """ + # get the 5-sigma limiting magnitudes for these bands + m5 = np.array([self._all_m5[band] for band in bands]) + # and the values for gamma + gamma = np.array([self.params.gamma[band] for band in bands]) + # and the number of visits per year + nVisYr = np.array([self.params.nVisYr[band] for band in bands]) + # and the scales + scale = np.array([self.params.scale[band] for band in bands]) + + # get the irreducible system NSR + if self.params.highSNR: + nsrSys = self.params.sigmaSys + else: + nsrSys = 10 ** (self.params.sigmaSys / 2.5) - 1 + + # calculate the random NSR + nsrRand = np.sqrt(nsr**2 - nsrSys**2) + + # rescale the NSR + nsrRand /= scale + + # rescale according to the area ratio + if majors is not None and minors is not None: + if self.params.extendedType == "auto": + A_ratio = self._get_area_ratio_auto(majors, minors, bands) + elif self.params.extendedType == "gaap": + A_ratio = self._get_area_ratio_gaap(majors, minors, bands) + else: + A_ratio = 1 # type: ignore + nsrRand = nsrRand / np.sqrt(A_ratio) + + # get the number of exposures + if coadded: + nStackedObs = nVisYr * self.params.nYrObs + else: + nStackedObs = 1 # type: ignore + + # calculate the NSR for a single image + nsrRandSingleExp = nsrRand * np.sqrt(nStackedObs) + + # calculate the value of x that corresponds to this NSR + # this is just the quadratic equation applied to Eq 5 from Ivezic 2019 + x = ( + (gamma - 0.04) + + np.sqrt((gamma - 0.04) ** 2 + 4 * gamma * nsrRandSingleExp**2) + ) / (2 * gamma) + + # convert x to magnitudes + mags = m5 + 2.5 * np.log10(x) + + return mags + def _get_obs_and_errs( self, mags: np.ndarray, majors: np.ndarray, minors: np.ndarray, bands: list, - sigLim: np.ndarray, rng: np.random.Generator, ) -> Tuple[np.ndarray, np.ndarray]: """Calculate the noise-to-signal ratio. @@ -253,8 +336,6 @@ def _get_obs_and_errs( The semi-minor axes of the galaxies in arcseconds bands : list The list of bands the galaxy is observed in - sigLim : np.ndarray - The n-sigma limits for non-detection rng : np.random.Generator A numpy random number generator @@ -264,7 +345,7 @@ def _get_obs_and_errs( The ratio of aperture size to PSF size for each band and galaxy. """ # get the NSR for all galaxies - nsr = self._get_NSR(mags, majors, minors, bands) + nsr = self._get_nsr_from_mags(mags, majors, minors, bands) if self.params.highSNR: # in the high SNR approximation, mag err ~ nsr, and we can model @@ -273,16 +354,6 @@ def _get_obs_and_errs( # calculate observed magnitudes obsMags = rng.normal(loc=mags, scale=nsr) - # if ndMode == sigLim, then clip all magnitudes at the n-sigma limit - if self.params.ndMode == "sigLim": - obsMags = np.clip(obsMags, None, sigLim) - - # if decorrelate, then we calculate new errors using the observed mags - if self.params.decorrelate: - nsr = self._get_NSR(obsMags, majors, minors, bands) - - obsMagErrs = nsr - else: # in the more accurate error model, we acknowledge err != nsr, # and we model errors as Gaussian in flux space @@ -295,14 +366,26 @@ def _get_obs_and_errs( with np.errstate(divide="ignore"): obsMags = -2.5 * np.log10(np.clip(obsFluxes, 0, None)) - # if ndMode == sigLim, then clip all magnitudes at the n-sigma limit - if self.params.ndMode == "sigLim": - obsMags = np.clip(obsMags, None, sigLim) + # if decorrelate, then we calculate new errors using the observed mags + if self.params.decorrelate: + nsr = self._get_nsr_from_mags(obsMags, majors, minors, bands) + + # if ndMode == sigLim, then clip at the n-sigma limit + if self.params.ndMode == "sigLim": + # determine the nsr limit + with np.errstate(divide="ignore"): + nsrLim = np.divide(1, self.params.sigLim) + + # calculate limiting magnitudes for each galaxy + magLim = self._get_mags_from_nsr(nsrLim, majors, minors, bands) - # if decorrelate, then we calculate new errors using the observed mags - if self.params.decorrelate: - nsr = self._get_NSR(obsMags, majors, minors, bands) + # clip mags and nsr's at this limit + nsr = np.clip(nsr, 0, nsrLim) + obsMags = np.clip(obsMags, None, magLim) + if self.params.highSNR: + obsMagErrs = nsr + else: obsMagErrs = 2.5 * np.log10(1 + nsr) return obsMags, obsMagErrs @@ -336,9 +419,6 @@ def __call__( # get the bands we will calculate errors for bands = [band for band in catalog.columns if band in self._bands] - # calculate the n-sigma limits - sigLim = np.array([self._sigLim[band] for band in bands]) - # get the numpy array of magnitudes mags = catalog[bands].to_numpy() @@ -351,13 +431,18 @@ def __call__( minors = None # get observed magnitudes and errors - obsMags, obsMagErrs = self._get_obs_and_errs( - mags, majors, minors, bands, sigLim, rng - ) + obsMags, obsMagErrs = self._get_obs_and_errs(mags, majors, minors, bands, rng) # flag all non-detections with the ndFlag if self.params.ndMode == "flag": - idx = (~np.isfinite(obsMags)) | (obsMags > sigLim) + # calculate SNR + if self.params.highSNR: + snr = 1 / obsMagErrs + else: + snr = 1 / (10 ** (obsMagErrs / 2.5) - 1) + + # flag non-finite mags and where SNR is below sigLim + idx = (~np.isfinite(obsMags)) | (snr < self.params.sigLim) obsMags[idx] = self.params.ndFlag obsMagErrs[idx] = self.params.ndFlag @@ -382,12 +467,10 @@ def __call__( return obsCatalog - def getLimitingMags(self, nSigma: float = 5, coadded: bool = True) -> dict: - """Return the limiting magnitudes for point sources in all bands. - - This method essentially reverse engineers the _get_NSR method so that we - calculate what magnitude results in NSR = 1/nSigma. - (NSR is noise-to-signal ratio; NSR = 1/SNR) + def getLimitingMags( + self, nSigma: float = 5, coadded: bool = True, aperture: float = 0 + ) -> dict: + """Return the limiting magnitudes in all bands. Parameters ---------- @@ -398,52 +481,30 @@ def getLimitingMags(self, nSigma: float = 5, coadded: bool = True) -> dict: coadded : bool; default=True If True, returns the limiting magnitudes for a coadded image. If False, returns the limiting magnitudes for a single visit. + aperture : float, default=0 + Radius of circular aperture in arcseconds. If zero, limiting magnitudes + are for point sources. If greater than zero, limiting magnitudes are + for objects of that size. Increasing the aperture decreases the + signal-to-noise ratio. This only has an impact if extendedType != "point. Returns ------- dict The dictionary of limiting magnitudes """ - bands = self._bands - - # if nSigma is zero, return infinite magnitude limits - if np.isclose(0, nSigma): - return {band: np.inf for band in bands} - - # get the 5-sigma limiting magnitudes for these bands - m5 = np.array([self._all_m5[band] for band in bands]) - # and the values for gamma - gamma = np.array([self.params.gamma[band] for band in bands]) - # and the number of visits per year - nVisYr = np.array([self.params.nVisYr[band] for band in bands]) - - # get the number of exposures - if coadded: - nStackedObs = nVisYr * self.params.nYrObs - else: - nStackedObs = 1 # type: ignore - - # get the irreducible system error - if self.params.highSNR: - nsrSys = self.params.sigmaSys - else: - nsrSys = 10 ** (self.params.sigmaSys / 2.5) - 1 - - # calculate the random NSR that a single exposure must have - nsrRandSingleExp = np.sqrt((1 / nSigma**2 - nsrSys**2) * nStackedObs) - - # calculate the value of x that corresponds to this NSR - # this is just the quadratic equation applied to Eq 5 from Ivezic 2019 - x = ( - (gamma - 0.04) - + np.sqrt((gamma - 0.04) ** 2 + 4 * gamma * nsrRandSingleExp**2) - ) / (2 * gamma) - - # convert x to a limiting magnitude - limiting_mags = m5 + 2.5 * np.log10(x) - # return as a dictionary - return dict(zip(bands, limiting_mags)) + return dict( + zip( + self._bands, + self._get_mags_from_nsr( + 1 / nSigma, # type: ignore + np.array([aperture]), + np.array([aperture]), + self._bands, + coadded, + ).flatten(), + ) + ) def __repr__(self) -> str: # pragma: no cover return "Photometric error model with parameters:\n\n" + str(self.params) diff --git a/photerr/params.py b/photerr/params.py index 4ef930c..279d385 100644 --- a/photerr/params.py +++ b/photerr/params.py @@ -83,6 +83,11 @@ then Eq. 5 from that paper is used to calculate (N/S)^2 in flux, and errors are Gaussian in flux space. If True, Eq. 5 is used to calculate the Gaussian variance in magnitude space. + scale : dict; default=dict() + A dictionary that rescales the error for the given bands. For example, if + scale = {"u": 2}, then all the errors for the u band are doubled. This allows + you to answer questions like "what happens to my science if the u band errors + are doubled." errLoc: str; default="after" Where to place the error columns in the output table. If "after", then the error columns will be placed immediately after the corresponding magnitude @@ -163,6 +168,7 @@ "minorCol": [False, [str], [], None], "decorrelate": [False, [bool], [], None], "highSNR": [False, [bool], [], None], + "scale": [True, [int, float], [], None], "errLoc": [False, [str], ["after", "end", "alone"], None], } @@ -201,10 +207,11 @@ class ErrorParams: decorrelate: bool = True highSNR: bool = False + scale: Dict[str, float] = field(default_factory=lambda: {}) + errLoc: str = "after" renameDict: InitVar[Dict[str, str]] = None - validate: InitVar[bool] = True def __post_init__( @@ -243,13 +250,16 @@ def _clean_dictionaries(self) -> None: Additionally, we will remove any band from all dictionaries that don't have an entry in all of nVisYr, gamma, and (m5 OR Cm + msky + theta + km). + + Finally, we will set scale=1 for all bands for which the scale isn't + explicitly set. """ # keep a running set of all the bands we will calculate errors for all_bands = set(self.nVisYr.keys()).intersection(self.gamma.keys()) # remove the m5 bands from all other parameter dictionaries, and remove # bands from all_bands for which we don't have m5 data for - ignore_dicts = ["m5", "nVisYr", "gamma"] + ignore_dicts = ["m5", "nVisYr", "gamma", "scale"] for key, param in self.__dict__.items(): # get the parameters that are dictionaries if isinstance(param, dict) and key not in ignore_dicts: @@ -264,8 +274,7 @@ def _clean_dictionaries(self) -> None: set(param.keys()).union(set(self.m5.keys())) ) - # finally, remove all of the data for bands that we will not be - # calculating errors for + # remove all data for bands we will not be calculating errors for for key, param in self.__dict__.items(): if isinstance(param, dict): self.__dict__[key] = { @@ -283,6 +292,12 @@ def _clean_dictionaries(self) -> None: "in all of Cm, msky, theta, and km." ) + # finally, fill out the rest of the scale dictionary + self.scale = { + band: float(self.scale[band]) if band in self.scale else 1.0 + for band in self.nVisYr + } + def rename_bands(self, renameDict: Dict[str, str]) -> None: """Rename the bands used in the error model. diff --git a/pyproject.toml b/pyproject.toml index 05dbd4c..8413834 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "photerr" -version = "1.0.1" +version = "1.1.0" description = "Photometric error model for astronomical imaging surveys" authors = ["John Franklin Crenshaw "] readme = "README.md" diff --git a/tests/test_models.py b/tests/test_models.py index 124ad11..d86454e 100644 --- a/tests/test_models.py +++ b/tests/test_models.py @@ -10,8 +10,8 @@ ErrorModel, EuclidErrorModel, LsstErrorModel, - RomanErrorModel, LsstErrorParams, + RomanErrorModel, ) @@ -33,6 +33,20 @@ def data() -> pd.DataFrame: return dataframe +@pytest.fixture() +def lsst_data() -> pd.DataFrame: + """Return 1000 random LSST galaxies for error model tests.""" + rng = np.random.default_rng(0) + array = rng.normal( + loc=[25, 25, 25, 25, 25, 25, 0, 0], + scale=[2, 2, 2, 2, 2, 2, 0.2, 0.2], + size=(10_000, 8), + ) + array[:, -2:] = np.abs(array[:, -2:]) + dataframe = pd.DataFrame(array, columns=list("ugrizy") + ["major", "minor"]) + return dataframe + + def test_random_state(data: pd.DataFrame) -> None: """Test that the random state behaves as expected.""" data = data[["g"]].iloc[:2] @@ -133,26 +147,125 @@ def test_ndFlag(data: pd.DataFrame) -> None: assert LsstErrorModel(ndFlag=-999, sigLim=1)(data, 0).iloc[-1, 0] == -999 -def test_sigLim(data: pd.DataFrame) -> None: - """Test that sigLim works correctly.""" - # test that everything beyond sigLim is flagged +def test_pointSource_sigLim(data: pd.DataFrame) -> None: + """Test that everything beyond sigLim is flagged.""" assert np.all(LsstErrorModel(sigLim=10)(data, 0)[["g", "g_err"]].iloc[1:] == np.inf) - # test that everything beyond sigLim is cut to sigLim when ndMode==sigLim - # first with highSNR=False - sigLimData = LsstErrorModel(sigLim=10, ndMode="sigLim", highSNR=False)(data, 0) - assert np.isclose(sigLimData["g"][1], sigLimData["g"][2]) - assert np.isclose(sigLimData["g_err"][1], sigLimData["g_err"][2]) - # now with highSNR=True - sigLimData = LsstErrorModel(sigLim=10, ndMode="sigLim", highSNR=True)(data, 0) + +@pytest.mark.parametrize("extendedType", ["auto", "gaap"]) +@pytest.mark.parametrize("ndMode", ["flag", "sigLim"]) +@pytest.mark.parametrize("highSNR", [True, False]) +@pytest.mark.parametrize("decorrelate", [True, False]) +def test_extended_sigLim( + extendedType: str, + ndMode: str, + highSNR: bool, + decorrelate: bool, + lsst_data: pd.DataFrame, +) -> None: + """Test that sigLim works with extended source errors.""" + # generate the photometric errors + sigLim = 3.0 + errM = LsstErrorModel( + sigLim=sigLim, + ndMode=ndMode, + extendedType=extendedType, + highSNR=highSNR, + decorrelate=decorrelate, + ) + data = errM(lsst_data, 0) + + # keep only the galaxies with finite magnitudes + data = data[np.isfinite(data).all(axis=1)] + + # calculate the SNR of every galaxy + if highSNR: + snr = 1 / data[[col for col in data.columns if "err" in col]] + else: + snr = 1 / ( + 10 ** (data[[col for col in data.columns if "err" in col]] / 2.5) - 1 + ) + snr.rename(columns=lambda name: name.replace("err", "snr"), inplace=True) + + # make sure the minimum SNR exceeds sigLim + assert np.all((snr.min() >= sigLim) | np.isclose(snr.min(), sigLim)) + + +@pytest.mark.parametrize("highSNR", [True, False]) +@pytest.mark.parametrize("decorrelate", [True, False]) +def test_pointSource_ndFlag_sigLim( + highSNR: bool, decorrelate: bool, data: pd.DataFrame +) -> None: + """Test that ndFlag=sigLim works for point sources. + + i.e. everything beyond sigLim is cut to sigLim when ndMode==sigLim + """ + sigLimData = LsstErrorModel( + sigLim=10, + ndMode="sigLim", + highSNR=highSNR, + decorrelate=decorrelate, + )(data, 0) assert np.isclose(sigLimData["g"][1], sigLimData["g"][2]) assert np.isclose(sigLimData["g_err"][1], sigLimData["g_err"][2]) -def test_absFlux(data: pd.DataFrame) -> None: - """Test that absFlux results in all finite magnitudes.""" - assert ~np.all(np.isfinite(LsstErrorModel(absFlux=False)(data, 5))) - assert np.all(np.isfinite(LsstErrorModel(absFlux=True)(data, 5))) +@pytest.mark.parametrize("highSNR", [False, True]) +def test_decorrelate_doesnt_change_siglim(highSNR: bool, data: pd.DataFrame) -> None: + """Test that decorrelate doesn't change the sigLim-ed outputs.""" + decorrData = LsstErrorModel( + sigLim=100, ndMode="sigLim", highSNR=highSNR, decorrelate=True + )(data, 0) + + corrData = LsstErrorModel( + sigLim=100, ndMode="sigLim", highSNR=highSNR, decorrelate=False + )(data, 0) + + assert np.allclose(decorrData[1:], corrData[1:]) + + +@pytest.mark.parametrize("decorrelate", [False, True]) +@pytest.mark.parametrize( + "sigLim,ndMode,absFlux", + [ + (5, "sigLim", False), + (0, "flag", True), + ], +) +def test_finite_values( + decorrelate: bool, + sigLim: float, + ndMode: str, + absFlux: bool, + lsst_data: pd.DataFrame, +) -> None: + """Test settings that should result in all finite values. + + These two settings should result in all finite values: + - absFlux=True with sigLim=0 + - ndMode="sigLim" with sigLim>0 + """ + # first let's make sure that our default has some non-finite elements + assert ~np.all( + np.isfinite( + LsstErrorModel( + decorrelate=decorrelate, + sigLim=sigLim, + )(lsst_data, 0) + ) + ) + + # now check that enabling these settings makes everything finite + assert np.all( + np.isfinite( + LsstErrorModel( + decorrelate=decorrelate, + sigLim=sigLim, + ndMode=ndMode, + absFlux=absFlux, + )(lsst_data, 0) + ) + ) def test_limitingMags() -> None: @@ -178,6 +291,29 @@ def test_limitingMags() -> None: assert all(np.isclose(errModel.params.m5[band], magLim[band]) for band in magLim) +@pytest.mark.parametrize("extendedType", ["auto", "gaap"]) +def test_limitMags_aperture(extendedType: str) -> None: + """Test that increasing the aperture decreasing depth.""" + errModel = LsstErrorModel(extendedType=extendedType) + + magLimAp0 = errModel.getLimitingMags(aperture=0) + magLimAp1 = errModel.getLimitingMags(aperture=1) + assert all(magLimAp0[band] > magLimAp1[band] for band in magLimAp0) + + +def test_limitingMags_scale() -> None: + """Test that increasing the error scale decreases the limiting mags.""" + magLimScale1 = LsstErrorModel( + scale={band: 1 for band in "ugrizy"} + ).getLimitingMags() + + magLimScale2 = LsstErrorModel( + scale={band: 2 for band in "ugrizy"} + ).getLimitingMags() + + assert all(magLimScale1[band] > magLimScale2[band] for band in magLimScale1) + + def test_errLoc(data: pd.DataFrame) -> None: """Test that errLoc works as expected.""" # the error column should come right after the magnitude column @@ -233,3 +369,96 @@ def test_alternate_instantiation() -> None: errM1 = ErrorModel(LsstErrorParams(), nYrObs=2) errM2 = ErrorModel(LsstErrorParams(nYrObs=2)) assert errM1.params == errM2.params + + +@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 +) -> None: + """Test that the mag->nsr and nsr->mag methods are inverses.""" + # get the arrays of data + bands = list("ugrizy") + mags = lsst_data[bands].to_numpy() + majors = lsst_data["major"].to_numpy() + minors = lsst_data["minor"].to_numpy() + + # create the error model + errM = LsstErrorModel(extendedType=extendedType, highSNR=highSNR) + + # compute NSR + nsr = errM._get_nsr_from_mags(mags, majors, minors, bands) + + # use these NSRs to re-compute input mags + mags_recomputed = errM._get_mags_from_nsr(nsr, majors, minors, bands) + + # make sure we got back what we put in + assert np.allclose(mags, mags_recomputed) + + +@pytest.mark.parametrize("sigLim", [0, 5]) +def test_both_not_finite(sigLim: float, lsst_data: pd.DataFrame) -> None: + """Test that the band and error are both finite at the same time.""" + data = LsstErrorModel(sigLim=sigLim)(lsst_data, 0) + + # loop through every band + for band in "ugrizy": + # where band is not finite, make sure error is not finite + assert ~np.isfinite(data[f"{band}_err"][~np.isfinite(data[band])]).any() + + # where error is not finite, make sure band is not finite + assert ~np.isfinite(data[band][~np.isfinite(data[f"{band}_err"])]).any() + + +@pytest.mark.parametrize("highSNR", [True, False]) +def test_scale(highSNR: bool) -> None: + """Test that scale linearly scales the final errors and changes limiting mags. + + Note we have to set decorrelate=False to make sure the errors are re-sampled, + and absFlux=True to make sure all the errors are finite. + """ + # some custom data for this test + # one super high SNR, and one very low SNR + data = pd.DataFrame( + np.array([[10], [30]]) * np.ones((2, 6)), + columns=list("ugrizy"), + ) + + # calculate errors with u scale = 1 and 2 + errsScale1 = LsstErrorModel( + errLoc="alone", + decorrelate=False, + absFlux=True, + highSNR=highSNR, + )(data, 0) + errsScale2 = LsstErrorModel( + errLoc="alone", + decorrelate=False, + absFlux=True, + highSNR=highSNR, + scale={"u": 2}, + )(data, 0) + + # convert to numpy array + errsScale1 = errsScale1.to_numpy() + errsScale2 = errsScale2.to_numpy() + + # calculate nsr + if highSNR: + nsrScale1 = errsScale1 + nsrScale2 = errsScale2 + else: + nsrScale1 = 10 ** (errsScale1 / 2.5) - 1 + nsrScale2 = 10 ** (errsScale2 / 2.5) - 1 + + # calculate the error ratio + ratio = nsrScale2 / nsrScale1 + + # check that super high SNR is not impacted + assert np.allclose(ratio[0, :], 1) + + # and the grizy bands of the low SNR + np.allclose(ratio[1, 1:], 1) + + # but u band of low SNR is doubled + np.allclose(ratio[1, 0], 2) diff --git a/tests/test_params.py b/tests/test_params.py index f3da75b..d5f99b5 100644 --- a/tests/test_params.py +++ b/tests/test_params.py @@ -115,9 +115,9 @@ def test_param_val_dict() -> None: ) def test_bad_params(params: dict, error: Exception) -> None: """Test that instantiation and updating fails with bad parameters.""" - with pytest.raises(error): + with pytest.raises(error): # type: ignore LsstErrorParams(**params) - with pytest.raises(error): + with pytest.raises(error): # type: ignore LsstErrorParams().update(**params)