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

Improve spectrum analysis #403

Merged
merged 22 commits into from Jan 15, 2016
Merged

Improve spectrum analysis #403

merged 22 commits into from Jan 15, 2016

Conversation

@joleroi
Copy link
Contributor

@joleroi joleroi commented Dec 14, 2015

This PR introduces a new background method to the spectral fitting utilities: OFF vectors from Background models.

The appropriate statistics in Sherpa (WStat vs. cash) is not yet available. Default stat is set to CASH

  • add method to compute OFF vec from background model
  • add debug plots
  • add functionality to exclude runs without sufficient number of reflected regions
  • add copy function to SpectrumAnalysis (to try different bkg methods)
@joleroi joleroi self-assigned this Dec 14, 2015
@joleroi joleroi added the feature label Dec 14, 2015
@joleroi joleroi added this to the 0.4 milestone Dec 14, 2015
@joleroi joleroi force-pushed the joleroi:off_from_background branch 2 times, most recently from efbb64c to c52fc71 Dec 14, 2015
@joleroi joleroi force-pushed the joleroi:off_from_background branch from b1491c4 to 4ed117e Dec 18, 2015
@joleroi
Copy link
Contributor Author

@joleroi joleroi commented Dec 18, 2015

I would like to merge this today. Comments welcome

@joleroi
Copy link
Contributor Author

@joleroi joleroi commented Dec 18, 2015

@cdeil : one question - in gammapy-extra/datasets/hess-crab4_pha I have some example PHA data. In order to use this for tests I set the reference to the RMF and ARF files in the PHA header relative to $GAMMAPY_EXTRA.
Now I get this error
https://travis-ci.org/gammapy/gammapy/jobs/97637674#L1431

Can I in principle use gammapy-extra this way or does it only work via gammapy.datasets.gammapy-extra. If the latter is the case, I don't know how to fix this test

@joleroi
Copy link
Contributor Author

@joleroi joleroi commented Dec 18, 2015

@cdeil: I found the error, it was a bug on my side

@cdeil
Copy link
Member

@cdeil cdeil commented Jan 7, 2016

@joleroi - Rebase on master, please.

@joleroi joleroi force-pushed the joleroi:off_from_background branch from b2f3b12 to de0ffa9 Jan 7, 2016
@joleroi
Copy link
Contributor Author

@joleroi joleroi commented Jan 7, 2016

Tests still fail because of #408 . I guess that used to happen only locally, but I guess now it's also a problem for travis (which is a good thing I guess). I have no idea how to attack this, though.

@joleroi
Copy link
Contributor Author

@joleroi joleroi commented Jan 8, 2016

Any comments before I merge?

@property
def area(self):
"""Circle Area
"""

This comment has been minimized.

@cdeil

cdeil Jan 8, 2016
Member

This formula isn't obvious.

Please leave a comment that it's the spherical cap area and a link to:
https://en.wikipedia.org/wiki/Solid_angle#Cone.2C_spherical_cap.2C_hemisphere
https://en.wikipedia.org/wiki/Versine#Definitions
Maybe even putting this in the docstring is better than as a code comment?

This should be equivalent and slightly simpler / faster, no?

val = 2 * np.pi * (1 - np.cos(self.radius))
@@ -50,3 +52,11 @@ def test_sky_to_pix_to_sky(self):
assert_allclose(sky.pos.l, sky2.pos.l)
assert_allclose(sky.pos.b, sky2.pos.b)
assert_allclose(sky.radius, sky2.radius)

def test_area(self):
pos = SkyCoord.from_name('crab')

This comment has been minimized.

@cdeil

cdeil Jan 8, 2016
Member

The SkyCoord.from_name method makes a query to the internet SESAME service.
This is better avoided in tests.

You could hard-code the Crab position or choose some random other position, no?

>>> from astropy.coordinates import SkyCoord
>>> SkyCoord.from_name('crab')
<SkyCoord (ICRS): (ra, dec) in deg
    (83.633083, 22.0145)>
pos = SkyCoord.from_name('crab')
rad = Angle('1.3451 deg')
reg = SkyCircleRegion(pos, rad)
desired = (rad ** 2 * np.pi).to('steradian')

