Merge branch 'main' into BMF_and_ionized_maps

Author: yonboyage

.codecov_trigger

@@ -0,0 +1 @@
+# This file triggers the CI workflow

.github/workflows/python-tests.yml

@@ -0,0 +1,67 @@
+name: Python Tests
+
+on:
+  push:
+    branches: [ main ]
+  pull_request:
+    branches: [ main ]
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+    
+    steps:
+    - uses: actions/checkout@v4
+    - name: Set up Python 3.11
+      uses: actions/setup-python@v4
+      with:
+        python-version: '3.11'
+    
+    - name: Install CLASS dependencies
+      run: |
+        sudo apt-get update
+        sudo apt-get install -y build-essential gfortran
+
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install pytest pytest-cov
+        pip install Cython
+        pip install -r requirements.txt
+        
+    - name: Install CLASS
+      run: |
+        git clone https://github.com/lesgourg/class_public.git
+        cd class_public
+        make
+        cd python
+        python setup.py install
+        cd ../../
+        # Keep the directory for data files
+        
+    - name: Install package
+      run: |
+        pip install -e .
+    
+    - name: Debug SFR_III function
+      env:
+        CLASSDIR: ${{ github.workspace }}/class_public
+      run: |
+        python -c "import zeus21; from zeus21.sfrd import SFR_III; import inspect; print('SFR_III parameters:', inspect.signature(SFR_III)); print('Parameter count:', len(inspect.signature(SFR_III).parameters))"
+        
+    - name: Run tests with coverage
+      env:
+        CLASSDIR: ${{ github.workspace }}/class_public
+      run: |
+        python -m pytest --cov=zeus21 --cov-report=xml --cov-report=term tests/
+    
+    - name: Upload coverage reports to Codecov
+      uses: codecov/codecov-action@v4
+      env:
+        CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
+      with:
+        files: ./coverage.xml
+        flags: unittests
+        name: codecov-umbrella
+        verbose: true
+        fail_ci_if_error: false

README.md

@@ -4,6 +4,9 @@
 
 # Zeus21: Lightning-fast simulations of cosmic dawn
 
