# HMcode
A pure-`Python` implementation of the `HMcode-2020` method ([Mead et al. 2021](https://arxiv.org/abs/2009.01858)) for computing accurate non-linear matter power spectra across a wide range of cosmological parameters for $w(a)$-CDM models including curvature and massive neutrinos.
## Installation
Either
```bash
> pip install hmcode
```
or, if you want to edit the code, use the demo notebook, or run the tests or consistency checks, then clone the repository, `cd` into the directory, and then
```bash
> poetry install
```
to create a virtual environment with everything you need to get going.
## Dependencies
- `numpy`
- `scipy`
- `camb`
- `pyhalomodel`
## Use
```
import numpy as np
import camb
import hmcode
# Ranges
k = np.logspace(-3, 1, 100) # Wavenumbers [h/Mpc]
zs = [3., 2., 1., 0.5, 0.] # Redshifts
T_AGN = 10**7.8 # Feedback temperature [K]
# Run CAMB
parameters = camb.CAMBparams(WantCls=False)
parameters.set_cosmology(H0=70.)
parameters.set_matter_power(redshifts=zs, kmax=100.) # kmax should be much larger than the wavenumber of interest
results = camb.get_results(parameters)
# HMcode
Pk = hmcode.power(k, zs, results, T_AGN)
```
## Note
To whom it may concern,
I coded this `Python` version of `HMcode-2020` ([Mead et al. 2021](https://arxiv.org/abs/2009.01858)) up quite quickly before leaving academia in February 2023. It is written in pure `Python` and doesn't use any of the original Fortran code whatsoever. There is something amazing/dispiriting about coding something up in 3 days that previously took 5 years. A tragic last hoorah! At least I switched to `Python` eventually...
You might also be interested in [`pyhalomodel`](https://pypi.org/project/pyhalomodel/), upon which this code depends, which implements a vanilla halo-model calculation for any desired large-scale-structure tracer. Alternatively, and very confusingly, you might be interested in this [`pyhmcode`](https://pypi.org/project/pyhmcode/), which provides a wrapper around the original `Fortran` `HMcode` implementation.
I compared it against the `CAMB-HMcode` version for 100 random sets of cosmological parameters ($k < 10 h\mathrm{Mpc}^{-1}$; $z < 3$). The level of agreement between the two codes is as follows:
- LCDM: Mean error: 0.10%; Standard deviation of error: 0.03%; Worst-case error; 0.28%
- k-LCDM: Mean error: 0.12%; Standard deviation of error: 0.04%; Worst-case error; 0.32%
- w-CDM: Mean error: 0.11%; Standard deviation of error: 0.03%; Worst-case error; 0.31%
- w(a)-CDM: Mean error: 0.14%; Standard deviation of error: 0.08%; Worst-case error; 0.77%
- nu-LCDM: Mean error: 0.45%; Standard deviation of error: 0.43%; Worst-case error; 1.97% (larger errors strongly correlated with neutrino mass)
- nu-k-w(a)-CDM: Mean error: 0.41%; Standard deviation of error: 0.42%; Worst-case error; 1.99% (larger errors strongly correlated with neutrino mass)
These comparisons can be reproduced using the `comparisons/CAMB.py` script.
The power-spectrum suppression caused by baryonic feedback has also been implemented, and the level of agreement between the *suppression* predicted by this code and the same implementation in `CAMB` is as follows:
- LCDM: Mean error: 0.02%; Standard deviation of error: <0.01%; Worst-case error; 0.03%
- nu-k-w(a)-CDM: Mean error: 0.06%; Standard deviation of error: 0.07%; Worst-case error; 0.49% (larger errors strongly correlated with neutrino mass)
The comparisons can be reproduced using the `comparisons/CAMB_feedback.py` script.
While the quoted accuracy of `HMcode-2020` relative to simulations is RMS ~2.5%, note that the accuracy is anti-correlated with neutrino masses (cf. Fig. 2 of [Mead et al. 2021](https://arxiv.org/abs/2009.01858)). The larger discrepancies between the codes for massive neutrinos (2% for ~1eV) may seem worrisome, but here are some reasons why I am not that worried:
- Here, neutrinos are treated as completely cold matter when calculating the linear growth factors, whereas in `CAMB-HMcode` the transition from behaving like radiation to behaving like matter is accounted for in the linear growth.
- Here the *cold* matter power spectrum is taken directly from `CAMB` whereas in `CAMB-HMcode` the *cold* spectrum is calculated approximately from the total matter power spectrum using approximations for the scale-dependent growth rate from [Eisenstein & Hu (1999)](https://arxiv.org/abs/astro-ph/9710252).
- If I resort to the old approximation for the *cold* matter power spectrum (and therefore the cold $\sigma(R)$) then the level of agreement between the codes for nu-k-w(a)-CDM improves to: Mean error: 0.15%; Std error: 0.11%; Worst error; 1.15%.
Using the actual cold matter spectrum is definitely an improvement from a theoretical perspective (and for speed), so I decided to keep that at the cost of the disagreement between this code at `CAMB-HMcode` for models with very high neutrino mass. If better agreement between this code and `CAMB-HMcode` is important to you then the old approximate cold spectrum can be used by changing the `sigma_cold_approx` flag at the top of `hmcode.py`. Note that while ignoring the actual energy-density scaling of massive neutrinos might seem to be a (small) problem, keep in mind the comments below regarding the linear growth factor.
I think any residual differences between codes must therefore stem from:
- The BAO de-wiggling process (different `k` grids)
- The $\sigma_\mathrm{v}$ numerical integration
- The $n_\mathrm{eff}$ calculation (numerical differentiation here; numerical integration in `CAMB-HMcode`)
- The $\sigma(R)$ numerical integration (using `CAMB` here; done internally in `CAMB-HMcode`)
- The linear growth ODE solution
- Root finding for the halo-collapse redshift and for $R_\mathrm{nl}$
But I didn't have time to investigate these differences more thoroughly. Note that there are accuracy parameters in `CAMB-HMcode` fixed at the $10^{-4}$ level, so you would never expect better than 0.01% agreement. Given that `HMcode` is only accurate at the ~2.5% level compared to simulated power spectra, the level of agreement between the codes seems okay to me, with the caveats above regarding very massive neutrinos.
While writing this code I had a few ideas for future improvements:
- ~~Add the `HMcode-2020` baryon-feedback model; this would not be too hard for the enthusiastic student/postdoc.~~ (Thanks Laila Linke!)
- The predictions are a bit sensitive to the smoothing $\sigma$ used for the dewiggling. This should probably be a fitted parameter.
- It's annoying having to calculate linear growth functions (all, LCDM), especially since the linear growth doesn't really exist. One should probably use the $P(k)$ amplitude evolution over time at some cleverly chosen scale instead, or instead the evolution of $\sigma(R)$ over time at some pertinent $R$ value. Note that the growth factors are *only* used to calculate the [Dolag et al. (2004)](https://arxiv.org/abs/astro-ph/0309771) correction and [Mead (2017)](https://arxiv.org/abs/1606.05345) $\delta_\mathrm{c}$, $\Delta_\mathrm{v}$.
- I never liked the halo bloating parameter, it's hard to understand the effect of modifying halo profiles in Fourier space. Someone should get rid of this (maybe modify the halo mass function instead?).
- Redshift 'infinity' for the Dolag correction is actually $z_\infty = 10$. `HMcode` predictions *are* sensitive to this (particularly w(a)CDM). Either the redshift 'infinity' should be fitted or the halo-concentration model for beyond-LCDM should be improved.
- The massive neutrino correction for the [Mead (2017)](https://arxiv.org/abs/1606.05345) $\delta_\mathrm{c}$, $\Delta_\mathrm{v}$ formulae (appendix A of [Mead et al. 2021](https://arxiv.org/abs/2009.01858)) is crude and should be improved. I guess using the intuition that hot neutrinos are ~smoothly distributed on halo scales. Currently neutrinos are treated as cold matter in the linear/accumulated growth calculation (used by [Mead 2017](https://arxiv.org/abs/1606.05345)), which seems a bit wrong.
- I haven't checked how fast this code is, but there are a couple of TODO in the code that might improve speed if necessary.
- The choices regarding how to account for massive neutrinos could usefully be revisited. This whole subject is a bit confusing and the code doesn't help to alleviate the confusion. Choices like what to use for: $\delta_\mathrm{c}$; $\Delta_\mathrm{v}$; $\sigma(R)$; $R_\mathrm{nl}$; $n_\mathrm{eff}$; $c(M)$.
- Refit model (including $\sigma$ for BAO smoothing and $z_\infty$ for [Dolag et al. 2004](https://arxiv.org/abs/astro-ph/0309771)) to new emulator(s) (e.g., [Mira Titan IV](https://arxiv.org/abs/2207.12345)).
- Don't be under any illusions that the `HMcode` parameters, or the forms of their dependence on the underlying power spectrum, are special in any particular way. A lot of experimentation went into finding these, but it was by no means exhaustive. However, please note that obviously these parameters should only depend on the underlying linear spectrum (rather than being random functions of $z$, $\Omega_\mathrm{m}$, $w$, or whatever).
Have fun,
Alexander Mead (2023/02/28)
Raw data
{
"_id": null,
"home_page": "https://github.com/alexander-mead/HMcode-python",
"name": "hmcode",
"maintainer": null,
"docs_url": null,
"requires_python": "<3.13,>=3.9",
"maintainer_email": null,
"keywords": null,
"author": "Alexander Mead",
"author_email": "alexander.j.mead@googlemail.com",
"download_url": "https://files.pythonhosted.org/packages/50/2c/4e099610e51230878f05f0b545629eedc1b01f6067b61a2c9ae87af95455/hmcode-1.1.1.tar.gz",
"platform": null,
"description": "# HMcode\n\n\n\nA pure-`Python` implementation of the `HMcode-2020` method ([Mead et al. 2021](https://arxiv.org/abs/2009.01858)) for computing accurate non-linear matter power spectra across a wide range of cosmological parameters for $w(a)$-CDM models including curvature and massive neutrinos.\n\n## Installation\n\nEither\n```bash\n> pip install hmcode\n```\nor, if you want to edit the code, use the demo notebook, or run the tests or consistency checks, then clone the repository, `cd` into the directory, and then\n```bash\n> poetry install\n```\nto create a virtual environment with everything you need to get going.\n\n## Dependencies\n\n- `numpy`\n- `scipy`\n- `camb`\n- `pyhalomodel`\n\n## Use\n\n```\nimport numpy as np\nimport camb\nimport hmcode\n\n# Ranges\nk = np.logspace(-3, 1, 100) # Wavenumbers [h/Mpc]\nzs = [3., 2., 1., 0.5, 0.] # Redshifts\nT_AGN = 10**7.8 # Feedback temperature [K]\n\n# Run CAMB\nparameters = camb.CAMBparams(WantCls=False)\nparameters.set_cosmology(H0=70.)\nparameters.set_matter_power(redshifts=zs, kmax=100.) # kmax should be much larger than the wavenumber of interest\nresults = camb.get_results(parameters)\n\n# HMcode\nPk = hmcode.power(k, zs, results, T_AGN)\n```\n\n## Note\n\nTo whom it may concern,\n\nI coded this `Python` version of `HMcode-2020` ([Mead et al. 2021](https://arxiv.org/abs/2009.01858)) up quite quickly before leaving academia in February 2023. It is written in pure `Python` and doesn't use any of the original Fortran code whatsoever. There is something amazing/dispiriting about coding something up in 3 days that previously took 5 years. A tragic last hoorah! At least I switched to `Python` eventually...\n\nYou might also be interested in [`pyhalomodel`](https://pypi.org/project/pyhalomodel/), upon which this code depends, which implements a vanilla halo-model calculation for any desired large-scale-structure tracer. Alternatively, and very confusingly, you might be interested in this [`pyhmcode`](https://pypi.org/project/pyhmcode/), which provides a wrapper around the original `Fortran` `HMcode` implementation.\n\nI compared it against the `CAMB-HMcode` version for 100 random sets of cosmological parameters ($k < 10 h\\mathrm{Mpc}^{-1}$; $z < 3$). The level of agreement between the two codes is as follows:\n- LCDM: Mean error: 0.10%; Standard deviation of error: 0.03%; Worst-case error; 0.28%\n- k-LCDM: Mean error: 0.12%; Standard deviation of error: 0.04%; Worst-case error; 0.32%\n- w-CDM: Mean error: 0.11%; Standard deviation of error: 0.03%; Worst-case error; 0.31%\n- w(a)-CDM: Mean error: 0.14%; Standard deviation of error: 0.08%; Worst-case error; 0.77%\n- nu-LCDM: Mean error: 0.45%; Standard deviation of error: 0.43%; Worst-case error; 1.97% (larger errors strongly correlated with neutrino mass)\n- nu-k-w(a)-CDM: Mean error: 0.41%; Standard deviation of error: 0.42%; Worst-case error; 1.99% (larger errors strongly correlated with neutrino mass)\n\nThese comparisons can be reproduced using the `comparisons/CAMB.py` script.\n\nThe power-spectrum suppression caused by baryonic feedback has also been implemented, and the level of agreement between the *suppression* predicted by this code and the same implementation in `CAMB` is as follows:\n- LCDM: Mean error: 0.02%; Standard deviation of error: <0.01%; Worst-case error; 0.03%\n- nu-k-w(a)-CDM: Mean error: 0.06%; Standard deviation of error: 0.07%; Worst-case error; 0.49% (larger errors strongly correlated with neutrino mass)\n\nThe comparisons can be reproduced using the `comparisons/CAMB_feedback.py` script.\n\nWhile the quoted accuracy of `HMcode-2020` relative to simulations is RMS ~2.5%, note that the accuracy is anti-correlated with neutrino masses (cf. Fig. 2 of [Mead et al. 2021](https://arxiv.org/abs/2009.01858)). The larger discrepancies between the codes for massive neutrinos (2% for ~1eV) may seem worrisome, but here are some reasons why I am not that worried:\n- Here, neutrinos are treated as completely cold matter when calculating the linear growth factors, whereas in `CAMB-HMcode` the transition from behaving like radiation to behaving like matter is accounted for in the linear growth.\n- Here the *cold* matter power spectrum is taken directly from `CAMB` whereas in `CAMB-HMcode` the *cold* spectrum is calculated approximately from the total matter power spectrum using approximations for the scale-dependent growth rate from [Eisenstein & Hu (1999)](https://arxiv.org/abs/astro-ph/9710252).\n- If I resort to the old approximation for the *cold* matter power spectrum (and therefore the cold $\\sigma(R)$) then the level of agreement between the codes for nu-k-w(a)-CDM improves to: Mean error: 0.15%; Std error: 0.11%; Worst error; 1.15%.\n\nUsing the actual cold matter spectrum is definitely an improvement from a theoretical perspective (and for speed), so I decided to keep that at the cost of the disagreement between this code at `CAMB-HMcode` for models with very high neutrino mass. If better agreement between this code and `CAMB-HMcode` is important to you then the old approximate cold spectrum can be used by changing the `sigma_cold_approx` flag at the top of `hmcode.py`. Note that while ignoring the actual energy-density scaling of massive neutrinos might seem to be a (small) problem, keep in mind the comments below regarding the linear growth factor.\n\nI think any residual differences between codes must therefore stem from:\n- The BAO de-wiggling process (different `k` grids)\n- The $\\sigma_\\mathrm{v}$ numerical integration\n- The $n_\\mathrm{eff}$ calculation (numerical differentiation here; numerical integration in `CAMB-HMcode`)\n- The $\\sigma(R)$ numerical integration (using `CAMB` here; done internally in `CAMB-HMcode`)\n- The linear growth ODE solution\n- Root finding for the halo-collapse redshift and for $R_\\mathrm{nl}$\n\nBut I didn't have time to investigate these differences more thoroughly. Note that there are accuracy parameters in `CAMB-HMcode` fixed at the $10^{-4}$ level, so you would never expect better than 0.01% agreement. Given that `HMcode` is only accurate at the ~2.5% level compared to simulated power spectra, the level of agreement between the codes seems okay to me, with the caveats above regarding very massive neutrinos.\n\nWhile writing this code I had a few ideas for future improvements:\n- ~~Add the `HMcode-2020` baryon-feedback model; this would not be too hard for the enthusiastic student/postdoc.~~ (Thanks Laila Linke!)\n- The predictions are a bit sensitive to the smoothing $\\sigma$ used for the dewiggling. This should probably be a fitted parameter.\n- It's annoying having to calculate linear growth functions (all, LCDM), especially since the linear growth doesn't really exist. One should probably use the $P(k)$ amplitude evolution over time at some cleverly chosen scale instead, or instead the evolution of $\\sigma(R)$ over time at some pertinent $R$ value. Note that the growth factors are *only* used to calculate the [Dolag et al. (2004)](https://arxiv.org/abs/astro-ph/0309771) correction and [Mead (2017)](https://arxiv.org/abs/1606.05345) $\\delta_\\mathrm{c}$, $\\Delta_\\mathrm{v}$.\n- I never liked the halo bloating parameter, it's hard to understand the effect of modifying halo profiles in Fourier space. Someone should get rid of this (maybe modify the halo mass function instead?).\n- Redshift 'infinity' for the Dolag correction is actually $z_\\infty = 10$. `HMcode` predictions *are* sensitive to this (particularly w(a)CDM). Either the redshift 'infinity' should be fitted or the halo-concentration model for beyond-LCDM should be improved.\n- The massive neutrino correction for the [Mead (2017)](https://arxiv.org/abs/1606.05345) $\\delta_\\mathrm{c}$, $\\Delta_\\mathrm{v}$ formulae (appendix A of [Mead et al. 2021](https://arxiv.org/abs/2009.01858)) is crude and should be improved. I guess using the intuition that hot neutrinos are ~smoothly distributed on halo scales. Currently neutrinos are treated as cold matter in the linear/accumulated growth calculation (used by [Mead 2017](https://arxiv.org/abs/1606.05345)), which seems a bit wrong.\n- I haven't checked how fast this code is, but there are a couple of TODO in the code that might improve speed if necessary.\n- The choices regarding how to account for massive neutrinos could usefully be revisited. This whole subject is a bit confusing and the code doesn't help to alleviate the confusion. Choices like what to use for: $\\delta_\\mathrm{c}$; $\\Delta_\\mathrm{v}$; $\\sigma(R)$; $R_\\mathrm{nl}$; $n_\\mathrm{eff}$; $c(M)$.\n- Refit model (including $\\sigma$ for BAO smoothing and $z_\\infty$ for [Dolag et al. 2004](https://arxiv.org/abs/astro-ph/0309771)) to new emulator(s) (e.g., [Mira Titan IV](https://arxiv.org/abs/2207.12345)).\n- Don't be under any illusions that the `HMcode` parameters, or the forms of their dependence on the underlying power spectrum, are special in any particular way. A lot of experimentation went into finding these, but it was by no means exhaustive. However, please note that obviously these parameters should only depend on the underlying linear spectrum (rather than being random functions of $z$, $\\Omega_\\mathrm{m}$, $w$, or whatever).\n\nHave fun,\n\nAlexander Mead (2023/02/28)\n",
"bugtrack_url": null,
"license": "MIT",
"summary": "Pure Python implementation of HMcode",
"version": "1.1.1",
"project_urls": {
"Homepage": "https://github.com/alexander-mead/HMcode-python",
"Repository": "https://github.com/alexander-mead/HMcode-python"
},
"split_keywords": [],
"urls": [
{
"comment_text": "",
"digests": {
"blake2b_256": "7e2ac05788f9590b5b3802d5949be6607ae365b8715b03bffe911aa375142482",
"md5": "c98d42e66b94783c8c846f44701c169c",
"sha256": "02c7e48572e801bae9609411b31161609c3ee67e4892f272ed696b7f226d3a46"
},
"downloads": -1,
"filename": "hmcode-1.1.1-py3-none-any.whl",
"has_sig": false,
"md5_digest": "c98d42e66b94783c8c846f44701c169c",
"packagetype": "bdist_wheel",
"python_version": "py3",
"requires_python": "<3.13,>=3.9",
"size": 18749,
"upload_time": "2024-04-26T21:39:40",
"upload_time_iso_8601": "2024-04-26T21:39:40.166674Z",
"url": "https://files.pythonhosted.org/packages/7e/2a/c05788f9590b5b3802d5949be6607ae365b8715b03bffe911aa375142482/hmcode-1.1.1-py3-none-any.whl",
"yanked": false,
"yanked_reason": null
},
{
"comment_text": "",
"digests": {
"blake2b_256": "502c4e099610e51230878f05f0b545629eedc1b01f6067b61a2c9ae87af95455",
"md5": "4a6802951dc0ece706790a34bb76952d",
"sha256": "1f745a1c400d669afab5cf8d32a4f0db64484d98e4f318ff028bb77b2922b3a8"
},
"downloads": -1,
"filename": "hmcode-1.1.1.tar.gz",
"has_sig": false,
"md5_digest": "4a6802951dc0ece706790a34bb76952d",
"packagetype": "sdist",
"python_version": "source",
"requires_python": "<3.13,>=3.9",
"size": 19477,
"upload_time": "2024-04-26T21:39:42",
"upload_time_iso_8601": "2024-04-26T21:39:42.009328Z",
"url": "https://files.pythonhosted.org/packages/50/2c/4e099610e51230878f05f0b545629eedc1b01f6067b61a2c9ae87af95455/hmcode-1.1.1.tar.gz",
"yanked": false,
"yanked_reason": null
}
],
"upload_time": "2024-04-26 21:39:42",
"github": true,
"gitlab": false,
"bitbucket": false,
"codeberg": false,
"github_user": "alexander-mead",
"github_project": "HMcode-python",
"travis_ci": false,
"coveralls": false,
"github_actions": false,
"lcname": "hmcode"
}
JSON