This comment has been minimized.

@cdeil

cdeil Jan 8, 2016
Member

Leave a comment that this is the small angle approximation, and the method is the spherical cap area.

from ..image import ExclusionMask
from astropy.wcs.utils import skycoord_to_pixel

from gammapy.extern import pathlib

This comment has been minimized.

@cdeil

cdeil Jan 8, 2016
Member

Let's use relative explicit imports everywhere. I.e. here:

from ..extern.pathlib import Path

I think this import is superfluous, no?

from gammapy.extern import pathlib

This comment has been minimized.

@cdeil

cdeil Jan 14, 2016
Member

@cdeil : Do you know how to teach PyCharm to make import relative when optimizing imports?

I don't think that's possible. It probably would have to be added as an option under Editor -> General -> Auto Import -> Python -> Preferred import style.
The only feature request I could find is this, which at a very quick glance looks different:
https://youtrack.jetbrains.com/issue/PY-9987

I'd also like to have that option ... can you take the time to file a feature request in the PyCharm tracker? Or should I?

@@ -40,6 +46,33 @@ def main(args=None):
run_spectral_fit_using_config(config)


def filter_reflected_regions_analysis(analysis, n_min):

This comment has been minimized.

@cdeil

cdeil Jan 8, 2016
Member

Make this a method on gammapy.spectrum.SpectrumAnalysis?

I know having more modularity and not aggregating all functionality in *Analysis is good,
but if the standalone function takes one single *Analysis object as input to do something,
there's no advantage over making it a method, no?
(I'd claim the disadvantage is that it's harder to find / understand for users if such things are split into separate functions in the gammapy.spectrum namespace.)

Returns
-------
mask : `np.array`

This comment has been minimized.

@cdeil

cdeil Jan 8, 2016
Member

I usually use mask for boolean mask arrays and idx or index for index arrays.
I think that's also the case in numpy / scipy / ....

So may I suggest you call this idx and call numpy.nonzero here instead of numpy.where?

