Incorporating the Gammapy wrapper in agnpy

[cosimoNigro]

I think that some people are using our tutorial notebook for fitting with the Gammapy wrapper.
https://agnpy.readthedocs.io/en/latest/tutorials/ssc_gammapy_fit.html

I don't think it is particularly pretty, especially because of the large cell of code defining the wrapper.
There was some discussion in one of the developer calls of Gammapy if we were to include the wrapper in agnpy, or in Gammapy, as currently done for the NaimaSpectralModel.

I think it is ok to have the wrapper within agnpy, we will be introducing Gammapy as an additional dependency but I think it's also fine. Now we will have a centralised wrapper, and hopefully adding a line of code we can solve the issue encountered by @guillaumegrolleron and @jepcor97 when fitting and simulating with the latest Gammapy versions (see the comments in the wrappers).

I thought about having a single wrapper per each radiative processes but this implies:

1. re-writing the same code several times (one per each radiative processes);
2. users having to link all the different parameters between the two (or more) radiative processes they are considering and then define yet another Gammapy.SpectralModel summing them.

So I implemented a wrapper per each scenario / source type, I have written:

1. one with synchrotron + SSC for BL Lacs;
2. one with synchrotron + SSC + external Compton (on BLR and DT) for FSRQs (the EC on BLR computation can be switched off as it is time consuming).

Parameters are read from an instance of the emission region and target.

Many things are missing, still I wanted to open the review to start to discuss the design and the organisation of the wrapper.
@jsitarek what do you think? Let me know what you think about the strategy, the naming scheme, the position of the new submodule.

I will also be happy to have the opinions of @QRemy and @AtreyeeS, from the Gammapy side - if they have time.

Follows a script to cross-check the normal computation against the new wrappers implemented (we can incorporate this in the tests).

import numpy as np
import astropy.units as u
from astropy.coordinates import Distance
from agnpy.emission_regions import Blob
from agnpy.targets import SphericalShellBLR, RingDustTorus
from agnpy.synchrotron import Synchrotron
from agnpy.compton import SynchrotronSelfCompton, ExternalCompton
from agnpy.wrappers import SpectralModelSSC, SpectralModelEC
from agnpy.utils.plot import sed_x_label, sed_y_label
import matplotlib.pyplot as plt

# emission region and targets
# - SSC blob
spectrum_norm = 1e48 * u.Unit("erg")
spectrum_dict = {
    "type": "PowerLaw",
    "parameters": {
        "p": 2.8,
        "gamma_min": 1e2,
        "gamma_max": 1e7
    }
}
R_b = 1e16 * u.cm
B = 1 * u.G
z = Distance(1e27, unit=u.cm).z
delta_D = 10
Gamma = 10
ssc_blob = Blob(R_b, z, delta_D, Gamma, B, spectrum_norm, spectrum_dict)

# - EC blob
spectrum_norm = 6e42 * u.erg
parameters = {
    "p1": 2.0,
    "p2": 3.5,
    "gamma_b": 1e4,
    "gamma_min": 20,
    "gamma_max": 5e7,
}
spectrum_dict = {"type": "BrokenPowerLaw", "parameters": parameters}
R_b = 1e16 * u.cm
B = 0.56 * u.G
z = 1
delta_D = 40
Gamma = 40
ec_blob = Blob(R_b, z, delta_D, Gamma, B, spectrum_norm, spectrum_dict)

# - BLR definition
xi_line = 0.024
L_disk = 1e46 * u.Unit("erg s-1")
R_line = 1e17 * u.cm
blr = SphericalShellBLR(L_disk, xi_line, "Lyalpha", R_line)
# - dust torus definition
T_dt = 1e3 * u.K
csi_dt = 0.1
dt = RingDustTorus(L_disk, csi_dt, T_dt)

# distance
r = 1e16 * u.cm

# agnpy radiative processes
synch = Synchrotron(ssc_blob)
ssc = SynchrotronSelfCompton(ssc_blob)

synch_ec = Synchrotron(ec_blob)
ssc_ec = SynchrotronSelfCompton(ec_blob)
ec_blr = ExternalCompton(ec_blob, blr, r)
ec_dt = ExternalCompton(ec_blob, dt, r)

# radiative processes wrapped by Gammapy
ssc_model = SpectralModelSSC(ssc_blob)
ec_model = SpectralModelEC(ec_blob, blr, dt, r, ec_blr=False)

# SED computation
nu = np.logspace(9, 29, 100) * u.Hz
E = nu.to("eV", equivalencies=u.spectral())

# SSC model
sed_ssc_agnpy = synch.sed_flux(nu) + ssc.sed_flux(nu)
sed_ssc_gammapy = (E**2 * ssc_model(E)).to("erg cm-2 s-1")

