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 light curve computation #950

Merged
merged 8 commits into from Mar 17, 2017
Merged

Add light curve computation #950

merged 8 commits into from Mar 17, 2017

Conversation

@jjlk
Copy link
Contributor

@jjlk jjlk commented Mar 14, 2017

Hi @cdeil,
As we discussed, here is a first attempt to deal with light curves. Implementation follows closely what has been done in the the H.E.S.S. publication: http://adsabs.harvard.edu/abs/2010A%26A...520A..83H.

This PR:

  • add a LightCurveFactory class
  • modify LightCurve.plot and LightCurve.simulate functions (I replaced seconds by MJD values)
  • add a small test (only checking the length of the output LightCurve object since values might change after checking that the results are okés)
    ++

PS: forgot to say that there is room for improvements, different time bining definitions (only run bye run, or user defined binning are implemented) and optimisation of the code, I guess I'm doing too much loops =)

@cdeil cdeil self-assigned this Mar 14, 2017
@cdeil cdeil added the feature label Mar 14, 2017
@cdeil cdeil added this to the 0.6 milestone Mar 14, 2017
@cdeil cdeil changed the title Light curve Add light curve computation Mar 14, 2017
Copy link
Member

@cdeil cdeil left a comment

@jjlk - Thank you!

I've left a bunch of inline comments.
Please decide for yourself what you want to implement here and defer to future PRs (or reject as a suggestion), and let me know when this is ready to merge from your side.

