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 Sherpa cube analysis prototype #484

Merged
merged 6 commits into from Mar 16, 2016

Conversation

Projects
None yet
2 participants
@adonath
Member

adonath commented Mar 7, 2016

This PR implements basic Sherpa classes for simple cube style analysis. An example notebook how to use it, can be found here. Any suggestion how to proceed with this? @cdeil @joleroi
E.g. add PSF and background handling and create a command-line tool? So, we can say cube style analysis is possible in gammapy?

@cdeil cdeil added the feature label Mar 7, 2016

@cdeil cdeil added this to the 0.4 milestone Mar 7, 2016

@cdeil cdeil self-assigned this Mar 7, 2016

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Mar 7, 2016

Member

I'll do a quick review now (but not try out the code / review implementation details).

What about the example files produced by
https://github.com/gammapy/gammapy-extra/blob/master/experiments/sherpa_cube_analysis_prototype/prepare_data_cubes.py ?

Could you save them as .fits.gz, and if they are no more than a few MB, commit them to the repo?
(just as a convenience / to have the inputs under version control / available for tests and examples?)

Member

cdeil commented Mar 7, 2016

I'll do a quick review now (but not try out the code / review implementation details).

What about the example files produced by
https://github.com/gammapy/gammapy-extra/blob/master/experiments/sherpa_cube_analysis_prototype/prepare_data_cubes.py ?

Could you save them as .fits.gz, and if they are no more than a few MB, commit them to the repo?
(just as a convenience / to have the inputs under version control / available for tests and examples?)

Show outdated Hide outdated gammapy/image/tests/test_utils.py Outdated
@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Mar 7, 2016

Member

As you're already noted at https://github.com/gammapy/gammapy-extra/blob/master/experiments/sherpa_cube_analysis_prototype/sherpa_cube_analysis_prototype.ipynb ... it would be good if you put that code somewhere in Gammapy ... these could be accessible as properties or via methods on SpectralCube, right?

# TODO: attach getting coordinate cubes to SpectralCube class
n_ebins = len(elo)
ra_cube = np.tile(ra, (n_ebins, 1, 1))
dec_cube = np.tile(dec, (n_ebins, 1, 1))
elo_cube = elo.reshape(n_ebins, 1, 1) * np.ones_like(ra)
ehi_cube = ehi.reshape(n_ebins, 1, 1) * np.ones_like(ra)

Even this could be spectral_cube.sherpa_data3d, no?
(with a delayed import of the Sherpa class inside that property)

# Internally sherpa uses flattened arrays, therefore call .ravel()
cube = Data3D('counts', elo_cube.ravel(), ehi_cube.ravel(), ra_cube.ravel(),
              dec_cube.ravel(), counts.data.value.ravel(), counts.data.value.shape)
Member

cdeil commented Mar 7, 2016

As you're already noted at https://github.com/gammapy/gammapy-extra/blob/master/experiments/sherpa_cube_analysis_prototype/sherpa_cube_analysis_prototype.ipynb ... it would be good if you put that code somewhere in Gammapy ... these could be accessible as properties or via methods on SpectralCube, right?

# TODO: attach getting coordinate cubes to SpectralCube class
n_ebins = len(elo)
ra_cube = np.tile(ra, (n_ebins, 1, 1))
dec_cube = np.tile(dec, (n_ebins, 1, 1))
elo_cube = elo.reshape(n_ebins, 1, 1) * np.ones_like(ra)
ehi_cube = ehi.reshape(n_ebins, 1, 1) * np.ones_like(ra)

Even this could be spectral_cube.sherpa_data3d, no?
(with a delayed import of the Sherpa class inside that property)

# Internally sherpa uses flattened arrays, therefore call .ravel()
cube = Data3D('counts', elo_cube.ravel(), ehi_cube.ravel(), ra_cube.ravel(),
              dec_cube.ravel(), counts.data.value.ravel(), counts.data.value.shape)
@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Mar 7, 2016

Member

Can you print the model at the end of
https://github.com/gammapy/gammapy-extra/blob/master/experiments/sherpa_cube_analysis_prototype/sherpa_cube_analysis_prototype.ipynb ?
So that the best-fit parameters and errors are shown nicely?

Concerning this:

