Merge main and quadratic SFRD

requirements.txt

@@ -1,10 +1,11 @@
 pytest
-numpy
+numpy>=2.0
 scipy
 mcfit
 numexpr
 astropy
 powerbox
 pyfftw
 sphinx
-myst_parser
+myst_parser
+tqdm

setup.py

@@ -14,7 +14,7 @@
           packages=['zeus21'],
           long_description=open('README.md').read(),
           install_requires=[
-           "numpy",
+           "numpy>=2.0",
            "scipy",
            "mcfit",
            "classy",

tests/test_UVLFs.py

@@ -138,4 +138,67 @@ def test_UVLF_binned():
     assert uvlf_nodust.shape == (3,)
 
     # Without dust, we expect different values than with dust
-    assert not np.array_equal(uvlf, uvlf_nodust)
+    assert not np.array_equal(uvlf, uvlf_nodust)
+
+
+def test_UVLF_binned_with_min_t_formation():
+    """Test that min_t_formation_Myr produces finite outputs and suppresses the bright end.
+
+    When sigmaUV is large, scatter can push small halos into unphysically bright bins.
+    Setting min_t_formation_Myr places a physical upper limit on each halo's SFR based on
+    its maximum stellar mass (all baryons converted to stars) and the minimum formation time.
+    This should suppress the very bright end of the UVLF without affecting the faint end.
+    """
+    UserParams = zeus21.User_Parameters()
+    CosmoParams_input = zeus21.Cosmo_Parameters_Input(kmax_CLASS=10., zmax_CLASS=20.)
+    ClassyCosmo = zeus21.runclass(CosmoParams_input)
+    CosmoParams = zeus21.Cosmo_Parameters(UserParams, CosmoParams_input, ClassyCosmo)
+    HMFintclass = zeus21.HMF_interpolator(UserParams, CosmoParams, ClassyCosmo)
+
+    # Use a large sigmaUV to create unphysical scatter into the bright end
+    large_sigmaUV = 2.0
+    min_t_Myr = 10.0
+
+    # AstroParams with the physicality cutoff applied
+    AstroParams_cut = zeus21.Astro_Parameters(
+        UserParams, CosmoParams,
+        sigmaUV=large_sigmaUV,
+        min_t_formation_Myr=min_t_Myr
+    )
+
+    # AstroParams without the cutoff (default None)
+    AstroParams_nocut = zeus21.Astro_Parameters(
+        UserParams, CosmoParams,
+        sigmaUV=large_sigmaUV
+    )
+
+    z_center = 6.0
+    z_width = 0.5
+    # Include a very bright bin (-25) where small-halo scatter is cut off,
+    # a typical bin (-20), and a faint bin (-15) that should be unaffected
+    MUV_centers = np.array([-25.0, -20.0, -15.0])
+    MUV_widths = np.full_like(MUV_centers, 1.0)
+
+    uvlf_cut = UVLF_binned(
+        AstroParams_cut, CosmoParams, HMFintclass,
+        z_center, z_width, MUV_centers, MUV_widths,
+        DUST_FLAG=False, RETURNBIAS=False
+    )
+    uvlf_nocut = UVLF_binned(
+        AstroParams_nocut, CosmoParams, HMFintclass,
+        z_center, z_width, MUV_centers, MUV_widths,
+        DUST_FLAG=False, RETURNBIAS=False
+    )
+
+    # Output must be finite (no NaNs or Infs) with the cutoff applied
+    assert np.all(np.isfinite(uvlf_cut)), "UVLF with min_t_formation_Myr cutoff contains NaN or Inf values"
+
+    # All values must be non-negative
+    assert np.all(uvlf_cut >= 0.0), "UVLF with min_t_formation_Myr cutoff contains negative values"
+
+    # The cutoff should suppress the very bright end: small halos that could not
+    # physically produce MUV=-25 galaxies (min_MUV~-18.7 for 1e8 Msun with t_min=10 Myr)
+    # no longer contribute via scatter, so the bright-end UVLF should be lower
+    assert uvlf_cut[0] < uvlf_nocut[0], (
+        "min_t_formation_Myr cutoff should suppress the very bright end (MUV=-25) of the UVLF"
+    )

tests/test_inputs.py

@@ -96,23 +96,23 @@ def test_inputs():
     #test Pop II Xray SED
     Energylisttest = np.logspace(2,np.log10(AstroParams.Emax_xray_norm),100)
     SEDXtab_test = AstroParams.SED_XRAY(Energylisttest, 2) #same in both models
-    normalization_XraySED = np.trapz(Energylisttest * SEDXtab_test,Energylisttest)
+    normalization_XraySED = np.trapezoid(Energylisttest * SEDXtab_test,Energylisttest)
     assert( normalization_XraySED == pytest.approx(1.0, 0.05) ) #5% is enough here
 
     #test Pop III Xray SED
     SEDXtab_test = AstroParams.SED_XRAY(Energylisttest, 3) #same in both models
-    normalization_XraySED = np.trapz(Energylisttest * SEDXtab_test,Energylisttest)
+    normalization_XraySED = np.trapezoid(Energylisttest * SEDXtab_test,Energylisttest)
     assert( normalization_XraySED == pytest.approx(1.0, 0.05) ) #5% is enough here
 
 
     #test Pop II LyA SED
     nulisttest = np.linspace(zeus21.constants.freqLyA, zeus21.constants.freqLyCont, 100)
     SEDLtab_test = AstroParams.SED_LyA(nulisttest, 2) #same in both models
-    normalization_LyASED = np.trapz(SEDLtab_test,nulisttest)
+    normalization_LyASED = np.trapezoid(SEDLtab_test,nulisttest)
     assert( normalization_LyASED == pytest.approx(1.0, 0.05) ) #5% is enough here
 
     #test Pop III LyA SED
     nulisttest = np.linspace(zeus21.constants.freqLyA, zeus21.constants.freqLyCont, 100)
     SEDLtab_test = AstroParams.SED_LyA(nulisttest, 3) #same in both models
-    normalization_LyASED = np.trapz(SEDLtab_test,nulisttest)
+    normalization_LyASED = np.trapezoid(SEDLtab_test,nulisttest)
     assert( normalization_LyASED == pytest.approx(1.0, 0.05) ) #5% is enough here

zeus21/UVLFs.py

@@ -80,9 +80,26 @@ def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwi
 
     xhi = np.subtract.outer(MUVcuthi, currMUV)/(np.sqrt(2) * sigmaUV)
     xlo = np.subtract.outer(MUVcutlo, currMUV )/(np.sqrt(2) * sigmaUV)
-    weights = (erf(xhi) - erf(xlo)).T/(2.0 * MUVwidths)
+
+    if (Astro_Parameters.min_t_formation_Myr == None):
+        min_MUV = -100.0 # essentially no cutoff, since the scatter is large at low masses and can cause numerical issues if we try to integrate over unphysically bright galaxies there. This is just a numerical cutoff, not a physical one, and the exact value doesn't matter much since the scatter is large there anyway.
+    else:
+        Mstarmax = HMF_interpolator.Mhtab * Cosmo_Parameters.OmegaB /Cosmo_Parameters.OmegaM #max stellar mass in each halo, if all baryons turned to stars
+        _tmaxSFR = Astro_Parameters.min_t_formation_Myr * 1e6 #arbitrary timescale to determine max SFR in yrs
+        SFRmax  = Mstarmax / (_tmaxSFR)
+        min_MUV = MUV_of_SFR(SFRmax, Astro_Parameters._kappaUV) #min MUV in each halo, if all baryons turned to stars at max SFR for 10 Myr. This is a very rough cutoff to avoid unphysically small MUVs (bright galaxies) at low masses, which can cause numerical issues since the scatter is large there. It's not a physical cutoff, just a numerical one. The exact value doesn't matter much since the scatter is large there anyway, but it prevents the code from trying to integrate over unphysically bright galaxies in low-mass halos.
+    x_min = (min_MUV - currMUV)/(np.sqrt(2) * sigmaUV)
+    xhi_cut = np.fmax(xhi, x_min)
+    xlo_cut = np.fmax(xlo, x_min)
+
+    weights_unnormalized = (erf(xhi_cut) - erf(xlo_cut)).T/(2.0 * MUVwidths)
+    weights = weights_unnormalized/ (0.5*(1-erf(x_min)+1e-6))[:,None] # Renormalize distributions based on the portion cut off by min_MUV
+
+    ### Standard as usual, no cuts:
+    # weights = (erf(xhi) - erf(xlo)).T/(2.0 * MUVwidths) #comment to myself, this 2 in denominator is correct here, nothing to do with the MUVwidths/2 a few lines above
 
-    UVLF_filtered = np.trapz(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1)
+    UVLF_filtered = np.trapezoid(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1)
+
 
     if(Astro_Parameters.USE_POPIII==False):
         return UVLF_filtered
@@ -98,7 +115,7 @@ def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwi
         xlo = np.subtract.outer(MUVcutlo, MUVbarlist_III)/(np.sqrt(2) * sigmaUV)
         weights = (erf(xhi) - erf(xlo)).T/(2.0 * MUVwidths)
 
-        UVLF_filtered_III = np.trapz(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1)
+        UVLF_filtered_III = np.trapezoid(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1)
 
         return UVLF_filtered, UVLF_filtered_III
 

zeus21/correlations.py

@@ -382,7 +382,7 @@ def __init__(self, User_Parameters, Cosmo_Parameters, Astro_Parameters, ClassCos
             #the z>zmax part of the integral we do aside. Assume Tk=Tadiabatic from CLASS.
             _zlisthighz_ = np.linspace(T21_coefficients.zintegral[-1], 99., 100) #beyond z=100 need to explictly tell CLASS to save growth
             _dgrowthhighz_ = cosmology.dgrowth_dz(Cosmo_Parameters, _zlisthighz_)
-            _hizintegral = np.trapz(cosmology.Tadiabatic(Cosmo_Parameters,_zlisthighz_)
+            _hizintegral = np.trapezoid(cosmology.Tadiabatic(Cosmo_Parameters,_zlisthighz_)
             /(1+_zlisthighz_)**2 * _dgrowthhighz_, _zlisthighz_)
 
         self._betaTad_ = -2./3. * _factor_adi_/self._lingrowthd * (np.cumsum(_integrand_adi[::-1])[::-1] + _hizintegral) #units of Tk_avg. Internal sum goes from high to low z (backwards), minus sign accounts for it properly so it's positive.
@@ -1261,7 +1261,7 @@ def get_Pk_from_xi(self, rsinput, xiinput):
 #
 #                 Probdtab = np.exp(-dtab**2/sigmaRRsq/2.0)
 #
-#                 norm = np.trapz(NionEPS * Probdtab, dtab)
+#                 norm = np.trapezoid(NionEPS * Probdtab, dtab)
 #                 NionEPS/=norm
 #
 #                 bindex = min(range(len(NionEPS)), key=lambda i: abs(NionEPS[i]-_invQbar))

zeus21/cosmology.py

@@ -75,13 +75,13 @@ def runclass(CosmologyIn):
         theta_b = velTransFunc['t_b']
         theta_c = velTransFunc['t_cdm']
 
-        sigma_vcb = np.sqrt(np.trapz(CosmologyIn.As * (kVel/0.05)**(CosmologyIn.ns-1) /kVel * (theta_b - theta_c)**2/kVel**2, kVel)) * constants.c_kms
+        sigma_vcb = np.sqrt(np.trapezoid(CosmologyIn.As * (kVel/0.05)**(CosmologyIn.ns-1) /kVel * (theta_b - theta_c)**2/kVel**2, kVel)) * constants.c_kms
         ClassCosmo.pars['sigma_vcb'] = sigma_vcb
 
         ###HAC: now computing average velocity assuming a Maxwell-Boltzmann distribution of velocities
         velArr = np.geomspace(0.01, constants.c_kms, 1000) #in km/s
         vavgIntegrand = (3 / (2 * np.pi * sigma_vcb**2))**(3/2) * 4 * np.pi * velArr**2 * np.exp(-3 * velArr**2 / (2 * sigma_vcb**2))
-        ClassCosmo.pars['v_avg'] = np.trapz(vavgIntegrand * velArr, velArr)
+        ClassCosmo.pars['v_avg'] = np.trapezoid(vavgIntegrand * velArr, velArr)
 
         ###HAC: Computing Vcb Power Spectrum
         ClassCosmo.pars['k_vcb'] = kVel
@@ -99,8 +99,8 @@ def runclass(CosmologyIn):
         j0bessel = lambda x: np.sin(x)/x
         j2bessel = lambda x: (3 / x**2 - 1) * np.sin(x)/x - 3*np.cos(x)/x**2
 
-        psi0 = 1 / 3 / (sigma_vcb/constants.c_kms)**2 * np.trapz(kVelIntp**2 / 2 / np.pi**2 * p_vcb_intp(np.log(kVelIntp)) * j0bessel(kVelIntp * np.transpose([rVelIntp])), kVelIntp, axis = 1)
-        psi2 = -2 / 3 / (sigma_vcb/constants.c_kms)**2 * np.trapz(kVelIntp**2 / 2 / np.pi**2 * p_vcb_intp(np.log(kVelIntp)) * j2bessel(kVelIntp * np.transpose([rVelIntp])), kVelIntp, axis = 1)
+        psi0 = 1 / 3 / (sigma_vcb/constants.c_kms)**2 * np.trapezoid(kVelIntp**2 / 2 / np.pi**2 * p_vcb_intp(np.log(kVelIntp)) * j0bessel(kVelIntp * np.transpose([rVelIntp])), kVelIntp, axis = 1)
+        psi2 = -2 / 3 / (sigma_vcb/constants.c_kms)**2 * np.trapezoid(kVelIntp**2 / 2 / np.pi**2 * p_vcb_intp(np.log(kVelIntp)) * j2bessel(kVelIntp * np.transpose([rVelIntp])), kVelIntp, axis = 1)
 
         k_eta, P_eta = mcfit.xi2P(rVelIntp, l=0, lowring = True)((6 * psi0**2 + 3 * psi2**2), extrap = False)
 

zeus21/inputs.py

@@ -270,7 +270,8 @@ def __init__(self, UserParams, Cosmo_Parameters,
                     A_vcb = 1.0,
                     beta_vcb = 1.8,
 
-                    quadratic_SFRD_lognormal = False # Sarah Libanore, use second order in lognormal
+                    quadratic_SFRD_lognormal = False, # Sarah Libanore, use second order in lognormal
+                    min_t_formation_Myr = None 
 
                 ):
 
@@ -412,7 +413,13 @@ def __init__(self, UserParams, Cosmo_Parameters,
                 if Cosmo_Parameters.Flag_emulate_21cmfast:
                     print('Quadratic SFRD not yet implemented when Flag_emulate_21cmfast = True; the code will use quadratic_SFRD_lognormal = False')
                 self.quadratic_SFRD_lognormal = False
-
+
+        if min_t_formation_Myr is not None:
+            if (not np.isscalar(min_t_formation_Myr)
+                    or not np.isfinite(min_t_formation_Myr)
+                    or min_t_formation_Myr <= 0):
+                raise ValueError("min_t_formation_Myr must be None or a strictly positive finite number.")
+        self.min_t_formation_Myr = min_t_formation_Myr #Minimum formation time of galaxies in Myr for UVLF, sets a minimum M*dot = M*/t_formation with fstar = 1
 
 
     def SED_XRAY(self, En, pop = 0): #pop set to zero as default, but it must be set to either 2 or 3