LSSTDESC / CCL

DESC Core Cosmology Library: cosmology routines with validated numerical accuracy
BSD 3-Clause "New" or "Revised" License
133 stars 64 forks source link

Log definitions in Zheng HOD #1143

Closed cheymans closed 6 months ago

cheymans commented 7 months ago

Hi CCL-crew,

In hod.py the Zheng HOD for centrals is defined using $\log_e$. Lines 386-390

    def _Nc(self, M, a):
        # Number of centrals
        Mmin = 10.**(self.log10Mmin_0 + self.log10Mmin_p * (a - self.a_pivot))
        siglnM = self.siglnM_0 + self.siglnM_p * (a - self.a_pivot)
        return 0.5 * (1 + erf(np.log(M/Mmin)/siglnM))

encode this:

Screenshot 2023-12-06 at 12 56 02

with some extra bells and whistles.

We've always used $\log_{10}$ here i.e:

Screenshot 2023-12-06 at 12 57 05

And going back to the original Zheng et al 2005 paper, I think this is correct as they have written just before eqn 3 "All logarithms in this paper are base-10."

This is how the two compare image

It's probably an academic discussion as few use the Zheng HOD any more - but as this mismatch lost me a couple of hours in a benchmarking comparison I thought I'd highlight it to you!

tilmantroester commented 6 months ago

The HOD model in CCL is benchmarked against https://github.com/damonge/pysz_gal, which uses the natural log here: https://github.com/damonge/pysz_gal/blob/master/pysz_gal/source/hod.f90#L13. pysz_gal says it's benchmarked against class_sz, but class_sz uses base-10 logs for the HOD: https://github.com/CLASS-SZ/class_sz/blob/master/source/class_sz.c#L21006 @damonge @ryumakiya

damonge commented 6 months ago

Thanks both. Yes, indeed, the key parameter here is sigma_logM, and whether it will operate on natural or base-10 logarithms. For historical reasons (e.g. the code is based on what was used in https://arxiv.org/abs/1912.08209), we followed the prescription of https://arxiv.org/pdf/1706.05422.pdf, which assumes natural logarithms. Mathematicians and theoretical physicists would argue (I think generally) that log = log_e, and otherwise you should specify the base, which is what we based ourselves on here. This is probably one of the most frustrating convention-related sources of confusion that exist.

For this reason, to minimise potential confusion, we called the parameter siglnM, and since it's a keyword-only argument, folks need to actually spell it out! However, we should probably add a more visible warning in the docs about this... I think we can do that in #1122.

Sorry that this caused so much of a hassle @cheymans !