diff --git a/.github/workflows/continuous_integration.yml b/.github/workflows/continuous_integration.yml index 7cf25b6c13..f15ec9990c 100644 --- a/.github/workflows/continuous_integration.yml +++ b/.github/workflows/continuous_integration.yml @@ -36,4 +36,10 @@ jobs: export PATH_TO_POSYDON=./ export PATH_TO_POSYDON_DATA=./posydon/unit_tests/_data/ export MESA_DIR=./ - python -m pytest posydon/unit_tests/ --cov=posydon.config --cov=posydon.utils --cov=posydon.grids --cov-branch --cov-report term-missing --cov-fail-under=100 + python -m pytest posydon/unit_tests/ --cov=posydon.utils \ + --cov=posydon.config \ + --cov=posydon.grids \ + --cov=posydon.popsyn.star_formation_history \ + --cov-branch \ + --cov-report term-missing \ + --cov-fail-under=100 diff --git a/posydon/popsyn/rate_calculation.py b/posydon/popsyn/rate_calculation.py index 3657403ffb..70f4deb531 100644 --- a/posydon/popsyn/rate_calculation.py +++ b/posydon/popsyn/rate_calculation.py @@ -14,12 +14,13 @@ from astropy import units as u -DEFAULT_MODEL = { +DEFAULT_SFH_MODEL = { "delta_t": 100, # Myr "SFR": "IllustrisTNG", "sigma_SFR": None, - "Z_max": 1.0, - "select_one_met": False, + "Z_max": None, # Zsun + "Z_min": None, # Zsun + "normalise": True, # normalise the SFR to 1 "dlogZ": None, # e.g, [np.log10(0.0142/2),np.log10(0.0142*2)] "Zsun": Zsun, } @@ -150,7 +151,6 @@ def get_redshift_from_cosmic_time(t_cosm): return trained_tz_interp(t_cosm) - def get_redshift_bin_edges(delta_t): """Compute the redshift bin edges. diff --git a/posydon/popsyn/star_formation_history.py b/posydon/popsyn/star_formation_history.py index c314a66663..196c79a2b5 100644 --- a/posydon/popsyn/star_formation_history.py +++ b/posydon/popsyn/star_formation_history.py @@ -12,6 +12,7 @@ import os import numpy as np import scipy as sp +import pandas as pd from scipy import stats from posydon.config import PATH_TO_POSYDON_DATA from posydon.utils.constants import age_of_universe @@ -22,7 +23,8 @@ ) from posydon.utils.constants import Zsun from posydon.utils.interpolators import interp1d -from astropy.cosmology import Planck15 as cosmology +from abc import ABC, abstractmethod +from posydon.utils.posydonwarning import Pwarn SFH_SCENARIOS = [ @@ -35,7 +37,1010 @@ ] -def get_formation_times(N_binaries, star_formation="constant", **kwargs): +class SFHBase(ABC): + """Abstract class for star formation history models""" + def __init__(self, SFH_MODEL): + """Initialise the SFH model + + Adds the model parameters as attributes. + + Parameters + ---------- + SFH_MODEL : dict + SFH model parameters. + - Z_max : float + The maximum metallicity in absolute units. + - Z_min : float + The minimum metallicity in absolute units. + - normalise : bool + Normalise the metallicity distribution to 1. + """ + self.Z_max = None + self.Z_min = None + self.normalise = False + + self.SFH_MODEL = SFH_MODEL + # Automatically attach all model parameters as attributes + for key, value in SFH_MODEL.items(): + setattr(self, key, value) + + # check if Z_max is not larger than 1 + if self.Z_max is not None: + if self.Z_max > 1: + raise ValueError("Z_max must be in absolute units! " + "It cannot be larger than 1!") + elif self.Z_max < 0: + raise ValueError("Z_max must be in absolute units! " + "It cannot be negative!") + if self.Z_min is not None: + if self.Z_min < 0: + raise ValueError("Z_min must be in absolute units! " + "It cannot be negative!") + elif self.Z_min > 1: + raise ValueError("Z_min must be in absolute units! " + "It cannot be larger than 1!") + if self.Z_min is not None and self.Z_max is not None: + if self.Z_min >= self.Z_max: + raise ValueError("Z_min must be smaller than Z_max!") + + + @abstractmethod + def CSFRD(self, z): # pragma: no cover + """Compute the cosmic star formation rate density. + + This is an abstract method that must be implemented by subclasses. + The implementation should calculate and return the cosmic star formation + rate density at the given redshift(s). + + Parameters + ---------- + z : float or array-like + Cosmological redshift. + + Returns + ------- + float or array-like + The cosmic star formation rate density at the given redshift(s). + + Raises + ------- + NotImplementedError: If the subclass does not implement this method. + """ + pass + + @abstractmethod + def mean_metallicity(self, z): # pragma: no cover + """Return the mean metallicity at redshift z. + + This is an abstract method that must be implemented by subclasses. + The implementation should calculate and return the mean metallicity + at the given redshift(s). + + Parameters + ---------- + z : float or array-like + Cosmological redshift. + + Returns + ------- + float or array-like + The mean metallicity at the given redshift(s). + + Raises + ------- + NotImplementedError: If the subclass does not implement this method. + """ + pass + + @abstractmethod + def fSFR(self, z, metallicity_bins): # pragma: no cover + """Compute the star formation rate fraction (fSFR) at a given redshift + using the specified metallicity bins. + + This is an abstract method that must be implemented by subclasses. + The implementation should calculate and return the fractional SFR per + metallicity bins at the provided redshift (z). + + Parameters + --------- + z : float or array-like + The redshift(s) at which to compute the star formation rate. + metallicity_bins : list or array-like + The metallicity bin boundaries or labels used in the computation to + account for different metallicity contributions. + + Returns + ------- + float or array-like + The calculated star formation rate at the given redshift(s) and + metallicity bins in Msun/yr. + + Raises: + NotImplementedError: If the subclass does not implement this method. + """ + pass + + def _distribute_cdf(self, cdf_func, metallicity_bins): + '''Distribute the SFR over the metallicity bins using the CDF + + Parameters + ---------- + cdf_func : function + The cumulative density function to use. + metallicity_bins : array + Metallicity bins edges in absolute metallicity. + + Returns + ------- + ndarray + Fraction of the SFR in the given metallicity bin at the given redshift. + ''' + # verify if the metallicity bins are sorted + if not np.all(np.diff(metallicity_bins) > 0): + raise ValueError("Metallicity bins must be sorted " + "in ascending order.") + + fSFR = (np.array(cdf_func(metallicity_bins[1:])) + - np.array(cdf_func(metallicity_bins[:-1]))) + + first_bin_index = 0 + last_bin_index = len(fSFR) - 1 + + # include material outside the metallicity bounds if requested + if self.Z_max is not None: + if self.Z_max >= metallicity_bins[-1]: + fSFR[-1] = cdf_func(self.Z_max) - cdf_func(metallicity_bins[-2]) + else: + Pwarn('Z_max is smaller than the highest metallicity bin.', + 'SFHModelWarning') + # find the index of the last bin that is smaller than Z_max + last_bin_index = np.searchsorted(metallicity_bins, self.Z_max) - 1 + fSFR[last_bin_index] = cdf_func(self.Z_max) - cdf_func(metallicity_bins[last_bin_index]) + fSFR[last_bin_index + 1:] = 0.0 + + if self.Z_min is not None: + if self.Z_min <= metallicity_bins[0]: + fSFR[0] = cdf_func(metallicity_bins[1]) - cdf_func(self.Z_min) + else: + Pwarn('Z_min is larger than the lowest metallicity bin.', + 'SFHModelWarning') + # find the index of the first bin that is larger than Z_min + first_bin_index = np.searchsorted(metallicity_bins, self.Z_min) -1 + fSFR[:first_bin_index] = 0.0 + fSFR[first_bin_index] = cdf_func(metallicity_bins[first_bin_index+1]) - cdf_func(self.Z_min) + + # Check if in the same bin + if self.Z_max is not None and self.Z_min is not None: + if first_bin_index == last_bin_index: + fSFR[first_bin_index] = cdf_func(self.Z_max) - cdf_func(self.Z_min) + fSFR[first_bin_index + 1:] = 0.0 + fSFR[:first_bin_index] = 0.0 + + if self.normalise: + fSFR /= np.sum(fSFR) + + return fSFR + + + def __call__(self, z, met_bins): + """Return the star formation history at a given redshift and metallicity + bins + + Parameters + ---------- + z : float or array-like + Cosmological redshift. + met_bins : array + Metallicity bins edges in absolute metallicity. + + Returns + ------- + ndarray + Star formation history per metallicity bin at the given redshift(s). + """ + return self.CSFRD(z)[:, np.newaxis] * self.fSFR(z, met_bins) + +class MadauBase(SFHBase): + """ + Base class for Madau-style star-formation history implementations. + This class implements common methods for CSFRD, mean metallicity, + and fractional SFR based on the chosen Madau parameterisation. + The specific parameters for CSFRD must be provided by subclasses. + """ + def __init__(self, SFH_MODEL): + """Initialise the MadauBase class + + Parameters + ---------- + SFH_MODEL : dict + SFH model parameters. MadauBase requires the following parameters: + - sigma : float or str + The standard deviation of the log-normal metallicity distribution. + Options are: + - Bavera+20 + - Neijssel+19 + - float + Additional SFH model parameters: + - Z_max : float + The maximum metallicity in absolute units. + - Z_min : float + The minimum metallicity in absolute units. + - normalise : bool + Normalise the metallicity distribution to 1. + - CSFRD_params: dict + Parameters for the cosmic star formation rate density (CSFRD) + - a, b, c, d : float + Follows the Madau & Dickinson (2014) CSFRD formula (Eq. 15) in [1]_ + + References + ---------- + .. [1] Madau, P., & Dickinson, M. (2014). Cosmic star formation history. + Annual Review of Astronomy and Astrophysics, 52, 415-486. + https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..415M/abstract + + """ + if "sigma" not in SFH_MODEL: + raise ValueError("sigma not given!") + self.CSFRD_params = None + super().__init__(SFH_MODEL) + + def CSFRD(self, z): + """The cosmic star formation rate density at a given redshift. + + Follows the Madau & Dickinson (2014) cosmic star formation rate + density formula. + + Parameters + ---------- + z : float or np.array + Cosmological redshift. + + Returns + ------- + float or array + The cosmic star formation rate density at the given redshift. + """ + p = self.CSFRD_params + return p["a"] * (1.0 + z) ** p["b"] / (1.0 + ((1.0 + z) / p["c"]) ** p["d"]) + + def std_log_metallicity_dist(self): + """Return the standard deviation of the log-normal metallicity + distribution + + Either recognised the strings "Bavera+20" (sigma=0.5) + or "Neijssel+19" (sigma=0.39) or a float value. + + Returns + ------- + float + The standard deviation of the log-normal metallicity distribution. + """ + sigma = self.sigma + if isinstance(sigma, str): + if sigma == "Bavera+20": + return 0.5 + elif sigma == "Neijssel+19": + return 0.39 + else: + raise ValueError("Unknown sigma choice!") + elif isinstance(sigma, (float, int)): + return np.float64(sigma) + else: + raise ValueError(f"Invalid sigma value {sigma}!") + + def mean_metallicity(self, z): + """The mean metallicity at a given redshift + + Follows Madau & Fragos (2017) mean metallicity evolution + + Parameters + ---------- + z : float or np.array + Cosmological redshift. + + Returns + ------- + float or array + The mean metallicity at the given redshift. + """ + return 10 ** (0.153 - 0.074 * z ** 1.34) * Zsun + + def fSFR(self, z, metallicity_bins): + """Fraction of the SFR at a given redshift z in a given metallicity + bin as described in Bavera et al. (2020). + + Parameters + ---------- + z : np.array + Cosmological redshift. + metallicity_bins : array + Metallicity bins edges in absolute metallicity. + + Returns + ------- + ndarray + Fraction of the SFR in the given metallicity bin at the given + redshift. + """ + sigma = self.std_log_metallicity_dist() + # Compute mu; if z is an array, mu will be an array. + mu = np.log10(self.mean_metallicity(z)) - sigma ** 2 * np.log(10) / 2.0 + # Ensure mu is an array for consistency + mu_array = np.atleast_1d(mu) + + fSFR = np.zeros((len(z), len(metallicity_bins) - 1)) + for i, mean in enumerate(mu_array): + cdf_func = lambda x: stats.norm.cdf(np.log10(x), mean, sigma) + fSFR[i, :] = self._distribute_cdf(cdf_func, metallicity_bins) + + return fSFR + +class MadauDickinson14(MadauBase): + """Madau & Dickinson (2014) [1]_ star formation history model using the + mean metallicity evolution of Madau & Fragos (2017) [2]_. + + References + ---------- + .. [1] Madau, P., & Dickinson, M. (2014). ARA&A, 52, 415-486. + https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..415M + .. [2] Madau, P., & Fragos, T. (2017). ApJ, 840(1), 39. + https://ui.adsabs.harvard.edu/abs/2017ApJ...840...39M + """ + + def __init__(self, SFH_MODEL): + """Initialise the Madau & Dickinson (2014) [1]_ SFH model with the + metallicity evolution of Madau & Fragos (2017) [2]_. + + Parameters + ---------- + SFH_MODEL : dict + SFH model parameters. Madau+14 requires the following parameters: + - sigma : float or str + The standard deviation of the log-normal metallicity + distribution. + Options are: + - Bavera+20 + - Neijssel+19 + - float + Additional SFH model parameters: + - Z_max : float + The maximum metallicity in absolute units. + - Z_min : float + The minimum metallicity in absolute units. + - normalise : bool + Normalise the metallicity distribution to 1. + - CSFRD_params: dict + Parameters for the cosmic star formation rate density (CSFRD) + - a, b, c, d : float + Follows the Madau & Dickinson (2014) CSFRD formula (Eq. 15) in [1]_ + + References + ---------- + .. [1] Madau, P., & Dickinson, M. (2014). ARA&A, 52, 415-486. + https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..415M + .. [2] Madau, P., & Fragos, T. (2017). ApJ, 840(1), 39. + https://ui.adsabs.harvard.edu/abs/2017ApJ...840...39M + """ + super().__init__(SFH_MODEL) + # Parameters for Madau+Dickinson14 CSFRD + self.CSFRD_params = { + "a": 0.015, + "b": 2.7, + "c": 2.9, + "d": 5.6, + } + +class MadauFragos17(MadauBase): + """The Madau & Fragos (2017) star formation history model with the + metallicity evolution of Madau & Fragos (2017) [1]_. + + + References + ---------- + .. [1] Madau, P., & Fragos, T. (2017). ApJ, 840(1), 39. + https://ui.adsabs.harvard.edu/abs/2017ApJ...840...39M/abstract + """ + + def __init__(self, SFH_MODEL): + """Initialise the Madau+17 model + + Parameters + ---------- + SFH_MODEL : dict + SFH model parameters. Madau+17 requires the following parameters: + - sigma : float or str + The standard deviation of the log-normal metallicity distribution. + Options are: + - Bavera+20 + - Neijssel+19 + - float + Additional SFH model parameters: + - Z_max : float + The maximum metallicity in absolute units. + - Z_min : float + The minimum metallicity in absolute units. + - normalise : bool + Normalise the metallicity distribution to 1. + """ + super().__init__(SFH_MODEL) + # Parameters for Madau+Fragos17 CSFRD + self.CSFRD_params = { + "a": 0.01, + "b": 2.6, + "c": 3.2, + "d": 6.2, + } + +class Neijssel19(MadauBase): + """The Neijssel et al. (2019) [1]_ star formation history model, which fits + the Madau & Dickinson (2014) [2]_ cosmic star formation rate density formula + with the BBH merger rate and uses a truncated log-normal distribution for + the mean metallicity distribution. + The mean metallicity evolution follows the Langer and Normal parameterisation + also fitted to the BBH merger rate. + + + References + ---------- + .. [1] Neijssel, C. J., et al. (2019). MNRAS, 490, 3740. + https://ui.adsabs.harvard.edu/abs/2019MNRAS.490.3740N + .. [2] Madau, P., & Fragos, T. (2017). ApJ, 840(1), 39. + https://ui.adsabs.harvard.edu/abs/2017ApJ...840...39M/ + """ + def __init__(self, SFH_MODEL): + """Initialise the Neijssel+19 model + + Parameters + ---------- + SFH_MODEL : dict + SFH model parameters. Neijssel+19 requires the following parameters: + - sigma : float or str + The standard deviation of the log-normal metallicity distribution. + Options are: + - Bavera+20 + - Neijssel+19 + - float + Additional SFH model parameters: + - Z_max : float + The maximum metallicity in absolute units. + - Z_min : float + The minimum metallicity in absolute units. + - normalise : bool + Normalise the metallicity distribution to 1. + """ + super().__init__(SFH_MODEL) + # Parameters for Neijssel+19 CSFRD + self.CSFRD_params = { + "a": 0.01, + "b": 2.77, + "c": 2.9, + "d": 4.7, + } + + # overwrite mean_metallicity method of MadauBase + def mean_metallicity(self, z): + """Calculate the mean metallicity at a given redshift + + Overwrites the mean_metallicity method of MadauBase class. + + Parameters + ---------- + z : float or array-like + Cosmological redshift. + + Returns + ------- + float or array-like + The mean metallicity at the given redshift(s). + """ + return 0.035 * 10 ** (-0.23 * z) + + def fSFR(self, z, metallicity_bins): + """Fraction of the SFR at a given redshift z in a given metallicity bin + as described in Neijssel et al. (2019). + + Overwrites the fSFR method of MadauBase class. + + Parameters + ---------- + z : np.array + Cosmological redshift. + metallicity_bins : array + Metallicity bins edges in absolute metallicity. + + Returns + ------- + ndarray + Fraction of the SFR in the given metallicity bins at the + given redshift. + """ + # assume a truncated ln-normal distribution of metallicities + sigma = self.std_log_metallicity_dist() + mu = np.log(self.mean_metallicity(z)) - sigma**2 / 2.0 + mu_array = np.atleast_1d(mu) + fSFR = np.zeros((len(z), len(metallicity_bins) - 1)) + for i, mean in enumerate(mu_array): + cdf_func = lambda x: stats.norm.cdf(np.log(x), mean, sigma) + fSFR[i,:] = self._distribute_cdf(cdf_func, metallicity_bins) + return fSFR + +class IllustrisTNG(SFHBase): + """The IllustrisTNG star formation history model. + + Uses the TNG100-1 model from the IllustrisTNG simulation. + + https://www.tng-project.org/ + """ + + def __init__(self, SFH_MODEL): + """Initialise the IllustrisTNG model + + Parameters + ---------- + SFH_MODEL : dict + Additional SFH model parameters: + - Z_max : float + The maximum metallicity in absolute units. + - Z_min : float + The minimum metallicity in absolute units. + - normalise : bool + Normalise the metallicity distribution to 1. + """ + super().__init__(SFH_MODEL) + # load the TNG data + illustris_data = self._get_illustrisTNG_data() + # the data is stored in reverse order high to low redshift + self.CSFRD_data = np.flip(illustris_data["SFR"]) + self.redshifts = np.flip(illustris_data["redshifts"]) + + self.Z = illustris_data["mets"] + self.M = np.flip(illustris_data["M"], axis=0) # Msun + + def _get_illustrisTNG_data(self, verbose=False): # pragma: no cover + """Load IllustrisTNG SFR dataset into the class. + + Parameters + ---------- + verbose : bool, optional + Print information about the data loading. + """ + if verbose: + print("Loading IllustrisTNG data...") + return np.load(os.path.join(PATH_TO_POSYDON_DATA, "SFR/IllustrisTNG.npz")) + + def CSFRD(self, z): + """The cosmic star formation rate density at a given redshift. + + Parameters + ---------- + z : float or np.array + Cosmological redshift. + + Returns + ------- + float or array + The cosmic star formation rate density at the given redshift(s). + """ + SFR_interp = interp1d(self.redshifts, self.CSFRD_data) + return SFR_interp(z) + + def mean_metallicity(self, z): + """Calculate the mean metallicity at a given redshift + + Parameters + ---------- + z : float or array-like + Cosmological redshift. + + Returns + ------- + float or array-like + The mean metallicity at the given redshift(s). + """ + out = np.zeros_like(self.redshifts) + for i in range(len(out)): + if np.sum(self.M[i, :]) == 0: + out[i] = np.nan + else: + out[i] = np.average(self.Z, weights=self.M[i, :]) + Z_interp = interp1d(self.redshifts, out) + return Z_interp(z) + + def fSFR(self, z, metallicity_bins): + """Calculate the fractional SFR as a function of redshift and + metallicity bins. + + Parameters + ---------- + z : float or array-like + Cosmological redshift. + metallicity_bins : array + Metallicity bins edges in absolute metallicity. + + Returns + ------- + ndarray + Fraction of the SFR in the given metallicity bin at the given redshift. + """ + redshift_indices = np.array([np.where(self.redshifts <= i)[0][-1] for i in z]) + Z_dist = self.M[redshift_indices] + fSFR = np.zeros((len(z), len(metallicity_bins) - 1)) + for i in range(len(z)): + if Z_dist[i].sum() == 0.0: + continue + else: + # At a specific redshift, the SFR is distributed over the metallicities + # according to the mass distribution + # Add a final point to the CDF and metallicities to ensure normalisation to 1 + Z_dist_cdf = np.cumsum(Z_dist[i]) / Z_dist[i].sum() + Z_dist_cdf = np.append(Z_dist_cdf, 1) + + Z_x_values = np.append(np.log10(self.Z), 0) + Z_dist_cdf_interp = interp1d(Z_x_values, Z_dist_cdf) + cdf_fun = lambda x: Z_dist_cdf_interp(np.log10(x)) + fSFR[i, :] = self._distribute_cdf(cdf_fun, metallicity_bins) + + return fSFR + +class Chruslinska21(SFHBase): + """The Chruślińska+21 star formation history model [1]_. + + + References + ---------- + .. [1] Chruślińska, M., et al. (2021). MNRAS, 508, 4994. + https://ui.adsabs.harvard.edu/abs/2021MNRAS.508.4994C + + Data source: + https://ftp.science.ru.nl/astro/mchruslinska/Chruslinska_et_al_2021/ + """ + def __init__(self, SFH_MODEL): + """Initialise the Chruslinska+21 model + + Parameters + ---------- + SFH_MODEL : dict + SFH model parameters. Chruslinska+21 requires the + following parameters: + - sub_model : str + The sub-model to use. + This is the name of the file containing the data. + - Z_solar_scaling : str + The scaling of the solar metallicity. Options are: + - Asplund09 + - AndersGrevesse89 + - GrevesseSauval98 + - Villante14 + Additional SFH model parameters: + - Z_max : float + The maximum metallicity in absolute units. + - Z_min : float + The minimum metallicity in absolute units. + - normalise : bool + Normalise the metallicity distribution to 1. + """ + if "sub_model" not in SFH_MODEL: + raise ValueError("Sub-model not given!") + if "Z_solar_scaling" not in SFH_MODEL: + raise ValueError("Z_solar_scaling not given!") + super().__init__(SFH_MODEL) + self._load_chruslinska_data() + + def _load_chruslinska_data(self, verbose=False): + """Load the data from the Chruslinska+21 models. + Transforms the data to the format used in the classes. + + Parameters + ---------- + verbose : bool, optional + Print information about the data loading. + """ + # oxygen to hydrogen abundance ratio ( FOH == 12 + log(O/H) ) + # as used in the calculations - do not change + # This is the metallicity bin edges used in the Chruslinska+21 calculations + FOH_min = 5.3 + FOH_max = 9.7 + self.FOH_bins = np.linspace(FOH_min, FOH_max, 200) + self.dFOH = self.FOH_bins[1] - self.FOH_bins[0] + # Need to use the Z_solar_scaling parameter to + # convert the FOH bins to absolute metallicity. + # Use solar metallicity as the reference point + self.Z = self._FOH_to_Z(self.FOH_bins) + + self._data_folder = os.path.join(PATH_TO_POSYDON_DATA, "SFR/Chruslinska+21") + _, self.redshifts, delta_T = self._load_redshift_data(verbose) + M = self._load_raw_data() + self.SFR_data = np.array( [M[ii]/(1e6*delta_T[ii]) for ii in range(len(delta_T))])/self.dFOH + + def _FOH_to_Z(self, FOH): + """Convert the oxygen to hydrogen abundance ratio to absolute metallicity + + Parameters + ---------- + FOH : float or array-like + The oxygen to hydrogen abundance ratio. + + Returns + ------- + float or array-like + The absolute metallicity. + """ + # scalings from Chruslinksa+21 + if self.Z_solar_scaling == "Asplund09": + Zsun = 0.0134 + FOHsun = 8.69 + elif self.Z_solar_scaling == "AndersGrevesse89": + Zsun = 0.017 + FOHsun = 8.83 + elif self.Z_solar_scaling == "GrevesseSauval98": + Zsun = 0.0201 + FOHsun = 8.93 + elif self.Z_solar_scaling == "Villante14": + Zsun = 0.019 + FOHsun = 8.85 + else: + valid_options = ["Asplund09", "AndersGrevesse89", + "GrevesseSauval98", "Villante14"] + raise ValueError(f"Invalid Z_solar_scaling " + f"'{self.Z_solar_scaling}'. " + f"Valid options: {valid_options}") + logZ = np.log10(Zsun) + FOH - FOHsun + return 10**logZ + + def mean_metallicity(self, z): + """Calculate the mean metallicity at a given redshift + + Parameters + ---------- + z : float or array-like + Cosmological redshift. + + Returns + ------- + float or array-like + The mean metallicity at the given redshift(s). + """ + mean_over_redshift = np.zeros_like(self.redshifts) + for i in range(len(mean_over_redshift)): + if np.sum(self.SFR_data[i]) == 0: + mean_over_redshift[i] = np.nan + else: + mean_over_redshift[i] = np.average(self.Z, weights=self.SFR_data[i,:]*self.dFOH) + + Z_interp = interp1d(self.redshifts, mean_over_redshift) + return Z_interp(z) + + def fSFR(self, z, metallicity_bins): + """Calculate the fractional SFR as a function of redshift and metallicity bins + + Parameters + ---------- + z : float or array-like + Cosmological redshift. + metallicity_bins : array + Metallicity bins edges in absolute metallicity. + + Returns + ------- + ndarray + Fraction of the SFR in the given metallicity bin at the given redshift. + """ + # only use data within the metallicity bounds (no lower bound) + redshift_indices = np.array([np.where(self.redshifts <= i)[0][0] for i in z]) + Z_dist = self.SFR_data[redshift_indices] + fSFR = np.zeros((len(z), len(metallicity_bins) - 1)) + + for i in range(len(z)): + if Z_dist[i].sum() == 0.0: + continue + else: + Z_dist_cdf = np.cumsum(Z_dist[i]) / Z_dist[i].sum() + Z_dist_cdf = np.append(Z_dist_cdf, 1) + Z_x_values = np.append(np.log10(self.Z), 0) + Z_dist_cdf_interp = interp1d(Z_x_values, Z_dist_cdf) + cdf_fun = lambda x: Z_dist_cdf_interp(np.log10(x)) + fSFR[i, :] = self._distribute_cdf(cdf_fun, metallicity_bins) + return fSFR + + def _load_redshift_data(self, verbose=False): # pragma: no cover + """Load the redshift data from a Chruslinsk+21 model file. + + Returns + ------- + time : array + the center of the time bins + redshift : array + the redshifts corresponding to the time bins + delt : array + the width of the time bins + """ + if verbose: + print("Loading redshift data...") + + time, redshift, delt = np.loadtxt( + os.path.join(self._data_folder, "Time_redshift_deltaT.dat"), unpack=True) + return time, redshift, delt + + def _load_raw_data(self): # pragma: no cover + """Read the sub-model data from the file + + The data structure is as follows: + - mass per unit (comoving) volume formed in each z (row) - FOH (column) bin + + Returns + ------- + array + Mass formed per unit volume in each redshift and FOH bin + """ + input_file = os.path.join(self._data_folder, f"{self.sub_model}.dat") + data = np.loadtxt(input_file) + return data + + def CSFRD(self, z): + """Interpolate the cosmic star formation rate density at the given redshift(s) + + Parameters + ---------- + z : float or array-like + Cosmological redshift. + + Returns + ------- + float or array-like + The cosmic star formation rate density at the given redshift(s). + """ + SFR_interp = interp1d(self.redshifts, np.sum(self.SFR_data*self.dFOH, axis=1)) + return SFR_interp(z) + +class Fujimoto24(MadauBase): + """The Fujimoto et al. (2024) star formation history model. + mean metallicity evolution of Madau & Fragos (2017). + + Fujimoto et al. (2024), ApJ SS, 275, 2, 36, 59 + https://ui.adsabs.harvard.edu/abs/2024ApJS..275...36F/abstract + """ + def __init__(self, SFH_MODEL): + """Initialise the Fujimoto+24 model + + Parameters + ---------- + SFH_MODEL : dict + SFH model parameters. Fujimoto+24 requires the following parameters: + - sigma : float or str + The standard deviation of the log-normal metallicity distribution. + Options are: + - Bavera+20 + - Neijssel+19 + - float + Additional SFH model parameters: + - Z_max : float + The maximum metallicity in absolute units. + - Z_min : float + The minimum metallicity in absolute units. + - normalise : bool + Normalise the metallicity distribution to 1. + """ + super().__init__(SFH_MODEL) + # Parameters for Fujimoto+24 CSFRD + self.CSFRD_params = { + "a": 0.010, + "b": 2.8, + "c": 3.3, + "d": 6.6, + } + +class Zavala21(MadauBase): + """The Zavala et al. (2021) star formation history model [1]_. + + The "min" and "max" models are based on the obscured and unobscured + star formation rate density models, respectively. + + References + ---------- + .. [1] Zavala, J., et al. (2021). ApJ, 909(2), 165. + https://ui.adsabs.harvard.edu/abs/2021ApJ...909..165Z/ + """ + def __init__(self, SFH_MODEL): + """Initialise the Zavala+21 model + + Parameters + ---------- + SFH_MODEL : dict + SFH model parameters. Zavala+21 requires the following parameters: + - sub_model : str + The sub-model to use. Either "min" or "max". + Additional SFH model parameters: + - Z_max : float + The maximum metallicity in absolute units. + - Z_min : float + The minimum metallicity in absolute units. + - normalise : bool + Normalise the metallicity distribution to 1. + """ + if "sub_model" not in SFH_MODEL: + raise ValueError("Sub-model not given!") + + super().__init__(SFH_MODEL) + self._load_zavala_data() + + def _load_zavala_data(self): # pragma: no cover + """Load the data from the Zavala+21 models + Transforms the data to the format used in the classes. + + """ + data_file = os.path.join(PATH_TO_POSYDON_DATA, "SFR/Zavala+21.txt") + tmp_data = pd.read_csv(data_file, + names=["redshift", "SFRD_min", "SFRD_max"], + skiprows=1, + sep=r"\s+") + self.redshifts = tmp_data["redshift"].values + # The min / max values originally come from their obscured + # and unobscured SFRD model. + if self.sub_model == "min": + self.SFR_data = tmp_data["SFRD_min"].values + elif self.sub_model == "max": + self.SFR_data = tmp_data["SFRD_max"].values + else: + raise ValueError("Invalid sub-model!") + + # overwrite the CSFRD method of MadauBase + def CSFRD(self, z): + SFR_interp = interp1d(self.redshifts, self.SFR_data) + return SFR_interp(z) + +def get_SFH_model(SFH_MODEL): + """Return the appropriate SFH model based on the given parameters + + Parameters + ---------- + SFH_MODEL : dict + SFH model parameters. + + Returns + ------- + a SFHBase instance or subclass + The SFH model instance. + """ + if SFH_MODEL["SFR"] == "Madau+Fragos17": + return MadauFragos17(SFH_MODEL) + elif SFH_MODEL["SFR"] == "Madau+Dickinson14": + return MadauDickinson14(SFH_MODEL) + elif SFH_MODEL["SFR"] == "Fujimoto+24": + return Fujimoto24(SFH_MODEL) + elif SFH_MODEL["SFR"] == "Neijssel+19": + return Neijssel19(SFH_MODEL) + elif SFH_MODEL["SFR"] == "IllustrisTNG": + return IllustrisTNG(SFH_MODEL) + elif SFH_MODEL["SFR"] == "Chruslinska+21": + return Chruslinska21(SFH_MODEL) + elif SFH_MODEL["SFR"] == "Zavala+21": + return Zavala21(SFH_MODEL) + else: + raise ValueError("Invalid SFR!") + +def SFR_per_met_at_z(z, met_bins, SFH_MODEL): + """Calculate the SFR per metallicity bin at a given redshift(s) + + Parameters + ---------- + z : float or array-like + Cosmological redshift. + met_bins : array + Metallicity bins edges in absolute metallicity. + SFH_MODEL : dict + SFH model parameters. + + Returns + ------- + SFH : 2D array + Star formation history per metallicity bin at the given redshift(s). + + """ + SFH = get_SFH_model(SFH_MODEL) + return SFH(z, met_bins) + +# TODO: No testing coverage for the following function, but should be added +def get_formation_times(N_binaries, star_formation="constant", **kwargs): # pragma: no cover """Get formation times of binaries in a population based on a SFH scenario. Parameters @@ -94,263 +1099,4 @@ def get_formation_times(N_binaries, star_formation="constant", **kwargs): "Unknown star formation scenario '{}' given. Valid options: {}".format( star_formation, ",".join(SFH_SCENARIOS) ) - ) - - -def get_illustrisTNG_data(verbose=False): - """Load IllustrisTNG SFR dataset.""" - if verbose: - print("Loading IllustrisTNG data...") - return np.load(os.path.join(PATH_TO_POSYDON_DATA, "SFR/IllustrisTNG.npz")) - - -def star_formation_rate(SFR, z): - """Star formation rate in M_sun yr^-1 Mpc^-3. - - Parameters - ---------- - SFR : string - Star formation rate assumption: - - Madau+Fragos17 see arXiv:1606.07887 - - Madau+Dickinson14 see arXiv:1403.0007 - - Neijssel+19 see arXiv:1906.08136 - - IllustrisTNG see see arXiv:1707.03395 - z : double - Cosmological redshift. - - Returns - ------- - double - The total mass of stars in M_sun formed per comoving volume Mpc^-3 - per year. - """ - if SFR == "Madau+Fragos17": - return ( - 0.01 * (1.0 + z) ** 2.6 / (1.0 + ((1.0 + z) / 3.2) ** 6.2) - ) # M_sun yr^-1 Mpc^-3 - elif SFR == "Madau+Dickinson14": - return ( - 0.015 * (1.0 + z) ** 2.7 / (1.0 + ((1.0 + z) / 2.9) ** 5.6) - ) # M_sun yr^-1 Mpc^-3 - elif SFR == "Neijssel+19": - return ( - 0.01 * (1.0 + z) ** 2.77 / (1.0 + ((1.0 + z) / 2.9) ** 4.7) - ) # M_sun yr^-1 Mpc^-3 - elif SFR == "IllustrisTNG": - illustris_data = get_illustrisTNG_data() - SFR = illustris_data["SFR"] # M_sun yr^-1 Mpc^-3 - redshifts = illustris_data["redshifts"] - SFR_interp = interp1d(redshifts, SFR) - return SFR_interp(z) - else: - raise ValueError("Invalid SFR!") - - -def mean_metallicity(SFR, z): - """Empiric mean metallicity function. - - Parameters - ---------- - SFR : string - Star formation rate assumption: - - Madau+Fragos17 see arXiv:1606.07887 - - Madau+Dickinson14 see arXiv:1403.0007 - - Neijssel+19 see arXiv:1906.08136 - z : double - Cosmological redshift. - - Returns - ------- - double - Mean metallicty of the universe at the given redhist. - - """ - - if SFR == "Madau+Fragos17" or SFR == "Madau+Dickinson14": - return 10 ** (0.153 - 0.074 * z**1.34) * Zsun - elif SFR == "Neijssel+19": - return 0.035 * 10 ** (-0.23 * z) - else: - raise ValueError("Invalid SFR!") - - -def std_log_metallicity_dist(sigma): - """Standard deviation of the log-metallicity distribution. - - Returns - ------- - double - Standard deviation of the adopted distribution. - - """ - if isinstance(sigma, str): - if sigma == "Bavera+20": - return 0.5 - elif sigma == "Neijssel+19": - return 0.39 - else: - raise ValueError("Uknown sigma choice!") - elif isinstance(sigma, float): - return sigma - else: - raise ValueError(f"Invalid sigma value {sigma}!") - - -def SFR_Z_fraction_at_given_redshift( - z, SFR, sigma, metallicity_bins, Z_max, select_one_met -): - """'Fraction of the SFR at a given redshift z in a given metallicity bin as in Eq. (B.8) of Bavera et al. (2020). - - Parameters - ---------- - z : np.array - Cosmological redshift. - SFR : string - Star formation rate assumption: - - Madau+Fragos17 see arXiv:1606.07887 - - Madau+Dickinson14 see arXiv:1403.0007 - - IllustrisTNG see see arXiv:1707.03395 - - Neijssel+19 see arXiv:1906.08136 - sigma : double / string - Standard deviation of the log-metallicity distribution. - If string, it can be 'Bavera+20' or 'Neijssel+19'. - metallicity_bins : array - Metallicity bins edges in absolute metallicity. - Z_max : double - Maximum metallicity in absolute metallicity. - select_one_met : bool - If True, the function returns the fraction of the SFR in the given metallicity bin. - If False, the function returns the fraction of the SFR in the given metallicity bin and the fraction of the SFR in the metallicity bin - - Returns - ------- - array - Fraction of the SFR in the given metallicity bin at the given redshift. - In absolute metallicity. - """ - - if SFR == "Madau+Fragos17" or SFR == "Madau+Dickinson14": - sigma = std_log_metallicity_dist(sigma) - mu = np.log10(mean_metallicity(SFR, z)) - sigma**2 * np.log(10) / 2.0 - # renormalisation constant. We can use mu[0], since we integrate over the whole metallicity range - norm = stats.norm.cdf(np.log10(Z_max), mu[0], sigma) - fSFR = np.empty((len(z), len(metallicity_bins) - 1)) - fSFR[:, :] = np.array( - [(stats.norm.cdf(np.log10(metallicity_bins[1:]), m, sigma) / norm - - stats.norm.cdf(np.log10(metallicity_bins[:-1]), m, sigma) / norm - ) for m in mu ] - ) - if not select_one_met: - fSFR[:, 0] = stats.norm.cdf(np.log10(metallicity_bins[1]), mu, sigma) / norm - fSFR[:,-1] = norm - stats.norm.cdf(np.log10(metallicity_bins[-1]), mu, sigma)/norm - - elif SFR == "Neijssel+19": - # assume a truncated ln-normal distribution of metallicities - sigma = std_log_metallicity_dist(sigma) - mu = np.log(mean_metallicity(SFR, z)) - sigma**2 / 2.0 - # renormalisation constant - norm = stats.norm.cdf(np.log(Z_max), mu[0], sigma) - fSFR = np.empty((len(z), len(metallicity_bins) - 1)) - fSFR[:, :] = np.array( - [ - ( - stats.norm.cdf(np.log(metallicity_bins[1:]), m, sigma) / norm - - stats.norm.cdf(np.log(metallicity_bins[:-1]), m, sigma) / norm - ) - for m in mu - ] - ) - if not select_one_met: - fSFR[:, 0] = stats.norm.cdf(np.log(metallicity_bins[1]), mu, sigma) / norm - fSFR[:,-1] = norm - stats.norm.cdf(np.log(metallicity_bins[-1]), mu, sigma)/norm - - elif SFR == "IllustrisTNG": - # numerically itegrate the IlluystrisTNG SFR(z,Z) - illustris_data = get_illustrisTNG_data() - redshifts = illustris_data["redshifts"] - Z = illustris_data["mets"] - M = illustris_data["M"] # Msun - # only use data within the metallicity bounds (no lower bound) - Z_max_mask = Z <= Z_max - redshift_indices = np.array([np.where(redshifts <= i)[0][0] for i in z]) - Z_dist = M[:, Z_max_mask][redshift_indices] - fSFR = np.zeros((len(z), len(metallicity_bins) - 1)) - - for i in range(len(z)): - if Z_dist[i].sum() == 0.0: - continue - else: - # We add a final point to the CDF and metallicities to ensure normalisation to 1 - Z_dist_cdf = np.cumsum(Z_dist[i]) / Z_dist[i].sum() - Z_dist_cdf = np.append(Z_dist_cdf, 1) - Z_x_values = np.append(np.log10(Z[Z_max_mask]), 0) - Z_dist_cdf_interp = interp1d(Z_x_values, Z_dist_cdf) - - fSFR[i, :] = (Z_dist_cdf_interp(np.log10(metallicity_bins[1:])) - - Z_dist_cdf_interp(np.log10(metallicity_bins[:-1]))) - - if not select_one_met: - # add the fraction of the SFR in the first and last bin - # or the only bin without selecting one metallicity - if len(metallicity_bins) == 2: - fSFR[i, 0] = 1 - else: - fSFR[i, 0] = Z_dist_cdf_interp(np.log10(metallicity_bins[1])) - fSFR[i, -1] = 1 - Z_dist_cdf_interp(np.log10(metallicity_bins[-1])) - else: - raise ValueError("Invalid SFR!") - - return fSFR - - -def integrated_SFRH_over_redshift(SFR, sigma, Z, Z_max): - """Integrated SFR history over z as in Eq. (B.10) of Bavera et al. (2020). - - Parameters - ---------- - SFR : string - Star formation rate assumption: - - Madau+Fragos17 see arXiv:1606.07887 - - Madau+Dickinson14 see arXiv:1403.0007 - - Neijssel+19 see arXiv:1906.08136 - Z : double - Metallicity. - - Returns - ------- - double - The total mass of stars formed per comoving volume at a given - metallicity Z. - - """ - - def E(z, Omega_m=cosmology.Om0): - Omega_L = 1.0 - Omega_m - return (Omega_m * (1.0 + z) ** 3 + Omega_L) ** (1.0 / 2.0) - - def f(z, Z): - if SFR == "Madau+Fragos17" or SFR == "Madau+Dickinson14": - sigma = std_log_metallicity_dist(sigma) - mu = np.log10(mean_metallicity(SFR, z)) - sigma**2 * np.log(10) / 2.0 - H_0 = cosmology.H0.to("1/yr").value # yr - # put a cutoff on metallicity at Z_max - norm = stats.norm.cdf(np.log10(Z_max), mu, sigma) - return ( - star_formation_rate(SFR, z) - * stats.norm.pdf(np.log10(Z), mu, sigma) - / norm - * (H_0 * (1.0 + z) * E(z)) ** (-1) - ) - elif SFR == "Neijssel+19": - sigma = std_log_metallicity_dist(sigma) - mu = np.log10(mean_metallicity(SFR, z)) - sigma**2 / 2.0 - H_0 = cosmology.H0.to("1/yr").value # yr - return ( - star_formation_rate(SFR, z) - * stats.norm.pdf(np.log(Z), mu, sigma) - * (H_0 * (1.0 + z) * E(z)) ** (-1) - ) - else: - raise ValueError("Invalid SFR!") - - return sp.integrate.quad(f, 1e-10, np.inf, args=(Z,))[0] # M_sun yr^-1 Mpc^-3 + ) \ No newline at end of file diff --git a/posydon/popsyn/synthetic_population.py b/posydon/popsyn/synthetic_population.py index bcc6fb1425..6d5e190ddb 100644 --- a/posydon/popsyn/synthetic_population.py +++ b/posydon/popsyn/synthetic_population.py @@ -52,15 +52,12 @@ get_comoving_distance_from_redshift, get_cosmic_time_from_redshift, redshift_from_cosmic_time_interpolator, - DEFAULT_MODEL, + DEFAULT_SFH_MODEL, get_redshift_bin_edges, get_redshift_bin_centers, ) -from posydon.popsyn.star_formation_history import ( - star_formation_rate, - SFR_Z_fraction_at_given_redshift, -) +from posydon.popsyn.star_formation_history import SFR_per_met_at_z from posydon.popsyn.binarypopulation import ( BinaryPopulation, @@ -2030,20 +2027,16 @@ def calculate_cosmic_weights(self, SFH_identifier, MODEL_in=None): Examples -------- >>> transient_population = TransientPopulation('filename.h5', 'transient_name') - >>> transient_population.calculate_cosmic_weights('IllustrisTNG', MODEL_in=DEFAULT_MODEL) + >>> transient_population.calculate_cosmic_weights('IllustrisTNG', MODEL_in=DEFAULT_SFH_MODEL) """ # Set model to DEFAULT or provided MODEL parameters # Allows for partial model specification if MODEL_in is None: - MODEL = DEFAULT_MODEL + MODEL = DEFAULT_SFH_MODEL else: - for key in MODEL_in: - if key not in DEFAULT_MODEL: - raise ValueError(key + " is not a valid parameter name!") - - # write the DEFAULT_MODEL with updates parameters to self.MODEL. - MODEL = DEFAULT_MODEL + # write the DEFAULT_SFH_MODEL with updated parameters to MODEL. + MODEL = DEFAULT_SFH_MODEL MODEL.update(MODEL_in) path_in_file = ( @@ -2083,23 +2076,12 @@ def calculate_cosmic_weights(self, SFH_identifier, MODEL_in=None): get_redshift_from_cosmic_time = redshift_from_cosmic_time_interpolator() indices = self.indices + met_edges = rates.edges_metallicity_bins # sample the SFH for only the events that are within the Hubble time # only need to sample the SFH at each metallicity and z_birth - # Not for every event! - SFR_at_z_birth = star_formation_rate(rates.MODEL["SFR"], z_birth) - # get metallicity bin edges - met_edges = rates.edges_metallicity_bins - - # get the fractional SFR at each metallicity and z_birth - fSFR = SFR_Z_fraction_at_given_redshift( - z_birth, - rates.MODEL["SFR"], - rates.MODEL["sigma_SFR"], - met_edges, - rates.MODEL["Z_max"], - rates.MODEL["select_one_met"], - ) - + # Not for every event! + SFR_per_met_at_z_birth = SFR_per_met_at_z(z_birth, met_edges, rates.MODEL) + # simulated mass per given metallicity corrected for the unmodeled # single and binary stellar mass M_model = rates.mass_per_metallicity.loc[rates.centers_metallicity_bins / Zsun][ @@ -2153,7 +2135,7 @@ def calculate_cosmic_weights(self, SFH_identifier, MODEL_in=None): ) weights = np.zeros((len(met_events), nr_of_birth_bins)) - for i, met in enumerate(rates.centers_metallicity_bins): + for j, met in enumerate(rates.centers_metallicity_bins): mask = met_events == met weights[mask, :] = ( 4.0 @@ -2161,8 +2143,8 @@ def calculate_cosmic_weights(self, SFH_identifier, MODEL_in=None): * c * D_c[mask] ** 2 * deltaT - * (fSFR[:, i] * SFR_at_z_birth) - / M_model[i] + * SFR_per_met_at_z_birth[:, j] + / M_model[j] ) # yr^-1 with pd.HDFStore(self.filename, mode="a") as store: diff --git a/posydon/unit_tests/popsyn/test_star_formation_history.py b/posydon/unit_tests/popsyn/test_star_formation_history.py new file mode 100644 index 0000000000..8c5514ef27 --- /dev/null +++ b/posydon/unit_tests/popsyn/test_star_formation_history.py @@ -0,0 +1,1069 @@ +"""Unit tests for posydon/popsyn/star_formation_history.py +""" + +__authors__ = [ + "Max Briel ", +] + +import numpy as np +import pandas as pd +import pytest +from posydon.utils.posydonwarning import SFHModelWarning +from posydon.popsyn.star_formation_history import SFHBase, MadauBase +from posydon.popsyn.star_formation_history import ( + MadauDickinson14, + MadauFragos17, + Neijssel19, + Fujimoto24, + IllustrisTNG, + Chruslinska21, + Zavala21, + get_SFH_model +) + + +class TestSFHBase: + + @pytest.fixture + def ConcreteSFH(self): + """Create a concrete subclass of SFHBase for testing.""" + class ConcreteSFH(SFHBase): + def CSFRD(self, z): + return z + + def mean_metallicity(self, z): + return z + + def fSFR(self, z, metallicity_bins): + return np.ones((len(z), len(metallicity_bins)-1)) + + return ConcreteSFH + + def test_init_attributes(self, ConcreteSFH): + """Test that the initialization sets attributes correctly.""" + model_dict = { + "test_param": 42, + "Z_max": 0.03, + "another_param": "test" + } + sfh = ConcreteSFH(model_dict) + + # Check that attributes are set correctly + for key, value in model_dict.items(): + assert getattr(sfh, key) == value + + # additional SFH_model set check + assert sfh.SFH_MODEL == model_dict + + + @pytest.mark.parametrize("model_dict, error_msg", [ + # Z_max + ({"Z_max": 1.5}, "Z_max must be in absolute units! It cannot be larger than 1!"), + ({"Z_max": -0.1}, "Z_max must be in absolute units! It cannot be negative!"), + # Z_min + ({"Z_min": -0.1}, "Z_min must be in absolute units! It cannot be negative!"), + ({"Z_min": 1.2}, "Z_min must be in absolute units! It cannot be larger than 1!"), + # Z_min > Z_max + ({"Z_max": 0.1, "Z_min": 0.2}, "Z_min must be smaller than Z_max!"), + ]) + def test_validation(self, ConcreteSFH, model_dict, error_msg): + with pytest.raises(ValueError) as excinfo: + sfh = ConcreteSFH(model_dict) + assert error_msg in str(excinfo.value) + + def test_abstract_methods(self): + """Test that abstract methods must be implemented.""" + # Create incomplete subclasses that don't implement all abstract methods + class IncompleteSFH1(SFHBase): + def CSFRD(self, z): + return z + + class IncompleteSFH2(SFHBase): + def CSFRD(self, z): + return z + def mean_metallicity(self, z): + return z + + model_dict = {"Z_max": 0.03} + with pytest.raises(TypeError) as excinfo: + IncompleteSFH1(model_dict) + assert ("Can't instantiate abstract class IncompleteSFH1 " + "with abstract methods fSFR, mean_metallicity") in str(excinfo.value) + + with pytest.raises(TypeError) as excinfo: + IncompleteSFH2(model_dict) + assert ("Can't instantiate abstract class IncompleteSFH2 " + "with abstract method fSFR") in str(excinfo.value) + + @pytest.mark.parametrize( + "model_dict, normalise, met_edges, expected, warning", + [ + # Simple CDF with Z_max=1, Z_min=0.0 + ({"Z_max": 1, "Z_min": 0.0}, False, np.array([0.0, 0.01, 0.02, 0.03]), + np.array([0.01, 0.01, 0.98]), None), + # No Z_min/Z_max set + ({}, False, np.array([0.0, 0.01, 0.02, 0.03]), + 0.01 * np.ones(3), None), + # Model dict warning + ({"Z_max": 0.02, "Z_min": 0.0}, False, np.array([0.0, 0.01, 0.02, 0.03]), + None, 'Z_max is smaller than the highest metallicity bin.'), + # Different model dicts + ({"Z_max": 1, "Z_min": 0.015}, False, np.array([0.3, 0.6, 0.9]), + np.array([0.585, 0.4]), None), + # With normalization + ({"Z_max": 1, "Z_min": 0.015}, True, np.array([0.3, 0.6, 0.9]), + None, None), + # Restrict upper bound + ({"Z_max": 0.95, "Z_min": 0.2}, False, np.array([0.3, 0.6, 0.9]), + np.array([0.4, 0.35]), None), + # With normalization + ({"Z_max": 0.95, "Z_min": 0.2}, True, np.array([0.3, 0.6, 0.9]), + None, None), + # Minimum in lowest bin + ({"Z_min": 0.25}, False, np.array([0.2, 0.3, 0.6, 0.9]), + np.array([0.05, 0.3, 0.3]), + 'Z_min is larger than the lowest metallicity bin.'), + # Minimum higher than minimum bin + ({"Z_min": 0.35}, False, np.array([0.2, 0.3, 0.6, 0.9]), + np.array([0.0, 0.25, 0.3]), + 'Z_min is larger than the lowest metallicity bin.'), + # Minimum in lowest bin and maximum + ({"Z_min": 0.25, "Z_max": 0.8}, False, np.array([0.2, 0.3, 0.6, 0.9]), + np.array([0.05, 0.3, 0.2]), + 'Z_max is smaller than the highest metallicity bin.'), + # Minimum higher than minimum bin, narrow range + ({"Z_min": 0.35, "Z_max": 0.4}, False, np.array([0.2, 0.3, 0.6, 0.9]), + np.array([0.0, 0.05, 0.0]), + 'Z_max is smaller than the highest metallicity bin.'), + # Minimum higher than minimum bin, medium range + ({"Z_min": 0.35, "Z_max": 0.65}, False, np.array([0.2, 0.3, 0.6, 0.9]), + np.array([0.0, 0.25, 0.05]), + 'Z_max is smaller than the highest metallicity bin.'), + ]) + def test_distribute_cdf(self, ConcreteSFH, model_dict, normalise, met_edges, expected, warning): + """Test the _distribute_cdf method with various scenarios.""" + # Create a simple CDF function + cdf_func = lambda x: x + + # Create the SFH instance + sfh = ConcreteSFH(model_dict) + sfh.normalise = normalise + + # Test execution with or without warning check + if warning is not None: + with pytest.warns(SFHModelWarning, match=warning): + result = sfh._distribute_cdf(cdf_func, met_edges) + else: + result = sfh._distribute_cdf(cdf_func, met_edges) + + # Check results + if normalise: + # For normalise=True, check sum is 1.0 + np.testing.assert_allclose(np.sum(result), 1.0) + elif expected is not None: + # For specific expected values + np.testing.assert_allclose(result, expected) + + def test_distribute_cdf_invalid(self, ConcreteSFH): + """Test the _distribute_cdf method with invalid inputs.""" + # Create a simple CDF function + cdf_func = lambda x: x + + # Create the SFH instance + model_dict = {} + sfh = ConcreteSFH(model_dict) + + # Test with invalid metallicity edges + met_edges = np.array([0.04, 0.01, 0.02]) + with pytest.raises(ValueError) as excinfo: + sfh._distribute_cdf(cdf_func, met_edges) + assert "Metallicity bins must be sorted" in str(excinfo.value) + + def test_call_method(self): + """Test the __call__ method.""" + class ConcreteSFH(SFHBase): + def CSFRD(self, z): + return z*2.0 + + # just a placeholder. Doesn't contribute + def mean_metallicity(self, z): + return z*0.01 + + def fSFR(self, z, metallicity_bins): + # Return a simple array for testing + delta = np.diff(metallicity_bins) + # normalise + delta /= delta.sum() + + out = np.zeros((len(z), len(delta))) + out[:, :] = delta + return out + + model_dict = {"Z_max": 0.03} + sfh = ConcreteSFH(model_dict) + + z = np.array([0.5, 1.0]) + met_edges = np.array([0.0, 0.02, 0.06]) + + result = sfh(z, met_edges) + + # Expected: CSFRD(z)[:, np.newaxis] * fSFR(z, met_edges) + expected = np.array([ + [1.0 * 1/3, 1.0 * 2/3], + [2.0 * 1/3, 2.0 * 2/3] + ]) + + np.testing.assert_allclose(result, expected) + + +class TestMadauBase: + """Test class for MadauBase""" + + @pytest.fixture + def e_CSFRD_params(self): + return {"a": 0.01, + "b": 2.6, + "c": 3.2, + "d": 6.2, + } + + @pytest.fixture + def ConcreteMadau(self, e_CSFRD_params): + class ConcreteMadau(MadauBase): + """Concrete subclass of MadauBase for testing""" + def __init__(self, MODEL): + super().__init__(MODEL) + self.CSFRD_params = e_CSFRD_params + return ConcreteMadau + + def test_init_requires_sigma(self, ConcreteMadau): + """Test that MadauBase requires a sigma parameter""" + model_dict = {"Z_max": 0.03} + with pytest.raises(ValueError) as excinfo: + ConcreteMadau(model_dict) + assert "sigma not given!" in str(excinfo.value) + + def test_init_sets_csfrd_params_to_none(self, ConcreteMadau, e_CSFRD_params): + """Test that CSFRD_params is not set to None initially""" + model_dict = {"sigma": 0.5} + madau = ConcreteMadau(model_dict) + assert madau.CSFRD_params is not None + assert madau.CSFRD_params["a"] == e_CSFRD_params["a"] + assert madau.CSFRD_params["b"] == e_CSFRD_params["b"] + assert madau.CSFRD_params["c"] == e_CSFRD_params["c"] + assert madau.CSFRD_params["d"] == e_CSFRD_params["d"] + + def test_std_log_metallicity_dist(self, ConcreteMadau): + """Test the std_log_metallicity_dist method with different sigma values""" + # Test Bavera+20 + model_dict = {"sigma": "Bavera+20"} + madau = ConcreteMadau(model_dict) + assert madau.std_log_metallicity_dist() == 0.5 + + # Test Neijssel+19 + model_dict = {"sigma": "Neijssel+19"} + madau = ConcreteMadau(model_dict) + assert madau.std_log_metallicity_dist() == 0.39 + + # Test float value + model_dict = {"sigma": 0.45} + madau = ConcreteMadau(model_dict) + assert madau.std_log_metallicity_dist() == 0.45 + + # Test unknown string + model_dict = {"sigma": "unknown"} + madau = ConcreteMadau(model_dict) + with pytest.raises(ValueError) as excinfo: + madau.std_log_metallicity_dist() + assert "Unknown sigma choice!" in str(excinfo.value) + + # Test integer type + model_dict = {"sigma": 1} + madau = ConcreteMadau(model_dict) + assert madau.std_log_metallicity_dist() == 1.0 + + # Test invalid type + model_dict = {"sigma": [0.5, 0.6]} + madau = ConcreteMadau(model_dict) + with pytest.raises(ValueError) as excinfo: + madau.std_log_metallicity_dist() + assert "Invalid sigma value" in str(excinfo.value) + + def test_csfrd(self, ConcreteMadau, e_CSFRD_params): + """Test the CSFRD method""" + model_dict = {"sigma": 0.5} + madau = ConcreteMadau(model_dict) + + def tmp_CSFRD(z): + return (e_CSFRD_params["a"] * ((1 + z)**e_CSFRD_params["b"]) + / (1 + ((1 + z)/e_CSFRD_params["c"])**e_CSFRD_params["d"])) + + # Test with single value + z = 0.0 + result = madau.CSFRD(z) + # a * (1+z)^b / (1 + ((1+z)/c)^d) with z=0 + expected = tmp_CSFRD(z) + np.testing.assert_allclose(result, expected) + + # Test with array of values + z_array = np.array([0.0, 1.0, 2.0]) + result = madau.CSFRD(z_array) + expected = np.array([ + tmp_CSFRD(z) for z in z_array + ]) + np.testing.assert_allclose(result, expected) + + def test_mean_metallicity(self, ConcreteMadau): + """Test the mean_metallicity method""" + from posydon.utils.constants import Zsun + + model_dict = {"sigma": 0.5, "Z_max": 0.03} + madau = ConcreteMadau(model_dict) + + # Test with single value + z = 0.0 + result = madau.mean_metallicity(z) + expected = 10**(0.153) * Zsun + np.testing.assert_allclose(result, expected) + + # Test with array of values + z_array = np.array([0.0, 1.0, 2.0]) + result = madau.mean_metallicity(z_array) + expected = 10**(0.153 - 0.074 * z_array**1.34) * Zsun + np.testing.assert_allclose(result, expected) + + def test_fsfr(self, ConcreteMadau): + """Test the fSFR method""" + model_dict = {"sigma": 0.5, "Z_max": 1} + madau = ConcreteMadau(model_dict) + + # Test with redshift array and metallicity bins + z = np.array([0.0, 1.0]) + met_bins = np.array([0.001, 0.01, 0.02, 0.03]) + + result = madau.fSFR(z, met_bins) + + # Shape check - should be (len(z), len(met_bins)-1) + assert result.shape == (2, 3) + + # THIS IS A VALIDATION TEST + from scipy.stats import norm + mean0, mean1 = (np.log10(madau.mean_metallicity(z)) + - model_dict['sigma']**2 * np.log(10) / 2) + + expected1 = np.array( + # integral from 0.001 to 0.01; Z_min = lowest bin edge + [(norm.cdf(np.log10(0.01), mean0, model_dict['sigma']) + - norm.cdf(np.log10(0.001), mean0, model_dict['sigma'])), + # integral from 0.01 to 0.02 + (norm.cdf(np.log10(0.02), mean0, model_dict['sigma']) + - norm.cdf(np.log10(0.01), mean0, model_dict['sigma'])), + # integral from 0.02 to 1 + (norm.cdf(np.log10(1), mean0, model_dict['sigma']) + - norm.cdf(np.log10(0.02), mean0, model_dict['sigma']))],) + + expected2 = np.array( + # integral from 0.001 to 0.01;Z_min = lowest bin edge + [(norm.cdf(np.log10(0.01), mean1, model_dict['sigma']) + - norm.cdf(np.log10(0.001), mean1, model_dict['sigma'])), + # integral from 0.01 to 0.02 + (norm.cdf(np.log10(0.02), mean1, model_dict['sigma']) + - norm.cdf(np.log10(0.01), mean1, model_dict['sigma'])), + # integral from 0.02 to 1 + (norm.cdf(np.log10(1), mean1, model_dict['sigma']) + - norm.cdf(np.log10(0.02), mean1, model_dict['sigma']))], + ) + expected = np.array([expected1, expected2]) + np.testing.assert_allclose(result, expected) + + # Change Z_min to a very small number to include the rest of the lowest mets + model_dict = {"sigma": 0.5, "Z_max": 0.3, "Z_min": 1e-11} + madau = ConcreteMadau(model_dict) + result = madau.fSFR(z, met_bins) + + expected1 = np.array( + [norm.cdf(np.log10(0.01), mean0, model_dict['sigma']), + (norm.cdf(np.log10(0.02), mean0, model_dict['sigma']) + - norm.cdf(np.log10(0.01), mean0, model_dict['sigma'])), + (norm.cdf(np.log10(0.3), mean0, model_dict['sigma']) + - norm.cdf(np.log10(0.02), mean0, model_dict['sigma']))],) + expected2 = np.array( + [norm.cdf(np.log10(0.01), mean1, model_dict['sigma']), + (norm.cdf(np.log10(0.02), mean1, model_dict['sigma']) + - norm.cdf(np.log10(0.01), mean1, model_dict['sigma'])), + (norm.cdf(np.log10(0.3), mean1, model_dict['sigma']) + - norm.cdf(np.log10(0.02), mean1, model_dict['sigma']))],) + expected = np.array([expected1, expected2]) + np.testing.assert_allclose(result, expected) + + # Test with normalise + model_dict = {"sigma": 0.5, "Z_max": 0.3, "Z_min": 1e-11, "normalise": True} + madau = ConcreteMadau(model_dict) + result = madau.fSFR(z, met_bins) + expected = np.ones(len(z)) + np.testing.assert_allclose(np.sum(result, axis=1), expected) + + # Test with Z_min > met_bins[1] + model_dict = {"sigma": 0.5, + "Z_max": 0.3, + "Z_min": 0.02, + "normalise": True} + madau = ConcreteMadau(model_dict) + warning_str = "Z_min is larger than the lowest metallicity bin." + with pytest.warns(SFHModelWarning, match=warning_str): + result = madau.fSFR(z, met_bins) + #result = madau.fSFR(z, met_bins) + expected = np.ones(len(z)) + np.testing.assert_allclose(np.sum(result, axis=1), expected) + +class TestIllustrisTNG: + """Tests for the IllustrisTNG SFH model with mocked data loading.""" + + @pytest.fixture + def mock_illustris_data(self): + """Create mock data for the IllustrisTNG class.""" + # Create mock data structure similar to the npz file + num_redshifts = 10 + num_metallicities = 5 + + mock_data = { + "SFR": np.linspace(0.1, 1.0, num_redshifts)[::-1], # SFR decreases with redshift + "redshifts": np.linspace(0.0, 9.0, num_redshifts)[::-1], # Redshifts from 0 to 9 + "mets": np.logspace(-4, -1, num_metallicities), # Metallicities from 1e-4 to 1e-1 + "M": np.ones((num_redshifts, num_metallicities)) # Equal mass in all bins for simplicity + } + + # Add some variation to mass distribution for testing mean_metallicity + for i in range(num_redshifts): + # Linear decrease in higher metallicities as redshift increases + scale = 1.0 - i / num_redshifts + mock_data["M"][i] = np.linspace(1.0, scale, num_metallicities) + + mock_data["M"] = np.flip(mock_data["M"], axis=0) # Reverse the mass array + return mock_data + + @pytest.fixture + def illustris_model(self, monkeypatch, mock_illustris_data): + """Create an IllustrisTNG model instance with mocked data.""" + # Create a function that returns the mock data + def mock_get_illustrisTNG_data(self, verbose=False): + return mock_illustris_data + + # Patch the _get_illustrisTNG_data method + monkeypatch.setattr(IllustrisTNG, "_get_illustrisTNG_data", mock_get_illustrisTNG_data) + + # Create and return the model + model_dict = {"Z_max": 0.3} + return IllustrisTNG(model_dict) + + def test_init_parameters(self, illustris_model, mock_illustris_data): + """Test that initialization sets the parameters correctly.""" + # Check that data was loaded correctly + np.testing.assert_array_equal(illustris_model.CSFRD_data, np.flip(mock_illustris_data["SFR"])) + np.testing.assert_array_equal(illustris_model.redshifts, np.flip(mock_illustris_data["redshifts"])) + np.testing.assert_array_equal(illustris_model.Z, mock_illustris_data["mets"]) + np.testing.assert_array_equal(illustris_model.M, np.flip(mock_illustris_data["M"], axis=0)) + + # Check that model parameters were set correctly + assert illustris_model.Z_max == 0.3 + + def test_csfrd_calculation(self, illustris_model, mock_illustris_data): + """Test the CSFRD method.""" + # Test at specific redshifts including boundary values + z_values = np.array([0.0, 4.5, 9.0]) + result = illustris_model.CSFRD(z_values) + + # Expected values come from interpolating flipped SFR data + flipped_sfr = np.flip(mock_illustris_data["SFR"]) + flipped_redshifts = np.flip(mock_illustris_data["redshifts"]) + expected = np.interp(z_values, flipped_redshifts, flipped_sfr) + + np.testing.assert_allclose(result, expected) + + def test_mean_metallicity(self, illustris_model, mock_illustris_data): + """Test the mean_metallicity method.""" + # Test at specific redshifts + z_values = np.array([0.0, 4.5, 9.0]) + result = illustris_model.mean_metallicity(z_values) + + # Calculate expected values manually + flipped_redshifts = np.flip(mock_illustris_data["redshifts"]) + flipped_masses = np.flip(mock_illustris_data["M"], axis=0) + metallicities = mock_illustris_data["mets"] + + # Calculate expected mean metallicities at each test redshift + out = np.zeros_like(flipped_redshifts) + for i, z in enumerate(out): + weights = flipped_masses[i, :] + if np.sum(weights) == 0: + out[i] = np.nan + else: + out[i] = np.average(metallicities, weights=weights) + Z_interp = np.interp(z_values, flipped_redshifts, out) + np.testing.assert_allclose(result, Z_interp) + + # Test empty mass array + illustris_model.M[0] = np.zeros_like(flipped_masses[0]) + with pytest.raises(AssertionError): + result = illustris_model.mean_metallicity(z_values) + + def test_fsfr_calculation(self, illustris_model): + """Test the fSFR method.""" + # Test with redshift array and metallicity bins + z = np.array([0.0, 4.5]) + met_bins = np.array([0.001, 0.01, 0.05, 0.1]) + + result = illustris_model.fSFR(z, met_bins) + + # Shape check - should be (len(z), len(met_bins)-1) + assert result.shape == (2, 3) + + # Test with normalise=True + illustris_model.normalise = True + result = illustris_model.fSFR(z, met_bins) + for row in result: + if np.sum(row) > 0: + np.testing.assert_allclose(np.sum(row), 1.0) + + # Test for Z_dist[i].sum = 0 + # Force the first mass array to be all zeros + illustris_model.M[0] = np.zeros_like(illustris_model.M[0]) + result = illustris_model.fSFR(z, met_bins) + np.testing.assert_allclose(result[0], np.zeros_like(result[0])) + +class TestMadauDickinson14: + """Tests for the MadauDickinson14 SFH model""" + + def test_init_parameters(self): + """Test that initialization sets the correct CSFRD parameters""" + model_dict = {"sigma": 0.5, "Z_max": 0.03} + madau = MadauDickinson14(model_dict) + + # Check that CSFRD_params were set correctly + assert madau.CSFRD_params["a"] == 0.015 + assert madau.CSFRD_params["b"] == 2.7 + assert madau.CSFRD_params["c"] == 2.9 + assert madau.CSFRD_params["d"] == 5.6 + + # Check that it inherits correctly from MadauBase + assert isinstance(madau, MadauBase) + + def test_csfrd_calculation(self): + """Test that CSFRD calculations match expected values""" + model_dict = {"sigma": 0.5, "Z_max": 0.03} + madau = MadauDickinson14(model_dict) + + # Test at specific redshifts + z_values = np.array([0.0, 1.0, 2.0, 6.0]) + result = madau.CSFRD(z_values) + + # Calculate expected values manually + p = madau.CSFRD_params + expected = (p["a"] * (1.0 + z_values) ** p["b"] + / (1.0 + ((1.0 + z_values) / p["c"]) ** p["d"])) + + np.testing.assert_allclose(result, expected) + +class TestMadauFragos17: + """Tests for the MadauFragos17 SFH model""" + + def test_init_parameters(self): + """Test that initialization sets the correct CSFRD parameters""" + model_dict = {"sigma": 0.5, "Z_max": 0.03} + madau = MadauFragos17(model_dict) + + # Check that CSFRD_params were set correctly + assert madau.CSFRD_params["a"] == 0.01 + assert madau.CSFRD_params["b"] == 2.6 + assert madau.CSFRD_params["c"] == 3.2 + assert madau.CSFRD_params["d"] == 6.2 + + # Check that it inherits correctly from MadauBase + assert isinstance(madau, MadauBase) + +class TestNeijssel19: + """Tests for the Neijssel19 SFH model""" + + def test_init_parameters(self): + """Test that initialization sets the correct CSFRD parameters""" + model_dict = {"sigma": 0.5, "Z_max": 0.03} + neijssel = Neijssel19(model_dict) + + # Check that CSFRD_params were set correctly + assert neijssel.CSFRD_params["a"] == 0.01 + assert neijssel.CSFRD_params["b"] == 2.77 + assert neijssel.CSFRD_params["c"] == 2.9 + assert neijssel.CSFRD_params["d"] == 4.7 + + # Check that it inherits correctly from MadauBase + assert isinstance(neijssel, MadauBase) + + def test_mean_metallicity(self): + """Test the overridden mean_metallicity method""" + model_dict = {"sigma": 0.5, "Z_max": 0.03} + neijssel = Neijssel19(model_dict) + + # Test at specific redshifts + z_values = np.array([0.0, 1.0, 2.0, 6.0]) + result = neijssel.mean_metallicity(z_values) + + # Calculate expected values based on Neijssel19's formula + expected = 0.035 * 10 ** (-0.23 * z_values) + + np.testing.assert_allclose(result, expected) + + def test_fsfr_with_lognormal(self): + """Test the overridden fSFR method which uses a ln-normal distribution""" + model_dict = {"sigma": 0.5, "Z_max": 0.3} + neijssel = Neijssel19(model_dict) + + # Test with redshift array and metallicity bins + z = np.array([0.0, 1.0]) + met_bins = np.array([0.001, 0.01, 0.02, 0.03]) + + result = neijssel.fSFR(z, met_bins) + + # Shape check - should be (len(z), len(met_bins)-1) + assert result.shape == (2, 3) + + # Test normalization with normalise=True + model_dict = {"sigma": 0.5, "Z_max": 0.3, "normalise": True} + neijssel = Neijssel19(model_dict) + result = neijssel.fSFR(z, met_bins) + expected = np.ones(len(z)) + np.testing.assert_allclose(np.sum(result, axis=1), expected) + +class TestFujimoto24: + """Tests for the Fujimoto24 SFH model""" + + def test_init_parameters(self): + """Test that initialization sets the correct CSFRD parameters""" + model_dict = {"sigma": 0.5, "Z_max": 0.03} + fujimoto = Fujimoto24(model_dict) + + # Check that CSFRD_params were set correctly + assert fujimoto.CSFRD_params["a"] == 0.010 + assert fujimoto.CSFRD_params["b"] == 2.8 + assert fujimoto.CSFRD_params["c"] == 3.3 + assert fujimoto.CSFRD_params["d"] == 6.6 + + # Check that it inherits correctly from MadauBase + assert isinstance(fujimoto, MadauBase) + +class TestChruslinska21: + """Tests for the Chruslinska21 SFH model with mocked data loading.""" + + @pytest.fixture + def mock_chruslinska_data(self, monkeypatch): + """Create mock data for the Chruslinska21 class.""" + # Create mock data for FOH bins + FOH_bins = np.linspace(5.3, 9.7, 200) + dFOH = FOH_bins[1] - FOH_bins[0] + redshifts = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0]) + delta_T = np.array([1e9, 1e9, 1e9, 1e9, 1e9, 1e9]) # Time bin widths + + # Mock SFR data - decreasing with redshift, varying with metallicity + SFR_data = np.zeros((len(redshifts), len(FOH_bins))) + for i in range(len(redshifts)): + # Simple pattern: peak at middle metallicity, decreasing with redshift + peak_idx = len(FOH_bins) // 2 + SFR_data[i] = np.exp(-0.5 * ((np.arange(len(FOH_bins)) - peak_idx) / 20)**2) + SFR_data[i] *= np.exp(-redshifts[i] / 2) # Decrease with redshift + + # Create method that returns tuple of time, redshift, deltaT + def mock_load_redshift_data(self, verbose=False): + time = np.array([1e9, 2e9, 3e9, 4e9, 5e9, 6e9]) # Fake times + return time, redshifts, delta_T + + # Create method that returns the mock SFR data + def mock_load_raw_data(self): + return SFR_data * 1e6 * delta_T[:, np.newaxis] + + # Patch the methods + monkeypatch.setattr(Chruslinska21, "_load_redshift_data", mock_load_redshift_data) + monkeypatch.setattr(Chruslinska21, "_load_raw_data", mock_load_raw_data) + + return { + "FOH_bins": FOH_bins, + "dFOH": dFOH, + "redshifts": redshifts, + "SFR_data": SFR_data + } + + @pytest.fixture + def chruslinska_model(self, mock_chruslinska_data): + """Create a Chruslinska21 model instance with mocked data.""" + model_dict = { + "sub_model": "test_model", + "Z_solar_scaling": "Asplund09", + "Z_max": 0.03, + "select_one_met": False + } + return Chruslinska21(model_dict) + + def test_init_parameters(self): + """Test that initialization validates required parameters.""" + # Test missing sub_model + with pytest.raises(ValueError) as excinfo: + Chruslinska21({"Z_solar_scaling": "Asplund09", "Z_max": 0.03, "select_one_met": False}) + assert "Sub-model not given!" in str(excinfo.value) + + # Test missing Z_solar_scaling + with pytest.raises(ValueError) as excinfo: + Chruslinska21({"sub_model": "test", "Z_max": 0.03, "select_one_met": False}) + assert "Z_solar_scaling not given!" in str(excinfo.value) + + def test_foh_to_z_conversion(self, chruslinska_model): + """Test the _FOH_to_Z method for all scaling options.""" + # Test Asplund09 scaling + FOH_test = np.array([7.0, 8.0, 8.69, 9.0]) + result = chruslinska_model._FOH_to_Z(FOH_test) + + # Expected: 10^(log10(0.0134) + FOH - 8.69) + expected = 10**(np.log10(0.0134) + FOH_test - 8.69) + np.testing.assert_allclose(result, expected) + + # Test other scaling options + model_dict = { + "sub_model": "test_model", + "Z_solar_scaling": "AndersGrevesse89", + "Z_max": 0.03, + "select_one_met": False + } + model = Chruslinska21(model_dict) + result = model._FOH_to_Z(FOH_test) + expected = 10**(np.log10(0.017) + FOH_test - 8.83) + np.testing.assert_allclose(result, expected) + + # Test GrevesseSauval98 scaling + model_dict["Z_solar_scaling"] = "GrevesseSauval98" + model = Chruslinska21(model_dict) + result = model._FOH_to_Z(FOH_test) + expected = 10**(np.log10(0.0201) + FOH_test - 8.93) + np.testing.assert_allclose(result, expected) + + # Test Villante14 scaling + model_dict["Z_solar_scaling"] = "Villante14" + model = Chruslinska21(model_dict) + result = model._FOH_to_Z(FOH_test) + expected = 10**(np.log10(0.019) + FOH_test - 8.85) + np.testing.assert_allclose(result, expected) + + # Test invalid scaling + model_dict["Z_solar_scaling"] = "InvalidScaling" + with pytest.raises(ValueError) as excinfo: + model = Chruslinska21(model_dict) + expected_str = ("Invalid Z_solar_scaling 'InvalidScaling'. " + "Valid options: ['Asplund09', 'AndersGrevesse89', " + "'GrevesseSauval98', 'Villante14']") + assert expected_str in str(excinfo.value) + + def test_mean_metallicity(self, chruslinska_model, mock_chruslinska_data): + """Test the mean_metallicity method.""" + # Test at specific redshifts + z_values = np.array([0.0, 2.0, 4.0]) + result = chruslinska_model.mean_metallicity(z_values) + print(result) + # Should be an array of the same length as z_values + assert len(result) == len(z_values) + # Should be the same value at all redshifts + assert np.isclose(result[0], result[1]) + assert np.isclose(result[0], result[2]) + assert np.isclose(result[0], 0.0014903210118641882) + + # Test with SFR_data == 0 + chruslinska_model.SFR_data = np.zeros_like(chruslinska_model.SFR_data) + with pytest.raises(AssertionError): + result = chruslinska_model.mean_metallicity(z_values) + + def test_csfrd_calculation(self, chruslinska_model, mock_chruslinska_data): + """Test the CSFRD method.""" + # Test at specific redshifts + z_values = np.array([0.0, 2.0, 4.0]) + result = chruslinska_model.CSFRD(z_values) + + # Should be an array of the same length as z_values + assert len(result) == len(z_values) + + # Should decrease with increasing redshift in the test case + assert result[0] > result[1] > result[2] + + def test_fsfr_calculation(self, chruslinska_model): + """Test the fSFR method.""" + # Test with redshift array and metallicity bins + z = np.array([0.0, 2.0]) + met_bins = np.array([0.001, 0.01, 0.02, 0.03]) + + result = chruslinska_model.fSFR(z, met_bins) + # Shape check - should be (len(z), len(met_bins)-1) + assert result.shape == (2, 3) + + # Test with normalization + chruslinska_model.normalise = True + result = chruslinska_model.fSFR(z, met_bins) + for row in result: + if np.sum(row) > 0: + np.testing.assert_allclose(np.sum(row), 1.0) + + # Test with Z_dist[i].sum = 0 + # Force the first mass array to be all zeros + chruslinska_model.SFR_data[0] = np.zeros_like(chruslinska_model.SFR_data[0]) + result = chruslinska_model.fSFR(z, met_bins) + np.testing.assert_allclose(result[0], np.zeros_like(result[0])) + + + +class TestZavala21: + """Tests for the Zavala21 SFH model with mocked data loading.""" + + @pytest.fixture + def mock_zavala_data(self, monkeypatch): + """Create mock data for the Zavala21 class.""" + # Create mock data - simple decreasing function with redshift + redshifts = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]) + SFRD_min = 0.1 * np.exp(-redshifts / 3.0) # Simple declining function + SFRD_max = 0.2 * np.exp(-redshifts / 3.0) # Double the min values + + def mock_read_csv(self, **kwargs): + return pd.DataFrame(data={ + "redshift": redshifts, + "SFRD_min": SFRD_min, + "SFRD_max": SFRD_max + }) + + monkeypatch.setattr(pd, "read_csv", mock_read_csv) + + def test_init_parameters(self, mock_zavala_data): + """Test that initialization validates and sets parameters correctly.""" + # Test missing sub_model + with pytest.raises(ValueError) as excinfo: + Zavala21({"Z_max": 0.03, "sigma": 0.5}) + assert "Sub-model not given!" in str(excinfo.value) + + # Test valid initialization with min model + model_dict = {"sub_model": "min", "Z_max": 0.03, "sigma": 0.5} + zavala_min = Zavala21(model_dict) + assert zavala_min.sub_model == "min" + assert zavala_min.Z_max == 0.03 + assert zavala_min.sigma == 0.5 + + # Test valid initialization with max model + model_dict = {"sub_model": "max", "Z_max": 0.03, "sigma": 0.5} + zavala_max = Zavala21(model_dict) + assert zavala_max.sub_model == "max" + + # Test invalid sub_model + model_dict = {"sub_model": "invalid", "Z_max": 0.03, "sigma": 0.5} + with pytest.raises(ValueError) as excinfo: + zavala_invalid = Zavala21(model_dict) + assert "Invalid sub-model!" in str(excinfo.value) + + def test_csfrd_min_model(self, mock_zavala_data): + """Test the CSFRD method with min sub-model.""" + model_dict = {"sub_model": "min", "Z_max": 0.03, "sigma": 0.5} + zavala = Zavala21(model_dict) + + # Test at specific redshifts + z_values = np.array([0.0, 2.0, 4.0, 6.0]) + result = zavala.CSFRD(z_values) + + # Expected values come from interpolating the mock data + expected = 0.1 * np.exp(-z_values / 3.0) + np.testing.assert_allclose(result, expected) + + def test_csfrd_max_model(self, mock_zavala_data): + """Test the CSFRD method with max sub-model.""" + model_dict = {"sub_model": "max", "Z_max": 0.03, "sigma": 0.5} + zavala = Zavala21(model_dict) + + # Test at specific redshifts + z_values = np.array([0.0, 2.0, 4.0, 6.0]) + result = zavala.CSFRD(z_values) + + # Expected values come from interpolating the mock data + expected = 0.2 * np.exp(-z_values / 3.0) + np.testing.assert_allclose(result, expected) + + def test_fsfr_calculation(self, mock_zavala_data): + """Test the fSFR method which is inherited from MadauBase.""" + model_dict = {"sub_model": "min", "Z_max": 0.03, "sigma": 0.5} + zavala = Zavala21(model_dict) + + # Test with redshift array and metallicity bins + z = np.array([0.0, 2.0]) + met_bins = np.array([0.001, 0.01, 0.02, 0.03]) + + result = zavala.fSFR(z, met_bins) + + # Shape check - should be (len(z), len(met_bins)-1) + assert result.shape == (2, 3) + + # Test normalization + zavala.normalise = True + result = zavala.fSFR(z, met_bins) + for row in result: + np.testing.assert_allclose(np.sum(row), 1.0) + +class TestGetSFHModel: + """Tests for the get_SFH_model function.""" + + def test_returns_correct_instance(self): + """Test that get_SFH_model returns the correct instance for each model.""" + # Test for MadauDickinson14 + model_dict = {"SFR": "Madau+Dickinson14", "sigma": 0.5, "Z_max": 0.03} + model = get_SFH_model(model_dict) + assert isinstance(model, MadauDickinson14) + + # Test for MadauFragos17 + model_dict = {"SFR": "Madau+Fragos17", "sigma": 0.5, "Z_max": 0.03} + model = get_SFH_model(model_dict) + assert isinstance(model, MadauFragos17) + + # Test for Neijssel19 + model_dict = {"SFR": "Neijssel+19", "sigma": 0.5, "Z_max": 0.03} + model = get_SFH_model(model_dict) + assert isinstance(model, Neijssel19) + + # Test for Fujimoto24 + model_dict = {"SFR": "Fujimoto+24", "sigma": 0.5, "Z_max": 0.03} + model = get_SFH_model(model_dict) + assert isinstance(model, Fujimoto24) + + def test_illustris_tng_model(self, monkeypatch): + """Test that get_SFH_model returns IllustrisTNG instance.""" + # Mock the data loading method + def mock_get_data(self, verbose=False): + # Return minimal mock data structure + return { + "SFR": np.array([0.1, 0.2, 0.3]), + "redshifts": np.array([0.0, 1.0, 2.0]), + "mets": np.array([0.001, 0.01, 0.02]), + "M": np.ones((3, 3)) + } + + # Patch the data loading method + monkeypatch.setattr(IllustrisTNG, "_get_illustrisTNG_data", mock_get_data) + + # Test the model creation + model_dict = {"SFR": "IllustrisTNG", "Z_max": 0.03} + model = get_SFH_model(model_dict) + assert isinstance(model, IllustrisTNG) + + def test_chruslinska_model(self, monkeypatch): + """Test that get_SFH_model returns Chruslinska21 instance.""" + # Mock the methods needed for initialization + def mock_load_data(self): + # Minimal setup to make initialization work + self.FOH_bins = np.linspace(5.3, 9.7, 10) + self.dFOH = self.FOH_bins[1] - self.FOH_bins[0] + self.Z = np.array([0.001, 0.01, 0.02]) + self.redshifts = np.array([0.0, 1.0, 2.0]) + self.SFR_data = np.ones((3, 10)) + + def mock_load_redshift(self, verbose=False): + # Return mock time, redshift, deltaT + return (np.array([1e9, 2e9, 3e9]), + np.array([0.0, 1.0, 2.0]), + np.array([1e9, 1e9, 1e9])) + + def mock_load_raw(self): + # Return mock data matrix + return np.ones((3, 10)) * 1e6 + + # Patch the methods + monkeypatch.setattr(Chruslinska21, "_load_chruslinska_data", mock_load_data) + monkeypatch.setattr(Chruslinska21, "_load_redshift_data", mock_load_redshift) + monkeypatch.setattr(Chruslinska21, "_load_raw_data", mock_load_raw) + + # Test the model creation + model_dict = { + "SFR": "Chruslinska+21", + "sub_model": "test", + "Z_solar_scaling": "Asplund09", + "Z_max": 0.03 + } + model = get_SFH_model(model_dict) + assert isinstance(model, Chruslinska21) + + def test_zavala_model(self, monkeypatch): + """Test that get_SFH_model returns Zavala21 instance.""" + # Mock the data loading method + def mock_load_data(self): + # Set required attributes directly + self.redshifts = np.array([0.0, 1.0, 2.0]) + if self.sub_model == "min": + self.SFR_data = np.array([0.1, 0.08, 0.06]) + else: + self.SFR_data = np.array([0.2, 0.16, 0.12]) + + # Patch the data loading method + monkeypatch.setattr(Zavala21, "_load_zavala_data", mock_load_data) + + # Test for min model + model_dict = { + "SFR": "Zavala+21", + "sub_model": "min", + "sigma": 0.5, + "Z_max": 0.03 + } + model = get_SFH_model(model_dict) + assert isinstance(model, Zavala21) + assert model.sub_model == "min" + + # Test for max model + model_dict = { + "SFR": "Zavala+21", + "sub_model": "max", + "sigma": 0.5, + "Z_max": 0.03 + } + model = get_SFH_model(model_dict) + assert isinstance(model, Zavala21) + assert model.sub_model == "max" + + def test_invalid_model(self): + """Test that get_SFH_model raises an error for an invalid model.""" + + model_dict = {"SFR": "InvalidModel"} + with pytest.raises(ValueError) as excinfo: + model = get_SFH_model(model_dict) + assert "Invalid SFR!" in str(excinfo.value) + +class TestSFR_per_met_at_z: + """Tests for SFR_per_met_at_z function.""" + + def test_SFR_per_met_at_z(self, monkeypatch): + """Test that SFR_per_met_at_z correctly calls the model.""" + # Create a mock model result + expected_result = np.array([[0.1, 0.2], [0.3, 0.4]]) + + # Mock SFH model class + class MockSFH: + def __call__(self, z, met_bins): + return expected_result + + mock_model = MockSFH() + + # Mock the get_SFH_model function to return our mock model + def mock_get_sfh_model(MODEL): + assert MODEL["SFR"] == "TestModel" # Verify correct model is requested + assert MODEL["param"] == "value" # Verify parameters are passed + return mock_model + + # Patch the function + monkeypatch.setattr( + "posydon.popsyn.star_formation_history.get_SFH_model", + mock_get_sfh_model + ) + + # Test the function + from posydon.popsyn.star_formation_history import SFR_per_met_at_z + + z = np.array([0.0, 1.0]) + met_bins = np.array([0.001, 0.01, 0.02]) + model_dict = {"SFR": "TestModel", "param": "value"} + + result = SFR_per_met_at_z(z, met_bins, model_dict) + + # Verify the result + np.testing.assert_array_equal(result, expected_result) diff --git a/posydon/unit_tests/utils/test_posydonwarning.py b/posydon/unit_tests/utils/test_posydonwarning.py index 3fc65c5f43..0641a15206 100644 --- a/posydon/unit_tests/utils/test_posydonwarning.py +++ b/posydon/unit_tests/utils/test_posydonwarning.py @@ -27,8 +27,9 @@ def test_dir(self): 'ClassificationWarning', 'EvolutionWarning',\ 'InappropriateValueWarning', 'IncompletenessWarning',\ 'InterpolationWarning', 'MissingFilesWarning',\ - 'NoPOSYDONWarnings', 'OverwriteWarning', 'POSYDONWarning',\ - 'Pwarn', 'ReplaceValueWarning', 'SetPOSYDONWarnings',\ + 'NoPOSYDONWarnings', 'OverwriteWarning', 'SFHModelWarning',\ + 'POSYDONWarning','Pwarn', 'ReplaceValueWarning',\ + 'SetPOSYDONWarnings',\ 'UnsupportedModelWarning', '_CAUGHT_POSYDON_WARNINGS',\ '_Caught_POSYDON_Warnings', '_POSYDONWarning_subclasses',\ '_POSYDON_WARNINGS_REGISTRY', '__authors__',\ @@ -101,6 +102,10 @@ def test_instance_UnsupportedModelWarning(self): assert isclass(totest.UnsupportedModelWarning) assert issubclass(totest.UnsupportedModelWarning,\ totest.POSYDONWarning) + + def test_instance_SFHModelWarning(self): + assert isclass(totest.SFHModelWarning) + assert issubclass(totest.SFHModelWarning, totest.POSYDONWarning) def test_instance_POSYDONWarning_subclasses(self): assert isinstance(totest._POSYDONWarning_subclasses, (dict)) diff --git a/posydon/utils/posydonwarning.py b/posydon/utils/posydonwarning.py index 738b4c86a2..78b347bad7 100644 --- a/posydon/utils/posydonwarning.py +++ b/posydon/utils/posydonwarning.py @@ -102,6 +102,11 @@ class UnsupportedModelWarning(POSYDONWarning): """Warnings related to selecting a model that is not supported.""" def __init__(self, message=''): super().__init__(message) + +class SFHModelWarning(POSYDONWarning): + """Warnings related to the SFH model.""" + def __init__(self, message=''): + super().__init__(message) # All POSYDON warnings subclasses should be defined beforehand