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

Allow different true and reco energy in map analysis #1846

Merged
merged 10 commits into from Oct 4, 2018

Conversation

3 participants
@AtreyeeS
Contributor

AtreyeeS commented Oct 3, 2018

This PR lets maps have different geometries for the true and reco energy.
The exposure and PSF map is in true energies, the counts and background map in reco energies.
Changes are made in map creation (MapMaker) and model evaluation and fitting.

@cdeil cdeil added the feature label Oct 3, 2018

@cdeil cdeil added this to the 0.9 milestone Oct 3, 2018

@cdeil cdeil added this to To do in Map analysis via automation Oct 3, 2018

exclusion_mask=exclusion_mask,
).run(selection)
# Stack observation maps to total
for name in selection:
data = maps_obs[name].quantity.to(self.maps[name].unit).value
self.maps[name].fill_by_coord(coords, data)
self.maps[name].fill_by_coord(coords_etrue, data) if name == "exposure" else self.maps[name].fill_by_coord(coords, data)

This comment has been minimized.

@cdeil

cdeil Oct 3, 2018

Member

Break the line, put an if, else

def compute_npred(self):
"""Evaluate model predicted counts.
"""
Evaluate model predicted counts.

This comment has been minimized.

@cdeil

cdeil Oct 3, 2018

Member

Always put an empty line before Returns section.
To check locally

python setup.py build_docs
make docs-show
Evaluate model predicted counts.
Returns
-------
npred.data : numpy array

This comment has been minimized.

@cdeil

cdeil Oct 3, 2018

Member

Write type like this: ~numpy.ndarray

data = np.dot(data, self.edisp.pdf_matrix)
return np.rollaxis(data, 2, 0)
data1 = np.moveaxis(data, -1, loc) #now dim is in ereco
npred1 = Map.from_geom(self.background.geom, unit="")

This comment has been minimized.

@cdeil

cdeil Oct 3, 2018

Member

Suggest to leave a comment:

# TODO: maybe pass counts geom to evaluator separately and use that here

This comment has been minimized.

@cdeil

cdeil Oct 3, 2018

Member

One idea here could be to make the energy axis from the EDISP.

if self.edisp is not None:
npred.data = self.apply_edisp(npred.data)
npred = self.apply_edisp(npred)
if self.background:
npred.data += self.background.data
return npred.data

This comment has been minimized.

@cdeil

cdeil Oct 3, 2018

Member

Suggest to return npred map here and adapt the MapFit

This comment has been minimized.

@AtreyeeS

AtreyeeS Oct 4, 2018

Contributor

This is breaking the simulation notebooks. The MapFit only needs a ~numpy.ndarray to compute the statistics, so not changing it now.

assert_allclose(out.data.sum(), 4.004864e+12, rtol=1e-5)
assert_allclose(out.data[0, 0, 0], 1.380614e+09, rtol=1e-5)
assert out.data.shape == (3, 4, 5)
assert_allclose(out.data.sum(), 1.106404e+12, rtol=1e-5)

This comment has been minimized.

@cdeil

cdeil Oct 3, 2018

Member

Why did this value change from 4e12 to 1.1e12? Is this expected?
Can you please put a TODO comment to improve the test cases in this file to have something where we know what to expect (but do that test case change in a follow-up PR).

@@ -1421,3 +1421,26 @@ def copy(self, **kwargs):
Copied map geometry.
"""
return self._init_copy(**kwargs)
def get_axes_dimension_by_name(self, name):

This comment has been minimized.

@cdeil

cdeil Oct 3, 2018

Member

Suggest to call this get_axis_index_by_name and put it next to get_axis_by_name.

For the implementation, I would suggest to do it very similar to the existing get_axis_by_name, i.e. put this:

axes = {axis.name.upper(): idx for idx, axis in enumerate(self.axes)}
return axes[name.upper()]

Add a unit test case that calls it once with a name that exists to show that it works.
And call it once with a test case for a name that doesn't exist, and check that an appropriate error is raised via

with pytest.raises(IndexError):
    ...

This comment has been minimized.

@AtreyeeS

AtreyeeS Oct 4, 2018

Contributor

I have put the test in test_wcs, with a raises(ValueError)

@cdeil

A few more small comments inline.

Otherwise, 👍 to merge.

def apply_edisp(self, npred):
"""Convolve map data with energy dispersion.
npred: npred map in e_true

This comment has been minimized.

@cdeil

cdeil Oct 3, 2018

Member

Fix docstring formatting

"""
loc = npred.geom.get_axes_dimension_by_name("energy")
data = np.moveaxis(npred.data, loc, -1)

This comment has been minimized.

@cdeil

cdeil Oct 3, 2018

Member

Use numpy rollaxis, we still support old Numpy versions.

data = np.dot(data, self.edisp.pdf_matrix)
return np.rollaxis(data, 2, 0)
data1 = np.moveaxis(data, -1, loc) #now dim is in ereco

This comment has been minimized.

@cdeil

cdeil Oct 3, 2018

Member

You are modifying data and still keep calling the variable data on the lines above.
Suggest to continue here and not start calling it data1.

if self.edisp is not None:
npred.data = self.apply_edisp(npred.data)
npred = self.apply_edisp(npred)

This comment has been minimized.

@cdeil

cdeil Oct 3, 2018

Member

I think that's a good change, to always return Map objects now.
If this is a significant inefficiency we can change again later, but probably it's nicer / cleaner like this.

@cdeil cdeil assigned cdeil and unassigned registerrier Oct 3, 2018

"""
axes = {axis.name.upper(): idx for idx, axis in enumerate(self.axes)}
return axes[name.upper()]

This comment has been minimized.

@registerrier

registerrier Oct 4, 2018

Contributor

you should return an index not the axis here, no?

names = [axis.name for axis in self.axes]
return names.index(name)

AtreyeeS and others added some commits Oct 4, 2018

@cdeil

cdeil approved these changes Oct 4, 2018

Looks great. Thank you!

@cdeil cdeil merged commit 9ca3dde into gammapy:master Oct 4, 2018

0 of 3 checks passed

Scrutinizer Running
Details
continuous-integration/appveyor/pr Waiting for AppVeyor build to complete
Details
continuous-integration/travis-ci/pr The Travis CI build is in progress
Details

Map analysis automation moved this from To do to Done Oct 4, 2018

@cdeil cdeil changed the title from Allow different energy binning for true and reconstructed energy in Maps to Allow different true and reco energy in man analysis Nov 23, 2018

@cdeil cdeil changed the title from Allow different true and reco energy in man analysis to Allow different true and reco energy in map analysis Nov 23, 2018

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