def total_alpha(self):
"""Averaged exposure ratio between ON and OFF regions
Alpha for all observations is calculated as :

This comment has been minimized.

@cdeil

cdeil Jan 8, 2016
Member

Did you want to add a formula here?

self.energy_dispersion.write(str(outdir/rmffile), energy_unit='keV',
clobber=clobber)

def plot_exclusion_mask(self, size=None, **kwargs):

This comment has been minimized.

@cdeil

cdeil Jan 8, 2016
Member

Pass ax so that it's possible to easily create this plot as one of several panels in a larger Figure?

This comment has been minimized.

@joleroi

joleroi Jan 14, 2016
Author Contributor

I had a discussion about this with Axel, the problem is that I need a WCS axis object, right?
This object has to know the transformation defined in the Exclusion mask header apparently due to some binning issue. So, how can I pass an axis and then give it the WCS info from the exclusion mask?

This comment has been minimized.

@cdeil

cdeil Jan 14, 2016
Member

I vaguely remember discussing this in the past as well and don't know if we found a good solution.
Maybe split this out in a separate issue and leave as-is for now?


self.pha = [make_path(f) for f in pha]
self._model = None
self._thres_lo = None
self._thres_hi = None
self.statistic = stat

This comment has been minimized.

@cdeil

cdeil Jan 8, 2016
Member

Call it consistently statistic or stat?
I think Sherpa uses stat in most places, no?


self._model = model

@property
def statistic(self):

This comment has been minimized.

@cdeil

cdeil Jan 8, 2016
Member

In __init__ you called the data member statistic.
Now you're creating a property statistic that accesses a data member _statistic.
Is this a bug?

This comment has been minimized.

@joleroi

joleroi Jan 14, 2016
Author Contributor

Yes

@@ -1,10 +1,17 @@
# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
from __future__ import absolute_import, division, print_function, \
unicode_literals

This comment has been minimized.

@cdeil

cdeil Jan 8, 2016
Member

Let's keep this on one line.
We have > 100 of those lines in Gammapy, and I put them all on one line.

from numpy.testing import assert_equal

from gammapy.datasets import gammapy_extra

This comment has been minimized.

@cdeil

cdeil Jan 8, 2016
Member

Use relative explicit imports consistently within Gammapy.

@@ -60,28 +62,28 @@ def test_spectral_fit(tmpdir):
fit.energy_threshold_low = '100 GeV'
fit.energy_threshold_high = '10 TeV'
fit.run(method='sherpa')
assert_allclose(fit.model.gamma.val, 2.0, rtol=1e-1)
assert_allclose(fit.model.gamma.val, 2.23, rtol=1e-1)

This comment has been minimized.

@cdeil

cdeil Jan 8, 2016
Member

It's great to now have first high-level tests passing on travis-ci:
https://travis-ci.org/gammapy/gammapy/jobs/101052536#L1449

For spectral index, I'd find atol more intuitive than rtol. Can you please change to atol here?

@cdeil
Copy link
Member

@cdeil cdeil commented Jan 8, 2016

I've left a few inline comments.

Maybe change the PR title from "OFF vector from cube background model" to "Improve spectrum analysis" since this PR now does several important improvements?

Also this comment in the PR description is cryptic:

The appropriate statistics in Sherpa (WStat vs. cash) is not yet available. Default stat is set to CASH

When using a background spectral template extracted from stacked off runs, using CASH is the right thing to do, and it's available, and you're using it. Maybe rephrase?

If you have time, please look over the coverage of the files you've worked on:
https://coveralls.io/builds/4663382
Maybe this could be improved a bit?
https://coveralls.io/builds/4663382/source?filename=gammapy%2Fscripts%2Fspectrum_pipe.py

Wouldn't it be better to split the SpectrumAnalysis.make_comparison_plot functionality out into a standalone function or class? I.e. have SpectrumAnalysis only be responsible for running one analysis, have a way to serialise the results, and then have separate functions / classes to compare analyses?
Of course it's OK to do this in future PRs ... I just saw that this here for the first time:
https://coveralls.io/builds/4663382/source?filename=gammapy%2Fscripts%2Fspectrum_pipe.py#L139

@joleroi
Copy link
Contributor Author

@joleroi joleroi commented Jan 11, 2016

Thank for the comments, I will go trough them today. About your general comments:

  • Concerning the statistics issue: At the moment ONLY cash is availabe for all spectral fits. The Version of sherpa I have installed (4.7) does not have WStat, and I don't manage to install the head version
checking for flex... flex
checking lex output file root... configure: error: cannot find output from flex; giving up
configure: error: ./configure failed for region-4.7

This looks related https://gist.github.com/cdeil/39e9404e91954cb67270. Suggestions?

  • Continued ...

When using a background spectral template extracted from stacked off runs, using CASH is the right thing to do

Yes and I don't know, the background model gives non-integer couns, sherpa only understands PHA (i.e. integer counts). I am happy for any suggestions what to do 😁

  • Make comparison plot

Wouldn't it be better to split the SpectrumAnalysis.make_comparison_plot functionality out into a standalone function or class?

Of course this is just something I hacked together. It will go away soon!

@cdeil
Copy link
Member

@cdeil cdeil commented Jan 14, 2016

Do you mean on the CI servers?

You have to rebase this branch:

screen shot 2016-01-14 at 09 46 00

If you mean local testing, I could come by later today and help debug it.

@joleroi
Copy link
Contributor Author

@joleroi joleroi commented Jan 14, 2016

I meant on the CI servers. Was it always the case that the need for rebase made the test not run? But anyhow, thanks!

@joleroi joleroi force-pushed the joleroi:off_from_background branch from 000d201 to 433476e Jan 14, 2016
@joleroi
Copy link
Contributor Author

@joleroi joleroi commented Jan 14, 2016

Can I merge?

@cdeil
Copy link
Member

@cdeil cdeil commented Jan 14, 2016

👍 to merge!

(I didn't review this in detail or try it out because I don't have time at the moment. So I might have comments in the future.)

joleroi added a commit that referenced this pull request Jan 15, 2016
Improve spectrum analysis
@joleroi joleroi merged commit da10a84 into gammapy:master Jan 15, 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
@joleroi joleroi deleted the joleroi:off_from_background branch Jan 15, 2016
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

2 participants
You can’t perform that action at this time.