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 SkyDiffuseCube model for 3D maps #1634

Merged
merged 4 commits into from Aug 6, 2018

Conversation

Projects
4 participants
@robertazanin
Contributor

robertazanin commented Aug 3, 2018

Add a class for SkyMap3d.
The test is not working because the energy interpolation is not properly working

@robertazanin robertazanin added the feature label Aug 3, 2018

@robertazanin robertazanin added this to the 0.8 milestone Aug 3, 2018

@robertazanin robertazanin added this to In progress in Modeling via automation Aug 3, 2018

@registerrier

Thanks. See some suggestions inline.

Show outdated Hide outdated gammapy/cube/models.py
Show outdated Hide outdated gammapy/cube/models.py
Show outdated Hide outdated gammapy/cube/models.py
Show outdated Hide outdated gammapy/cube/models.py
Show outdated Hide outdated gammapy/cube/models.py
@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Aug 4, 2018

Member

@robertazanin - I figures out why the map inter_by_coord wasn't giving what we wanted yesterday.

First of all, for diffuse model cubes, we need to create an energy MapAxis with node_type='center', either by passing that as an option to MapAxis or by calling Map.from_nodes. Then the map values will represent the diffuse model fluxes at the given energies.

Secondly, by default interp_by_coord on Map does nearest neighbor interpolation. Here we want to pass interp='linear' so that interpolation in energy ("linear" in log-log, i.e. actually piecewise power-law interpolation) takes place.

Here's an example one can use to understand how it works. A cube with 2 pixels in energy, with value 10 at energy 1 TeV and value 20 at energy 100 TeV:

import numpy as np
from gammapy.maps import Map, MapAxis
axis = MapAxis([1, 100], node_type='center', name="energy", unit='TeV', interp='log')
m = Map.create(npix=(4, 3), binsz=2, axes=[axis])
m.data += [[[10]], [[20]]]
coord = {'lon': 0, 'lat': 0, 'energy': 10}
pix = m.geom.coord_to_pix(coord)
val = m.interp_by_coord(coord, fill_value=0, interp='linear')

gives

>>> print(pix, val)
(array([1.5]), array([1.]), array([0.5])) [15.]

So the action items here are:

  • Change the test case to such an energy axis, using the example given here.
  • Add an evaluate for a sky position outside the map range
  • Add an evaluate for an energy above or below the first or last node.
  • Change both the 2D and the 3D class to use interp='linear'

In __init__, we could check if m.geom.axis[0].node_type is "center", and if it's "edge" error out with a good message.

A second test case should be added where you actually read a diffuse model cube file.
I tried this, and it does come with node_type="center" as it should, so you can use that one:

from gammapy.maps import Map
m = Map.read('test_datasets/unbundled/fermi/gll_iem_v02_cutout.fits')
axis = m.geom.axes[0]
print(axis)

@robertazanin - I'll leave this PR to you to finish up, or just come by on Monday and we do it together.

Member

cdeil commented Aug 4, 2018

@robertazanin - I figures out why the map inter_by_coord wasn't giving what we wanted yesterday.

First of all, for diffuse model cubes, we need to create an energy MapAxis with node_type='center', either by passing that as an option to MapAxis or by calling Map.from_nodes. Then the map values will represent the diffuse model fluxes at the given energies.

Secondly, by default interp_by_coord on Map does nearest neighbor interpolation. Here we want to pass interp='linear' so that interpolation in energy ("linear" in log-log, i.e. actually piecewise power-law interpolation) takes place.

Here's an example one can use to understand how it works. A cube with 2 pixels in energy, with value 10 at energy 1 TeV and value 20 at energy 100 TeV:

import numpy as np
from gammapy.maps import Map, MapAxis
axis = MapAxis([1, 100], node_type='center', name="energy", unit='TeV', interp='log')
m = Map.create(npix=(4, 3), binsz=2, axes=[axis])
m.data += [[[10]], [[20]]]
coord = {'lon': 0, 'lat': 0, 'energy': 10}
pix = m.geom.coord_to_pix(coord)
val = m.interp_by_coord(coord, fill_value=0, interp='linear')

gives

>>> print(pix, val)
(array([1.5]), array([1.]), array([0.5])) [15.]

So the action items here are:

  • Change the test case to such an energy axis, using the example given here.
  • Add an evaluate for a sky position outside the map range
  • Add an evaluate for an energy above or below the first or last node.
  • Change both the 2D and the 3D class to use interp='linear'

In __init__, we could check if m.geom.axis[0].node_type is "center", and if it's "edge" error out with a good message.

A second test case should be added where you actually read a diffuse model cube file.
I tried this, and it does come with node_type="center" as it should, so you can use that one:

