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

Simplify and extend background model handling #2334

Merged
merged 8 commits into from Sep 4, 2019

Conversation

@adonath
Copy link
Member

commented Sep 2, 2019

This PR improve the background model handling on the tutorial and implements related changes:

  • Add the background models as SkyModel and not separate BackgroundModel components in the Fermi-LAT tutorial
  • Add tilt and reference parameter to SkyDiffuseCube
  • Add tilt and reference parameter to TableModel
  • Add dtype argument for Map.from_geom(), which was previously missing
  • Set SkyDiffuseConstant.evaluation_radius to None, which switches of the evaluation radius handling
@adonath adonath self-assigned this Sep 2, 2019
@adonath adonath added this to the 0.14 milestone Sep 2, 2019
@adonath adonath force-pushed the adonath:remove_fermi_background_model branch from 1e6df8b to 3028b6d Sep 2, 2019
Copy link
Member

left a comment

@adonath - Two inline comments.

For the PR title, suggest to drop the "Fermi-LAT tutorial" part.
This PR is really changing a few different things in Gammapy, the title shouldn't suggest that it's just a docs tutorial change.

@@ -226,7 +226,7 @@ def from_geom(geom, meta=None, data=None, map_type="auto", unit=""):
raise ValueError("Unrecognized geom type.")

cls_out = Map._get_map_cls(map_type)
return cls_out(geom, data=data, meta=meta, unit=unit)
return cls_out(geom, data=data, meta=meta, unit=unit, dtype=dtype)

This comment has been minimized.

Copy link
@cdeil

cdeil Sep 3, 2019

Member

This calls Map.__init__? That doesn't seem to have a separate dtype argument. Do we really need the option to do a data dtype change here (instead of calling astype explicitly)? Does this work? Add test?

This comment has been minimized.

Copy link
@adonath

adonath Sep 3, 2019

Author Member

Both WcsNDMap.init and HpxNDMap.init have a dtype argument and I guess it's needed, because you want to create the default data given a geometry. I can add a separate test for this...

This comment has been minimized.

Copy link
@cdeil

cdeil Sep 3, 2019

Member

Ah, OK. I only looked at Map.init which doesn't have this.

I wonder if we can simplify the whole from_geom and copy logic for maps. It's really complex, calling up and down the class hierarchy.
One thought: would changing to Map.from_geom to geom.make_map help? This hard-coded if / else on two pre-defined types means it's non-extensible to other sub-classes, and I think generally is a bad pattern with class inheritance:

def from_geom(geom, meta=None, data=None, map_type="auto", unit=""):
"""Generate an empty map from a `MapGeom` instance.
Parameters
----------
geom : `MapGeom`
Map geometry.
data : `numpy.ndarray`
data array
meta : `~collections.OrderedDict`
Dictionary to store meta data.
map_type : {'wcs', 'wcs-sparse', 'hpx', 'hpx-sparse', 'auto'}
Map type. Selects the class that will be used to
instantiate the map. The map type should be consistent
with the geometry. If map_type is 'auto' then an
appropriate map type will be inferred from type of ``geom``.
unit : str or `~astropy.units.Unit`
Data unit.
Returns
-------
map_out : `Map`
Map object
"""
if map_type == "auto":
from .hpx import HpxGeom
from .wcs import WcsGeom
if isinstance(geom, HpxGeom):
map_type = "hpx"
elif isinstance(geom, WcsGeom):
map_type = "wcs"
else:
raise ValueError("Unrecognized geom type.")
cls_out = Map._get_map_cls(map_type)
return cls_out(geom, data=data, meta=meta, unit=unit)

This comment has been minimized.

Copy link
@adonath

adonath Sep 3, 2019

Author Member

I think API-wise the pattern is not too bad. It follows the logic of creating the map geometry
first and then creating a map by combining it with the data or crating default data. I guess Map.from_geom() was introduced to have a factory function, that is agnostic of the map type and I think you definitely need something like this because in the code you want to create maps, without knowing the map-type in advance. The problem of the hard-coded map sub-types you also fave for Map.create().

My proposal would be to stick with the .from_geom(), but I agree we should change to a registry pattern for sub-map types and get rid of the hard-coded classes mappings.

This comment has been minimized.

Copy link
@cdeil

cdeil Sep 3, 2019

Member

It doesn't solve all problems, but moving Map.from_geom to geom.make_make does solve this inheritance issue, allowing us to remove this if / else / ininstance check, because the geom is the concrete subclass and knows it's type. But OK, it's not clear if it's worth making the change. Let's drop this discussion for now, and then, if we have time in a few weeks, consider improving gammapy.maps.

def eval_radius_gauss(sigma, **kwargs):
return 8 * sigma
gauss.evaluation_radius = eval_radius_gauss

This comment has been minimized.

Copy link
@cdeil

cdeil Sep 3, 2019

Member

Why did you change this to staticmethod _evaluation_radius and dispatch from evaluation_radius? Why not keep as-is. Seems simpler and all use cases you have here can still be done.

