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

Stack EDISP for a set of observations #772

Merged
merged 12 commits into from Nov 16, 2016
Merged

Conversation

@JouvinLea
Copy link

@JouvinLea JouvinLea commented Nov 14, 2016

@cdeil
This is a PR for stacking the RMF for a set of observation and a given source position. This is coded in the same philosophy than the mean PSF. You will need it for 3D spectral analysis if you want to convert your exposure in true energy to reco energy.

Copy link
Member

@cdeil cdeil left a comment

@joleroi - Thanks!

Why does the PR title say "stack RMF for images"? There's no image involved in this PR, no?

Left two inline comments.

@@ -711,3 +713,59 @@ def make_psf(self, position, energy=None, theta=None):
psf_tot = EnergyDependentTablePSF(energy=energy, offset=theta, exposure=exposure,
psf_value=psf_value.T)
return psf_tot

def make_mean_rmf(self, position, e_true, e_reco, low_reco_threshold=None, high_reco_threshold=None):

This comment has been minimized.

@cdeil

cdeil Nov 14, 2016
Member

You call the method "make_mean_rmf", but return an "EnergyDispersion" object.

My understanding is that "rmf" is "Redistribution Matrix File", a file format:
http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html#tth_sEc3

So I'd be +1 to use edisp for "energy dispersion" in Gammapy, and to only use "RMF" in docs or on file I/O,
i.e. call this method "make_mean_edisp".

This comment has been minimized.

@joleroi

joleroi Nov 14, 2016
Contributor

👍

Parameters
----------
position : `~astropy.coordinates.SkyCoord`
Position at which to compute the PSF

This comment has been minimized.

@cdeil

cdeil Nov 14, 2016
Member

Copy & paste error: "PSF"-> "EDISP" (or put something shorter, less specific in both docstrings).


assert len(rmf.e_true.nodes) == 80
assert len(rmf.e_reco.nodes) == 15
assert_quantity_allclose(rmf.data[53, 8], 0.0559785805550798)

This comment has been minimized.

@cdeil

cdeil Nov 14, 2016
Member

Assert on thresholds and / or behaviour near threshold for the stacked EDISP?

@cdeil cdeil added the feature label Nov 14, 2016
@cdeil cdeil added this to the 0.5 milestone Nov 14, 2016
@JouvinLea JouvinLea changed the title Stack RMF for images Stack RMF for a set of observations Nov 14, 2016
@joleroi
Copy link
Contributor

@joleroi joleroi commented Nov 14, 2016

@cdeil I don't know why you adressed me in the first comment, I don't have anything to do with this PR

@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Nov 14, 2016

@cdeil
I took your comments into account and add a test for the thresholds! Tell me if it's ok.

Copy link
Member

@cdeil cdeil left a comment

I left some more inline comments. That's all from my side, @joleroi I'll leave final review and merging to you.

@@ -150,3 +151,38 @@ def test_make_psf(pars, result):
assert_quantity_allclose(psf.energy[10], result["psf_energy"])
assert_quantity_allclose(psf.exposure[10], result["psf_exposure"])
assert_quantity_allclose(psf.psf_value[10, 50], result["psf_value"])


def test_make_meanedisp(tmpdir):

This comment has been minimized.

@cdeil

cdeil Nov 15, 2016
Member

You have to mark this test

@requires_data('gammapy-extra')

https://travis-ci.org/gammapy/gammapy/jobs/175732342#L1722

And if it uses scipy or some other optional dependency, also

@requires_dependency('scipy')

Also, please rename: test_make_meanedisp -> test_make_mean_edisp

low_reco_threshold=Energy(1, "TeV"), high_reco_threshold=Energy(60, "TeV"))
# Test that the energy dispersion is equal to zero for the reco energy bin inferior to low_reco_threshold and
# superior to high_reco_threshold
i2 = np.where(rmf2.data[:, 13] != 0)[0]

This comment has been minimized.

@cdeil

cdeil Nov 15, 2016
Member

I find this hard to read. I don't know what index 13 is by looking at the test.
Isn't there an evaluate method or something like that where you can pass an energy?
That would make the test higher-level, i.e. it would be clear that you evaluate EDISP below or above the threshold from the code (what you now have in the comment).

