Beta esc #538

zhaotingchen · 2025-07-20T22:10:32Z

Attempt to fix #537

zhaotingchen · 2025-07-20T22:17:22Z

A naive test showing that it is affecting the ioniosation efficiency:

import numpy as np
import matplotlib.pyplot as plt
import py21cmfast as p21c
from multiprocessing import cpu_count

max_num_thread = cpu_count()
sim_op = {
    "HII_DIM":32, 
    "BOX_LEN":250,
    "N_THREADS":max_num_thread,
}
matter_op = {
    "PERTURB_ON_HIGH_RES":False, # True for production
    "MINIMIZE_MEMORY":True, # False for production
    "USE_HALO_FIELD":True, # True for production
    "HALO_STOCHASTICITY":False, # True for production
}
astro_op = {
    "USE_TS_FLUCT":False, # True for production
    "SUBCELL_RSD":False, # True for production
    "INHOMO_RECO":False, # True for production
    "USE_MINI_HALOS":False,
    "USE_MASS_DEPENDENT_ZETA":True,
    "PHOTON_CONS_TYPE":"no-photoncons", #v4 # change to something else for production
    #"PHOTON_CONS":True,
}
astro_pars = p21c.AstroParams()
z_min = 5.0
z_max = 10.0
cosmo_pars = p21c.CosmoParams()
cache = p21c.OutputCache(cache_dir)

inputs = p21c.InputParameters(
    random_seed=12345,
    cosmo_params=cosmo_pars,
    matter_options=matter_op,
    astro_options=astro_op,
    astro_params=astro_pars,
    simulation_options=sim_op,
)

coevals = []
for coeval,in_outputs in p21c.generate_coeval(
    inputs=inputs,
    out_redshifts=[8, 9],
    cache=cache,
):
    if in_outputs:
        coevals.append(coeval)

plt.scatter(np.log10(coevals[0].halo_mass).ravel(),np.log10(coevals[0].n_ion.ravel()),)

The plot will show a straight line due to the scaling between f_esc and halo_mass.

In default, BETA_ESC is 0 and this scaling does not change with redshift:

print(coevals[0].redshift,np.polyfit(np.log10(coevals[0].halo_mass).ravel(),np.log10(coevals[0].n_ion.ravel()),1))
print(coevals[1].redshift,np.polyfit(np.log10(coevals[1].halo_mass).ravel(),np.log10(coevals[1].n_ion.ravel()),1))

gives

9.0 [ 1.66731392 -6.39808893]
8.0 [ 1.66509232 -6.36885117]

changing the input to astro_pars = p21c.AstroParams(BETA_ESC = -1.5) and repeat gives

9.0 [ 1.66731389 -6.54345369]
8.0 [ 1.66509231 -6.44557988]

now there is a noticeable change in the intercept. Larger z have smaller n_ion due to beta<0.

BETA_ESC = 1.5 gives

9.0 [ 1.66731392 -6.25272388]
8.0 [ 1.66509231 -6.29212233]

Where the relation flips as expected (but the change is too small?)

daviesje

I'm assuming from your issue that you have already spoken with @andreimesinger or @steven-murray about the pros/cons of adding new parameters to the escape fraction so I'm just going to review the code here.

As for the "photon-conservation" model, I don't immediately see a problem with running our calibration simulation with a BETA_ESC in any of the cases, but maybe you found something wrong with the f-photoncons model here? In any case I wouldn't want these modules getting in the way of new features so in the worst case we can simply throw parameter errors when a user enables both

daviesje · 2025-07-20T23:41:00Z

src/py21cmfast/src/hmf.c

+    double Zesc = log(pow((1 + p.redshift) / 8.0, p.beta_esc));

-    return exp(Fstar + Fesc - p.Mturn_acg / exp(lnM) + lnM);
+    return exp(Fstar + Fesc + Zesc - p.Mturn_acg / exp(lnM) + lnM);


These are very hot functions. So saving log/pow calls here is very important. You can see how the above function log_scaling_PL_limit does this for the other power-laws and you could use two less here. However more importantly, this term is only calculated once per redshift, so it can be saved as a constant in the scaling parameters and simply loaded in.

thanks, I see that I actually double-counted this factor twice which is wrong (I already added this factor in Nion_ConditionalM). Sorry the code structure is too complicated for me to understand rn so this is probably stupid: Should I add it at prefactor_nion level in set_fixed_grids and get_box_averages? Or should I add this in set_scaling_constants level by directly changing f_esc10? Or somewhere else?

So I would implement this in the following way:

Add a new double field to the ScalingConstants struct (in scaling_relations.h) which contains the entire power-law term from BETA_ESC

