Skip to content
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 energy dispersion for 3D spectral analysis #827

Merged
merged 7 commits into from Jan 10, 2017

Conversation

@JouvinLea
Copy link

@JouvinLea JouvinLea commented Dec 19, 2016

@adonath
Hi,
This PR is for taking into account the true energy for the spectral analysis.
For taking into account this true energy you have to reshape your model in the good shape to have your spectral and spatial model in true energy and then convolve by the energy dispersion. When you use the to_sherpa_data3d in the class SkyCube, you can not store the true energy since the elo and eli must have the same shape than the data. The data are the counts cube in reconstructed energy. So when you want to take into account the true energy you have to do it directly in the method calc of the object CombinedModel3DIntConvolveEdisp.
Therefore, in order to have the same philosophy in CombinedModel3DIntConvolveEdispand in CombinedModel3DInt, I changed a little bit the class CombinedModel3DInt. More over you gain a little bit in time to not compute the spectral model for all the (x,y) and to not compute the spatial_model for the energies.

I hope this is clear.
I don't know what I can add for test since we don't have a model that take xlo,xhi,ylo,yhi,elo,ehi and that I am using the mine locally from your Gauss2dInt.
But I test it with my script on the crab for HAP-FR data and I have the same result with 1D analysis and 3D analysis if I take into account the energy resolution.
cheers,
Lea

@cdeil cdeil requested a review from adonath Dec 19, 2016
@cdeil cdeil added the feature label Dec 19, 2016
@cdeil cdeil added this to the 0.6 milestone Dec 19, 2016
@cdeil
Copy link
Member

@cdeil cdeil commented Dec 19, 2016

I don't know what I can add for test since we don't have a model that take xlo,xhi,ylo,yhi,elo,ehi and that I am using the mine locally from your Gauss2dInt.

If there's no test, it's very hard to get the code working for others and improving it in the future.

If the issue is that a model is missing, why not add the model?
(either in the Gammapy library, or at the top of the test file if it's not appropriate for Gammapy for some reason).

I won't have time to review this before Christmas. I'm unsubscribing from this issue now.
@adonath - Can you please have a look?

@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Dec 20, 2016

@adonath
Hi,
Maybe Christoph have a good point I can add at the top of the test file the model and like this use Combined3dInt and Combined3dIntETrue...I will send you on slack my file for the Norm2dGaussInt model.

1)I don't understand what is the pb on Travis....
2)For the test since I will have the Cube analysis now I can use the outputs of the test test_cube_pipe.py. Do you think I have to put the output in the repository GAMMAPY_EXTRAand then run a test in test_sherpa_.py from these outputs or I re run the cube analysis in the test_sherpa_.py and then do a test on the class for spectral 3D analysis.

Cheers,
Lea

@adonath
Copy link
Member

@adonath adonath commented Dec 20, 2016

@JouvinLea I'd much prefer the solution to put the test files into https://github.com/gammapy/gammapy-extra/tree/master/test_datasets/cube (which has to be created...), and then read it for the test.

The Travis-CI fail seems to be unrelated.

@JouvinLea JouvinLea force-pushed the JouvinLea:AddTrueEnergy_For3DSpectralAnalysis branch 2 times, most recently from 7e90b95 to 0728a05 Dec 22, 2016
@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Dec 22, 2016

@adonath @cdeil
I add a test for this PR with cube I add in $GAMMAPY_EXTRA. As I said, the fit didn't converge to reasonable value due to the psf in the dummy data... If you use the counts, bkg, exposure and rmf in the $GAMMAPY_EXTRA/test_datasets/cube with a mean psf cube of these 4 crab runs with real data from PA or HAP-FR it converge for the position on right values.
The test will only work if the PR in GAMMAPY_EXTRA is merge since I change de energy binning.

Cheers,
LEa

@cdeil
Copy link
Member

@cdeil cdeil commented Jan 6, 2017

I've merged gammapy/gammapy-extra#47 just now, and I'm re-starting the continuous integration tests now.

@adonath - Do you have time to have another look at this? Or should I?

@cdeil
Copy link
Member

@cdeil cdeil commented Jan 6, 2017

Tests are currently failing if the optional dependency sherpa is not available:
https://ci.appveyor.com/project/cdeil/gammapy/build/1.0.1820/job/o8cswojtq165e6um#L1349

This happens because you added module-level (i.e., at the top) sherpa imports to gammapy/cube/tests/test_sherpa_.py, and the pytest test collection imports all files starting with test_.

