From 875748200d6cd0e86f6354356f6ddb6fe4b6b83c Mon Sep 17 00:00:00 2001 From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com> Date: Wed, 6 Jul 2022 14:52:20 +0100 Subject: [PATCH 1/6] Ad functions and tests for the redshift-dependent amplitude --- skypy/galaxies/stellar_mass.py | 66 +++++++++++++++++++++++ skypy/galaxies/tests/test_stellar_mass.py | 29 ++++++++++ 2 files changed, 95 insertions(+) diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py index 0fd3cd80..373772b9 100644 --- a/skypy/galaxies/stellar_mass.py +++ b/skypy/galaxies/stellar_mass.py @@ -2,6 +2,7 @@ """ import numpy as np +from scipy.integrate import trapz from ..utils.random import schechter from ..utils import dependent_argument @@ -9,6 +10,7 @@ __all__ = [ 'schechter_smf_mass', + 'schechter_smf_phi_active_redshift', ] @@ -78,3 +80,67 @@ def schechter_smf_mass(redshift, alpha, m_star, m_min, m_max, size=None, m = schechter(alpha, x_min, x_max, resolution, size=size, scale=m_star) return m + + +def schechter_smf_phi_active_redshift(redshift, + phi_active_today, + probability_satellites, + cosmology, + redshift_initial=10, + resolution=100): + r'''Time-dependent Schechter amplitude of active galaxies. + This function returns the time-dependent Schechter mass function amplitude + for the total active population based on equation (17) + in de la Bella et al. 2021 [1]_. + + Parameters + ---------- + redshift: + Values of redshift at which to evaluate the amplitude. + + phi_active_today: array_like + Schechter mass function amplitude for the entire active + sample of galaxies in the past when density and mass was very low. + + probability_satellites: func + Probability of active galaxies becoming satellites + as a function of redshift. + + cosmology: astropy.cosmology.Cosmology + Cosmology object providing methods for the evolution history + of omega_matter and omega_lambda with redshift. + + redshift_initial: float + Value of redshift in the past when density and mass was very low. + Default is 10. + + resolution: float + Resolution of the integral. Default is 100. + + Returns + ------- + amplitude: array_like + Amplitude of the Schechter mass function. + + References + ---------- + .. [1] de la Bella et al. 2021, Quenching and Galaxy Demographics, + arXiv 2112.11110. + + ''' + + def integrand(z): + return probability_satellites(z) / cosmology.H(z).value / (1 + z) + + # Calculate the amplitude in the past + redshift_today = np.linspace(0, redshift_initial, resolution) + integrand_today = integrand(redshift_today) + integral_today = trapz(integrand_today) + B = phi_active_today * np.exp(- integral_today) / (1 - integral_today) + + # Calculate the amplitude for the given redshift + redshift_array = np.linspace(redshift, redshift_initial, resolution) + integrand_redshift = integrand(redshift_array) + integral = trapz(integrand_redshift) + + return B * (1 - integral) * np.exp(integral) diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py index b8809235..ee6d5242 100644 --- a/skypy/galaxies/tests/test_stellar_mass.py +++ b/skypy/galaxies/tests/test_stellar_mass.py @@ -4,6 +4,8 @@ from scipy.special import gammaln import pytest from astropy.modeling.models import Exponential1D +from hypothesis import given +from hypothesis.strategies import integers from skypy.galaxies import stellar_mass from skypy.utils import special @@ -79,3 +81,30 @@ def calc_cdf(m): size=1000, resolution=100) p_value = scipy.stats.kstest(sample, calc_cdf)[1] assert p_value >= 0.01 + + +@given(integers(1, 100)) +def test_schechter_smf_phi_active_redshift(redshift): + # Amplitude today and cosmology + from astropy.cosmology import Planck15 + phi_today = 10**-2.423 + + # Simple probability of satellites + def prob_sat(z, a=0.4): + return a * np.exp(- z) + + # Test that today the output of the active amplitude equals the input + output_today = stellar_mass.schechter_smf_phi_active_redshift(0, phi_today, prob_sat, Planck15) + assert output_today == phi_today + + # Test that at any redshift for a null probability of creating satellite galaxies + # the output of the active amplitude equals the input + def prob_null(z): + return prob_sat(z, 0) + + output_null = stellar_mass.schechter_smf_phi_active_redshift(redshift, + phi_today, + prob_null, + Planck15 + ) + assert output_null == phi_today From 3afa6f44a293fbf133b9a566a3aab1f3402d26ac Mon Sep 17 00:00:00 2001 From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com> Date: Wed, 6 Jul 2022 14:56:33 +0100 Subject: [PATCH 2/6] Add function to docs --- docs/galaxies.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/galaxies.rst b/docs/galaxies.rst index bfddeeec..e7efb4c7 100644 --- a/docs/galaxies.rst +++ b/docs/galaxies.rst @@ -89,6 +89,7 @@ The following models are found in the `skypy.galaxies.stellar_mass` package. :nosignatures: schechter_smf_mass + schechter_smf_phi_active_redshift Reference/API From 2acb259f0947a2700443e73ca5ce3619c3a22d3e Mon Sep 17 00:00:00 2001 From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com> Date: Wed, 6 Jul 2022 15:01:02 +0100 Subject: [PATCH 3/6] Fix title of docstring --- skypy/galaxies/stellar_mass.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py index 373772b9..0ebe6aa0 100644 --- a/skypy/galaxies/stellar_mass.py +++ b/skypy/galaxies/stellar_mass.py @@ -88,7 +88,7 @@ def schechter_smf_phi_active_redshift(redshift, cosmology, redshift_initial=10, resolution=100): - r'''Time-dependent Schechter amplitude of active galaxies. + r'''Redshift-dependent Schechter amplitude of active galaxies. This function returns the time-dependent Schechter mass function amplitude for the total active population based on equation (17) in de la Bella et al. 2021 [1]_. From 5f09795b6d2839bfd6533aa93e643e9184636538 Mon Sep 17 00:00:00 2001 From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com> Date: Wed, 13 Jul 2022 09:51:06 +0100 Subject: [PATCH 4/6] Fix test --- skypy/galaxies/tests/test_stellar_mass.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py index ee6d5242..0ab0bb3a 100644 --- a/skypy/galaxies/tests/test_stellar_mass.py +++ b/skypy/galaxies/tests/test_stellar_mass.py @@ -95,7 +95,7 @@ def prob_sat(z, a=0.4): # Test that today the output of the active amplitude equals the input output_today = stellar_mass.schechter_smf_phi_active_redshift(0, phi_today, prob_sat, Planck15) - assert output_today == phi_today + assert np.testing.assert_allclose(output_today, phi_today, rtol=1e-5, atol=0) # Test that at any redshift for a null probability of creating satellite galaxies # the output of the active amplitude equals the input @@ -107,4 +107,4 @@ def prob_null(z): prob_null, Planck15 ) - assert output_null == phi_today + assert np.testing.assert_allclose(output_null, phi_today, rtol=1e-5, atol=0) From bde5aaac8f77f8b24a03745694d4388af9283dda Mon Sep 17 00:00:00 2001 From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com> Date: Wed, 13 Jul 2022 09:56:38 +0100 Subject: [PATCH 5/6] Fix test with pytest.approx --- skypy/galaxies/tests/test_stellar_mass.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py index 0ab0bb3a..f6a6e0a3 100644 --- a/skypy/galaxies/tests/test_stellar_mass.py +++ b/skypy/galaxies/tests/test_stellar_mass.py @@ -95,7 +95,7 @@ def prob_sat(z, a=0.4): # Test that today the output of the active amplitude equals the input output_today = stellar_mass.schechter_smf_phi_active_redshift(0, phi_today, prob_sat, Planck15) - assert np.testing.assert_allclose(output_today, phi_today, rtol=1e-5, atol=0) + assert output_today == approx(phi_today) # Test that at any redshift for a null probability of creating satellite galaxies # the output of the active amplitude equals the input @@ -107,4 +107,4 @@ def prob_null(z): prob_null, Planck15 ) - assert np.testing.assert_allclose(output_null, phi_today, rtol=1e-5, atol=0) + assert output_null == approx(phi_today) From 83b7389a4834f0624459b64244459b98f2057b07 Mon Sep 17 00:00:00 2001 From: "Lucia F. de la Bella" <55983939+Lucia-Fonseca@users.noreply.github.com> Date: Wed, 13 Jul 2022 10:01:29 +0100 Subject: [PATCH 6/6] Fix flake8 errors --- skypy/galaxies/tests/test_stellar_mass.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py index f6a6e0a3..b5b38277 100644 --- a/skypy/galaxies/tests/test_stellar_mass.py +++ b/skypy/galaxies/tests/test_stellar_mass.py @@ -95,7 +95,7 @@ def prob_sat(z, a=0.4): # Test that today the output of the active amplitude equals the input output_today = stellar_mass.schechter_smf_phi_active_redshift(0, phi_today, prob_sat, Planck15) - assert output_today == approx(phi_today) + assert output_today == pytest.approx(phi_today) # Test that at any redshift for a null probability of creating satellite galaxies # the output of the active amplitude equals the input @@ -107,4 +107,4 @@ def prob_null(z): prob_null, Planck15 ) - assert output_null == approx(phi_today) + assert output_null == pytest.approx(phi_today)