Haloprop Correlation

[daviesje]

Note: I've written this description in a hurry, I intend to add to it after digging through my notes more, I'm also happy to answer questions
This PR serves as documentation and status for where I left off with the source model upgrade project. This is still unfinished in a few aspects but hopefully it's mostly just bugfixes, performance and cleanup.

The Idea: Self-Consistent Star-formation histories.

Determining SFR from stellar mass results in an inconsistency where integrating the SFR over time does not necessarily result in the stellar mass of the halo. this is exacerbated by the random sampling introduced in v4.

So instead of sampling the stellar mass, sfr and x-ray luminosity semi-independently, we sample the following three quantities:

• The SFR in the last 10 Myr
• The SFR in the last 100 Myr
• The SFR in the last snapshot

All three of these quantities are sampled from a single spectral-density, (Carvajal-Bohorquez et al. 2025). Covariance matrices are computed analytically at the start of each snapshot and then sampled for each halo (see https://arxiv.org/html/2507.14746v1#S2). We get a 6x6 full covariance matrix since we condition each of the three quantities on their descendant.

After sampling the quantities, the 10Myr and 100Myr SFR can be used to compute emissivities correlated with shorter and longer timescales. At the moment 10Myr is used for X-ray, and 100Myr is used for Gamma 12. The snapshot-SFRs are added together, resulting in the stellar mass.

Status

• All the matrix math in correlated_sfh.c is complete and tested:
• Below is an exaggerated example (high sigma setting) of a halo SFH history, imagine an RNG value of +1 as 1 sigma above mean etc:

- Results are reasonable, where 10Myr and snapshot (~30Myr max) are very similar, 100Myr is much more slowly varying, and closer to the mean. - The implementation into the rest of the code (HaloBox and afterwards) is unfinished. While all the components are there, it has not been run through an entire simulation, and likely contains many bugs - There is test code still present, much of the old source model (including parameters) has not been deleted - Since I needed access to multiple redshifts in the C code, I've created a structure similar to the parameters which holds the node redshifts in the backend. We don't need to pass the redshift directly to the backend anymore

TODOS

• Get the model running through HaloBox, use the PerturbHalofield to check halo emissivity histories.
• Check the inputs to Spintemp and reionisation algorithms, these should not have changed much.
• rewrite all the halo sampler and property tests to accomodate the new model
• clean up all the old code and parameters
• We need to consider performance more thoroughly, speed seems similar at the moment, we need to hang on to the previous halo fields a little longer though. There is one large allocation in HaloBox.c summing progenitor properties which I was meaning to find a way around.

Known Issues

• I believe there is a bug in the PSD timescale SFH_TAU, When I change it the SFHs don't respond the way I expect
• The first (dummy) halofield does not have sfh_computed set properly, causing a crash at the first call of ComputeHaloBox

Summary by Sourcery

Introduce a self-consistent, temporally correlated star-formation history model and propagate it through halo sampling, halo gridding, and radiative transfer pipelines, using globally broadcast snapshot metadata instead of per-call redshift plumbing.

New Features:

• Add a correlated SFH model based on a PSD parametrization, including covariance construction, sampling utilities, and backend test hooks.
• Extend astro- and input-parameter handling to support SFH PSD parameters and node redshift metadata, exposing a Python test helper for SFH correlations.

Bug Fixes:

• Ensure halo-box and global-mean calculations use consistent Hubble time and star-formation timescales derived from cosmology and astro-parameters.
• Fix several redshift-handling inconsistencies by deriving current, previous, and descendant redshifts from a shared snapshot index instead of function arguments.

Enhancements:

• Refactor halo property definitions to work with multi-timescale SFR (10 Myr, 100 Myr, snapshot-integrated) and track progenitor contributions for self-consistent stellar mass growth.
• Replace legacy stochastic RNG-based halo property assignment with SFH-based correlated sampling throughout the halo catalog, progenitor sampler, and mapping to grids.
• Unify and simplify C API signatures for core compute kernels (perturbed fields, halo catalogs, spin temperature, ionization, brightness temperature, halo boxes) to use internal redshift state instead of explicit arguments.
• Improve logging and debug output for halo and box averages, including new SFH-related quantities.
• Extend struct and FFI layers (C and Python) to support new halo catalog fields, 64-bit indices, and SFH bookkeeping, while cleaning up obsolete scatter/median-handling paths.

Documentation:

• Update Python-level driver and wrapper expectations around required halo catalog inputs, especially previous catalogs for discrete halo models and SFH-enabled runs.

Tests:

• Add backend SFH correlation tests that validate PSD-based correlation tables and covariance construction against analytic expectations.

Chores:

• Remove unused or redundant scaling-relation scatter plumbing and related helper functions now superseded by the SFH model.

[sourcery-ai]

Reviewer's Guide

Introduces a self-consistent, temporally correlated star-formation history (SFH) model driven by a new correlated_sfh backend, rewires halo/IGM pipelines to use node-based redshift broadcasting and SFH-based halo properties instead of independent RNG draws, and updates Python wrappers and structs accordingly while deprecating several redshift arguments and legacy stochastic parameters.

File-Level Changes

Change	Details	Files
Replace independent halo property RNGs with a correlated SFH sampler and propagate new SFR/stellar-mass fields through HaloCatalog, HaloBox, and mapping routines.	Add correlated_sfh.{c,h} implementing PSD-based SFH correlation functions, covariance construction, and Cholesky-based sampling across 10 Myr, 100 Myr, and snapshot scales, plus test hooks. Change HaloCatalog fields from star_rng/sfr_rng/xray_rng to sfr_10/sfr_100/stellar_mass plus descendant_index and sfh_computed, and adapt stochastic sampling (Stochasticity.c) to fill these using sample_correlated_sfh for both first-generation and progenitor halos. Refactor scaling_relations: replace get_halo_stellarmass/get_halo_sfr with get_halo_sfh that integrates mass growth per snapshot using correlated SFH draws, adds new mean-correction arrays in ScalingConstants, and adjust metallicity/X-ray calculations to use SFH outputs. Update map_mass move_halo_galprops/move_grid_galprops and HaloBox.set_halo_properties to consume the new SFH fields, compute progenitor-summed masses, and store SFR10/SFR100 and stellar masses in HaloBox and halos instead of legacy halo_sfr/sfr_mini.	src/py21cmfast/src/correlated_sfh.c src/py21cmfast/src/correlated_sfh.h src/py21cmfast/src/Stochasticity.c src/py21cmfast/src/Stochasticity.h src/py21cmfast/src/HaloCatalog.c src/py21cmfast/src/HaloCatalog.h src/py21cmfast/src/HaloBox.c src/py21cmfast/src/HaloBox.h src/py21cmfast/src/map_mass.c src/py21cmfast/src/map_mass.h src/py21cmfast/src/scaling_relations.c src/py21cmfast/src/scaling_relations.h src/py21cmfast/src/hmf.c src/py21cmfast/src/_outputstructs_wrapper.h
Introduce node-based snapshot redshift broadcasting and remove explicit redshift arguments from core C compute entrypoints and Python wrappers.	Define NodeRedshifts struct and Broadcast_snapshot_info/get_current_redshift/get_previous_redshift/get_descendant_redshift helpers in InputParameters.{c,h} and wire them through cfuncs.broadcast_input_struct. Remove redshift parameters from ComputePerturbedField, ComputeHaloCatalog, ComputePerturbedHaloCatalog, ComputeHaloBox, ComputeTsBox, ComputeIonizedBox, and ComputeBrightnessTemp prototypes and implementations, using get_*redshift accessors instead. Update Python OutputStruct wrappers and drivers (outputs.py, single_field.py, _param_config.py, cfuncs.py) to drop now-redundant redshift arguments, pass previous halo catalogs where needed, and ensure snapshot info is broadcast per call.	src/py21cmfast/src/InputParameters.c src/py21cmfast/src/InputParameters.h src/py21cmfast/src/_inputparams_wrapper.h src/py21cmfast/src/_functionprototypes_wrapper.h src/py21cmfast/src/PerturbedField.c src/py21cmfast/src/PerturbedField.h src/py21cmfast/src/HaloCatalog.c src/py21cmfast/src/HaloCatalog.h src/py21cmfast/src/PerturbedHaloCatalog.c src/py21cmfast/src/PerturbedHaloCatalog.h src/py21cmfast/src/SpinTemperatureBox.c src/py21cmfast/src/SpinTemperatureBox.h src/py21cmfast/src/IonisationBox.c src/py21cmfast/src/IonisationBox.h src/py21cmfast/src/BrightnessTemperatureBox.c src/py21cmfast/src/BrightnessTemperatureBox.h src/py21cmfast/wrapper/cfuncs.py src/py21cmfast/wrapper/outputs.py src/py21cmfast/drivers/single_field.py src/py21cmfast/drivers/_param_config.py
Adjust scaling relations, interpolation tables, and photon conservation to work with SFH-driven stochasticity and revised ScalingConstants.	Augment ScalingConstants with sampled_mean_correction/integral_mean_correction arrays, initialise them via correlated_sfh.get_current_vars in set_scaling_constants, and strip out scatter-related fields (sigma_star, sigma_sfr_*, sigma_xray) and mimic_scatter_in_consts. Update SFRD/Nion/X-ray interpolation and evaluation routines (interp_tables.c, photoncons.c, hmf.c, HaloBox.c) to call evolve_scaling_constants_to_redshift without photoncons flags, rely on sc->t_h/t_star for time-scales, and use the new correction factors when building grids and integrals. Modify get_halo_xray to use the SFH-based sampled_mean_correction and per-halo RNG rather than sigma_xray-based lognormal scatter. Ensure UHMF and halobox average calculations now map to sfr_10/sfr_10_mcg instead of legacy halo_sfr/sfr_mini fields and update debug logging accordingly.	src/py21cmfast/src/scaling_relations.c src/py21cmfast/src/scaling_relations.h src/py21cmfast/src/interp_tables.c src/py21cmfast/src/photoncons.c src/py21cmfast/src/hmf.c src/py21cmfast/src/HaloBox.c src/py21cmfast/src/map_mass.c
Extend cosmology utilities with a robust time_between_z integrator and reuse it in SFH time-scale handling.	Add time_between_z(z_low,z_high) in cosmology.c using gsl_integration of dtdz, and export it from cosmology.h. Switch various time-scale computations to use time_between_z and double-precision dtdz; adjust comments and TODOs accordingly.	src/py21cmfast/src/cosmology.c src/py21cmfast/src/cosmology.h
Expose new SFH and halo-catalog configuration parameters and C test hooks to Python, and update struct type mappings.	Extend AstroParams with SFH_TAU and SFH_INDEX fields in both C and Python (inputs.py, _inputparams_wrapper.h) and use them in the SFH PSD and correlation functions. Update StructWrapper._TYPEMAP to support int64 arrays, enabling descendant_index on HaloCatalog; wire descendant_index and sfh_computed through _outputstructs_wrapper.h and wrapper/outputs.HaloCatalog. Add cfuncs.test_sfh_print wrapper for backend test_sfh_corr, and modify testing helpers (integral_wrappers.c, Stochasticity.c single_test_sample) to use the new stoc_set_consts_z signature and SFH initialisation. Clean up or update various debug and log messages to reflect new SFR10/SFR100 and snapshot semantics.	src/py21cmfast/wrapper/inputs.py src/py21cmfast/src/_inputparams_wrapper.h src/py21cmfast/wrapper/structs.py src/py21cmfast/src/_outputstructs_wrapper.h src/py21cmfast/wrapper/cfuncs.py src/py21cmfast/src/integral_wrappers.c src/py21cmfast/src/Stochasticity.c

---

Tips and commands

Interacting with Sourcery

• Trigger a new review: Comment @sourcery-ai review on the pull request.
• Continue discussions: Reply directly to Sourcery's review comments.
• Generate a GitHub issue from a review comment: Ask Sourcery to create an
  issue from a review comment by replying to it. You can also reply to a
  review comment with @sourcery-ai issue to create an issue from it.
• Generate a pull request title: Write @sourcery-ai anywhere in the pull
  request title to generate a title at any time. You can also comment
  @sourcery-ai title on the pull request to (re-)generate the title at any time.
• Generate a pull request summary: Write @sourcery-ai summary anywhere in
  the pull request body to generate a PR summary at any time exactly where you
  want it. You can also comment @sourcery-ai summary on the pull request to
  (re-)generate the summary at any time.
• Generate reviewer's guide: Comment @sourcery-ai guide on the pull
  request to (re-)generate the reviewer's guide at any time.
• Resolve all Sourcery comments: Comment @sourcery-ai resolve on the
  pull request to resolve all Sourcery comments. Useful if you've already
  addressed all the comments and don't want to see them anymore.
• Dismiss all Sourcery reviews: Comment @sourcery-ai dismiss on the pull
  request to dismiss all existing Sourcery reviews. Especially useful if you
  want to start fresh with a new review - don't forget to comment
  @sourcery-ai review to trigger a new review!

Customizing Your Experience

Access your dashboard to:

• Enable or disable review features such as the Sourcery-generated pull request
  summary, the reviewer's guide, and others.
• Change the review language.
• Add, remove or edit custom review instructions.
• Adjust other review settings.

Getting Help

• Contact our support team for questions or feedback.
• Visit our documentation for detailed guides and information.
• Keep in touch with the Sourcery team by following us on X/Twitter, LinkedIn or GitHub.

[codecov]

⚠️ JUnit XML file not found

The CLI was unable to find any JUnit XML files to upload.
For more help, visit our troubleshooting guide.

[steven-murray]

@daviesje thanks for this!

One question -- you mention that the reionization module should stay pretty similar. However, as you may see in #516, there is currently a problem with how we calculate reionization, due to the way cumulative recombinations are handled on short time-steps. We thought that perhaps your new source model might allow computing this based on the differential ionisation rate, dnion/dt rather than the cumulative nion. Do you have any thoughts on this?

[daviesje]

@daviesje thanks for this!

One question -- you mention that the reionization module should stay pretty similar. However, as you may see in #516, there is currently a problem with how we calculate reionization, due to the way cumulative recombinations are handled on short time-steps. We thought that perhaps your new source model might allow computing this based on the differential ionisation rate, dnion/dt rather than the cumulative nion. Do you have any thoughts on this?

I don't think this model specifically "allows" that change (none of the techincal details really overlap), but having more consistent SF histories might avoid some pitfalls that could occur imnplementing it

[jordanflitter]

Thanks @daviesje, this sounds really cool!

Could you share with us your notes (with equations)? Having them would make it easier to understand better the idea of this PR and review it. Plus, they could be the basis of any future paper that will come out of this work.

[daviesje]

Thanks @daviesje, this sounds really cool!

Could you share with us your notes (with equations)? Having them would make it easier to understand better the idea of this PR and review it. Plus, they could be the basis of any future paper that will come out of this work.

https://www.wolframcloud.com/obj/2d233119-1277-40e1-a3b2-a8e561a58851

[steven-murray]

Hi @daviesje -- thanks for this work again. I just wanted to make one general comment: you've added a new global node redshift array, but we're actively trying to diminish the number of globals we are dealing with. Passing in the full array of node redshifts rather than just a single redshift is great (and we were planning on doing this already!), but can we make it so that we pass these in as a parameter, rather than a global?

[daviesje]

@steven-murray of course! Effectively I've added another parameter structure, minus some of the class logic. I guess when you're looking to reduce globals it should follow the whatever pattern you implement for the rest of the InputParameters classes.