physical min_MUV in P(MUV|Mh)#35
Conversation
There was a problem hiding this comment.
Pull request overview
This PR adds support for enforcing a physical bright-end cutoff (minMUV) in the conditional magnitude distribution used to compute binned UV luminosity functions, and introduces an option to return the per-halo bin weights directly.
Changes:
- Extend
UVLF_binnedwith a halo-mass-dependentminMUVcutoff and renormalize the truncated distribution. - Add
RETURNWEIGHTSoption to return the computed MUV-bin weights. - Update module header edit attribution/date.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| # Cut distributions based on minMUV (user-input, halo mass depenent) | ||
| #MUVmin can be set to be MUV_of_SFR(max_SFR, Astro_Parameters._kappaUV), with max_SFR = Mstar/min_t & Mstar = fb*Mh |
There was a problem hiding this comment.
The comment suggests minMUV can be computed via MUV_of_SFR(max_SFR, ...), which is an intrinsic (pre-dust) magnitude, but the cutoff is applied in terms of currMUV (which is dust-iterated when DUST_FLAG=True). Please clarify whether minMUV is expected to be in intrinsic or dust-attenuated (observed) magnitude space, and adjust either the implementation or the comment accordingly so the cutoff is applied in the intended domain.
| # Cut distributions based on minMUV (user-input, halo mass depenent) | |
| #MUVmin can be set to be MUV_of_SFR(max_SFR, Astro_Parameters._kappaUV), with max_SFR = Mstar/min_t & Mstar = fb*Mh | |
| # Cut distributions based on minMUV (user-input, halo mass dependent). | |
| # NOTE: minMUV is expected to be in the same magnitude space as currMUV, i.e. dust-attenuated/observed | |
| # magnitudes when DUST_FLAG == True. The helper MUV_of_SFR(max_SFR, Astro_Parameters._kappaUV) returns | |
| # intrinsic (pre-dust) magnitudes; if you use it to set minMUV you should either (i) work without dust | |
| # attenuation (DUST_FLAG == False), or (ii) convert those intrinsic magnitudes to observed ones using the | |
| # same dust model that defines currMUV. A typical choice is max_SFR = Mstar/min_t with Mstar = fb * Mh. |
| #MUVmin can be set to be MUV_of_SFR(max_SFR, Astro_Parameters._kappaUV), with max_SFR = Mstar/min_t & Mstar = fb*Mh | ||
| if minMUV is None: | ||
| minMUV = np.full_like(HMF_interpolator.Mhtab, -100) | ||
| x_min = (minMUV - 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)))[:,None] # Renormalize distributions based on the portion cut off by minMUV | ||
|
|
There was a problem hiding this comment.
weights = weights_unnormalized / (0.5*(1-erf(x_min))) can divide by ~0 when minMUV is much fainter than currMUV (large positive x_min), producing inf/NaN weights and UVLFs. Please guard against a zero tail probability (e.g., clip the denominator to an epsilon and/or set weights to 0 where the surviving probability is 0), and consider validating/broadcasting minMUV to a 1D array of length len(Mhtab) to avoid shape surprises.
| #MUVmin can be set to be MUV_of_SFR(max_SFR, Astro_Parameters._kappaUV), with max_SFR = Mstar/min_t & Mstar = fb*Mh | |
| if minMUV is None: | |
| minMUV = np.full_like(HMF_interpolator.Mhtab, -100) | |
| x_min = (minMUV - 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)))[:,None] # Renormalize distributions based on the portion cut off by minMUV | |
| # MUVmin can be set to be MUV_of_SFR(max_SFR, Astro_Parameters._kappaUV), | |
| # with max_SFR = Mstar/min_t & Mstar = fb*Mh | |
| if minMUV is None: | |
| # Default to a very faint cutoff for all halo masses | |
| minMUV = np.full_like(HMF_interpolator.Mhtab, -100.0) | |
| else: | |
| # Ensure minMUV is a 1D array aligned with the halo-mass grid | |
| minMUV_arr = np.asarray(minMUV) | |
| if minMUV_arr.ndim == 0: | |
| # Broadcast scalar to all halo masses | |
| minMUV = np.full_like(HMF_interpolator.Mhtab, float(minMUV_arr)) | |
| else: | |
| if minMUV_arr.shape != HMF_interpolator.Mhtab.shape: | |
| raise ValueError( | |
| "minMUV must be None, a scalar, or have the same shape as " | |
| "HMF_interpolator.Mhtab (got shape %s vs %s)" | |
| % (minMUV_arr.shape, HMF_interpolator.Mhtab.shape,) | |
| ) | |
| minMUV = minMUV_arr | |
| x_min = (minMUV - 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) | |
| # Renormalize distributions based on the portion cut off by minMUV. | |
| # Guard against zero/negative tail probabilities to avoid inf/NaN. | |
| tail_prob = 0.5 * (1.0 - erf(x_min)) | |
| # Define where the surviving probability is effectively zero or invalid | |
| zero_tail_mask = tail_prob <= 0.0 | |
| # Use a safe denominator for valid entries; dummy value (1.0) where masked | |
| safe_tail_prob = np.where(zero_tail_mask, 1.0, tail_prob) | |
| weights = weights_unnormalized / safe_tail_prob[:, None] | |
| # Set weights to zero where there is no surviving probability | |
| if np.any(zero_tail_mask): | |
| weights[zero_tail_mask, :] = 0.0 |
| # Cut distributions based on minMUV (user-input, halo mass depenent) | ||
| #MUVmin can be set to be MUV_of_SFR(max_SFR, Astro_Parameters._kappaUV), with max_SFR = Mstar/min_t & Mstar = fb*Mh | ||
| if minMUV is None: | ||
| minMUV = np.full_like(HMF_interpolator.Mhtab, -100) |
There was a problem hiding this comment.
Using -100 as the “no cutoff” default for minMUV is a magic number. Using -np.inf (or a named constant) would make the intent explicit and avoid having to reason about whether -100 is always safely beyond the bright-end of the model.
| minMUV = np.full_like(HMF_interpolator.Mhtab, -100) | |
| minMUV = np.full_like(HMF_interpolator.Mhtab, -np.inf) |
| # Cut distributions based on minMUV (user-input, halo mass depenent) | ||
| #MUVmin can be set to be MUV_of_SFR(max_SFR, Astro_Parameters._kappaUV), with max_SFR = Mstar/min_t & Mstar = fb*Mh | ||
| if minMUV is None: | ||
| minMUV = np.full_like(HMF_interpolator.Mhtab, -100) | ||
| x_min = (minMUV - 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)))[:,None] # Renormalize distributions based on the portion cut off by minMUV | ||
|
|
||
| if RETURNWEIGHTS: | ||
| return weights |
There was a problem hiding this comment.
New behavior (minMUV truncation + renormalization and the RETURNWEIGHTS early-return path) isn’t covered by the existing tests/test_UVLFs.py::test_UVLF_binned. Please add tests that (1) verify bins brighter than minMUV get zero weight/number density and (2) verify RETURNWEIGHTS=True returns correctly-shaped weights that integrate to 1 over MUV for each halo mass after truncation.
|
|
||
| #and combine to get UVLF: | ||
| def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwidth, MUVcenters, MUVwidths, DUST_FLAG=True, RETURNBIAS = False): | ||
| def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwidth, MUVcenters, MUVwidths, minMUV = None, DUST_FLAG=True, RETURNBIAS = False, RETURNWEIGHTS = False): |
There was a problem hiding this comment.
Adding minMUV before DUST_FLAG changes the positional-argument order of UVLF_binned, which can silently break downstream callers that pass DUST_FLAG/RETURNBIAS positionally. To keep backward compatibility, consider moving minMUV (and RETURNWEIGHTS) to the end of the signature or making the new args keyword-only. Also note minMUV/RETURNWEIGHTS currently only affect the Pop II branch; when USE_POPIII is enabled the Pop III weights are computed without this cutoff/return path, so the API behavior is inconsistent unless that’s explicitly intended/documented.
| def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwidth, MUVcenters, MUVwidths, minMUV = None, DUST_FLAG=True, RETURNBIAS = False, RETURNWEIGHTS = False): | ||
| 'Binned UVLF in units of 1/Mpc^3/mag, for bins at <zcenter> with a Gaussian width zwidth, centered at MUV centers with tophat width MUVwidths. z width only in HMF since that varies the most rapidly. If flag RETURNBIAS set to true it returns number-avgd bias instead of UVLF, still have to divide by UVLF' |
There was a problem hiding this comment.
The function docstring still describes only DUST_FLAG / RETURNBIAS behavior, but minMUV and RETURNWEIGHTS materially change outputs (including an early return). Please update the docstring to document the meaning/expected shape/units of minMUV and what is returned when RETURNWEIGHTS=True (and how that interacts with USE_POPIII).
|
|
||
|
|
||
| if(RETURNBIAS==True): # weight by bias | ||
| if(RETURNBIAS==True): # weight by bias) |
There was a problem hiding this comment.
There’s an extra ) in this comment (# weight by bias)) which looks like a typo and may confuse readers/skimmers. Please remove the stray parenthesis.
| if(RETURNBIAS==True): # weight by bias) | |
| if(RETURNBIAS==True): # weight by bias |
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
No description provided.