Merge branch 'main' into galaxyq

Author: itrharrison

.zenodo.json

@@ -43,6 +43,11 @@
             "affiliation": "University of Manchester",
             "name": "Juan Pablo Cordero"
         },
+        {
+            "orcid": "0000-0003-1560-7959",
+            "affiliation": "University of Portsmouth",
+            "name": "Fox Davidson"
+        },
         {
             "orcid": "0000-0002-8218-563X",
             "affiliation": "University of Portsmouth",

docs/utils/index.rst

@@ -143,6 +143,8 @@ create their own named `speclite.filters.FilterResponse`.
    luminosity_in_band
    mag_ab
    SpectrumTemplates
+   magnitude_error_rykoff
+   logistic_completeness_function
 
 
 Random sampling (`skypy.utils.random`)

skypy/galaxies/__init__.py

@@ -5,6 +5,7 @@
 
 __all__ = [
     'schechter_lf',
+    'schechter_smf',
 ]
 
 from . import luminosity  # noqa F401,F403

skypy/galaxies/spectrum.py

@@ -143,6 +143,28 @@ def stellar_mass(self, coefficients, magnitudes, filter):
         Mt = self.absolute_magnitudes(coefficients, filter)
         return np.power(10, 0.4*(Mt-magnitudes))
 
+    def stellar_mass_remain(self, coefficients, magnitudes, filter):
+        r'''Compute total surviving stellar mass from absolute magnitudes.
+
+        Parameters
+        ----------
+        coefficients : (ng, nt) array_like
+            Array of template coefficients.
+        magnitudes : (ng,) array_like
+            The magnitudes to match in the reference bandpass.
+        filter : str
+            A single reference bandpass filter specification for
+            `~speclite.filters.load_filters`.
+
+        Returns
+        -------
+        stellar_mass : (ng,) array_like
+            Total surviving stellar mass of each galaxy in template units.
+        '''
+        m = self.stellar_mass(coefficients, magnitudes, filter)
+        m *= np.dot(coefficients, self.mremain)
+        return m
+
     def metallicity(self, coefficients):
         r'''Galaxy metallicities from kcorrect templates.
 

skypy/galaxies/stellar_mass.py

@@ -4,6 +4,7 @@
 import numpy as np
 
 from ..utils.random import schechter
+from ..utils import dependent_argument
 
 
 __all__ = [
@@ -15,6 +16,8 @@
 ]
 
 
+@dependent_argument('m_star', 'redshift')
+@dependent_argument('alpha', 'redshift')
 def schechter_smf_mass(redshift, alpha, m_star, m_min, m_max, size=None,
                        resolution=1000):
     r""" Stellar masses following the Schechter mass function [1]_.
@@ -23,9 +26,10 @@ def schechter_smf_mass(redshift, alpha, m_star, m_min, m_max, size=None,
     ----------
     redshift : array_like
         Galaxy redshifts for which to sample magnitudes.
-    alpha : float
-        The alpha parameter in the Schechter stellar mass function.
-    m_star : (nm,) array-like
+    alpha : float or function
+        The alpha parameter in the Schechter stellar mass function. If function,
+        it must return a scalar value.
+    m_star : (nm,) array-like or function
         Characteristic stellar mass m_*.
     size: int, optional
          Output shape of stellar mass samples. If size is None and m_star

skypy/galaxies/tests/test_spectrum.py

@@ -158,6 +158,39 @@ def test_kcorrect_stellar_mass():
     np.testing.assert_allclose(stellar_mass, truth)
 
 
+@pytest.mark.skipif(not HAS_SPECLITE, reason='test requires speclite')
+def test_kcorrect_stellar_mass_remain():
+
+    from astropy import units
+    from skypy.galaxies.spectrum import kcorrect
+    from speclite.filters import FilterResponse
+
+    # Gaussian bandpass
+    filt_lam = np.logspace(3, 4, 1000) * units.AA
+    filt_mean = 5000 * units.AA
+    filt_width = 100 * units.AA
+    filt_tx = np.exp(-((filt_lam-filt_mean)/filt_width)**2)
+    filt_tx[[0, -1]] = 0
+    FilterResponse(wavelength=filt_lam, response=filt_tx,
+                   meta=dict(group_name='test', band_name='filt'))
+
+    # Using the identity matrix for the coefficients yields trivial test cases
+    coeff = np.eye(5)
+    Mt = kcorrect.absolute_magnitudes(coeff, 'test-filt')
+
+    # Using the absolute magnitudes of the templates as reference magnitudes
+    # should return one solar mass for each template.
+    stellar_mass = kcorrect.stellar_mass_remain(coeff, Mt, 'test-filt')
+    truth = kcorrect.mremain
+    np.testing.assert_allclose(stellar_mass, truth)
+
+    # Solution for given magnitudes without template mixing
+    Mb = np.array([10, 20, 30, 40, 50])
+    stellar_mass = kcorrect.stellar_mass_remain(coeff, Mb, 'test-filt')
+    truth = np.power(10, -0.4*(Mb-Mt)) * kcorrect.mremain
+    np.testing.assert_allclose(stellar_mass, truth)
+
+
 def test_kcorrect_metallicity():
 
     from skypy.galaxies.spectrum import kcorrect

skypy/galaxies/tests/test_stellar_mass.py

@@ -3,8 +3,10 @@
 import scipy.integrate
 from scipy.special import gammaln
 import pytest
+
 from hypothesis import given
 from hypothesis.strategies import integers
+from astropy.modeling.models import Exponential1D
 
 from skypy.galaxies import stellar_mass
 from skypy.utils import special
@@ -27,12 +29,30 @@ def test_stellar_masses():
     with pytest.raises(ValueError):
         stellar_mass.schechter_smf_mass(0., -1.4, np.array([1e10, 2e10]), 1.e7, 1.e13, size=3)
 
-    # Test that an array with the sme shape as m_star is returned if m_star is
+    # Test that an array with the same shape as m_star is returned if m_star is
     # an array and size = None
     m_star = np.array([1e10, 2e10])
     sample = stellar_mass.schechter_smf_mass(0., -1.4, m_star, 1.e7, 1.e13, size=None)
     assert m_star.shape == sample.shape
 
+    # Test m_star can be a function that returns an array of values for each redshift
+    redshift = np.linspace(0, 2, 100)
+    alpha = 0.357
+    m_min = 10 ** 7
+    m_max = 10 ** 14
+    m_star_function = Exponential1D(10**10.626, np.log(10)/0.095)
+    sample = stellar_mass.schechter_smf_mass(redshift, alpha, m_star_function, m_min, m_max)
+    assert sample.shape == redshift.shape
+
+    # Test alpha can be a function returning a scalar value.
+    sample = stellar_mass.schechter_smf_mass(redshift, lambda z: alpha, 10**10.626, m_min, m_max)
+    assert sample.shape == redshift.shape
+
+    # Sampling with an array for alpha is not implemented
+    alpha_array = np.full_like(redshift, alpha)
+    with pytest.raises(NotImplementedError, match='only scalar alpha is supported'):
+        sample = stellar_mass.schechter_smf_mass(redshift, alpha_array, m_star, m_min, m_max)
+
     # Test that sampling corresponds to sampling from the right pdf.
     # For this, we sample an array of luminosities for redshift z = 1.0 and we
     # compare it to the corresponding cdf.

skypy/utils/photometry.py

@@ -5,6 +5,7 @@
 
 from abc import ABCMeta, abstractmethod
 import numpy as np
+import scipy.special
 
 
 __all__ = [
@@ -13,6 +14,8 @@
     'luminosity_in_band',
     'mag_ab',
     'SpectrumTemplates',
+    'magnitude_error_rykoff',
+    'logistic_completeness_function',
 ]
 
 try:
@@ -304,3 +307,141 @@ def absolute_magnitude_from_luminosity(luminosity, zeropoint=None):
         zeropoint = -luminosity_in_band[zeropoint]
 
     return -2.5*np.log10(luminosity) - zeropoint
+
+
+def magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b, error_limit=np.inf):
+    r"""Magnitude error according to the model from Rykoff et al. (2015).
+
+    Given an apparent magnitude calculate the magnitude error that is introduced
+    by the survey specifications and follows the model described in Rykoff et al. (2015).
+
+    Parameters
+    ----------
+    magnitude: array_like
+        Apparent magnitude. This and the other array_like parameters must
+        be broadcastable to the same shape.
+    magnitude_limit: array_like
+        :math:`10\sigma` limiting magnitude of the survey. This and the other
+        array_like parameters must be broadcastable to the same shape.
+    magnitude_zp: array_like
+        Zero-point magnitude of the survey. This and the other array_like parameters must
+        be broadcastable to the same shape.
+    a,b: array_like
+        Model parameters: a is the intercept and
+        b is the slope of the logarithmic effective time.
+        These and the other array_like parameters must be broadcastable to the same shape.
+    error_limit: float, optional
+        Upper limit of the returned error. If given, all values larger than this value
+        will be set to error_limit. Default is None.
+
+    Returns
+    -------
+    error: ndarray
+        The apparent magnitude error in the Rykoff et al. (2015) model. This is a scalar
+        if magnitude, magnitude_limit, magnitude_zp, a and b are scalars.
+
+    Notes
+    -----
+    Rykoff et al. (2015) (see [1]_) describe the error of the apparent magnitude :math:`m` as
+
+    .. math::
+
+        \sigma_m(m;m_{\mathrm{lim}}, t_{\mathrm{eff}}) &=
+            \sigma_m(F(m);F_{\mathrm{noise}}(m_{\mathrm{lim}}), t_{\mathrm{eff}}) \\
+        &= \frac{2.5}{\ln(10)} \left[ \frac{1}{Ft_{\mathrm{eff}}}
+            \left( 1 + \frac{F_{\mathrm{noise}}}{F} \right) \right]^{1/2} \;,
+
+    where
+
+    .. math::
+
+        F=10^{-0.4(m - m_{\mathrm{ZP}})}
+
+    is the source's flux,
+
+    .. math::
+
+        F_\mathrm{noise} = \frac{F_{\mathrm{lim}}^2 t_{\mathrm{eff}}}{10^2} - F_{\mathrm{lim}}
+
+    is the effective noise flux and :math:`t_\mathrm{eff}` is the effective exposure time
+    (we absorbed the normalisation constant :math:`k` in the definition of
+    :math:`t_\mathrm{eff}`).
+    Furthermore, :math:`m_\mathrm{ZP}` is the zero-point magnitude of the survey and
+    :math:`F_\mathrm{lim}` is the :math:`10\sigma` limiting flux.
+    Accordingly, :math:`m_\mathrm{lim}` is the :math:`10\sigma` limiting magnitud
+    associated with :math:`F_\mathrm{lim}`.
+
+    The effective exposure time is described by
+
+    .. math::
+
+        \ln{t_\mathrm{eff}} = a + b(m_\mathrm{lim} - 21)\;,
+
+    where :math:`a` and :math:`b` are free parameters.
+
+    Further note that the model was originally used for SDSS galaxy photometry.
+
+    References
+    ----------
+    .. [1] Rykoff E. S., Rozo E., Keisler R., 2015, eprint arXiv:1509.00870
+
+    """
+
+    flux = luminosity_from_absolute_magnitude(magnitude, -magnitude_zp)
+    flux_limit = luminosity_from_absolute_magnitude(magnitude_limit, -magnitude_zp)
+    t_eff = np.exp(a + b * np.subtract(magnitude_limit, 21.0))
+    flux_noise = np.square(flux_limit / 10) * t_eff - flux_limit
+    error = 2.5 / np.log(10) * np.sqrt((1 + flux_noise / flux) / (flux * t_eff))
+
+    return np.minimum(error, error_limit)
+
+
+def logistic_completeness_function(magnitude, magnitude_95, magnitude_50):
+    r'''Logistic completeness function.
+
+    This function calculates the logistic completeness function (based on eq. (7) in
+    López-Sanjuan, C. et al. (2017) [1]_.
+
+    .. math::
+
+        p(m) = \frac{1}{1 + \exp[\kappa (m - m_{50})]}\;,
+
+    which describes the probability :math:`p(m)` that an object of magnitude :math:`m` is detected
+    in a specific band and with
+
+    .. math::
+
+        \kappa = \frac{\ln(\frac{1}{19})}{m_{95} - m_{50}}\;.
+
+    Here, :math:`m_{95}` and :math:`m_{50}` are the 95% and 50% completeness
+    magnitudes, respectively.
+
+    Parameters
+    ----------
+    magnitude : array_like
+        Magnitudes. Can be multidimensional for computing with multiple filter bands.
+    magnitude_95 : scalar or 1-D array_like
+        95% completeness magnitude.
+        If `magnitude_50` is 1-D array it has to be scalar or 1-D array of the same shape.
+    magnitude_50 : scalar or 1-D array_like
+        50% completeness magnitude.
+        If `magnitude_95` is 1-D array it has to be scalar or 1-D array of the same shape.
+
+    Returns
+    -------
+    probability : scalar or array_like
+        Probability of detecting an object with magnitude :math:`m`.
+        Returns array_like of the same shape as magnitude.
+        Exemption: If magnitude is scalar and `magnitude_95` or `magnitude_50`
+        is array_like of shape (nb, ) it returns array_like of shape (nb, ).
+
+    References
+    -----------
+    .. [1] López-Sanjuan, C. et al., `2017A&A…599A..62L`_
+    .. _2017A&A…599A..62L: https://ui.adsabs.harvard.edu/abs/2017A%26A...599A..62L
+
+    '''
+
+    kappa = np.log(1. / 19) / np.subtract(magnitude_95, magnitude_50)
+    arg = kappa * np.subtract(magnitude, magnitude_50)
+    return scipy.special.expit(-arg)

skypy/utils/tests/test_photometry.py

@@ -175,3 +175,58 @@ def test_speclite_not_installed():
     filter = 'bessell-B'
     with pytest.raises(ImportError):
         mag_ab(wavelength, spectrum, filter)
+
+
+def test_magnitude_error_rykoff():
+    from skypy.utils.photometry import magnitude_error_rykoff
+
+    # Test broadcasting to same shape given array for each parameter and
+    # test for correct result.
+    magnitude = np.full((2, 1, 1, 1, 1), 21)
+    magnitude_limit = np.full((3, 1, 1, 1), 21)
+    magnitude_zp = np.full((5, 1, 1), 21)
+    a = np.full((7, 1), np.log(200))
+    b = np.zeros(11)
+    error = magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b)
+    # test result
+    assert np.allclose(error, 0.25 / np.log(10))
+    # test shape
+    assert error.shape == (2, 3, 5, 7, 11)
+
+    # second test for result
+    magnitude = 20
+    magnitude_limit = 22.5
+    magnitude_zp = 25
+    b = 2
+    a = np.log(10) - 1.5 * b
+    error = magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b)
+    assert np.isclose(error, 0.25 / np.log(10) / np.sqrt(10))
+
+    # test that error limit is returned if error is larger than error_limit
+    # The following set-up would give a value larger than 10
+    magnitude = 30
+    magnitude_limit = 25
+    magnitude_zp = 30
+    a = 0.5
+    b = 1.0
+    error_limit = 1
+    error = magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b, error_limit)
+    assert error == error_limit
+
+
+def test_logistic_completeness_function():
+    from skypy.utils.photometry import logistic_completeness_function
+
+    # Test that arguments broadcast correctly
+    m = np.full((2, 1, 1), 21)
+    m95 = np.full((3, 1), 22)
+    m50 = np.full(5, 23)
+    p = logistic_completeness_function(m, m95, m50)
+    assert p.shape == np.broadcast(m, m95, m50).shape
+
+    # Test result of completeness function for different given magnitudes
+    m95 = 24
+    m50 = 25
+    m = [np.finfo(np.float64).min, m95, m50, 2*m50-m95, np.finfo(np.float64).max]
+    p = logistic_completeness_function(m, m95, m50)
+    assert np.allclose(p, [1, 0.95, 0.5, 0.05, 0])