from gammapy.maps import Map
m = Map.read('test_datasets/unbundled/fermi/gll_iem_v02_cutout.fits')
axis = m.geom.axes[0]
print(axis)

@robertazanin - I'll leave this PR to you to finish up, or just come by on Monday and we do it together.

@adonath

This comment has been minimized.

Show comment
Hide comment
@adonath

adonath Aug 6, 2018

Member

Just a few general comments maybe worth thinking about:

Why is it needed to introduce a new class for this? Can't the existing SkyDiffuseMap easily be modified to work in arbitrary dimensions? This is what I would expect it to do from it's name...

For efficiency I would also strongly suggest to not interpolate in the .evaluate() method, but rather pass an geom object on init, which does a reprojection of the diffuse model to the data geometry and just returns the diffuse map data on evaluate. The repeated interpolation is not needed as map geometries will never change during a fit...

Member

adonath commented Aug 6, 2018

Just a few general comments maybe worth thinking about:

Why is it needed to introduce a new class for this? Can't the existing SkyDiffuseMap easily be modified to work in arbitrary dimensions? This is what I would expect it to do from it's name...

For efficiency I would also strongly suggest to not interpolate in the .evaluate() method, but rather pass an geom object on init, which does a reprojection of the diffuse model to the data geometry and just returns the diffuse map data on evaluate. The repeated interpolation is not needed as map geometries will never change during a fit...

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Aug 6, 2018

Member

@adonath - What you suggest is not so easy. Note that the goal is to finish this today, I plan to work on this next. Basically I would suggest to postpone your suggestions to v0.9, unless you have a suggestion how to address these difficulties:

Why is it needed to introduce a new class for this? Can't the existing SkyDiffuseMap easily be modified to work in arbitrary dimensions? This is what I would expect it to do from it's name...

It probably could, but the evaluate argument handling might get complex. I think arbitrary extra axes in the model evaluation scheme will need much more thought and discussion and work.

For efficiency I would also strongly suggest to not interpolate in the .evaluate() method, but rather pass an geom object on init, which does a reprojection of the diffuse model to the data geometry and just returns the diffuse map data on evaluate. The repeated interpolation is not needed as map geometries will never change during a fit...

Is fixing the geom on which one can evaluate on model init a good idea? It would mean special-casing model init for this one model to already pass a geom. So far our logic is that model evaluate is just on coords, not passed a geom.

How about we keep this model as-is, and then special-case it in the MapEvaluator. In my mind combining model and geom in a clever way it that class's responsibility. Although I have to admit I'm not sure it's clean and easy to do the caching of the interpolator or even pre-evaluated cube on a given geom there, i.e. to store that state on MapEvaluator.

Member

cdeil commented Aug 6, 2018

@adonath - What you suggest is not so easy. Note that the goal is to finish this today, I plan to work on this next. Basically I would suggest to postpone your suggestions to v0.9, unless you have a suggestion how to address these difficulties:

Why is it needed to introduce a new class for this? Can't the existing SkyDiffuseMap easily be modified to work in arbitrary dimensions? This is what I would expect it to do from it's name...

It probably could, but the evaluate argument handling might get complex. I think arbitrary extra axes in the model evaluation scheme will need much more thought and discussion and work.

For efficiency I would also strongly suggest to not interpolate in the .evaluate() method, but rather pass an geom object on init, which does a reprojection of the diffuse model to the data geometry and just returns the diffuse map data on evaluate. The repeated interpolation is not needed as map geometries will never change during a fit...

Is fixing the geom on which one can evaluate on model init a good idea? It would mean special-casing model init for this one model to already pass a geom. So far our logic is that model evaluate is just on coords, not passed a geom.

How about we keep this model as-is, and then special-case it in the MapEvaluator. In my mind combining model and geom in a clever way it that class's responsibility. Although I have to admit I'm not sure it's clean and easy to do the caching of the interpolator or even pre-evaluated cube on a given geom there, i.e. to store that state on MapEvaluator.

@adonath

This comment has been minimized.

Show comment
Hide comment
@adonath

adonath Aug 6, 2018

Member

@cdeil With taking a Map object the model init for the diffuse model is already special cased. In fact one could even ask the user to do the reprojection beforehand and pass in a map with the right geometry already.

The two suggestion I made above go along with each other and would greatly simplify the implementation. For now I would simplify the SkyDiffuseMap to just take a map with the same geometry as the data and have only one global norm parameter...interpolation (if needed at all) could be added later as an option.

Member

adonath commented Aug 6, 2018

@cdeil With taking a Map object the model init for the diffuse model is already special cased. In fact one could even ask the user to do the reprojection beforehand and pass in a map with the right geometry already.

