Update mcmc to be compatible with v1.2

[facero]

This PR updates the recipe to match with v1.2 and in particular to make use of the new Priors .
With the priors, there is no longer need for the sampling.py external file.

A comparison of two MCMC sampling algorithms is shown.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[bkhelifi]

Line #2.        "/Users/facero/Documents/Work/Program/gammapy-0.19/gammapy-datasets/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits"

The correct file is $GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits

Reply via ReviewNB

[bkhelifi]

The outputs of the notebook should be cleared for the final PR version.

Reply via ReviewNB

[Astro-Kirsty]

An easy way to do this is with the following commands:

from utils import notebook_black, notebook_strip
notebook_black(path_to_notebook)
notebook_strip(path_to_notebook)

[bkhelifi]

A trick to pu maybe in the dev documentation;)
Thanks Kirsty

[facero]

@Astro-Kirsty : I've tried the notebook_black thing and could not find it.
It's not in gammapy.utils
Which utils are you referring to ?

[Astro-Kirsty]

Apologies, it is in this repo under scripts

[bkhelifi]

Line #3.        # Update model parameters factors inplace

Can you add a description of the parameters?

Reply via ReviewNB

[Astro-Kirsty]

Agreed! And if it really is factor, then I would suggest to rename the pars to factors for consistency.

Adding a docstring to this would be really helpful.

Then you could mention all the comments inside the docstring instead

[bkhelifi]

Line #4.        for factor, par in zip(pars, dataset.models.parameters.free_parameters):`

pars is a factor in fact?

Reply via ReviewNB

[facero]

par was badly chosen, in fact par is a gammapy parameter object and pars is a list of factors.
But it more conventional to have lnprob() function taking list of pars as input
and the factor thing is an internal gammapy thing.
Not sure what to do
for now I'm changing to . What do you think ?

def lnprob(pars, dataset):
    """
    Estimate the likelihood of a model including prior on parameters.
    Input : 
    pars : a list of parameters
    dataset: a gammapy dataset 
    """
    # Update model parameters factors inplace
    # The MCMC sampler will evaluate the likelihood of the model given
    # a set of parameters. We need to update the model parameters value before 
    # evaluating the new likelihood value.
    for factor, parameter in zip(pars, dataset.models.parameters.free_parameters):
        parameter.factor = factor # What the sampler will need is the factor (e.g. 1) instead of the value (1e-12)
        
    # dataset.likelihood returns Cash statistics values
    # emcee will maximisise the LogLikelihood so we need -dataset.likelihood
    total_lnprob = -dataset.stat_sum() #stat_sum now includes stat + stat_priors

    return total_lnprob

[bkhelifi]

thanks

[bkhelifi]

Line #20.        nwalkers, ndim, lnprob, args=[dataset], #pool=pool

lnprob trakes two arguments. What are they?

Reply via ReviewNB

[facero]

the list of parameters (the free params in the model) and the gammapy dataset.
I've added a description in the lnprob() function definitoin (see previsou comment).

[bkhelifi]

@facero, thank you a lot for the revival of this recipes. Once finished, it would be great to mention it in the central documentation!
Otherwise, do you understand why the tests are failing? @cgalelli, any idea?

[facero]

I've updated the env yaml file as the previous failure seemed to be related to a python 3.12 not compatible .
but now the error message is even stranger :

Warning: Future Miniforge releases will NOT build Mambaforge installers. We advise you switch to Miniforge at your earliest convenience. More details at [https://conda-forge.org/news/2024/07/29/sunsetting-mambaforge/.](https://conda-forge.org/news/2024/07/29/sunsetting-mambaforge/)

We'll need to undertand how to fix this Mambaforge issue
Then :

  Warning: ERROR: executing pre_install.sh failed
  
  ERROR: executing pre_install.sh failed
Error: The process '/usr/bin/bash' failed with exit code 1

[Astro-Kirsty]

I would rather make this a markdown cell, to make it more clear.

Then for the code you could utilise the python code functionality of markdown?

Reply via ReviewNB

[Astro-Kirsty]

Change the first sentence because we do define priors now!

Reply via ReviewNB

[Astro-Kirsty]

Thanks for your effort here @facero!
Left some inline comments. Was there a reason you removed the last part PeVatrons in CTA? ?

edit: can you correct the mcmc plot for the amplitude. The amplitude above the plot is showing 0.00+-0.00 due to the scientific notation!

[cgalelli]

Otherwise, do you understand why the tests are failing? @cgalelli, any idea?

The full test fail is related to this i suspect

[facero]

Thanks for your effort here @facero! Left some inline comments. Was there a reason you removed the last part PeVatrons in CTA? ?

edit: can you correct the mcmc plot for the amplitude. The amplitude above the plot is showing 0.00+-0.00 due to the scientific notation!

This is related to the automatic corner plots generation not in my "simple" control.
Maybe there is a way to change that, for now I'm leaving it like that.

About the PeVatrons, No good reasons for removing it, I can add it back indeed.

[facero]

One question that we need to think about is whether we want the sampling to be done in terms of factor or values.
I would have thought that working on close to unity values would have been better but I compared convergence and speed of both MCMC and it was the similar (i.e sampling 1e-12 or sampling 1).

By default I was working on factor but I realized that the model I define by default has scale = 1 for the amplitude :
Parameter(name='amplitude', value=5e-12, factor=5e-12, scale=1

So the factor and the value were the same. In order to scale the amplitude I had to use :
dataset.models.parameters.autoscale()

But this brings some complications for example in the trace or corner plots where the samples are shown without rescaling (amplitude close to 1 and lambda_ close to 1 even if true lamba is 0.01). The corner plot would show lambda = 5 whereas it really means 0.05. So it can be a bit confusing.

Do we leave everything in terms of factor :
p0 = [free_par.factor for free_par in dataset.models.parameters.free_parameters]

But just not rescale or do I change everything from free_par.factor to free_par.values which would simplify and avoid confusion ?
Opinions @Astro-Kirsty , @bkhelifi , @registerrier ?

[facero]

After an in person discussion we decided to keep the free_par.values approach over the factor approach.
Using values has the advantage of printing/plotting the real values (1e-12) on the plots over the factor value (e.g. 1).

[facero]

I clean the notebook and applied formatting on it.
I think it is now ready to merge. Let me know.

[bkhelifi]

Thank you a lot @facero, We have now a very nice updated recipes;)

[registerrier]

Thanks @facero . Looks good to me too!