__all__ = [
'LightCurve',
'LightCurveFactory',

This comment has been minimized.

@cdeil

cdeil Mar 14, 2017
Member

I would suggest the name "LightCurveEstimator" or "LightCurveMaker" instead of "LightCurveFactory".

My reasoning is that we have several other classes called "Estimator" or "Maker" already in Gammapy, and shouldn't invent new terms for the same thing. There's also the fact that "factory" functions or classes to many Python programmers have a specific meaning, it's a term for things that create other functions or classes:
https://en.wikipedia.org/wiki/Factory_(object-oriented_programming)

These aren't strong reasons, and @jjlk - if you would like to keep the name "LightCurveFactory", then let's do that (at least until we have a review of all names and code organisation in Gammapy before doing a 1.0 release later this year).

This comment has been minimized.

@jjlk

jjlk Mar 16, 2017
Author Contributor

Done.

Parameters
----------
spectral_extraction: `~gammapy.spectrum.SpectrumExtraction`

This comment has been minimized.

@cdeil

cdeil Mar 14, 2017
Member

Suggest to call the variable spectrum_extraction for consistency with the class name SpectrumExtraction.
Or spec_extract or something shorter if you like, but sometimes calling it "spectral" and sometimes "spectrum" makes it harder to remember.

This comment has been minimized.

@jjlk

jjlk Mar 16, 2017
Author Contributor

Done.

self.on_evt_list = self._get_on_evt_list(spectral_extraction)


def _get_off_evt_list(self, spectral_extraction):

This comment has been minimized.

@cdeil

cdeil Mar 14, 2017
Member

This _get_off_evt_list and _get_on_evt_list is a bit weird, no?
You pass in self and spectrum_extraction, but then never access self!?

Suggest to move it as a helper method on the SpectrumExtraction class if you think it might be useful to have from other places as well, not just the lightcurve computation.
If you think it's really just something needed here, then I suggest you make it a @staticmethod (see other example in Gammapy) and don't pass in self.

This comment has been minimized.

@jjlk

jjlk Mar 16, 2017
Author Contributor

Done.

return on_evt_list


def _get_time_intervals(self, t_start, t_stop, time_binning_type):

This comment has been minimized.

@cdeil

cdeil Mar 14, 2017
Member

At least at the moment the time interval computation is decoupled from the lightcurve computation, just a step one has to execute before. My guess is that it will stay that way.

So my suggestion would be to not take t_start, t_stop, t_binning_type in the lightcurve method,
and to instead somehow expose a function or method to compute time binnings.
Then in the future, there could be serveral such methods to compute time binnings (or really often this will be something a user just does in their own script for their specific use case by themselves), and the lightcurve method would always be simple and not grow with other parameters concerning time binning.

That said, I'm not sure how to best do this.
One question is which things in this class to take as input in __init__ and which in the lightcurve or "run" method.
Maybe it doesn't matter much?

Another question is how to represent the time binning.
In your current method you don't have a Returns section in the docstring, but from the code it looks like it's a Python list of scalar Time objects.
That's not so bad, but maybe using
http://docs.gammapy.org/en/latest/api/gammapy.data.GTI.html
for that purpose would be better (i.e. an Astropy table o TSTART / TSTOP values, exposed as Time attributes)?
Then for that we have I/O already and we need to develop methods anyways to e.g. merge and clip GTIs from runs, so in the future there could be conveniences to work with GTIs and that could be used for people preparing a lightcurve time binning?

The last question would be where the code to compute the time binning using the method you have here should live. I'm not sure. It could live on the LightCurve class, on the GTI class, on the ObservationList class, or in a new function or class in this file where the LightCurveFactory is.

@jjlk - Thoughts?
Do you want to implement any of these suggestions in this PR, or postpone them to possible future PRs?

This comment has been minimized.

@jjlk

jjlk Mar 16, 2017
Author Contributor

So my suggestion would be to not take t_start, t_stop, t_binning_type in the lightcurve method,
and to instead somehow expose a function or method to compute time binnings.

Done, I added a time_intervals argument.

One question is which things in this class to take as input in init and which in the lightcurve or "run" method.
Maybe it doesn't matter much?

I keep the spectral_model and energy_range arguments since they are needed to compute the light curve. My goal was to be able to compute different time-binned light curves with the light_curve function for a given spectral extraction. With this we are able to easily compute light curves with different energy threshold (e.g. >100 GeV, >1 TeV), for a given model. Maybe someone might be interested to also test different spectral model (maybe less frequent that the energy threshold use case).

Another question is how to represent the time binning.
In your current method you don't have a Returns section in the docstring, but from the code it looks like it's a Python list of scalar Time objects.
That's not so bad, but maybe using
http://docs.gammapy.org/en/latest/api/gammapy.data.GTI.html
for that purpose would be better (i.e. an Astropy table o TSTART / TSTOP values, exposed as Time attributes)?
Then for that we have I/O already and we need to develop methods anyways to e.g. merge and clip GTIs from runs, so in the future there could be conveniences to work with GTIs and that could be used for people preparing a lightcurve time binning?

For now I put aside this question, i stick with astropy.time.core.Time objects.

The last question would be where the code to compute the time binning using the method you have here should live. I'm not sure. It could live on the LightCurve class, on the GTI class, on the ObservationList class, or in a new function or class in this file where the LightCurveFactory is.

Same here, I removed the _get_time_intervals function. We'll do something fancy later.

return intervals


def light_curve(self, t_start, t_stop, t_binning_type,

This comment has been minimized.

@cdeil

cdeil Mar 14, 2017
Member

Like I said, I'd suggest to simplify the arguments of this lightcurve method to compute lightcurves, and only take user-defined time binning here.

Also, this method is waaaaaaaay too long. I'll try to make suggestions how to split it into parts that can be read and tested separately.

# Gaussian errors, ToDo: should be improved
interval_flux_error *= np.sqrt(interval_on + interval_alpha_mean**2 * interval_off)

print('### start:{}, stop:{}'.format(interval_tmin, interval_tmax))

This comment has been minimized.

@cdeil

cdeil Mar 14, 2017
Member

No printing from Gammapy allowed.

Either just remove it, or, if you think it will be useful to people, then change this class to the following pattern of storing the relevant info and exposing a print_report or print_summary method.

est = LightCurveEstimator(...)
est.run(...) # what you call `lightcurve` now
est.lightcurve # The main result, stored on the estimator class for access after it runs
est.print_report() # print out summary or debugging info after the algorithm runs

This comment has been minimized.

@jjlk

jjlk Mar 16, 2017
Author Contributor

Done. I removed the outputs, even if we are measured in debugging phase, and added interesting values in the LightCurve output.

@@ -28,3 +35,61 @@ def test_lightcurve_fvar():
def test_lightcurve_plot():
lc = LightCurve.simulate_example()
lc.plot()


@requires_data('gammapy-extra')

This comment has been minimized.

@cdeil

cdeil Mar 14, 2017
Member

If I remember correctly, then @requires_data and @requires_dependency doesn't work on fixture functions.
They only work on test functions.

So remove the line @requires_data('gammapy-extra') from the fixture function here.

This comment has been minimized.

@jjlk

jjlk Mar 16, 2017
Author Contributor

Done

@pytest.fixture(scope='session')
def spec_extraction():
data_store = DataStore.from_dir('$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2/')
obs_ids = [23523, 23526, 23559, 23592]

This comment has been minimized.

@cdeil

cdeil Mar 14, 2017
Member

How long does this test run?
Maybe restrict to two or three runs if it takes long?

This comment has been minimized.

@jjlk

jjlk Mar 16, 2017
Author Contributor

Done

obs_ids = [23523, 23526, 23559, 23592]
obs_list = data_store.obs_list(obs_ids)

target_position = SkyCoord(ra=83.63308,

This comment has been minimized.

@cdeil

cdeil Mar 14, 2017
Member

You can just use the fact that frame='icrs' is the default (we do in many places in Gammapy) and put this on one line:

target_position = SkyCoord(ra=83.63308, dec=22.01450, unit='deg')

This comment has been minimized.

@jjlk

jjlk Mar 16, 2017
Author Contributor

Done

spectral_model=model,
energy_range=energy_range)

# ToDo, add real test when results are checked

This comment has been minimized.

@cdeil

cdeil Mar 14, 2017
Member

Please put asserts on what results you get now:

  • Time binning used
  • Flux and flux error obtained for at least one bin, maybe even the first and last.

Leave the comment that numbers need to be checked as you work on the implementation.

But it's useful to have an example under test that asserts on current results.
There might be people (like me) that want to refactor code and things in Gammapy, without being really familiar with the algorithms. Just try to improve the code structure and implementation, but results should remain the same.
So for that or anyone working on the LC code in the future really, it's good if there is a test asserting on some results.

This comment has been minimized.

@jjlk

jjlk Mar 16, 2017
Author Contributor

Done.

@cdeil
cdeil approved these changes Mar 16, 2017
jlk added 2 commits Mar 16, 2017
jlk
@cdeil
Copy link
Member

@cdeil cdeil commented Mar 17, 2017

Thanks!

@cdeil cdeil merged commit 2359586 into gammapy:master Mar 17, 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
@cdeil
Copy link
Member

@cdeil cdeil commented Mar 17, 2017

I've split out the computation for one time bin into a separate method in 3cd7d76
and done some minor cleanup, like renaming "on" to "n_on" (trying to go towards more consistency within Gammapy, and following the names chosen here.

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

2 participants