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 more SkySpatialModel subclasses #1377

Merged
merged 5 commits into from Apr 16, 2018

Conversation

Projects
None yet
3 participants
@joleroi
Contributor

joleroi commented Apr 16, 2018

  • Add SkyPointSource
  • Add SkyDisk2D
  • Add SkyShell2D
  • Add SkyTemplate2D
  • Make Spatial model base class abstract

Future PRs will

  • Remove existing model code
  • Implement some kind of model integration

joleroi added some commits Apr 16, 2018

@joleroi joleroi added the feature label Apr 16, 2018

@joleroi joleroi added this to the 0.8 milestone Apr 16, 2018

@joleroi joleroi self-assigned this Apr 16, 2018

"""Evaluate the model (static function)."""
tolerance = 1 * u.arcsec
sep = angular_separation(lon, lat, lon_0, lat_0)
val = 1 if sep < tolerance else 0

This comment has been minimized.

@adonath

adonath Apr 16, 2018

Member

This implementation is not numerically stable, even with the tolerance parameter. Currently it requires a pixel size < 2 arcsec, that it evaluates to one at all. There is also the problem, that the point source region is circular, while the coordinate grid is rectangular. I think the minimal required change here would be to make the tolerance adjustable so it can be adapted to the grid size or change the implementation to just set the nearest pixel to one.

For the record: The are two ways to get this right with sub-pixel accuracy. Either to integrate over the pixels (this I have done here, but it uses a rectangular pixel approximation and only works for arrays of coordinates) or to directly evaluate the PSF.

This comment has been minimized.

@joleroi

joleroi Apr 16, 2018

Contributor

I discussed this with @cdeil and he suggested to leave the proper treatment (whatever that means) of the delta function to the function that evaluated the spatial models on the skycube. I guess the current implementation of SpatialModel is meant to be evaluated at a given coordinate, so everything involving pixel intergration schemes is out of scope. Do you agree?

This comment has been minimized.

@cdeil

cdeil Apr 16, 2018

Member

Yes, the point source needs to be special-cases in the model evaluator. It requires knowledge of the pixels (WCS or HPX), and probably it's not possible or at least not a good idea to attempt to do it here, i.e. couple this in the model code.

So this spatial model is mostly there to do book-keeping of the parameters.
I'm not sure yet what this measn for the evaluate method.
Leave as-is? Remove possible? Always return 1?
Maybe postpone this discussion until tomorrow, when we have a first draft map model evaluator and see more clearly how it could work for this case?

This comment has been minimized.

@adonath

adonath Apr 16, 2018

Member

When it is special cased in the model evaluator anyway, I would either do the 0th order approximation and just set the nearest pixel center to one or raise a NotImplementedError.

This comment has been minimized.

@joleroi

joleroi Apr 16, 2018

Contributor

I guess the 'nearest pixel center' is not really defined if you call evaluate with one coordinate, is it?
I'm fine with raising NotImplementedError or always returning 1, or with leaving it as is (which is what gammalib does)

This comment has been minimized.

@cdeil

cdeil Apr 16, 2018

Member

Suggest to just leave as-is, and then continue the discussion tomorrow when we handle this in the model evaluator.

joleroi added some commits Apr 16, 2018

@joleroi joleroi changed the title from Add SkyPointSource and SkyDisk2D spatial models to Add more SkySpatialModel subclasses Apr 16, 2018

return norm * val
class SkyTemplate2D(SkySpatialModel):

This comment has been minimized.

@joleroi

joleroi Apr 16, 2018

Contributor

Note that this is a copy&paste (+ adaptions) from the old model classes, I'm nore sure if we want to keep it in this form

@joleroi

This comment has been minimized.

Contributor

joleroi commented Apr 16, 2018

Ready from my side

term2 = term1 - np.sqrt(r_i**2-sep**2)
norm = 3 / (2 * np.pi * (r_o ** 3 - r_i ** 3))
if sep < r_o:

This comment has been minimized.

@cdeil

cdeil Apr 16, 2018

Member

Usually spatial models will be evaluated on 2D coordinate grids, i.e. lon and lat are two-dimensional arrays.
Then you have to use numpy to evaluate the different cases, using if / else will error out.
Please change (probably to the implementation in the old classes), and add a test evaluating the models on input arrays.
What I sometimes do is to just add a test for the array evaluate and put in arrays that are constant and the same as for the scalar case, and only check that the shape of the output array and a single pixel is OK.

@cdeil

@joleroi - Thanks!

Do you have time to continue a bit here?
Making evaluate work with arrays would be the first thing.
changing the template model from SkyImage over to Map would be another task.
Both could also be done in follow-up PRs if you're done for today.

@cdeil cdeil assigned cdeil and unassigned joleroi Apr 16, 2018

@cdeil

cdeil approved these changes Apr 16, 2018

There was one fail:
https://travis-ci.org/gammapy/gammapy/jobs/367153942#L1655
Should be fixed by 4631f61, which also does a few docstring formatting fixes.

@cdeil cdeil merged commit b0788a6 into gammapy:master Apr 16, 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

@joleroi joleroi deleted the joleroi:spatialmodels branch Apr 17, 2018

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