list_arf = list()
list_edisp = list()
list_livetime = list()
if not low_reco_threshold:

This comment has been minimized.

@cdeil

cdeil Nov 15, 2016
Member

We looked at this in detail and decided that it's OK / best to use quantities as default keyword arguments.
So here you don't have to put None and then this if statement, please put low_reco_threshold=Energy(0.002, "TeV") directly in the function signature.

@@ -711,3 +713,59 @@ def make_psf(self, position, energy=None, theta=None):
psf_tot = EnergyDependentTablePSF(energy=energy, offset=theta, exposure=exposure,
psf_value=psf_value.T)
return psf_tot

def make_mean_edisp(self, position, e_true, e_reco, low_reco_threshold=None, high_reco_threshold=None):
r"""Make mean rmf for a given position and a set of observations.

This comment has been minimized.

@cdeil

cdeil Nov 15, 2016
Member

Docstring still says "rmf" three times where IMO it would be better to say "edisp".

This comment has been minimized.

@joleroi

joleroi Nov 15, 2016
Contributor

the same argument actually hold for arf vs. aeff. arf stands for auxilliary response file and is not really a synonym for effective area.

low_reco_threshold=Energy(1, "TeV"), high_reco_threshold=Energy(60, "TeV"))
# Test that the energy dispersion is equal to zero for the reco energy bin inferior to low_reco_threshold and
# superior to high_reco_threshold
i2 = np.where(rmf2.data[:, 13] != 0)[0]

This comment has been minimized.

@JouvinLea

JouvinLea Nov 15, 2016
Author

@cdeil @joleroi
I agree this will be more clear to evaluate at a reco energy but there is no method evaluate on edisp right?

This comment has been minimized.

@joleroi

joleroi Nov 15, 2016
Contributor

Yes, there is. It's inherited from NDDataArray. It shows up when you hit ? in IPython but not in the API docs. @cdeil can we change this somehow?

This comment has been minimized.

@cdeil

cdeil Nov 15, 2016
Member

It does show up here:
http://docs.gammapy.org/en/latest/api/gammapy.utils.nddata.NDDataArray.html#gammapy.utils.nddata.NDDataArray.evaluate

This can be configured globally or per-class in Sphinx, whether inherited attributes are shown for sub-classes. At the moment the Gammapy (and I think Astropy) docs don't show inherited methods, you have to click on the link to the base class at the top.

I'm not sure if changing the default to always show inherited attributes / methods, or to do it on a case-by-case basis to get "optimal" API docs is an improvement / maintainable.

@joleroi or @JouvinLea - If you think this could / should be improved, please open a separate issue stating your preference.

This comment has been minimized.

@joleroi

joleroi Nov 15, 2016
Contributor

I'm not sure if changing the default to always show inherited attributes / methods, or to do it on a case-by-case basis to get "optimal" API docs is an improvement / maintainable.

I agree. @JouvinLea if you think we should show the inherited methods for the IRF classes please open an issue.

Copy link
Contributor

@joleroi joleroi left a comment

Some minor comments

def make_mean_edisp(self, position, e_true, e_reco, low_reco_threshold=None, high_reco_threshold=None):
r"""Make mean rmf for a given position and a set of observations.
Compute the mean rmf of a set of observations j at a given position

This comment has been minimized.

@joleroi

joleroi Nov 15, 2016
Contributor

