Should the Theory class know about the data bins already at the construction step, besides redshift bins?

[andreufont]

While reviewing lya_theory.get_Px_AA function, I noticed these lines around here:

        # compute input k, theta to emulator in Mpc
        
        Nk = 0
        Ntheta = 0
        if Nz > 1:
            for iz in range(Nz):
                if len(k_AA[iz]) > Nk:
                    Nk = len(k_AA[iz])
                if len(theta_deg[iz]) > Ntheta:
                    Ntheta = len(theta_deg[iz])
        else:
            
            if len(k_AA) == 1:
                k_AA = k_AA[0]
            if len(theta_deg) == 1:
                theta_deg = theta_deg[0]
            Nk = len(k_AA)
            k_AA = [k_AA]
            # if theta is just 1 float
            if np.isscalar(theta_deg):
                Ntheta = 1
            else:
                Ntheta = len(theta_deg)
            theta_deg = [theta_deg]

        kin_Mpc = np.zeros((Nz, Nk))
        theta_in_Mpc = np.zeros((Nz, Ntheta))
        for iz in range(Nz):
            kin_Mpc[iz, : Nk] = k_AA[iz] * M_AA_of_z[iz]
            theta_in_Mpc[iz, : Ntheta] = theta_deg[iz] / M_tdeg_of_z[iz]

It makes sense that the coordinate transformations have to be computed every time, but only when varying cosmology... otherwise, we are repeating these type of operations (and I believe similar ones in other places of the code) every time we evaluate a model.

But even if we do vary cosmology, the grid of theta and k bins, both fine and coarse, are constant properties of the data object that the theory could store (just like it stores the redshifts zs). One could consider setting all of these (in observed coordinates) at the moment we setup the theory class.

[marlena6]

To me it makes more sense that the theory object doesn’t know about the data object until its functions are called, so that one can set up the theory and ask for predictions for various k’s and theta’s using the same object with the same background cosmo, nuisance parameters, etc. For example, I sometimes use the same theory object to evaluate the likelihood for multiple different mocks with different setups. I think the likelihood could store these conversions right away, though, and then avoid repeated unit conversions as you mention . What do you think?

[andreufont]

There are probably several different ways of doing this correctly, we just need to think them through the entire pipeline.

Note, however, that you are already passing data-related information (like the redshift bins that you'll need) to the theory object, so it seems ok to provide more bin-related information.

In terms of storing the conversions (dkms_dMpc and such), as long as we never vary the cosmology it might be possible to store it in the likelihood object (instead of in the theory object), but in my experience is always better to store cache results right where you compute them, to avoid inconsistencies. So if it is the theory object that computes these (using camb_model), it should be the theory object that stores them.

An alternative option would be to have a function in camb_model that does something like:

def get_dkms_dMpc(self, z, like_params): 
    if not self.fixed_background(like_params)
        return self.compute_dkms_dMpc(z, like_params)
    else:
        # see if this was already computed before
        if z not in self.cached_dkms_dMpc:
            self.cached_dkms_dMpc[z] = self.compute_dkms_dMpc(z, like_params)
        return self.cached_dkms_dMpc[z]

[andreufont]

As discussed in the new appendix about the "fully contaminated model" in the draft, the units conversion is a bit more complicated once you want to model metals as well, so we should keep this in mind in this discussion in mind (writing it here before I forget)