# EC model
sed_ec_agnpy = (
    synch_ec.sed_flux(nu) + ssc_ec.sed_flux(nu) + ec_dt.sed_flux(nu)
)
sed_ec_gammapy = (E**2 * ec_model(E)).to("erg cm-2 s-1")

# SSC comparison
plt.loglog(nu, sed_ssc_agnpy, ls="-", label="SSC, agnpy")
plt.loglog(nu, sed_ssc_gammapy, ls="--", label="SSC, Gammapy wrapper")
plt.xlabel(sed_x_label)
plt.ylabel(sed_y_label)
plt.ylim([1e-13, 1e-9])
plt.legend()
plt.show()

# EC comparison
plt.loglog(nu, sed_ec_agnpy, ls="-", label="EC, agnpy")
plt.loglog(nu, sed_ec_gammapy, ls="--", label="EC, Gammapy wrapper")
plt.xlabel(sed_x_label)
plt.ylabel(sed_y_label)
plt.ylim([1e-18, 1e-14])
plt.legend()
plt.show()

[QRemy]

Maybe the wrapped models can follow the naming convention of gammapy models as :
SynchrotronSelfComptonSpectralModel, ExternalComptonSpectralModel

[cosimoNigro]

Done, renamed following the Gammapy naming scheme.

[QRemy]

gammapy models parameters are usually defined at class level and come with default values, but I guess you prefer to avoid this because it's hard to define good defaults here.

[cosimoNigro]

Yes, it is a bit difficult to establish defaults, in general. Additionally each emission region comes with its own electron distribution of which I do not know the parameters. I added functions to fetch these parameters from the physical objects and set some reasonable ranges of values.

[QRemy]

maybe adding something like :
self.emission_region_parameters
self.targets_parameters
could help to simplify the sorting of the parameters later in evaluate
(and could also be convenient for users that want to look only at one group of parameters)

[cosimoNigro]

I grouped the parameters following your suggestion, and so now the user can access spectral_parameters and emission_region_parameters separately.
I don't know though how to re-use this grouping in the evaluate function. I see that the NaimaSpectralModel in Gammapy does something similar, i.e. uses in the evaluate some attribute of the class, but then this _update_naima_parameters is needed https://github.com/gammapy/gammapy/blob/master/gammapy/modeling/models/spectral.py#L2107.
Let me know if you have some suggestion here.

[AtreyeeS]

This should be quite useful @cosimoNigro . One general comment here... Currently, gammapy FluxPointEstimator needs either an amplitude or norm parameter to be defined on the model. (I imagine that it would be interesting to fit reduced datasets directly with physical models and extract flux points). This is not the case here.

It is not a blocker, as the user can still multiply SpectralModelSSC by a normed model to the specific case of extracting flux points, but it is probably important to be correctly conveyed

[AtreyeeS]