+[![Tests](https://github.com/JulianBMunoz/Zeus21/actions/workflows/python-tests.yml/badge.svg)](https://github.com/JulianBMunoz/Zeus21/actions/workflows/python-tests.yml)
+[![codecov](https://codecov.io/gh/JulianBMunoz/Zeus21/branch/main/graph/badge.svg)](https://codecov.io/gh/JulianBMunoz/Zeus21)
+
 Zeus21 encodes the effective model for the 21-cm power spectrum and global signal from [Muñoz 2023a](https://arxiv.org/abs/2302.08506). The goal is to capture all the nonlocal and nonlinear physics of cosmic dawn in a light and fully Pythonic code. Zeus21 takes advantage of the approximate log-normality of the star-formation rate density (SFRD) during cosmic dawn to compute the 21-cm power spectrum analytically. It agrees with more expensive semi-numerical simulations to roughly 10% precision, but has comparably negligible computational cost (~ s) and memory requirements. Now Zeus21 can also predict galaxy UV luminosity functions (UVLFs) and their linear clustering (galaxy bias) at any z, see [Muñoz et al. 2023b](https://arxiv.org/abs/2306.09403) for the implementation and application to JWST. Now Zeus21 includes Population III stars with inhomogeneous Lyman-Werner feedback and relative velocities (with their fluctuations), as described in [Cruz et al. 2024](https://arxiv.org/abs/2407.18294).
 
 Zeus21 (Zippy Early-Universe Solver for 21-cm) pairs well with data from [HERA](https://reionization.org/), but can be used for any 21-cm inference or prediction. Current capabilities include finding the 21-cm power spectrum (at a broad range of k and z), the global signal, IGM temperatures (Tk, Ts, Tcolor), neutral fraction xHI, Lyman-alpha fluxes, and the evolution of the SFRD; all across cosmic dawn z=5-35. Zeus21 can use three different astrophysical models, one of which emulates 21cmFAST, and can vary the cosmology through CLASS. 

codecov.yml

@@ -0,0 +1,28 @@
+codecov:
+  require_ci_to_pass: yes
+
+coverage:
+  precision: 2
+  round: down
+  range: "70...100"
+
+  status:
+    project:
+      default:
+        target: auto
+        threshold: 5%
+    patch: yes
+    changes: no
+
+parsers:
+  gcov:
+    branch_detection:
+      conditional: yes
+      loop: yes
+      method: no
+      macro: no
+
+comment:
+  layout: "reach,diff,flags,files,footer"
+  behavior: default
+  require_changes: no

docs/Tutorial_Zeus21_UVLFs.ipynb

@@ -226,14 +226,27 @@
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "#and the newly added pop3 UVLF - note they'll normally be very faint, so you'll need to go to fainter MUVs\n",
+    "AstroParams_popIII = zeus21.Astro_Parameters(UserParams, CosmoParams, accretion_model=0, USE_POPIII=True) "
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "UVLFs_pop2,UVLFs_pop3= zeus21.UVLFs.UVLF_binned(AstroParams_popIII,CosmoParams,HMFintclass,z,dz,MUVcenters,MUVwidths)\n",
+    "\n",
+    "\n",
+    "plt.semilogy(MUVcenters,UVLFs_pop2,'k-')\n",
+    "plt.semilogy(MUVcenters,UVLFs_pop3,'r--')\n",
+    "plt.xlim(-22,-17)\n",
+    "plt.ylim(1e-6,1e-1)\n",
+    "plt.xlabel(r'$M_{\\rm UV}$');\n",
+    "plt.ylabel(r'$\\Phi_{\\rm UV}\\,\\rm [Mpc^{-3}\\,mag^{-1}]$');"
+   ]
   }
  ],
  "metadata": {

requirements.txt

@@ -1,11 +1,12 @@
 pytest
 numpy
 scipy
-sphinx
-myst_parser
+mcfit
 numexpr
 astropy
 tqdm
 matplotlib
+pyfftw
 powerbox
-pyfftw
+sphinx
+myst_parser

tests/test_UVLFs.py

@@ -0,0 +1,141 @@
+"""
+
+Test UV luminosity functions for Zeus21
+
+Author: Claude AI
+April 2025
+
+"""
+
+import pytest
+import zeus21
+import numpy as np
+
+from zeus21.UVLFs import UVLF_binned, MUV_of_SFR, AUV, beta
+
+def test_MUV_of_SFR():
+    """Test the conversion from SFR to UV magnitudes"""
+    # Test a range of SFR values
+    SFR_test = np.logspace(-3, 2, 10)  # M_sun/yr
+    kappaUV_test = 1.15e-28  # Typical value
+    
+    # Calculate MUV
+    MUV_result = MUV_of_SFR(SFR_test, kappaUV_test)
+    
+    # Check that increasing SFR leads to brighter (more negative) MUV
+    assert np.all(np.diff(MUV_result) < 0)
+    
+    # Check specific value based on the formula M_UV = 51.63 - 2.5*log10(SFR/kappaUV)
+    # For SFR = 1 M_sun/yr with kappaUV = 1.15e-28
+    expected_MUV = 51.63 - 2.5 * np.log10(1.0/1.15e-28)
+    assert MUV_of_SFR(np.array([1.0]), kappaUV_test)[0] == pytest.approx(expected_MUV)
+    
+    # Test different kappaUV values
+    kappaUV_test2 = 2.0e-28
+    MUV_result2 = MUV_of_SFR(SFR_test, kappaUV_test2)
+    
+    # Higher kappaUV should result in fainter magnitudes (more positive)
+    assert np.all(MUV_result2 > MUV_result)
+
+def test_beta_function():
+    """Test the beta (UV slope) calculation"""
+    # Test a single redshift and magnitude but use arrays as the function expects
+    z_test = np.array([5.0])
+    MUV_test = np.array([-20.0])
+    
+    # Calculate beta value
+    beta_value = beta(z_test, MUV_test)
+    
+    # Check that beta value is reasonable (typical range is -3 to -1)
+    assert beta_value > -3.0
+    assert beta_value < -1.0
+    
+    # Test at pivot point
+    MUV_pivot = np.array([-19.5])  # The pivot point defined in the code
+    beta_at_pivot = beta(z_test, MUV_pivot)
+    
+    # Check that a value is returned
+    assert isinstance(beta_at_pivot, np.ndarray)
+
+def test_AUV_function():
+    """Test the dust attenuation calculation"""
+    # Set up parameters
+    UserParams = zeus21.User_Parameters()
+    CosmoParams_input = zeus21.Cosmo_Parameters_Input()
+    ClassyCosmo = zeus21.runclass(CosmoParams_input)
+    CosmoParams = zeus21.Cosmo_Parameters(UserParams, CosmoParams_input, ClassyCosmo)
+    AstroParams = zeus21.Astro_Parameters(UserParams, CosmoParams)
+    
+    # Test with arrays as the function expects
+    z_test = np.array([5.0])
+    MUV_test = np.array([-20.0])
+    
+    # Calculate dust attenuation
+    A_UV = AUV(AstroParams, z_test, MUV_test)
+    
+    # Check that attenuation is non-negative
+    assert np.all(A_UV >= 0.0)
+    
+    # Test the HIGH_Z_DUST flag behavior
+    z_high = np.array([9.0])  # High redshift above _zmaxdata
+    _zmaxdata = 8.0
+    
+    # Test with HIGH_Z_DUST=True (dust applied at high z)
+    A_UV_high = AUV(AstroParams, z_high, MUV_test, HIGH_Z_DUST=True)
+    
+    # Test with HIGH_Z_DUST=False (no dust above _zmaxdata)
+    A_UV_no_highz = AUV(AstroParams, z_high, MUV_test, HIGH_Z_DUST=False, _zmaxdata=_zmaxdata)
+    
+    # HIGH_Z_DUST=False should give zero attenuation for z > _zmaxdata
+    assert np.all(A_UV_no_highz == 0.0)
+
+def test_UVLF_binned():
+    """Test the binned UV luminosity function calculation"""
+    # Set up parameters
+    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)
+    AstroParams = zeus21.Astro_Parameters(UserParams, CosmoParams)
+    HMFintclass = zeus21.HMF_interpolator(UserParams, CosmoParams, ClassyCosmo)
+    
+    # Test data
+    z_center = 6.0
+    z_width = 0.5
+    MUV_centers = np.array([-22.0, -20.0, -18.0])
+    MUV_widths = np.full_like(MUV_centers, 1.0)
+    
+    # Calculate UVLF
+    uvlf = UVLF_binned(AstroParams, CosmoParams, HMFintclass, z_center, z_width, 
+                       MUV_centers, MUV_widths, DUST_FLAG=True, RETURNBIAS=False)
+    
+    # Check dimensions
+    assert uvlf.shape == (3,)
+    
+    # Check that values are positive
+    assert np.all(uvlf >= 0.0)
+    
+    # Test that fainter (more positive MUV) bins typically have higher number densities
+    # This is a general trend for LFs, but not strictly required
+    # We'll do a weak test that they're not all identical
+    assert len(np.unique(uvlf)) > 1
+    
+    # Test RETURNBIAS flag
+    bias_values = UVLF_binned(AstroParams, CosmoParams, HMFintclass, z_center, z_width, 
+                             MUV_centers, MUV_widths, DUST_FLAG=True, RETURNBIAS=True)
+    
+    # Check dimensions
+    assert bias_values.shape == (3,)
+    
+    # Check that biases are positive
+    assert np.all(bias_values >= 0.0)
+    
+    # Test without dust correction
+    uvlf_nodust = UVLF_binned(AstroParams, CosmoParams, HMFintclass, z_center, z_width, 
+                             MUV_centers, MUV_widths, DUST_FLAG=False, RETURNBIAS=False)
+    
+    # Check dimensions
+    assert uvlf_nodust.shape == (3,)
+    
+    # Without dust, we expect different values than with dust
+    assert not np.array_equal(uvlf, uvlf_nodust)

tests/test_astrophysics.py

@@ -57,7 +57,7 @@ def test_background():
     assert( (0 <= sSFR).all()) #positive
     assert( (sSFR/zeus21.cosmology.Hubinvyr(CosmoParams,ztest) <= 1).all()) #make sure sSFR/H < 1 (not all mass forms stars in a Hubble time)
 
-    sSFR3 = SFR_III(AstroParams, CosmoParams, ClassyCosmo, HMFintclass, HMFintclass.Mhtab, Coeffs_popIII.J21LW_interp_conv_avg, ztest, ztest, ClassyCosmo.pars['v_avg'])/HMFintclass.Mhtab
+    sSFR3 = SFR_III(AstroParams, CosmoParams, HMFintclass, HMFintclass.Mhtab, Coeffs_popIII.J21LW_interp_conv_avg, ztest, ztest, ClassyCosmo.pars['v_avg'])/HMFintclass.Mhtab
     assert( (0 <= sSFR3).all()) #positive
     assert( (sSFR3/zeus21.cosmology.Hubinvyr(CosmoParams,ztest) <= 1).all()) #make sure sSFR3/H < 1 (not all mass forms stars in a Hubble time)
 
@@ -67,7 +67,7 @@ def test_background():
     assert( (0 <= sSFR_exp).all())
     assert( (sSFR_exp/zeus21.cosmology.Hubinvyr(CosmoParams,ztest) <= 1).all())
 
-    sSFR_exp3 = SFR_III(AstroParams_expacc, CosmoParams, ClassyCosmo, HMFintclass, HMFintclass.Mhtab, Coeffs_popIII.J21LW_interp_conv_avg, ztest, ztest, ClassyCosmo.pars['v_avg'])/HMFintclass.Mhtab
+    sSFR_exp3 = SFR_III(AstroParams_expacc, CosmoParams, HMFintclass, HMFintclass.Mhtab, Coeffs_popIII.J21LW_interp_conv_avg, ztest, ztest, ClassyCosmo.pars['v_avg'])/HMFintclass.Mhtab
     assert( (0 <= sSFR_exp3).all())
     assert( (sSFR_exp3/zeus21.cosmology.Hubinvyr(CosmoParams,ztest) <= 1).all())
 
@@ -77,7 +77,7 @@ def test_background():
     assert( (0 <= sSFR_21cmfast).all())
     assert( (sSFR_21cmfast/zeus21.cosmology.Hubinvyr(CosmoParams_21cmfast,ztest) <= 1).all())
 
-    sSFR_21cmfast3 = SFR_III(AstroParams_expacc, CosmoParams_21cmfast, ClassyCosmo, HMFintclass, HMFintclass.Mhtab, Coeffs_popIII.J21LW_interp_conv_avg, ztest, ztest, ClassyCosmo.pars['v_avg'])/HMFintclass.Mhtab
+    sSFR_21cmfast3 = SFR_III(AstroParams_expacc, CosmoParams_21cmfast, HMFintclass, HMFintclass.Mhtab, Coeffs_popIII.J21LW_interp_conv_avg, ztest, ztest, ClassyCosmo.pars['v_avg'])/HMFintclass.Mhtab
     assert( (0 <= sSFR_21cmfast3).all())
     assert( (sSFR_21cmfast3/zeus21.cosmology.Hubinvyr(CosmoParams_21cmfast,ztest) <= 1).all())