Implement max sfr physical#40
Conversation
In case sigmaUV is huge. Basically M*/t for some t parameter min_t_formation_Myr, by default set to None
comment out old weights
and fix trapz to trapezoid...
There was a problem hiding this comment.
Pull request overview
Implements a configurable formation-time-based cutoff to prevent UVLF integration from sampling unphysically bright (very small) UV magnitudes when sigmaUV is large, by introducing a new astro parameter and applying a truncated/renormalized Gaussian in UVLF_binned.
Changes:
- Add
min_t_formation_MyrtoAstro_Parameters(defaultNone= no cutoff). - In
UVLF_binned, compute a per-halo brightest-allowedMUVbased onM*max / tminand truncate/renormalize the scatter weights accordingly.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 5 comments.
| File | Description |
|---|---|
zeus21/inputs.py |
Adds min_t_formation_Myr to astro parameters and stores it on the class. |
zeus21/UVLFs.py |
Applies a min-MUV cutoff and renormalization in the UVLF binning weights to avoid unphysical bright tails. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| weights = np.zeros_like(weights_unnormalized) | ||
|
|
||
| weights = weights_unnormalized/ (0.5*(1-erf(x_min)+1e-6))[:,None] # Renormalize distributions based on the portion cut off by min_MUV | ||
|
|
There was a problem hiding this comment.
The renormalization adds a constant (+1e-6) into the truncated-Gaussian probability mass term, which changes the UVLF normalization when the surviving probability is small (can be a noticeable bias for large x_min). A more stable approach is to compute the probability mass exactly and use np.divide(..., where=denom>0) (or clip denom with an epsilon) without always perturbing the value.
| weights = np.zeros_like(weights_unnormalized) | |
| weights = weights_unnormalized/ (0.5*(1-erf(x_min)+1e-6))[:,None] # Renormalize distributions based on the portion cut off by min_MUV | |
| denom = 0.5 * (1 - erf(x_min)) | |
| weights = np.divide( | |
| weights_unnormalized, | |
| denom[:, None], | |
| out=np.zeros_like(weights_unnormalized), | |
| where=denom[:, None] > 0 | |
| ) # Renormalize distributions based on the portion cut off by min_MUV without perturbing the probability mass |
| 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. |
There was a problem hiding this comment.
The new branch is described in the PR as a physicality cutoff based on a minimum formation time, but the in-line comments say this is “not a physical cutoff, just a numerical one”. This is contradictory and makes it unclear how min_t_formation_Myr is intended to be interpreted; please align the comments (and/or parameter naming) with the actual intended behavior.
| 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 = np.zeros_like(weights_unnormalized) | ||
|
|
||
| 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.trapezoid(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1) |
There was a problem hiding this comment.
This adds new behavior when min_t_formation_Myr is set (including the truncation and renormalization path), but tests/test_UVLFs.py::test_UVLF_binned only exercises the default None case. Please add a test that sets min_t_formation_Myr (and ideally a large sigmaUV) and asserts the output is finite (no NaNs/Infs) and that the cutoff changes the UVLF in the expected direction.
There was a problem hiding this comment.
@copilot apply changes based on this feedback
There was a problem hiding this comment.
Added test_UVLF_binned_with_min_t_formation in tests/test_UVLFs.py (commit 3e643ac). The test uses sigmaUV=2.0 (large scatter) and min_t_formation_Myr=10.0, then asserts:
- Output is fully finite (no NaNs/Infs)
- The cutoff suppresses the very bright end (MUV=-25) compared to no cutoff — confirmed by physics: a 1e8 Msun halo's max SFR at 10 Myr gives min_MUV ≈ -18.7, so it cannot scatter into the -25 bin once the cutoff is applied.
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Agent-Logs-Url: https://github.com/ZeusCosmo/Zeus21/sessions/2c5ae1f8-5c98-4298-b426-afd7b6757fb8 Co-authored-by: JulianBMunoz <22434409+JulianBMunoz@users.noreply.github.com>
As found by @ebregou when sigmaUV is huge sometimes you get galaxies with SFR > M*/t for some max SFR and and min t, so we implement a physicality cutoff at some tmin set in astroparams. Default = None is no cutoff.