eelregit / pmwd

Differentiable Cosmological Forward Model
BSD 3-Clause "New" or "Revised" License
69 stars 16 forks source link

A_s -\sigma_8 utiliy function #6

Closed modichirag closed 1 year ago

modichirag commented 1 year ago

Should we include a utility function to switch back and forth between $A_s$ and $\sigma_8$, maybe at the configuration level?

Often it is the case that we know the cosmology in terms of $\sigma_8$ and not $A_s$, in which case it might be a good option to have a way of letting the user specify config in terms of $\sigma_8$ while the conversion happens under the hood.

In the similar vein, once we generate samples (during inference) in terms of $A_s$, have a way to convert them to $\sigma_8$.

eelregit commented 1 year ago

I wonder if fftlog can be used to compute this oscillatory integral. See the last example of https://github.com/eelregit/mcfit

Actually , we need to cache $\sigma^2(R)$ in cosmo in the sto project, which includes $\sigma_8$. For that fftlog is necessary.

BTW, @modichirag do you have some spline utility in the growth_mlp branch that can be used for interpolation?

modichirag commented 1 year ago

This might be irrelevant now...but yes https://github.com/eelregit/pmwd/blob/growth_mlp/pmwd/growth_integrals.py#L53

This is taken from Francois' code https://github.com/EiffL/NeuralBSpline

I was using mcfit to do the integral in python, but in jax, I have implemented it in the sigma8 branch using romberg integral which seems to be accurate enough (tested against mcfit) https://github.com/eelregit/pmwd/blob/sigma8/pmwd/boltzmann.py#L333 This is also how jax_cosmo does it apparently.

eelregit commented 1 year ago

Now with #20, sigma8 can be computed by

from pmwd import boltzmann, varlin_integ, varlin

cosmo = boltzmann(cosmo, conf)  # or cosmo = varlin_integ(cosmo, conf)
R, a = 8, 1
sigma8 = jnp.sqrt(varlin(R, a, cosmo, conf))  # if conf.L is 1Mpc/h
eelregit commented 1 year ago
cosmo = boltzmann(cosmo, conf)  # or cosmo = varlin_integ(cosmo, conf)
cosmo.sigma8
eelregit commented 1 year ago

Added Cosmology.from_sigma8 constructor:

print(SimpleLCDM(conf).sigma8)  # 0.8121139683320571

print(Cosmology.from_sigma8(conf, 0.8121139683320571, 0.96, 0.3, 0.05, 0.7).A_s)  # 2e-9