We could configure the test runner to treat that file in a special way.
But I think maybe the simplest solution is that you move the NormGauss2DInt class and all the top-level Sherpa imports to gammapy/cube/sherpa_.py.
@JouvinLea - Can you try that, please?

@JouvinLea JouvinLea force-pushed the JouvinLea:AddTrueEnergy_For3DSpectralAnalysis branch from 0728a05 to 0c86225 Jan 9, 2017
@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Jan 9, 2017

@cdeil
Happy new year!
I tried what you proposed but it is still failing due to the import sherpa but this time in the class sherpa_.py

Copy link
Member

@adonath adonath left a comment

@JouvinLea I've left a few inline comments.

I think in general the code structure could be improved a lot. Having all computations of PSF convolution and energy dispersion in a single calc() method of one model, including a lot of ravel() and .reshape() calls, makes the code very confusing. Just by looking at the code I can't tell if does the right thing at all.

But judging from the tests, the code seems to work, so I'm fine with merging the code as is, after the minor comment have been addressed. Thanks for this contribution!

@@ -1,13 +1,68 @@
# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
import numpy as np
from ..utils.energy import EnergyBounds
import sherpa.astro.ui as sau

This comment has been minimized.

@adonath

adonath Jan 9, 2017
Member

Can you remove this import? We don't make use of sherpa's high level API...

from ..utils.energy import EnergyBounds
import sherpa.astro.ui as sau
from sherpa.astro.ui import erf
import astropy.wcs as pywcs

This comment has been minimized.

@adonath

adonath Jan 9, 2017
Member

Can you just explicitly import WCS here? I've never seen that convention of using pywcs

- erf.calc.calc([1, p[1], sigma_erf], ylo)))


sau.add_model(NormGauss2DInt)

This comment has been minimized.

@adonath

adonath Jan 9, 2017
Member

This is not needed, as we register the model with sherpa's high level API

self.n_ebins = None
ArithmeticModel.__init__(self, name, (self.xpos, self.ypos, self.ampl, self.fwhm))

def set_wcs(self, wcs):

This comment has been minimized.

@adonath

adonath Jan 9, 2017
Member

I don't see, that this method is being used somewhere (especially in calc()...), is that intended?

ArithmeticModel.__init__(self, name, pars)

def calc(self, pars, elo, xlo, ylo, ehi, xhi, yhi):
from scipy import signal

This comment has been minimized.

@adonath

adonath Jan 9, 2017
Member

I don't think this is the right place for the import statement, because it's executed in every iteration of the fit, can you move to __init__() and store the fftconvolve as e.g. self._fftconvolve

spatial_model=spatial_model, spectral_model=spectral_model)

# Set starting values
center = SkyCoord.from_name("Crab").galactic

This comment has been minimized.

@adonath

adonath Jan 9, 2017
Member

Can you rather set the position manually, instead of querying a web database using from_name(). Otherwise this will fail, if users don't have internet connection (and might cause random fails on Travis-CI...)

@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Jan 9, 2017

@adonath
Thank for your comments!! I agree that the code is not readable easily. Maybe be if the test pass we can merge just a way for me to continue my work with 3D analysis but we can discuss how to improve it for later if you have ideas.

@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Jan 9, 2017

@adonath
I don't understand the problem with the test for the ImportError: No module named 'sherpa. You were already calling some package from sherpa before so why it is failing here?

@JouvinLea JouvinLea force-pushed the JouvinLea:AddTrueEnergy_For3DSpectralAnalysis branch from cf7a9d1 to 6543804 Jan 9, 2017
from ...utils.testing import requires_dependency, requires_data
from .. import SkyCube
from ..sherpa_ import NormGauss2DInt

This comment has been minimized.

@adonath

adonath Jan 10, 2017
Member

@JouvinLea I think this line causes the test fails on Travis-CI. The import should be moved to testCombinedModel3DInt() and testCombinedModel3DIntConvolveEdisp(), as you did for the other imports.

Lea Jouvin
@adonath adonath merged commit 455d1da into gammapy:master Jan 10, 2017
2 checks passed
2 checks passed
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
@adonath
Copy link
Member

@adonath adonath commented Jan 10, 2017

Merged. Thanks @leajouvin!

@cdeil cdeil changed the title Add energy true for 3D spectral analysis Add energy dispersion for 3D spectral analysis Jan 10, 2017
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

3 participants