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 make_images method in MapMaker #1681

Merged
merged 7 commits into from Aug 12, 2018

Conversation

2 participants
@registerrier
Copy link
Contributor

registerrier commented Aug 8, 2018

This is a first commit for the code to produce 2D images in MapMaker. It relies on a new function in cube/exposure.py that compute a weighted exposure cube in reco energy, using an input energy dispersion in the form of EnergyDispersion.
This should provide an exposure in m2 s so that dividing excesses by this map should yield integral fluxes in the considered energy bins.

Test are missing. Will add them soon.

Opinions @cdeil, @adonath, @leajouvin ?

@registerrier registerrier self-assigned this Aug 8, 2018

@registerrier registerrier requested a review from cdeil Aug 8, 2018

@registerrier registerrier added the feature label Aug 8, 2018

@registerrier registerrier added this to To do in Map analysis via automation Aug 8, 2018

@registerrier registerrier added this to the 0.8 milestone Aug 8, 2018

@cdeil
Copy link
Member

cdeil left a comment

@registerrier - Thanks!

The function does two things: weight with a spectrum and apply an EDISP. Is it required to couple that computation? Would it make sense to do it separately?

if name == 'exposure':
expo_reco = make_map_exposure_reco_energy(map, spectrum, edisp)
images[name] = expo_reco.sum_over_axes()
images[name] = map.sum_over_axes()

This comment has been minimized.

@cdeil

cdeil Aug 8, 2018

Member

Move this line into an else block?

I think otherwise, you're overwriting the exposure image with just the sum image like for counts and background?

img *= weights[idx].value

if edisp is None:
edisp = EnergyDispersion.from_diagonal_matrix(etrue_edges, etrue_edges)

This comment has been minimized.

@cdeil

cdeil Aug 8, 2018

Member

Isn't applying the diagonal matrix a no-op?
If so, you could just now make one and do EDISP stuff if no EDISP is given.
Faster and clearer.

registerrier added some commits Aug 8, 2018

@registerrier

This comment has been minimized.

Copy link
Contributor Author

registerrier commented Aug 8, 2018

Since a priori, counts and true energy exposure maps can have different energy binning I tend to consider that convolution with EnergyDispersion is the natural mechanism to go from one to the other.
As such weighting the true energy exposure with spectral weights only does not mean much, expect that you assume a diagonal energy dispersion and identical true and reco energy bins. I don't think it needs duplication, no?

@cdeil

This comment has been minimized.

Copy link
Member

cdeil commented Aug 8, 2018

Shouldn't we (optionally) make a combined exposure / edisp cube, that can be used by people in likelihood fitting if they want something fast, so not apply EDISP in the 3D likelihood fit, but a bit more precise than not using EDISP at all?

And then make_images would just be sum_over_axes also in that case?

Would be really good to have a very simple test case (either power-law or constant spectrum) with some non-diagonal edisp (e.g. map 4 true to 2 reco bins) for discussion. If it's a bit hard to develop as automated tests directly, what we did often was put a script in examples for a PR / code review.


# Here we do the convolution with the edisp matrix
a = np.rollaxis(expo_map.data, 0, 3)
b = np.dot(a, edisp.pdf_matrix)

This comment has been minimized.

@cdeil

cdeil Aug 8, 2018

Member

We shouldn't duplicate this np.rollaxis, np.dot, np.rollaxis stuff, no?
Already figured out and tested here:
http://docs.gammapy.org/dev/_modules/gammapy/cube/fit.html#MapEvaluator.apply_edisp

Of course, we could also call from MapEvaluator into a function, but it shouldn't be duplicated IMO.

This comment has been minimized.

@registerrier

registerrier Aug 9, 2018

Author Contributor

We will have other usages of an apply_edisp taking a Map and an EnergyDispersion as input; for instance, if we want a PSFKernel in reco_energy for morphological analysis.

Therefore, I propose to create a stand-alone function apply_edisp(map, edisp). It can be based on a private function _map_dot_product(array, array) that will just to the axis rotation and dot product on the data. The latter can be also used in MapEvaluator.apply_edisp().

Is this OK @cdeil, @adonath ?

This comment has been minimized.

@cdeil

cdeil Aug 9, 2018

Member

I don't have the big picture yet how this function and the MapEvaluator relate.

This is on the line between MapMaker and MapEvaluator, could go either way.

This could also become a method on MapEvaluator, because it's doing Spectrum / npred related computations. Note that you can have @staticmethod just located on MapEvalutor, to avoid 10 standalone functions.

If you want to pursue this as a standalone function, you could also just make MapEvaluator.apply_edisp a `@staticmethod and call it here from now. Then the code duplication with the rollaxis is avoided. Or move it to a separate function and call it from both places. But I'm not super excited about ending up with ~ 10 functions that each just have 3 lines, like the ones we now have on the MapEvaluator.

@registerrier - ?

@registerrier

This comment has been minimized.

Copy link
Contributor Author

registerrier commented Aug 11, 2018

I have removed the usage of edisp for this PR. I have renamed the function to weighted_exposure_image to make it clear this is not folding the energy dispersion.

I have added a test assuming constant exposure.

@cdeil cdeil assigned cdeil and unassigned registerrier Aug 12, 2018

@cdeil

This comment has been minimized.

Copy link
Member

cdeil commented Aug 12, 2018

@registerrier - Thanks!

In cd20d72 I did a few small changes: add tests for make_images; rename weighted_exposure_image to _map_spectrum_weight to make it private for now, and restructure it a bit to compute the weights vector first and apply it second, in case we want to split that in two functions later. I've also left the sub_over_axes to the caller, seems more flexible. Overall the exposure image that comes out of MapMaker.make_images is exactly the same as before, just the helper function was changed slightly.

Merging this now.

@cdeil cdeil merged commit a43d786 into gammapy:master Aug 12, 2018

0 of 2 checks passed

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 Aug 12, 2018

@cdeil cdeil changed the title Make image mapmaker Add make_images method in MapMaker Aug 15, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.