bug when simulation of gammapy datasets with an agnpy spectrum model about agnpy HOT 9 CLOSED

Thanks a lot for reporting this @guillaumegrolleron.

I am surprised that the epsilon array has dimension (30, 10)...
It is converted from the nu that should be a one-dimensional array.
https://github.com/cosimoNigro/agnpy/blob/master/agnpy/synchrotron/synchrotron.py#L174

Do you encounter this problem only for the synchrotron SED?
This is the same problem you encountered in the fitting notebook, right?

from agnpy.

Thanks for your answer @cosimoNigro
In fact this bug is not specific for synchrotron because I have the same with inverse Compton.
This problem doesn't appear when I fit but when I try to simulate fake datasets, indeed when I use the fake() method, npred_signal() is called, which calls compute_flux_spectral() and then integral() and integrate_spectrum().
In this function the evaluator of the SpectralModel is computed with an energy defined as energy = np.geomspace(energy_min, energy_max, num=int(num), axis=-1) https://github.com/gammapy/gammapy/blob/ddbc5dc8a4291f125996a9ce0e9b18ecdcc4c6e5/gammapy/modeling/models/spectral.py#L97
energy_min and energy_max are the bounds of my energy_axis of size 30, and in this context num is equal to 10.
To sum up the problem is caused by the incompatibility between integrate_spectrum() (https://github.com/gammapy/gammapy/blob/ddbc5dc8a4291f125996a9ce0e9b18ecdcc4c6e5/gammapy/modeling/models/spectral.py#L78) which will call the evaluator of Gammapy SpectralModel on a multidimensional energy array and the fact that the evaluate_sed_flux() method (

from agnpy.

Hello @guillaumegrolleron,

first of all I am very sorry for such a slow reaction.
Finally I incorporated the sherpa and gammapy wrappers in the source code.
Could you retry now and tell me if you have the same problem?

If that's the case, in the source code of the gammapy wrapper I left an inline comment
https://github.com/cosimoNigro/agnpy/blob/master/agnpy/fit/gammapy_wrapper.py#L178
pointing to this issue, with a tentative solution.

Let me know, and sorry again that it took a while.

from agnpy.

Just to complement the discussion:
Following this tutorial
https://agnpy.readthedocs.io/en/latest/tutorials/ssc_gammapy_mcmc_fit.html
I was trying to use instead of just FluxPointsDataset, a MapDataset (Fermi data) and a SpectrumDatasetOnOff (VERITAS data)
and I was getting the same error, and found this solution. However, I would get the same error in the SSC calculation, not only synchrotron. I solved by changing

    emissivity = integrator(integrand, gamma, axis=0)

to

    emissivity = integrator(integrand, gamma, axis=0).reshape(epsilon.shape)

as suggested, but also in

    integral_epsilon = integrator(integral_gamma, epsilon, axis=0)

to

    integral_epsilon = integrator(integral_gamma, epsilon, axis=0).reshape(epsilon_s.shape)

But since I am not using SSA and/or EC, I am afraid that the same has to be done there, but I haven't checked.

from agnpy.

Hello @pivosb, thanks for reporting this.

Are you using the latest version of agnpy with the wrappers incorporated?
If so, did you try the solution I suggested here in the code
https://github.com/cosimoNigro/agnpy/blob/master/agnpy/fit/gammapy_wrapper.py#L178

This problem never occurred to me while fitting flux points, it seems it is occurring either when using a MapDataset or doing a simulation, so I guess in those cases some array with energy and some other axis is passed to the functions evaluating the SED.

Would you mind checking if my solution at the wrapper level works?
Thank you.

from agnpy.

Hi @cosimoNigro.

Yes, I am. And yes, I first tried the solution suggested in
https://github.com/cosimoNigro/agnpy/blob/master/agnpy/fit/gammapy_wrapper.py#L178

And you are correct, the problem is not with the flux points. It is specifically with any other dataset, and the error occurs before L178 in the wrapper. It happens when the code calls

    sed_synch = Synchrotron.evaluate_sed_flux(
        nu, z, d_L, delta_D, B, R_b, self._n_e, *args, ssa=self.ssa
    )
    sed_ssc = SynchrotronSelfCompton.evaluate_sed_flux(
        nu, z, d_L, delta_D, B, R_b, self._n_e, *args, ssa=self.ssa
    )

and then the evaluate is called in agnpy/agnpy/synchrotron/synchrotron.py and agnpy/agnpy/compton/synchrotron_self_compton.py.

So yes, it is related to using a MapDataset and/or a SpectrumDatasetOnOff (either a real one or a simulation)

from agnpy.

Ok, thanks, so it seems that my suggested fix will not help as the error occurs at the evaluate_sed_flux function calls.

@pivosb would you mind fixing it in a PR?

from agnpy.

@cosimoNigro

Yes, I can try! I will open a PR soon

from agnpy.

Should be solved by @pivosb PR #131.

from agnpy.