# Adding this constant background components the fit works with cash statistics as well
#spatial_model_bkg = Const2D('spatial-model-bkg')

My preference would be that you add this for now then and show the Cash stat as an example.
We probably don't want to spend time looking at numbers from the chi2 fit, so a working cash fit is a better starting point (even if the numbers should be more off for now), in my opinion.

Member

cdeil commented Mar 7, 2016

Can you print the model at the end of
https://github.com/gammapy/gammapy-extra/blob/master/experiments/sherpa_cube_analysis_prototype/sherpa_cube_analysis_prototype.ipynb ?
So that the best-fit parameters and errors are shown nicely?

Concerning this:

# Adding this constant background components the fit works with cash statistics as well
#spatial_model_bkg = Const2D('spatial-model-bkg')

My preference would be that you add this for now then and show the Cash stat as an example.
We probably don't want to spend time looking at numbers from the chi2 fit, so a working cash fit is a better starting point (even if the numbers should be more off for now), in my opinion.

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Mar 7, 2016

Member

Concerning priorities, my suggestion would be:

  1. implement the little things I pointed out above
  2. add test that runs the high-level example and asserts some fitted numbers
  3. add link to https://github.com/gammapy/gammapy-extra/blob/master/experiments/sherpa_cube_analysis_prototype/sherpa_cube_analysis_prototype.ipynb to the docs as minimal documentation
  4. merge pull request

Then it's less clear what should have priority, maybe like this?
5. Wrap the functionality into a CubeAnalysis class that holds counts, exposure, background, psf and model and fit result as data members and provides some convenience methods e.g. for debug plots. The most important debug plots are residual image for a given energy band (full energy band by default) and observed and predicted counts spectra for a given region (full field of view by default) ... see examples here.
6. Implement PSF convolution
7. Apply to 2FHL Crab example and see if you get correct results
8. Add IACT background model

Only then would I say we think about / how to expose this via a command line tool.

Member

cdeil commented Mar 7, 2016

Concerning priorities, my suggestion would be:

  1. implement the little things I pointed out above
  2. add test that runs the high-level example and asserts some fitted numbers
  3. add link to https://github.com/gammapy/gammapy-extra/blob/master/experiments/sherpa_cube_analysis_prototype/sherpa_cube_analysis_prototype.ipynb to the docs as minimal documentation
  4. merge pull request

Then it's less clear what should have priority, maybe like this?
5. Wrap the functionality into a CubeAnalysis class that holds counts, exposure, background, psf and model and fit result as data members and provides some convenience methods e.g. for debug plots. The most important debug plots are residual image for a given energy band (full energy band by default) and observed and predicted counts spectra for a given region (full field of view by default) ... see examples here.
6. Implement PSF convolution
7. Apply to 2FHL Crab example and see if you get correct results
8. Add IACT background model

Only then would I say we think about / how to expose this via a command line tool.

@cdeil cdeil changed the title from Sherpa cube analysis prototype to Add Sherpa cube analysis prototype Mar 14, 2016

@cdeil cdeil referenced this pull request Mar 14, 2016

Closed

Add Sherpa 3D modeling prototype #365

0 of 7 tasks complete
@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Mar 14, 2016

Member

Rebase, please.

Member

cdeil commented Mar 14, 2016

Rebase, please.

adonath added a commit that referenced this pull request Mar 16, 2016

@adonath adonath merged commit af3245f into gammapy:master Mar 16, 2016

2 checks passed

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

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Mar 16, 2016

Member

🎉

If I want to try this out, what do I do?

How about putting the equivalent of gammapy/cube/tests/test_sherpa_.py into a examples/example_cube_fit.py script and then including it on a subpage here:
https://gammapy.readthedocs.org/en/latest/tutorials/index.html
or directly here in a getting started section?
https://gammapy.readthedocs.org/en/latest/cube/index.html

Member

cdeil commented Mar 16, 2016

🎉

If I want to try this out, what do I do?

How about putting the equivalent of gammapy/cube/tests/test_sherpa_.py into a examples/example_cube_fit.py script and then including it on a subpage here:
https://gammapy.readthedocs.org/en/latest/tutorials/index.html
or directly here in a getting started section?
https://gammapy.readthedocs.org/en/latest/cube/index.html

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