clarify role of F_model_emcee in lya_theory

[andreufont]

We didn't have this before... is this necessary? If so, let's at least change its name to something more clear

[jchavesmontero]

We vary the parameters of X_model_emcee during the fits (minimizer or sampler), but we hold fixed the values of the X_model_fid parameters for comparison. Any suggestion for the name? X_model_optimize maybe?

[andreufont]

But why do we need to store as a variable of the LyaTheory class any other model than the fiducial?

I honestly don't remember if I changed this myself, but in the past we had a fiducial model that was a member variable and that would never change, and then we would create new model objects on the fly (using get_igm_models function) whenever we needed to make predictions as a function of likelihood parameters. Was this inefficient?

If for computational efficiency you wanted to have two sets of IGM models (the fiducial one and "another one") it would be good to clarify very well the difference and have proper documentation everywhere any of these are used. Would the other model stored be the current best guess (from the minimizer?) or would it be whatever was the last point evaluated? This seems dangerous.

It might be easier to discuss this in person!

[jchavesmontero]

Exactly, we decided to do this because creating new objects on the fly takes much more time than updating the ones we already had. I will clarify this. In the minimizer we also update the models. It is not dangerous since these are like dummy variables, every time we want to plot or store data we evaluate these again with the new parameters. We could simply could these X_dummy

[andreufont]

Multiprocessing is not a problem? In C++ this would be a problem when using OpenMP, but I don't have a good intuition for Python multiprocessing.

I will try to find time to investigate whether we need these objects at all (we could use the fiducial models objects and evaluate them with external likelihood parameters), but I will open a separate and more clear issue to discuss this.

Meanwhile, yes having very clear names for the different objects and clear documentation when we use them should be enough. Possible names:

- X_model_dummy
- X_model_free
- X_model_vary
- X_model

[jchavesmontero]

I have no idea about multiprocessing, but it is not a problem for mpi. Each thread deals with its own stuff so there are no conflicts.

[andreufont]

Sure, MPI should not be a problem. I thought we were using multiprocessing (at least Chris I think was), but if we are not using this then this should not be a problem.

[jchavesmontero]

We were using multiprocessing before, but it was not working as well as expected in NERSC (I asked around and nobody gave me a good solution) so I decided to bite the bullet and implement MPI

[andreufont]

I thought multiprocessing or OpenMP would be good options to avoid needing to load 128 archives in a single node at NERSC, but if this is not an issue we can continue with pure MPI.

[jchavesmontero]

This is what I meant about performance. I did a scaling test of the performance as a function of the number of threads in OpenMP and it was quite bad. On the other hand, the performance goes as expected for MPI.

As for the loading, this is not how the code works. We only have a rank reading the data, and this rank redistributes the data among the others for the sampling. Therefore, the number of requested ranks only changes how many of these go into the sampling, nothing else. Maybe now you understand why I was so reluctant to keep a copy of the whole archive in the emulator class.

[andreufont]

I was reviewing this again, and I'm not sure why we need the fiducial models in the current code. When are they used? And why do we need the get_IGM_models function now?

[andreufont]

Looking at it again, we should review the role of "models" and "parameters" in cup1d.

For instance, now we do something like (pseudo code):

mean_flux_model = MeanFlux(fid_mean_flux=fid_mean_flux, params=params)
mf = mean_flux_model.get_mean_flux(z=z)

but we could do instead something like:

mean_flux_model = MeanFlux(fid_mean_flux=fid_mean_flux)
mf = mean_flux_model.get_mean_flux(z=z, params=params)

The problem with the latter is that LyaTheory is also asking all these modules for their list of free parameters, and this would not be trivial in the latter setup. Unless we did something like this:

mean_flux_model = MeanFlux(fid_mean_flux=fid_mean_flux, free_params=free_params)
mf = mean_flux_model.get_mean_flux(z=z, param_values=param_values)

where free_params and param_values could actually be both a list of LikelihoodParameter objects, but the former would only be used to set the free parameters in the model, not their values, and therefore a call like this one: mf = mean_flux_model.get_mean_flux(z=z, param_values=None) or mf = mean_flux_model.get_mean_flux(z=z, param_values=[]) would return the fiducial model (and not the model corresponding to the free_params parameters passed in the constructor)

[jchavesmontero]

Sounds good, I could implement that

[jchavesmontero]

I implemented what you suggested in #70 for the mean flux, temperature, and pressure models. It still needs to be done for metals and HCD contamination, but I prefer to wait until you are back

[jchavesmontero]

By the way, the XXX_emcee models are gone as a result