Sorry my previous comment was not totally correct. As per the latest developments, users need to specify is_norm for scalable parameters (see discussions in gammapy/gammapy#3690 (comment))

So one option (to simplify things for the users) here maybe be to add a norm parameter which must be kept frozen, but can be used internally during flux point estimation...

[codecov]

Codecov Report

Merging #119 (d9e3a79) into master (86a28d4) will increase coverage by 0.43%.
The diff coverage is 98.70%.

@@            Coverage Diff             @@
##           master     #119      +/-   ##
==========================================
+ Coverage   96.99%   97.43%   +0.43%     
==========================================
  Files          30       38       +8     
  Lines        2233     2997     +764     
==========================================
+ Hits         2166     2920     +754     
- Misses         67       77      +10     

Flag	Coverage Δ
unittests	97.43% <98.70%> (+0.43%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
setup.py	0.00% <ø> (ø)
agnpy/fit/models.py	87.50% <87.50%> (ø)
agnpy/fit/gammapy_wrapper.py	95.39% <95.39%> (ø)
agnpy/tests/test_wrappers.py	99.20% <99.20%> (ø)
agnpy/__init__.py	100.00% <100.00%> (ø)
agnpy/emission_regions/blob.py	94.28% <100.00%> (+0.05%)	⬆️
agnpy/fit/__init__.py	100.00% <100.00%> (ø)
agnpy/fit/core.py	100.00% <100.00%> (ø)
agnpy/fit/data.py	100.00% <100.00%> (ø)
agnpy/fit/sherpa_wrapper.py	100.00% <100.00%> (ø)
... and 4 more

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

[cosimoNigro]

Thanks for your comments @AtreyeeS, now k_e, the normalisation of the electron distribution has is_norm=True when converted to a gammapy.modeling.Parameter. The flux or SED should scale linearly with the normalisation of the radiating electrons, as it is just a multiplicative factor in all the SED calculations.

[cosimoNigro]

Will add some tests later.

[AtreyeeS]

The flux or SED should scale linearly with the normalisation of the radiating electrons, as it is just a multiplicative factor in all the SED calculations.

I think that while the synchrotron and EC flux goes as k_e, the SSC goes as k_e^2 (because its the same population of electrons, so one factor of k_e from the electrons, and one from the sync photons).
You can see the same in the fig below:

[cosimoNigro]

Yes, you are totally right! Thank you!
I will then introduce a general amplitude multiplying the SED and that would be the norm of the model.
Would that be ok?

[AtreyeeS]

@cosimoNigro : In gammapy, we tend to use amplitude as a dimensional parameter, and norm as a scaling, so maybe norm is better? And document it well that the norm must be frozen during the fit. I don't know if it is possible to hide it from the users in any way, but if you are separating parameters as blob_parameters and jet_parameters as @QRemy suggested, then at least the norm should not show up in either of those

[cosimoNigro]

Ok, so I added the norm, and I put it a the bottom of the parameter list. It does appear if you call all the parameters together (don't know how to hide it).

from agnpy.emission_regions import Blob
from agnpy.wrappers import SynchrotronSelfComptonSpectralModel

blob = Blob()
ssc_model = SynchrotronSelfComptonSpectralModel(blob)

print(ssc_model.parameters.to_table())

  type      name     value    unit   error      min       max    frozen is_norm link
-------- --------- ---------- ---- --------- --------- --------- ------ ------- ----
spectral       k_e 9.2748e+06 cm-3 0.000e+00 1.000e+00 1.000e+09  False   False     
spectral         p 2.8000e+00      0.000e+00 1.000e+00 5.000e+00  False   False     
spectral gamma_min 1.0000e+02      0.000e+00 1.000e+00 1.000e+03   True   False     
spectral gamma_max 1.0000e+07      0.000e+00 1.000e+04 1.000e+08   True   False     
spectral         z 6.9000e-02      0.000e+00 1.000e-03 1.000e+01   True   False     
spectral       d_L 9.9203e+26   cm 0.000e+00 1.368e+25 3.271e+29   True   False     
spectral   delta_D 1.0000e+01      0.000e+00 1.000e+00 1.000e+02  False   False     
spectral         B 1.0000e+00    G 0.000e+00 1.000e-04 1.000e+03  False   False     
spectral       R_b 1.0000e+16   cm 0.000e+00 1.000e+12 1.000e+18  False   False     
spectral      norm 1.0000e+00      0.000e+00 1.000e-01 1.000e+01   True    True   

As @AtreyeeS said, it does not appear in the "grouped" parameters

print(ssc_model.spectral_parameters.to_table())

  type      name     value    unit   error      min       max    frozen is_norm link
-------- --------- ---------- ---- --------- --------- --------- ------ ------- ----
spectral       k_e 9.2748e+06 cm-3 0.000e+00 1.000e+00 1.000e+09  False   False     
spectral         p 2.8000e+00      0.000e+00 1.000e+00 5.000e+00  False   False     
spectral gamma_min 1.0000e+02      0.000e+00 1.000e+00 1.000e+03   True   False     
spectral gamma_max 1.0000e+07      0.000e+00 1.000e+04 1.000e+08   True   False

print(ssc_model.spectral_parameters.to_table())

  type     name    value    unit   error      min       max    frozen is_norm link
-------- ------- ---------- ---- --------- --------- --------- ------ ------- ----
spectral       z 6.9000e-02      0.000e+00 1.000e-03 1.000e+01   True   False     
spectral     d_L 9.9203e+26   cm 0.000e+00 1.368e+25 3.271e+29   True   False     
spectral delta_D 1.0000e+01      0.000e+00 1.000e+00 1.000e+02  False   False     
spectral       B 1.0000e+00    G 0.000e+00 1.000e-04 1.000e+03  False   False     
spectral     R_b 1.0000e+16   cm 0.000e+00 1.000e+12 1.000e+18  False   False  

I integrated the computation of the SED with the Gammapy wrappers in the tests, I find quite some discrepancy at the highest energies. I expect exactly identical results as I am using the same functions. Will try to dig deeper later.

[cosimoNigro]

Tests fail because the dev version of Gammapy is needed. Will fix it.

[cosimoNigro]

Ok, now it's complaining about the astropy version... 😮‍💨

[cosimoNigro]

Ok, very well, the tests were failing because Gammapy - that I am now installing - uses astropy 5.0 that has dropped support for python 3.7. Just removed python 3.7 from the tests.

[cosimoNigro]

Will take time to make a final couple of improvements in the EC wrapper and then merge.
In particular, the emission of the thermal components (disk and DT) is missing, also I would like to improve how the target for scattering is defined.

[cosimoNigro]

Ok, I am done with this.
I have introduced wrappers both for Gammapy and sherpa models.
I have also added some helper functions to load the SED data in a sherpa.data.Data1D object or in a list of gammapy.datasets.FluxPointsDatasets.
Proper tests have been added for the new classes / functions added.
I can now obtain the same results from the two wrappers for the two sources considered, Mrk421 and PKS1510-089.

I will make a new release after I merge.