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

Implement MapEventSampler class #2238

Merged
merged 21 commits into from Jul 8, 2019
Merged

Conversation

@fabiopintore
Copy link
Contributor

@fabiopintore fabiopintore commented Jun 14, 2019

This PR introduces the MapEventSampler class which samples a set of events from a given probability distribution function. The class is fundamental in the development of a Gammapy simulator, as described in the pig-009 (docs/development/pigs/pig-009.rst).

The MapEventSampler requires an MapEvaluator object, a start and stop time, and also a light-curves to take into account temporal variability.
The main methods are: a calculator of the predicted counts from the given input map (MapEvaluator.npred_total()), an event sampler that provides the coordinates and energies (MapEvaluator.sample_npred), and the times of arrival (MapEvaluator.sample_timepred) of the events. These are then stored into an astropy Table object (MapEvaluator.sample_events).

The MapEventSampler class is included in inverse_cdf.py in gammapy.utils.random. Some tests have also been added.

I highlight that some more work is necessary to improve the sampling of the temporal variability.

@adonath adonath self-assigned this Jun 15, 2019
@adonath adonath added this to the 0.13 milestone Jun 15, 2019
@cdeil
Copy link
Member

@cdeil cdeil commented Jul 5, 2019

Could you please add a line __all__ = ["InverseCDFSampler"] at the top of the file, below the imports?
(suggestion from here)

cdeil
Copy link
Member

cdeil commented on 99ff10c Jul 5, 2019

@fabiopintore - Thanks!

For code formatting, you can run

pip install black
black gammapy/utils/random/inverse_cdf.py gammapy/utils/random/tests/test_inverse_cdf.py
git commit -m 'black code formatting'

Usually formatting should be always done before committing edits, but in cases where it was forgotten, the way to do it is just via a follow-up commit like that, that just contains code formatting changes, not other functional changes mixed into the same commit.

Copy link
Member

@adonath adonath left a comment

Thanks @fabiopintore! I've left some more comments, once those are addressed the PR should be ready to merge.

Parameters
----------
Copy link
Member

@adonath adonath Jul 5, 2019

Remove empty line...

Copy link
Contributor Author

@fabiopintore fabiopintore Jul 5, 2019

Done.

"""

def __init__(
self, npred_map, random_state=0, lc=None, phase_lc=None, tmin=0, tmax=3600
Copy link
Member

@adonath adonath Jul 5, 2019

Do you think its possible to join the lc and lc_phase arguments to a single temporal_model?

Copy link
Member

@adonath adonath Jul 5, 2019

I think t_min and t_max should be required arguments. In case no light- or phase curve is given, it will be sampled from a uniform distribution if the temporal model is present, the it is sampled from the model between t_min and t_max.

Copy link
Contributor Author

@fabiopintore fabiopintore Jul 5, 2019

Done.

"""

n_events = self.n_events
if (self.lc is not None) and (self.tmax > 0):
Copy link
Member

@adonath adonath Jul 5, 2019

The code for the lc and lc_phase is to a large degree the same, except that in one case you call _interpolator(t, ext=3) and .evaluate_norm_at_time() in the other. I think it's probably the better solution
Maybe just add the ext argument to LightCurveTableModel.evaluate_norm_at_time(). Then everything becomes simpler...

Copy link
Contributor Author

@fabiopintore fabiopintore Jul 5, 2019

I've modified the ~gammapy.time.model.LightCurveTableModel.evaluate_norm_at_time() including the ext argument. By default, this will be equal to 3. I've added also the doc string for this argument.

Input phase-curve model of the source, given with columns labelled
as "Phase" and "normalization" (arbitrary units): the bin time
HAS to be costant.
tmin : float
Copy link
Member

@adonath adonath Jul 5, 2019

I think we should make t_min and t_max an astropy.time.Time object. This allows users to attach units as well as handling different reference times.

Copy link
Contributor Author

@fabiopintore fabiopintore Jul 6, 2019

I've changed the code as that tmin and tmax may be either astropy.units.Quantity or astropy.time.TIME objects. In the second case, the times are then all referred to a reference time, hardcoded in the MapEventSampler class. These are the changes:

from astropy.time import Time
tmin = ...
tmax = ...
if isinstance(tmin,Time) and isinstance(tmax,Time):
t_ref = Time('2019-01-01 00:00:00', scale='utc') #hardcoded reference epoch
self.tmin = ((tmin.decimalyear - t_ref.decimalyear) * u.year).to(u.s)
self.tmax = ((tmax.decimalyear - t_ref.decimalyear) * u.year).to(u.s)

else:
self.tmin = tmin.to(u.s)
self.tmax = tmax.to(u.s)

Returns
-------
coords : `~numpy.ndarray`
Copy link
Member

@adonath adonath Jul 5, 2019

I think this is not true, coords should be a MapCoord object...

Copy link
Contributor Author

@fabiopintore fabiopintore Jul 5, 2019

Done.

self.tmin = tmin
self.tmax = tmax

def npred_total(self):
Copy link
Member

@adonath adonath Jul 5, 2019

I think we should make npred_total a lazyproperty. This will compute the total number of events just once and the store it. See http://docs.astropy.org/en/stable/api/astropy.utils.decorators.lazyproperty.html.

Copy link
Contributor Author

@fabiopintore fabiopintore Jul 5, 2019

Done.

Array with coordinates and energies of the sampled events.
"""

self.n_events = self.npred_total()
Copy link
Member

@adonath adonath Jul 5, 2019

Once npred_total is a lazyproperty you can just access self.npred_total and don't have to set self.n_events.

Copy link
Contributor Author

@fabiopintore fabiopintore Jul 5, 2019

Done.

n_events = self.n_events
if (self.lc is not None) and (self.tmax > 0):
n_tbin = self.tmax - self.tmin
t = np.linspace(self.tmin, self.tmax, n_tbin)
Copy link
Member

@adonath adonath Jul 5, 2019

I think we should add t_delta argument, that user can choose the bin size for the time array.

Copy link
Contributor Author

@fabiopintore fabiopintore Jul 5, 2019

Done.

events["lon_true"] = coords[0] * u.deg
events["lat_true"] = coords[1] * u.deg
events["e_true"] = coords[2] * u.TeV
if ((self.lc is not None) and (self.tmax > 0)) or (self.tmax > 0):
Copy link
Member

@adonath adonath Jul 5, 2019

Once t_min and t_max are required arguments this becomes non-optional...

Copy link
Contributor Author

@fabiopintore fabiopintore Jul 5, 2019

Done.

@adonath adonath merged commit 369b335 into gammapy:master Jul 8, 2019
0 of 9 checks passed
@adonath adonath changed the title Introduce the MapEventSampler class Implement MapEventSampler class Jul 25, 2019
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