Is it a good idea to teach people to monkey-patch evaluation_radius here in the docstring?
I mean, yes, experts can use monkey-patching throughout Gammapy to customise the analysis, but genereally it's not something we advertise in the docs, and I don't see why this evaluation radius should be something they need to monkey-patch.

Having a docs page that explains how to make custom models, usually from scratch or by sub-classing, very rarely by monkey-patching the existing ones, could be a nice addition.

This comment has been minimized.

Copy link
@adonath

adonath Sep 3, 2019

Author Member

Yes, maybe we can find a better solution for this. My reasoning was as follows:

  • The model evaluation currently uses bounding boxes and this we definitely want to keep
  • The SkyDiffuseConstant is not spatially confined, so I needed a way to switch off the bounding box evaluation for individual model components. I decided to allow for evaluation_radius=None to control, whether the bounding boxes are used or not. At this point we could stop and just implement fixed default values for all models (that's how it was, evaluation_radius was a property and not settable)
  • But for some optimisers (e.g. sherpa's MC or similar) and models (e.g. SkyDiffuseConstant), it is useful to switch of the bounding box evaluation and I decided it would be a useful addition for to be able to define / switch-off the evaluation radius per model component
  • I implemented a setter for evaluation_radius. But this is where the problem starts, because typically the evaluation radius is a function of some model parameters. And that's where the monkey patching is not avoidable, except if you pass a function that determines the evaluation radius on model init

I guess there are a few options:

  • Make the evaluation_radius handling not configurable by the user
  • Only allow for fixed evaluation radii, that do not depend on model parameters
  • Invent a mini-language, so that users can pass evaluation_radius="5 * {sigma}"?
  • Add an evaluation_radius argument for model init?
  • Anything else...?

This comment has been minimized.

Copy link
@cdeil

cdeil Sep 3, 2019

Member

I think having None as declarative value for the case of no evaluation radius is useful.

With the callable, you could just have kept everything the same, and very advanced users could have monkey-patched just by attaching a function with self as first argument. That's the common way to monkey-patch methods, no?

This comment has been minimized.

Copy link
@adonath

adonath Sep 3, 2019

Author Member

Yes, I was struggling with the fact, that evaluation_radius was defined as a property and properties of instances can't be easily monkey-patched, I guess. But it's probably OK for power-users to do it on a class and not instance level?

This comment has been minimized.

Copy link
@adonath

adonath Sep 3, 2019

Author Member

Like so:

from gammapy.image.models import SkyGaussian

def eval_radius(self):
    return 8 * self.sigma.quantity

SkyGaussian.evaluation_radius = eval_radius

# or to switch off
SkyGaussian.evaluation_radius = None

This comment has been minimized.

Copy link
@cdeil

cdeil Sep 3, 2019

Member

I don't recall if I ever monkey-patched a property at the instance level.

We have 1000s of properties in Gammapy ... is there a real use case to monkey-patch this one?

This comment has been minimized.

Copy link
@adonath

adonath Sep 3, 2019

Author Member

To be honest, I don't want to spend much more time on this. I'll go back to the previous solution to have a global evaluation_mode argument on the MapDataset and remove the configuration per model component. The class-level monkey patching I've shown above is still possible...

@adonath adonath changed the title Simplify background model handling in Fermi-LAT tutorial Simplify and extend background model handling Sep 3, 2019
@adonath adonath force-pushed the adonath:remove_fermi_background_model branch from 3028b6d to eb78852 Sep 3, 2019
@adonath adonath force-pushed the adonath:remove_fermi_background_model branch from 82af9a5 to 94cd075 Sep 4, 2019
@adonath adonath force-pushed the adonath:remove_fermi_background_model branch from 94cd075 to c0d0ece Sep 4, 2019
@adonath

This comment has been minimized.

Copy link
Member Author

commented Sep 4, 2019

@QRemy Please note in this PR I had to slightly modify the serialisation tests: the model component gc09 is outside the image range of the gc dataset, which currently causes the model evaluation to raise an error. This behaviour is currently intended, because a model component outside of the data range (and "outside" means not even overlapping after PSF convolution) does not contribute to the likelihood of the dataset. So in gammapy/gammapy-extra@c5da4a1 I changed the datasets YAML file to not include the gc09 component in the gc dataset and vice versa.

@adonath adonath merged commit 46f541f into gammapy:master Sep 4, 2019
9 checks passed
9 checks passed
Codacy/PR Quality Review Up to standards. A positive pull request.
Details
Scrutinizer Analysis: No new issues – Tests: passed
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
gammapy.gammapy Build #20190904.5 succeeded
Details
gammapy.gammapy (DevDocs) DevDocs succeeded
Details
gammapy.gammapy (Lint) Lint succeeded
Details
gammapy.gammapy (Test Python36) Test Python36 succeeded
Details
gammapy.gammapy (Test Windows36) Test Windows36 succeeded
Details
gammapy.gammapy (Test Windows37) Test Windows37 succeeded
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
2 participants
You can’t perform that action at this time.