Instead of repeating the definiton here I would just add a link to the IRFStacker (because that is where the stacking is actually implemented). You can add links to a function like this
:func:~gammapy.irf.IRFStacker.stack_edisp`

low_reco_threshold = Energy(0.002, "TeV")
if not high_reco_threshold:
high_reco_threshold = Energy(150, "TeV")
list_low_threshold = [low_reco_threshold] * len(self)

This comment has been minimized.

@joleroi

joleroi Nov 15, 2016
Contributor

Why do you hard code the same energy threshold for all observations?

high_reco_threshold = Energy(150, "TeV")
list_low_threshold = [low_reco_threshold] * len(self)
list_high_threshold = [high_reco_threshold] * len(self)
for obs in self:

This comment has been minimized.

@joleroi

joleroi Nov 15, 2016
Contributor

In the following paragraph there are many lines that are longer than the recommended length of 79 characters. This makes the code really hard to read. Of course its no problem if a line is sometimes a little bit too long, but if you can avoid it, please do so. Actually, since you are using pycharm the editor should show you line that are too long and sometimes even give you the option to wrap them automatically (but I don't really know pycharm). In this case you can simply split one long line into several statements, e.g. instead of

list_arf.append(obs.aeff.to_effective_area_table(offset, energy=e_true))
aeff = obs.aeff.to_effective_area_table(offset=offset,
                                        energy=e_true)
list_arf.append(aeff)

For a general style guide for python see https://www.python.org/dev/peps/pep-0008/#introduction, or check out the pep8 command line tool

@@ -711,3 +713,59 @@ def make_psf(self, position, energy=None, theta=None):
psf_tot = EnergyDependentTablePSF(energy=energy, offset=theta, exposure=exposure,
psf_value=psf_value.T)
return psf_tot

def make_mean_edisp(self, position, e_true, e_reco, low_reco_threshold=None, high_reco_threshold=None):
r"""Make mean rmf for a given position and a set of observations.

This comment has been minimized.

@joleroi

joleroi Nov 15, 2016
Contributor

the same argument actually hold for arf vs. aeff. arf stands for auxilliary response file and is not really a synonym for effective area.

@joleroi
Copy link
Contributor

@joleroi joleroi commented Nov 15, 2016

@JouvinLea You pushed a new commit while I was adding some comments, so my comments are shown as outdated, sorry for that.

@cdeil cdeil changed the title Stack RMF for a set of observations Stack EDISP for a set of observations Nov 15, 2016
list_aeff = list()
list_edisp = list()
list_livetime = list()
list_low_threshold = [low_reco_threshold] * len(self)

This comment has been minimized.

@JouvinLea

JouvinLea Nov 15, 2016
Author

@joleroi
needed in the method stack_edisp since you can have e_reco thresholds different for each observation in your spectral extraction.
Here for the moment I don't know what we have to use for thresholds so if a user want to put one he can!!

This comment has been minimized.

@joleroi

joleroi Nov 15, 2016
Contributor

Ok, I don't know enough about image analysis, it seemed weird to have the same threshold everywhere. Thanks for clarifying!

list_high_threshold = [high_reco_threshold] * len(self)
for obs in self:
offset = position.separation(obs.pointing_radec)
list_aeff.append(obs.aeff.to_effective_area_table(offset, energy=e_true))

This comment has been minimized.

@JouvinLea

JouvinLea Nov 15, 2016
Author

@joleroi
For me this is ok on pycharm for the line format it's not too long....

This comment has been minimized.

@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Nov 15, 2016

@joleroi
I think I took most of your comments into account!
I don't understand the failing error for travis

@joleroi
Copy link
Contributor

@joleroi joleroi commented Nov 15, 2016

There is some error in the docstring of make_mean_edisp
https://travis-ci.org/gammapy/gammapy/jobs/176002042#L1665
I think you're just missing a ``` in line 730

@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Nov 15, 2016

@joleroi
It's failing on travis but me locally the docs is ok!

@joleroi
Copy link
Contributor

@joleroi joleroi commented Nov 15, 2016

When I run the sphinx build locally I also get this failure
https://travis-ci.org/gammapy/gammapy/jobs/176076832#L1665
You have to run

python setup.py build_sphinx -w 

in order to raise an error if there is a sphinx warning. Otherwise it just says

build succeeded, 1 warning.
@joleroi
Copy link
Contributor

@joleroi joleroi commented Nov 15, 2016

I don't have time today, but if you don't succeed I can have a look tomorrow. Sphinx errors are sometimes hard to track down.

theta : `~astropy.coordinates.Angle`
1-dim offset array for the output PSF.
If none is given, the energy array of the PSF from the first observation is used.
If none is given, the energy array of the PSF from the first
observation is used.

This comment has been minimized.

@cdeil

cdeil Nov 15, 2016
Member

Maybe this line is causing the Sphinx / RST error?
Try removing a space here, i.e. indent these two lines exactly the same:

If none is given, the energy array of the PSF from the first
 observation is used.
@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Nov 16, 2016

@cdeil
You were right it fixed the issue the line line you mentioned!!
Can we merge?

@joleroi joleroi merged commit 48acf76 into gammapy:master Nov 16, 2016
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
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