New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add spectrum fit serialisation code #421

Merged
merged 16 commits into from Feb 2, 2016

Conversation

Projects
None yet
3 participants
@joleroi
Contributor

joleroi commented Jan 20, 2016

Add functionality to serialise fit results from gammapy.spectrum

  • Add gammapy.spectrum.results.SpectrumFitResult with
    • read_sherpa method
    • read (FitSpectrum) JSON method
    • read method
    • write method
    • plot function
    • plot residuals
    • plot fluxpoints
  • Add gammapy.spectrum.results.SpectrumFitResultDict for comparing analyses with
    • to_table method
    • plot_method

TODO

  • Do not hardcode flux point units!
Fit results as return by :func:`~sherpa.astro.datastack.covar`
"""
model = SpectralModelPowerLaw()
model.meta = Bunch(covar.__dict__)

This comment has been minimized.

@joleroi

joleroi Jan 20, 2016

Contributor

@cdeil: in the from_json method, the json file is more or less written as it is into one attribute meta of the SpectralModel. I can do the same with the sherpa fit output, but this will lead to problems sooner or later I guess.
Do you have a proposal how to make this coherent (ASAP, not a final perfect solution)

This comment has been minimized.

@cdeil

cdeil Jan 21, 2016

Member

I need to look at the code and think about it a bit.
If you know a way how to add it that works for you and want something ASAP, go ahead.
Otherwise I can have a look after lunch.

This comment has been minimized.

@joleroi

joleroi Jan 21, 2016

Contributor

I also thought about it, my suggestion would be

  • Leave the meta attribute as it is in order to save all available information
  • Introduce attributes for the most important parameters and gradually add the ones I need. For now this would probably be
model.gamma
model.gamma_err
model.norm
model.norm_err

These would need to be filled in the from_json and from_sherpa methods. I guess that would be a start.

What do you think?

This comment has been minimized.

@cdeil

cdeil Jan 21, 2016

Member

@joleroi - I don't understand, sorry.

How are you interfacing with Sherpa models at the moment?
Are you sub-classing them to extend them, or do you use factory or converter functions to go back and forth with Gammapy classes?

As I said ... if you have a solution to compute what you need for now, I suggest you run with it, and then I'll read over and meditate about it in the coming days ...

@joleroi

This comment has been minimized.

Contributor

joleroi commented Jan 20, 2016

@cdeil I would like to have the functionality to serialize fit results ASAP (to finally be able to show some results). I have a question - see inline comment

@joleroi

This comment has been minimized.

Contributor

joleroi commented Jan 22, 2016

This is a proposal for the SpectrumFitResult class.

it write to a format different thatn the JSON format of FitSpectrum, because this format is tailored very much to FitSpectrum, e.g. it has a 'fitoptions' keys, which stores HAP options. So I think a slimmer new format, that concentrates only on FitResults is nicer to handle. But we will see....

@joleroi joleroi self-assigned this Jan 22, 2016

@joleroi joleroi added the feature label Jan 22, 2016

@joleroi joleroi force-pushed the joleroi:serialisation branch from d23308c to 5959390 Jan 25, 2016

joleroi added some commits Jan 26, 2016

@joleroi

This comment has been minimized.

Contributor

joleroi commented Jan 28, 2016

joleroi added some commits Jan 28, 2016

if pname == 'gamma':
name = 'Gamma'
unit = Unit('')
elif pname == 'ampl':

This comment has been minimized.

@tibaldo

tibaldo Jan 29, 2016

Contributor

Didn't we discuss about assigning parameter names for serialized results based on
gamma-astro-data-formats/source/source-models/spectral_models.yaml
?

This comment has been minimized.

@joleroi

joleroi Jan 29, 2016

Contributor

Yes, I forgot, I will change it

@tibaldo

This comment has been minimized.

Contributor

tibaldo commented Jan 29, 2016

@joleroi Is the results serialization part complete so that we can merge with master and I can add a read
_from_3ML method?

@joleroi

This comment has been minimized.

Contributor

joleroi commented Jan 29, 2016

I will work on this this afternoon, and maybe @cdeil wants to have a look. Whould you be fine with merging it monday morning?

joleroi added some commits Jan 29, 2016

@tibaldo

This comment has been minimized.

Contributor

tibaldo commented Jan 29, 2016

@joleroi absolutely, even better because I can bug you in person on Monday if there is something I don not understand :)

joleroi added some commits Jan 29, 2016

@cdeil cdeil added this to the 0.4 milestone Feb 1, 2016

from astropy.table import Table, Column, QTable
from astropy.units import Unit
from gammapy.extern.bunch import Bunch

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

Use explicit relative imports within Gammapy.

@requires_dependency('scipy')
@requires_dependency('sherpa')
@pytest.mark.skipif(not HAS_WSTAT, reason="Wstat only in sherpa head version")

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

"in head version" -> "in version >= 4.8"

Fitted parameters, for allowed names see :ref:`gadf:source-models`
parameter_errors : dict
Parameter errors, for allowed names see :ref:`gadf:source-models`
energy_range : `~gammapy.utils.energy.EnergyBounds

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

Missing backtick at the end of the line ...

Several models can be stored in the JSON file, the appropriate model
to read can be chosen with the ``model`` parameter.
TODO: Should this be a function?

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

Remove this TODO? Or do you think a function is better?

This comment has been minimized.

@joleroi

joleroi Feb 1, 2016

Contributor

The problem I saw is that this function won't be used by non-HESS users (since they don't have FitSpectrum). So I didn't want to have it so exposed. But this is a temporary solution anyways (in the long run FitSpectrum should directly write the 'correct' file)

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

Ah, OK.
Maybe change the TODO from "Should this be a function" to "Remove method when HESS FitSpectrum exporter is updated" then?
Some iteration on the serialisation format is needed for the coming months ... in the meantime it's fully OK to have such converter and convenience functions / methods in Gammapy.

def from_sherpa(cls, covar, filter, model, flux_graph = None):
"""Create SpectrumFitResult from sherpa objects
TODO: Should this be a function?

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