Add a line to set_scaling_constants (in scaling_relations.c) which initialises this term based on the redshift, this calculates the term once per redshift and stores it for fast multiplication later (no more pow` calls)

Multiply by this factor in two places: in HaloBox.c where the escape fraction is set in set_halo_properties, and in IonisationBox.c where the constant ion_eff_factor_gl is set in set_ionbox_constants.

Regarding the last step, the reason you have to alter two parts is just due to the way the discrete halo / density field models are split. It's not ideal but that's the way it's organised at the moment. Happy to answer any further questions

thanks for the detailed instruction, please let me know if the latest commit is doing things correctly. My only addition to your comment is that I think the prefactor_nion in set_fixed_grids and get_box_averages needs to be modified as well, otherwise the sampling below the mass limit when AVG_BELOW_SAMPLER will not have the correct scaling

daviesje · 2025-07-20T23:43:25Z

src/py21cmfast/src/hmf.c

-    return IntegratedNdM(lnM_Min, lnM_Max, params, &u_nion_integrand, 0);
+    // z factor can be taken out of the integral
+    double Zesc = pow((1 + z) / 8.0, sc->beta_esc);
+    return Zesc * IntegratedNdM(lnM_Min, lnM_Max, params, &u_nion_integrand, 0);


See above comment, redshift-based escape fractions should be pre-computed once per snapshot

daviesje · 2025-07-20T23:43:45Z

src/py21cmfast/src/hmf.c

    if (params.HMF != 0 && params.HMF != 1 && params.HMF != 4) params.HMF = 0;
-    return IntegratedNdM(lnM1, lnM2, params, &c_nion_integrand, method);
+    // z factor can be taken out of the integral
+    return Zesc * IntegratedNdM(lnM1, lnM2, params, &c_nion_integrand, method);


Here as well.

daviesje · 2025-07-20T23:44:35Z

src/py21cmfast/wrapper/inputs.py

            raise ValueError(
                "USE_MINI_HALOS is not compatible with the redshift-based"
-                " photon conservation corrections (PHOTON_CONS_TYPE=='z_photoncons')! "
+                " photon conservation corrections (PHOTON_CONS_TYPE=='z-photoncons')! "


Thanks for fixing this!

daviesje · 2025-07-20T23:47:58Z

src/py21cmfast/wrapper/inputs.py

+                    "This changes the escape fraction so it is not consistent with the manual setting of scaling."
+                    "Set PHOTON_CONS_TYPE to 'no-photoncons' or 'alpha-photoncons' if you want the scaling to be exact.",
+                    stacklevel=2,
+                )


Someone else may have an opinion on this, but if there is a real problem using these two models together then perhaps we should throw an error instead of a warning

I think an error will be better. Warnings are too easy to ignore.

I don't understand is what is happening when we turn on use_photoncons, but I see in the code f_esc10 gets modified per scaling relation if PHOTON_CONS_TYPE == 3. I will come back to it later when I understand the code better.

steven-murray · 2025-07-21T20:02:00Z

src/py21cmfast/wrapper/inputs.py

+    BETA_ESC = field(
+        default=0.0,
+        converter=float,
+    )


This will need to have a type!

Suggested change

BETA_ESC = field(

default=0.0,

converter=float,

)

BETA_ESC: float = field(

default=0.0,

converter=float,

)

daviesje

This is looking a lot better, but I still have two questions:

I know the Qin+25 model did not include minihalos at all, but would you like to include them? Either with the same scaling or another parameter? At the moment there is no scaling applied to them at all, and that might be confusing. If you wish not to include them, make it clear in the docstring for BETA_ESC
Did you try any runs with PHOTON_CONS=z-photoncons? I can't see anything that could go wrong here, since BETA_ESC should be applied to both the calibration and followup simulations. The warning/error check might not be necessary

steven-murray · 2026-01-13T21:53:40Z

@zhaotingchen do you have updates to @daviesje's comments above?

zhaotingchen added 4 commits July 19, 2025 05:16

feat: add dummy BETA_ESC to AstroParams

976e221

feat: add dummy BETAESC to input struct and print in debug

06b2112

feat: add BETAESC to halo

5ce8bea

fix: add warning for photocon and correct some typo

51a5b19

daviesje reviewed Jul 20, 2025

View reviewed changes

steven-murray reviewed Jul 21, 2025

View reviewed changes

zhaotingchen added 2 commits July 23, 2025 06:31

Merge branch 'main' into beta_esc

e5715bc

feat: refactor zesc into scaling relation

9b22eab

daviesje reviewed Aug 5, 2025

View reviewed changes

zhaotingchen added 4 commits August 30, 2025 19:47

feat: add beat_esc_mini

d84ba36

Merge branch 'master' into beta_esc

cec161f

Merge branch 'main' into beta_esc

53338d2

fix: debug logging and allow cheating single node z

32aecef

Beta esc #538

Are you sure you want to change the base?

Beta esc #538

Uh oh!

Conversation

zhaotingchen commented Jul 20, 2025

Uh oh!

zhaotingchen commented Jul 20, 2025

Uh oh!

daviesje left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

daviesje left a comment

Choose a reason for hiding this comment

Uh oh!

steven-murray commented Jan 13, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants