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 sensitivity computation #1053

Merged
merged 5 commits into from Jul 10, 2017

Conversation

Projects
None yet
5 participants
@bkhelifi
Copy link
Contributor

bkhelifi commented May 31, 2017

@cdeil
My first PR thanks to @jjlk: the sensitivity script and its test script.
Any comment or suggestion?
Thanks

@cdeil cdeil self-assigned this May 31, 2017

@cdeil cdeil added the feature label May 31, 2017

@cdeil cdeil added this to the 0.7 milestone May 31, 2017

@cdeil
Copy link
Member

cdeil left a comment

@bkhelifi @jjlk - Thanks for the addition!

I left inline comments on coding style stuff. I think most of them are straightforward, but if you have any questions or disagree, let me know.

One more thing is that as-is, this class will not show up in the HTML docs, i.e. user's won't find it. To get it to show up you need to do two things:

  1. list the class name in an __all__ list at the top of the gammapy/scripts/cta_sensitivity.py file.
  2. Add from .cta_sensitivity import * to the gammapy/scripts/__init__.py file.

If you want, you can then run this locally to check how your class appears in the Gammapy docs:

python setup.py build_docs # takes a two or three minutes, unfortunatly
open docs/_build/html/index.html

My other non-trivial comment would be that the get_excess method is too long. If possible it should be re-factored into smaller parts that can be read and tested more easily. But that could be left for another round of review, or, depending on your time for coding in the coming days, a future pull request.

@@ -0,0 +1,250 @@
#-*- coding: utf-8 -*-

This comment has been minimized.

@cdeil

cdeil May 31, 2017

Member

Remove line with #-*- coding: utf-8 -*-, add the two lines that we have at every file in Gammapy (and Astropy):

# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
#-*- coding: utf-8 -*-

"""
Compute the CTA sensitivity for a fixed observation config

This comment has been minimized.

@cdeil

cdeil May 31, 2017

Member

Suggest to remove this module-level docstring, and merge it's content into the docstring of the class.
(to a certain degree it's already a duplication, i.e. hard to maintain, but there's some things that should be copied over)

We don't put author / date / version information in files anywhere in Gammapy.
We use the git commit history, as well as
http://docs.gammapy.org/en/latest/about.html#contributors
and
http://docs.gammapy.org/en/latest/changelog.html#changelog
(I update it from time to time, and always before making a release).
Maintaining author / date / version info in files (we have > 100 in Gammapy) would be very difficult and a lot of work to keep accurate.

"""

# Inputs :

This comment has been minimized.

@cdeil

cdeil May 31, 2017

Member

Remove these comments?
They seem to be another variant of the class documentation.
Let's only have one version, in the class docstring.

import astropy.units as u
import numpy as np
import matplotlib.pyplot as plt
import sys

This comment has been minimized.

@cdeil

cdeil May 31, 2017

Member

This import sys can be removed, it's unused.


import astropy.units as u
import numpy as np
import matplotlib.pyplot as plt

This comment has been minimized.

@cdeil

cdeil May 31, 2017

Member

Matplotlib is an optional dependency. It should be possible to compute a sensitivity curve without matplotlib.
For this to work, you have to put the import matplotlib.pyplot as plt as a delayed import, into the plot method below.



if __name__ == '__main__':
cta_sensitivity_main()

This comment has been minimized.

@cdeil

cdeil May 31, 2017

Member

Add newline at end of file.


@requires_dependency('scipy')
@requires_data('gammapy-extra')

This comment has been minimized.

@cdeil

cdeil May 31, 2017

Member

Remove empty line here.

from astropy.tests.helper import assert_quantity_allclose
from numpy.testing import assert_allclose
from astropy.units import Quantity
from gammapy.utils.testing import requires_data, requires_dependency

This comment has been minimized.

@cdeil

cdeil May 31, 2017

Member

Use relative imports within Gammapy (even from the tests) for consistency with the rest of the codebase.

def test_cta_sensitivity():
"""Test that CTA sensitivity can be well evaluated."""

irffile = '$GAMMAPY_EXTRA/datasets/cta/CTA-Performance-North-20170327/CTA-Performance-North-20deg-average-30m_20170327.fits'

This comment has been minimized.

@cdeil

cdeil May 31, 2017

Member

This file doesn't exist in gammapy-extra. Can you pick an existing file or add this file?

sens.run()

true_value =6.3596014545e-12
actual = sens.diff_sens[np.where(np.abs(sens.energy.value - 1.) < 0.2)][0].value

This comment has been minimized.

@cdeil

cdeil May 31, 2017

Member

This line where you pick a value is very hard to read.
I suggest you access the results as a table (as I suggested in the implementation), and then just assert on the values of a given row or two of the table.

@bkhelifi

This comment has been minimized.

Copy link
Owner

bkhelifi commented on 5f16435 Jun 9, 2017

Modification of the code based on the @cdeil comments

@cdeil

This comment has been minimized.

Copy link
Member

cdeil commented Jun 20, 2017

@bkhelifi - Sorry for not responding here for a few weeks!

I've done some code cleanup in 0fea216 and pushed that commit here to your branch. I.e. if you pull from that branch you should get that version and be able to continue with it.

There's at least one change that will require an update of gammapy/gammapy-extra#80 : I changed __init__ to only take a livetime quantity, no longer a float without unit information.

The main remaining thing before this can go in is to fix the test. At least for me, it fails like this, complaining about missing SPECRESP:
https://gist.github.com/cdeil/908864a9aa1f226dd3b2536606225b76

@bkhelifi or @jjlk - Can you fix this, please? (I don't care how: either add a new data file that has that extension, or change the code here or in CTAPerf read to make it optional). Would be nice to have this functionality and example notebook in for the Gammapy presentation at the CTA ASWG bootcamp on Thursday.

@mackaiver
Copy link
Contributor

mackaiver left a comment

Nice work. I added a few questions here in there.

return bkg.value * u.Unit('')

def get_excess(self, bkg_counts):
"""Compute the gamma-ray excess for each energy bin.

This comment has been minimized.

@mackaiver

mackaiver Jun 21, 2017

Contributor

Could you provide a short (high level) description of how it does what it does? The code is not very obvious to me.

__all__ = ['SensiEstimator']


class SensiEstimator(object):

This comment has been minimized.

@mackaiver

mackaiver Jun 21, 2017

Contributor

Can we call it SensitivityEstimator? Not that much more to type.

if log.getEffectiveLevel() == 10:
self._ref_diff_sensi.pprint()

def plot(self, ax=None):

This comment has been minimized.

@mackaiver

mackaiver Jun 21, 2017

Contributor

Don't know how it is handled elsewhere in gammapy but maybe it makes sense to return the axis object here?
Also might it be better to let the user call plt.show() ? I might want to save it using savefig or something.

Return an array of background rate (bins in reconstructed energy)
"""
bkg = bkg_rate.data.data * self.livetime
return bkg.value * u.Unit('')

This comment has been minimized.

@mackaiver

mackaiver Jun 21, 2017

Contributor

This is returning a rate according to the docstring. Should this be unitless then?

return flux

def run(self):
"""Do the computation."""

This comment has been minimized.

@mackaiver

mackaiver Jun 21, 2017

Contributor

Maybe expand the docstring here a bit?

@mackaiver

This comment has been minimized.

Copy link
Contributor

mackaiver commented Jun 21, 2017

Also I noticed there are some functions in gammapy.scripts handling cta irf related stuff. It seems like there is some code duplication here. e.g. the plotting of the sensitvity curve.

BTW I was confused why things like cta_irf nad cta_utils appear in a 'scripts' folder. They don't seem to be scripts.

@jjlk

This comment has been minimized.

Copy link
Contributor

jjlk commented Jun 23, 2017

Hi @mackaiver ,

Also I noticed there are some functions in gammapy.scripts handling cta irf related stuff. It seems like there is some code duplication here. e.g. the plotting of the sensitvity curve.

BTW I was confused why things like cta_irf nad cta_utils appear in a 'scripts' folder. They don't seem to be scripts.

For now we are still waiting that the simulation group in CTA releases point-like IRFs in FITS format (they are available in ROOT format, the IRF FITS files contains only full-containement responses). We wrote a small script to export IRF from ROOT to FITS with an hybrid format.

All the classes in cta_utils and cta_irf are just a convenient and a temporary way to use the IRFs and extract some ``science'' for point-like analysis. As soon as a data format for point-like IRF will be defined and the IRF released we'll use the standard things implemented in Gammapy.

@cdeil

cdeil approved these changes Jun 29, 2017

@@ -330,7 +336,7 @@ def plot(self, ax=None, energy=None, **kwargs):
ax.errorbar(energy.value, values.value, xerr=xerr, fmt='o', **kwargs)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('Energy [{}]'.format(self.energy.unit))
ax.set_xlabel('Reco Energy [{}]'.format(self.energy.unit))

This comment has been minimized.

@mackaiver

mackaiver Jul 4, 2017

Contributor

Can we call it Reconstructed Energy ?

self.diff_sens = None

def get_bkg(self, bkg_rate):
"""Return the Background rate for each energy bin

This comment has been minimized.

@mackaiver

mackaiver Jul 4, 2017

Contributor

does this return a rate? should just be counts right?

self.diff_sens = diff_flux

@property
def diff_sensi_table(self):

This comment has been minimized.

@mackaiver

mackaiver Jul 4, 2017

Contributor

rename to differential_sensitivity_table?

@bkhelifi bkhelifi merged commit 2cc3258 into gammapy:master Jul 10, 2017

1 of 2 checks passed

continuous-integration/travis-ci/pr The Travis CI build could not complete due to an error
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
@MaxNoe

This comment has been minimized.

Copy link

MaxNoe commented Jul 24, 2017

I'm confused. It says this Pull Request has been merged into master, but I cannot find the files and I also cannot find the merge commit in the history of the master.

@jjlk

This comment has been minimized.

Copy link
Contributor

jjlk commented Jul 24, 2017

Hi @MaxNoe ,
Yeap. I have no idea of what is going on, it should have been merged two weeks ago.

@cdeil , @bkhelifi : any ideas of what's going wrong?

@bkhelifi

This comment has been minimized.

Copy link
Contributor

bkhelifi commented Aug 1, 2017

Dear all,
So, is it merged or not? How can I check it?
@MaxNoe @cdeil @jjlk

@cdeil cdeil changed the title Add the sensitivity script Add sensitivity computation Aug 31, 2017

@cdeil cdeil referenced this pull request Sep 11, 2017

Merged

Re-add SensitivityEstimator #1097

@cdeil

This comment has been minimized.

Copy link
Member

cdeil commented Sep 11, 2017

@MaxNoe @jjlk @bkhelifi - The SensitivityEstimator was re-added by @bkhelifi in #1097 . Unfortunately I'm not sure exactly what happened, see #1097 (comment) . Not to excuse the git screw-up on my part, but just to note: it's a one-time mistake in several years of development, and seems to be fixed now, usually we're doing OK with git / Github for development.

@MaxNoe

This comment has been minimized.

Copy link

MaxNoe commented Sep 11, 2017

No problem :) Nobody was questioning using git here !

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