Classmethod is nice, no? Remove this TODO?

Document parameters? (or leave a TODO to make it easier to remember ...)

This comment has been minimized.

@joleroi

joleroi Feb 1, 2016

Contributor

same problem, no one is going to use this function

"""
el, eh = float(filter.split(':')[0]), float(filter.split(':')[1])
energy_range = EnergyBounds((el, eh), 'keV')
if model.name.split('.')[0] == 'powlaw1d':

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

This looks very hacky:

model.name.split('.')[0] == 'powlaw1d'

There must be a better way to do this check, no?
Can you interactively explore the model object to see if you can get at the powerlaw1d string in a simpler way?
Or replace by isinstance check?

spectral_model=spectral_model, flux_points=flux_points,
fluxes=fluxes, flux_errors=flux_errors)
def to_yaml(self):

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

I still haven't found a good pattern for this. But combining the filling / computation of the dict with the yaml serialisation is not very good, because then you can't simply write to_json.
Maybe split out the dict filling part into a to_dict and call that from to_yaml?

@classmethod
def from_yaml(cls, val):
"""Create `~gammapy.spectrum.results.SpectrumResult` from YAML dict

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

This method doesn't do YAML de-serialisation, it just takes a dict, no?
If so, call the method from_dict and make the docstring correct` (it currently says "from YAML dict").

"""
import matplotlib.pyplot as plt
if 'fmt' not in kwargs:

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

These two lines are equivalent to this, no?

kwargs.setdefault('fmt', 'o')

This comment has been minimized.

@joleroi

joleroi Feb 1, 2016

Contributor

I knew there must be a better way to do this, thanks!

import matplotlib.pyplot as plt
ax = plt.gca() if ax is None else ax
if 'fmt' not in kwargs:

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

Same here ... use dict.setdefault to shorten the 2 lines to 1 line.

@@ -8,6 +8,7 @@
from astropy.extern import six
from astropy.wcs.utils import skycoord_to_pixel
from gammapy.extern.pathlib import Path

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

Use implicit relative imports within Gammapy

covar = datastack.get_covar_results()
efilter = datastack.get_filter()
# First go on calculation flux points following

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

Flux point computation is a big project on it's own, with some existing code in Gammapy.
For now, can you split this out into a separate function or method from _run_sherpa_fit and leave a TODO that this should be improved?

@@ -6,6 +6,8 @@
from astropy.utils.compat import NUMPY_LT_1_9
from numpy.testing import assert_allclose
from astropy.coordinates import SkyCoord, Angle
from gammapy.spectrum.results import SpectrumFitResult

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

Use explicit relative imports within Gammapy

@@ -19,6 +21,12 @@
SpectrumFit,
)
try:

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

This try ... except code that checks for WStat is present twice.
Can you please implement a SHERPA_GT_4_8 in gammapy.utils.testing.
(I haven't looked how from astropy.utils.compat import NUMPY_LT_1_9 is implemented ... can you check and do the same?)

@cdeil cdeil changed the title from Add functionality to serialise fit results from gammapy.spectrum to Add spectrum fit serialisation code Feb 1, 2016

@cdeil

This comment has been minimized.

Member

cdeil commented Feb 1, 2016

I've left minor inline comments ... +1 to merge once those are addressed.
(I also added the 0.4 milestone and shortened the PR title a bit. You could add the changelog entry before merging ... that's what I do when a PR is ready to go, to try to limit working directly in master.)

joleroi added some commits Feb 1, 2016

@joleroi

This comment has been minimized.

Contributor

joleroi commented Feb 1, 2016

@cdeil I cannot reproduce the Sphinx build error locally. What command is run exactly on Travis again?
I tried

python setup.py build_sphinx -w
python setup.py build_sphinx -l
@cdeil

This comment has been minimized.

Member

cdeil commented Feb 1, 2016

@joleroi - The Sphinx build error is present in gammapy master and is unrelated to this PR:
#428 (comment)
I don't know why it showed up now ... this docs file hasn't been touched in a long time ... probably Sphinx or some code it uses was updated and it's stricter now.

Feel free to merge this PR now. Or let it sit until tomorrow and I'll re-start the travis-ci build for this PR once #428 is merged.

@cdeil

This comment has been minimized.

Member

cdeil commented Feb 1, 2016

Now that #428 is merged, the Sphinx build should work. I've restarted travis-ci just now.

@joleroi - Can you remove the Long term -> abandon way spectral model is handled entry from the task list for this PR? Either make a new issue listing the next steps, or just remove it if having such a reminder issue with a task list is not helpful to you.

]
# Todo : Change to 4_8 once there is a sherpa release

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

@joleroi - I think this TODO can be resolved now!?

The Sherpa version we support and test at the moment is 4.8:
https://travis-ci.org/gammapy/gammapy/jobs/106282386#L1026

$ python -c 'import sherpa; print(sherpa.__version__)'
4.8.0

If you use the git master version it will be 4.8.something:

$ python2.7 -c 'import sherpa; print(sherpa.__version__)'
4.8.0+82.gbfcc195.dirty

Testing for '4.7+642' seems unnecessary / brittle to me.

This comment has been minimized.

@joleroi

joleroi Feb 2, 2016

Contributor

I did not know that there was a sherpa 4.8 release. The head version I was using was 4.7+642 and I wanted the tests to run locally.

def overplot(self):
"""Overplot spectra"""
pass

This comment has been minimized.

@cdeil

cdeil Feb 1, 2016

Member

I usually put raise NotImplementedError for functions / methods that aren't implemented yet.
If you put pass the user gets no feedback that something is wrong.

This comment has been minimized.

@joleroi

joleroi Feb 2, 2016

Contributor

poor user 😀

@cdeil

This comment has been minimized.

Member

cdeil commented Feb 1, 2016

@joleroi - All travis-ci builds passed for this PR.

@joleroi

This comment has been minimized.

Contributor

joleroi commented Feb 2, 2016

@tibaldo I am merging this, so feel free to add whatever you like!

joleroi added a commit that referenced this pull request Feb 2, 2016

Merge pull request #421 from joleroi/serialisation
Add spectrum fit serialisation code

@joleroi joleroi merged commit 79502a0 into gammapy:master Feb 2, 2016

2 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details

@joleroi joleroi deleted the joleroi:serialisation branch Feb 2, 2016

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment