Clean up spectrum code about gammapy HOT 9 CLOSED

@joleroi @adonath - I think we should make spectral models a priority in Gammapy ASAP.

There's #591 which roughly re-implements gammapy.spectrum.Model.integral and gammapy.spectrum.PowerLaw.integral.

That is very old code that was written before Astropy existed. It doesn't have tests and I think those models aren't used or usable with the other spectrum code. But I do think spectral model objects are better than standalone functions.

As we've discussed for a long time, moving forward on this is difficult because ideally we'd like to share model implementations with fermipy and naima (they both have their own at the moment) and have the models work with Sherpa and maybe be Astropy models as well. That's too much, we're stuck and we have to give up on some of those goals and implement something now and clean up the rest of the existing code.

@joleroi - Can you extend
http://nbviewer.jupyter.org/github/gammapy/gammapy-extra/blob/master/notebooks/hess_spectrum_analysis.ipynb or make a new Notebook that explains how spectral models in Gammapy with SpectrumFit work at the moment?

• Is there a way to use a built-in Sherpa spectral model?
• Is there a way for a user to define their own spectral model and fit it without having to edit Gammapy files?

from gammapy.

Is there a way to use a built-in Sherpa spectral model?

Yes, at the moment you can only fit built-in sherpa models with Gammapy
https://github.com/gammapy/gammapy/blob/master/gammapy/spectrum/spectrum_fit.py#L134

Is there a way for a user to define their own spectral model and fit it without having to edit Gammapy files?

See above

from gammapy.

Update. The Fit returns a SpectrumFitResult object, with the best fit parameters and errors. This object can be evaluated by calling into static methods of astropy models
https://github.com/gammapy/gammapy/blob/master/gammapy/spectrum/results.py#L397
I chose this as a temporary solution, because astropy models cannot be instantiated with Quantities. In the long run it would of course be nice to have some model classes that can be instantiated from the SpectrumFitResult. But as you said there was so much discussion that I didn't want to waste time with implementing something fancy.

from gammapy.

There has been a lot of progress regarding the clean up recently.

Notably we have the spectrum.models module which works quite well I think.

Is there a way for a user to define their own spectral model and fit it without having to edit Gammapy files?

Yes and no. You can't build a gammapy.spectrum.models.SpectralModel youself (but IMO that's out of scope). But you can simply pass a Sherpa model to the SpectrumFit class.

IMO What's left to do here is

• > @joleroi - Can you extend
  http://nbviewer.jupyter.org/github/gammapy/gammapy-extra/blob/master/notebooks/hess_spectrum_analysis.ipynb or make a new Notebook that explains how spectral models in Gammapy with SpectrumFit work at the moment?

from gammapy.

You can't build a gammapy.spectrum.models.SpectralModel youself (but IMO that's out of scope). But you can simply pass a Sherpa model to the SpectrumFit class.

Why not?
That's a very common use case.
Sherpa has explicit methods and docs for this, see e.g.
http://cxc.harvard.edu/sherpa/threads/user_model/

Why can't I subclass MySpectralModel(SpectralModel) in my script and then fit it at the moment?

Is it because there's a hard-coded list of supported models in a method?

If so, this hard-coded list should be replaced with a global dict spectral_model_registry where users can register their model.
(or just take the model as a parameter, why does there have to be a hard-coded list at all?)

from gammapy.

In Python using a global dict or other mutable data structure as a registry is the simple, low-tech solution (I'm using it in gammapy.catalogs already).

There's a slightly nicer, but also more complicated method using metaclasses:
https://github.com/faif/python-patterns/blob/master/registry.py

My advice would be to do the global dict for now and put the metaclass registry as an open issue to discuss and possibly implement later.

from gammapy.

Why can't I subclass MySpectralModel(SpectralModel) in my script and then fit it at the moment?

You would need to implement a to_sherpa method

So there's almost no benefit with regard to using sherpa models directly

from gammapy.

You would need to implement a to_sherpa method

So it is possible. Great!

So there's almost no benefit with regard to using sherpa models directly

Well, you get functionality to compute integral fluxes, energy fluxes, butterfly and have error propagation via uncertainties from the SpectralModel base class. That's something many people will want and need for their publication.

Maybe someday there could be a notebook showing how to define a user-defined spectral model and fit it using either or both ways?
It could just implement MyPowerLaw or do an ECPL, no need to implement think of and implement something new, just to show that it works and how to do it.

from gammapy.

I'm closing this issue now.

It just says "clean up spectrum code" in the title and that's not very helpful.
A lot of changes / cleanup have happened already, and more will.

@joleroi - If you feel like there is something useful here (like the suggestion to have an example of how to write a user-defined spectral model class and use it with FitSpectrum in the docs), please open a new more specific issue.

from gammapy.