Add function to determine pivot energy for any spectral model

[Astro-Kirsty]

Description
This PR is addressing the feature request from issue #4016.
I have added a pivot_energy function to the SpectralModel class that finds the energy at which the error butterfly is minimum. This function is then 'overridden' in the PowerLawSpectralModel where there is an analytical equation available.

Please note: this must be performed in log space (as I understand) because the butterfly minimum occurs in the log space plot.

Dear reviewer
I am certainly open to suggestions/feedback on this issue.

[MRegeard]

Thanks @Astro-Kirsty,
Really helpful feature ! See inline comments.

[bkhelifi]

Thanks a lot @Astro-Kirsty . With this PR, I think that it would be great to have some tutorials with this new feature.... e.g. the spectral analysis or 3D detailed analysis

[AtreyeeS]

Thanks @Astro-Kirsty !
Is is sufficient to see the thinnest band among the range of energies supplied by the user? The decorrelation energy may not even lie in the given range.

If the energy array is too coarsely binned, the estimate will be incorrect. Finding the minima of the error might be better, no?

[codecov]

Codecov Report

Merging #4635 (b9bf6a3) into main (74d0de2) will decrease coverage by 0.04%.
The diff coverage is 100.00%.

❗ Current head b9bf6a3 differs from pull request most recent head 98e9bef. Consider uploading reports for the commit 98e9bef to get more accurate results

@@            Coverage Diff             @@
##             main    #4635      +/-   ##
==========================================
- Coverage   95.09%   95.06%   -0.04%     
==========================================
  Files         222      221       -1     
  Lines       31750    31561     -189     
==========================================
- Hits        30193    30002     -191     
- Misses       1557     1559       +2     

Files Changed	Coverage Δ
gammapy/modeling/models/spectral.py	97.29% <100.00%> (+0.01%)	⬆️
gammapy/modeling/models/tests/test_spectral.py	100.00% <100.00%> (ø)

... and 21 files with indirect coverage changes

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[AtreyeeS]

Thanks @Astro-Kirsty ! This looks better, but it is still an issue that this does not guarantee you are returning the decorrelation energy. This is just the minimum error point in the given range (without letting users know if the actual pivot lies outside). I tried with fsolve, and the results look stable. Maybe this is not the best optimiser, but something along the lines of

from scipy.optimize import fsolve
from gammapy.modeling.models import *

lp = LogParabolaSpectralModel(
        amplitude="5.82442e-11 cm-2 s-1 GeV-1",
        reference="17.337 GeV",
        alpha="1.9134",
        beta="0.2217",
    )
lp.amplitude.error = "2.8804e-12 cm-2 s-1 GeV-1"
lp.alpha.error = "0.1126"
lp.beta.error = "0.0670"

def f1(x, unit=u.TeV):
    x=np.exp(x)
    dnde, dnde_error = lp.evaluate_error(x*unit)
    ene = dnde_error / dnde
    return ene

print("pivot = ", np.exp(fsolve(f1, x0=0))*u.TeV)

maybe an option? What do you think?

[Astro-Kirsty]

I agree @AtreyeeS. Using fsolve would be a nice way so that the use does not have to define any energy and therefore always gives the decorrelation energy for that SpectralModel

[Astro-Kirsty]

Implemented the code using fsolve, I tried the function on multiple spectra to find a reasonable initial guess.

[registerrier]

Thanks @Astro-Kirsty .

fsolve is a root finder, not a minimizer. Is it really adapted here? Shouldn't one use something like minimize_scalar instead?

Those methods require an interval to search for the minimum. We have to see if it is possible to keep pivot_energy a simple property in this context.

One should probably check that a proper minimum is found and yield an error if not.

[Astro-Kirsty]

An option would be to implement:
minimizer = scipy.optimize.minimize_scalar(min_energy, bounds=bounds)

where the bounds are predefined by default as the reference value from the SpectralModel plus or minus some small value. But making this an optional argument, so that the user can specific their own bounds?

And as you mention, to add that if the minimizer is not successful, then it should alert the user

[Astro-Kirsty]

Updated to make sure there is no issue with units.

Added an if statement for the std. The scipy minimize functions will always minimize a function, even if there is no minimum (i.e. a constant array of values). Hence providing the user with an incorrect, nonsensical pivot energy. Therefore, adding this if statement takes care of this.

I also adapted the test functions to ensure this was included.

[registerrier]

Thanks @Astro-Kirsty . There is a test failing, could you increase the tolerance? I suppose the method won't be precise below the percent.

[Astro-Kirsty]

Thanks @registerrier. It was not failing on my laptop which is odd. I have increased the tolerance

[registerrier]

Thanks @Astro-Kirsty . This looks good.

The only missing thing to me is a test with some correlation. e.g. use an exp-cutoff with the same parameters and correlation matrix as the PWL test and check you get the same result.

[QRemy]

Thanks @Astro-Kirsty this is a nice feature, I have just one question :

[QRemy]

does the choice of the number of steps in linespace can affect the precision ?
If yes it should be configurable.

[Astro-Kirsty]

It does not

[AtreyeeS]

Thanks! no more comments from me

[bkhelifi]

Thanks @Astro-Kirsty, all good for me.
I am still wondering if one could not have an optional argument about the energy range of computation. But this can be added later if needed.

[MRegeard]

Thanks @Astro-Kirsty,
Seems good to me !