The two suggestion I made above go along with each other and would greatly simplify the implementation. For now I would simplify the SkyDiffuseMap to just take a map with the same geometry as the data and have only one global norm parameter...interpolation (if needed at all) could be added later as an option.

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Aug 6, 2018

Member

So if we leave it up to the user or init to interpolate on a given geom, the evaluate becomes this, no?

def evaluate(self, *args, **kwargs):
    # Ignore any `args` or `kwargs`
    return self.parameters['norm'] * self._map.quantity

I'm still wondering if that is a good idea. I was assuming such model eval caching for a given geom is done (for many models) by the MapEvaluator and MapFit, not the model class itself.

I'll come by later to discuss ...

Member

cdeil commented Aug 6, 2018

So if we leave it up to the user or init to interpolate on a given geom, the evaluate becomes this, no?

def evaluate(self, *args, **kwargs):
    # Ignore any `args` or `kwargs`
    return self.parameters['norm'] * self._map.quantity

I'm still wondering if that is a good idea. I was assuming such model eval caching for a given geom is done (for many models) by the MapEvaluator and MapFit, not the model class itself.

I'll come by later to discuss ...

@adonath

This comment has been minimized.

Show comment
Hide comment
@adonath

adonath Aug 6, 2018

Member

I discussed with @cdel offline and I agree it makes sense to introduce the model evaluation caching in a general way in the MapEvaluator or MapFit classes. So no further discussion from my side needed in this PR...

Member

adonath commented Aug 6, 2018

I discussed with @cdel offline and I agree it makes sense to introduce the model evaluation caching in a general way in the MapEvaluator or MapFit classes. So no further discussion from my side needed in this PR...

@cdeil cdeil referenced this pull request Aug 6, 2018

Closed

Add SkyDiffuseCube model #1629

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Aug 6, 2018

Member

@registerrier or @adonath - Please review.

Member

cdeil commented Aug 6, 2018

@registerrier or @adonath - Please review.

@cdeil cdeil changed the title from add SkyMap3d to Add SkyDiffuseCube model for 3D maps Aug 6, 2018

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Aug 6, 2018

Member

I've fixed the coord handling in 824ef28 to match what the MapEvaluator passes.

This model now works with MapEvaluator. For the maps from analysis_3d.ipynb I find that the Galactic diffuse model produces 3567 counts. I'm not sure yet if that is correct, using

for obs in obs_list:
    events = obs.events.table
    events = events[events['ENERGY'] > 0.1]
    print(obs.obs_id, np.sum(events['MC_ID'] == 2))

I find

110380 4339
111140 4019
111159 4442

i.e. 4000 IEM events per run, so for 3 runs there should be 12,000 IEM events, a factor 3 more.

So maybe something is off, but I'd be inclined to merge now, to allow everyone using this and continue debugging in master.

screen shot 2018-08-06 at 18 34 17

Member

cdeil commented Aug 6, 2018

I've fixed the coord handling in 824ef28 to match what the MapEvaluator passes.

This model now works with MapEvaluator. For the maps from analysis_3d.ipynb I find that the Galactic diffuse model produces 3567 counts. I'm not sure yet if that is correct, using

for obs in obs_list:
    events = obs.events.table
    events = events[events['ENERGY'] > 0.1]
    print(obs.obs_id, np.sum(events['MC_ID'] == 2))

I find

110380 4339
111140 4019
111159 4442

i.e. 4000 IEM events per run, so for 3 runs there should be 12,000 IEM events, a factor 3 more.

So maybe something is off, but I'd be inclined to merge now, to allow everyone using this and continue debugging in master.

screen shot 2018-08-06 at 18 34 17

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Aug 6, 2018

Member

I'm merging this in now.

@robertazanin and all - Please try it out and report any issues you find, or open a new issue or PR with suggestions to make this class better.

@adonath - You were right that it's slow. Especially for the 1DC diffuse model cube, which is 800 MB, an evaluate takes almost a second, and I presume that most of that time is spent in interpolator init, i.e. could be avoided by caching that. Making a cutout of the diffuse model is one way to improve speed for now.

Member

cdeil commented Aug 6, 2018

I'm merging this in now.

@robertazanin and all - Please try it out and report any issues you find, or open a new issue or PR with suggestions to make this class better.

@adonath - You were right that it's slow. Especially for the 1DC diffuse model cube, which is 800 MB, an evaluate takes almost a second, and I presume that most of that time is spent in interpolator init, i.e. could be avoided by caching that. Making a cutout of the diffuse model is one way to improve speed for now.

@cdeil cdeil merged commit df3330a into gammapy:master Aug 6, 2018

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

Modeling automation moved this from In progress to Done Aug 6, 2018

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