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 gammapy.region and reflected region computation #379

Merged
merged 14 commits into from Nov 14, 2015

Conversation

@joleroi
Copy link
Contributor

@joleroi joleroi commented Nov 4, 2015

This PR rewrites the reflected regions algorithm in Gammapy. It is now a function in gammapy/regions/reflected.py
It also starts a regions library containing only circles so far (started by @cdeil)

It continues #372
TODO

  • Reflected regions
  • Basic exclusion region handling
  • Tests
  • basic region library (only circles)
  • Cleanup of config handling
  • Documentation with Getting Started section
  • move examples/make_reflected_regions to docs
@cdeil cdeil added the feature label Nov 4, 2015
@cdeil cdeil added this to the 0.4 milestone Nov 4, 2015
@cdeil cdeil self-assigned this Nov 4, 2015
@cdeil
Copy link
Member

@cdeil cdeil commented Nov 4, 2015

I'll have a quick look and leave some inline comments.

@@ -0,0 +1,18 @@
# Licensed under a 3-clause BSD style license - see LICENSE.rst

This comment has been minimized.

@cdeil

cdeil Nov 4, 2015
Member

Move examples/make_reflected_regions.py to the docs?

This comment has been minimized.

@cdeil

cdeil Nov 4, 2015
Member

(and add a nice image showing the reflected regions outside the mask?)

"""
xc, yc = self.pos
major = self.radius
patch = mpatches.Ellipse((xc, yc), 2*major, 2*major, angle=0, **kwargs)

This comment has been minimized.

@cdeil

cdeil Nov 4, 2015
Member

Make the factor 2 an argument to the method (with default value 2)?

This comment has been minimized.

@joleroi

joleroi Nov 4, 2015
Author Contributor

Why do you think this is useful?
Shouldn't the MPL patch reflect the actual region and not a scaled version?

This comment has been minimized.

@cdeil

cdeil Nov 4, 2015
Member

Shouldn't the MPL patch reflect the actual region and not a scaled version?

Yes.

(I thought the factor 2 was a scaled version to make it bigger, but now I see that Ellipse takes diameter as parameter.)

Why not use Circle from matplotlib to represent a circle? 😄

This comment has been minimized.

@joleroi

joleroi Nov 5, 2015
Author Contributor

I copied that code from astropy.regions without questioning it

wcs : `~astropy.wcs.WCS` instance
The world coordinate system transformation to assume
mode : str

This comment has been minimized.

@cdeil

cdeil Nov 4, 2015
Member

I think we have no intention of supporting "affine" and "full" mode, right?
Then remove the "mode" and "tolerance" arguments here.

"""
raise NotImplementedError("")

def to_shapely(self):

This comment has been minimized.

@cdeil

cdeil Nov 4, 2015
Member

Again, we probably won't support shapely, so I'd say remove for now.
(they do have some nice stuff, e.g. for overlap testing, but if we really want that we can always add it back)

wcs : `~astropy.wcs.WCS` instance
The world coordinate system transformation to assume
mode : str

This comment has been minimized.

@cdeil

cdeil Nov 4, 2015
Member

Again, remove "mode" and "tolerance"

----------
filename : str
Name of file to write
format : str

This comment has been minimized.

@cdeil

cdeil Nov 4, 2015
Member

Multiple-choice options are documented in docstrings like this:

parameter : {"choice1", "choice2"}
    A parameter

Also, it's better to have a method that generates string and then the method that writes to file just calls that one.
(I hate APIs where you have to write to temp files or in-memory file-like buffers even if you just want a string)

from __future__ import absolute_import, division, print_function, unicode_literals
import numpy as np
from astropy.wcs import WCS
from ..image import exclusion_distance

This comment has been minimized.

@cdeil

cdeil Nov 4, 2015
Member

In some other place you put the ..image input in a method.
Should gammapy.image be allowed to import gammapy.region?
Or the other way around?
Should ExclusionMask be located in gammapy.image?

(I haven't thought about it much, but at the moment ExclusionMask doesn't use region stuff, so maybe move it to gammapy.image and at least for now say that gammapy.region can import from gammapy.image, but we'll try to avoid the other way around?)

Rotation point
exclusion_mask : `~gammapy.region.ExclusionMask`
Exlusion mask
angle_increment : float

This comment has been minimized.

@cdeil

cdeil Nov 4, 2015
Member

Even if it's a bit more typing, I think we should require Angle objects as inputs and not use float with assumed unit "deg".

This comment has been minimized.

@joleroi

joleroi Nov 4, 2015
Author Contributor

You are right, proof: I was assuming radians

reflected_regions = reflected_regions_pix.to_sky(wcs)
return reflected_regions

def _compute_xy(pix_center, offset, angle):

This comment has been minimized.

@cdeil

cdeil Nov 4, 2015
Member

I'm sure this is available already in astropy.coordinates ... if you can't find it by browsing the docs / code or playing around in IPython let me know and I'll look around.

This comment has been minimized.

@joleroi

joleroi Nov 4, 2015
Author Contributor

I am not sure what you mean, I thought astropy coordinates was only for sky coordinates

This comment has been minimized.

@cdeil

cdeil Nov 4, 2015
Member

Astropy coordinates is also for 2D and 3D cartesian coordinates.
But if you're not familiar with astropy.coordinates (I only know it a little) it takes time to figure out how to to something like this.
So maybe just put a comment for now?

# TODO: replace by calculation using `astropy.coordinates`
sky = SkyCircleRegion(pos=pos1, radius=radius)
pix = sky.to_pixel(self.wcs)
sky2 = pix.to_sky(self.wcs)
assert sky.pos.l == sky2.pos.l

This comment has been minimized.

@cdeil

cdeil Nov 4, 2015
Member

This works?
You should never do bit-by-bit comparisons on computed floats, but always call np.testing.assert_allclose (or the equivalent assert_quantity if you have quantities).
Same in the tests above.

@joleroi joleroi force-pushed the joleroi:reflected branch from fbc9ea2 to 4c5f5dd Nov 4, 2015
@joleroi joleroi force-pushed the joleroi:reflected branch from 4c5f5dd to 3d5048d Nov 4, 2015
joleroi added 3 commits Nov 4, 2015
@joleroi joleroi mentioned this pull request Nov 5, 2015
3 of 7 tasks complete
@joleroi joleroi force-pushed the joleroi:reflected branch 5 times, most recently from 8c7467c to 3ea9e7a Nov 5, 2015
@joleroi joleroi force-pushed the joleroi:reflected branch 4 times, most recently from 83e7949 to 627c812 Nov 6, 2015
@joleroi joleroi force-pushed the joleroi:reflected branch from 627c812 to 1c3bbe0 Nov 6, 2015
joleroi added 3 commits Nov 6, 2015
----------
model : str
Sherpa model
thres_high : `~gammapy.spectrum.Energy`

This comment has been minimized.

@cdeil

cdeil Nov 12, 2015
Member

Combine thres_high and thres_low into a energy_range tuple?
If you don't like that, put low before high here in the docstring, to be consistent with the order of arguments in the function.

This comment has been minimized.

@joleroi

joleroi Nov 12, 2015
Author Contributor

shouldn't it be possible to give only one of the 2 {thres_high, thres_low}?

This comment has been minimized.

@cdeil

cdeil Nov 12, 2015
Member

You're right, this is a bit tricky.
If you want these to be optional, you could:

  • stick with two arguments and add default values of None to mean "no extra restriction.
  • combine into energy_band and allow a tuple (Energy, Energy) or (Energy, None).
  • combine into energy_band as one length-2 Energy and use np.nan as "no extra restriction" and default value.

Do whatever you like best for now .... we can review API and do a consistent pattern once we've figured out what we like.

@cdeil
Copy link
Member

@cdeil cdeil commented Nov 12, 2015

@kingj90 – you've been busy!

Can you wrap this up tomorrow or Saturday latest?
I need to cut the Gammapy 0.4 release before the workshop.

@cdeil
Copy link
Member

@cdeil cdeil commented Nov 12, 2015

I briefly browsed through this, and the code looks good, but of course this is very superficial.
The most important thing is to have a working example for the tutorial (using the one Crab run from gammapy-extra?), then it's not clear to me, probably try to increase the test coverage a bit or add some high-level docs?

Let me know if you want me to review this or try it out tomorrow.

@joleroi
Copy link
Contributor Author

@joleroi joleroi commented Nov 12, 2015

@cdeil : I will wrap this up today or tomorrow and then add some high level docs and examples for the tutorial over the weekend (in a separate PR, or does this have to be in the release?).
I have not had the time to prepare an example dataset, the problem is that gammapy-spectrum is based on the DataStore class so we cannot use the one crab run as it is. Who is going to take care of that (you, me or someone else). I guess the same problem holds for the image-pipe tutorial.

@cdeil
Copy link
Member

@cdeil cdeil commented Nov 12, 2015

Are you at MPIK tomorrow? Then we could figure out a solution given the very limited time we have.
I think I can put image-pipe together quickly, but I also have to prepare my tutorial for Tuesday and various other stuff...

joleroi added 2 commits Nov 14, 2015
@joleroi
Copy link
Contributor Author

@joleroi joleroi commented Nov 14, 2015

Merging this for the 0.4 release, I will copy the remaining TODOs to a new PR

joleroi added a commit that referenced this pull request Nov 14, 2015
Reflected Regions
@joleroi joleroi merged commit 007bd7a into gammapy:master Nov 14, 2015
1 check passed
1 check passed
continuous-integration/travis-ci/pr The Travis CI build passed
Details
@joleroi joleroi deleted the joleroi:reflected branch Nov 14, 2015
@cdeil
Copy link
Member

@cdeil cdeil commented Nov 15, 2015

I like the random reflected region image in the docs:
https://gammapy.readthedocs.org/en/latest/region/index.html#reflected-regions
Will have to bookmark that page and check back for the new one from time to time ... :-)

I get one test error now:

__________________________________________________________ test_spectrum_analysis_from_configfile ___________________________________________________________

tmpdir = local('/private/var/folders/sb/4qv5j4m90pz1rw7m70rj1b1r0000gn/T/pytest-deil/pytest-100/test_spectrum_analysis_from_co0')

    @requires_dependency('yaml')
    @requires_data('gammapy-extra')
    @requires_data('hess')
    def test_spectrum_analysis_from_configfile(tmpdir):
        import yaml

        configfile = gammapy_extra.filename('test_datasets/spectrum/spectrum_analysis_example_ring.yaml')

        import yaml
        with open(configfile) as fh:
            config = yaml.safe_load(fh)

        config['general']['outdir']=str(tmpdir)

>       fit = run_spectrum_analysis_using_config(config)

gammapy/spectrum/tests/test_spectrum_analysis.py:25: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
gammapy/spectrum/spectrum_analysis.py:483: in run_spectrum_analysis_using_config
    exclusion = ExclusionMask.from_fits(excl_file)
gammapy/image/mask.py:81: in from_fits
    hdulist = fits.open(excl_file)
/Users/deil/Library/Python/3.4/lib/python/site-packages/astropy-1.1.dev13977-py3.4-macosx-10.11-x86_64.egg/astropy/io/fits/hdu/hdulist.py:138: in fitsopen
    return HDUList.fromfile(name, mode, memmap, save_backup, cache, **kwargs)
/Users/deil/Library/Python/3.4/lib/python/site-packages/astropy-1.1.dev13977-py3.4-macosx-10.11-x86_64.egg/astropy/io/fits/hdu/hdulist.py:280: in fromfile
    save_backup=save_backup, cache=cache, **kwargs)
/Users/deil/Library/Python/3.4/lib/python/site-packages/astropy-1.1.dev13977-py3.4-macosx-10.11-x86_64.egg/astropy/io/fits/hdu/hdulist.py:801: in _readfrom
    ffo = _File(fileobj, mode=mode, memmap=memmap, cache=cache)
/Users/deil/Library/Python/3.4/lib/python/site-packages/astropy-1.1.dev13977-py3.4-macosx-10.11-x86_64.egg/astropy/io/fits/file.py:141: in __init__
    self._open_filename(fileobj, mode, clobber)
/Users/deil/Library/Python/3.4/lib/python/site-packages/astropy-1.1.dev13977-py3.4-macosx-10.11-x86_64.egg/astropy/io/fits/file.py:493: in _open_filename
    self._file = fileobj_open(self.name, PYFITS_MODES[mode])
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

filename = '/home/kingj/Software/gammapy-extra/test_datasets/spectrum/dummy_exclusion.fits', mode = 'rb'

    def fileobj_open(filename, mode):
        """
            A wrapper around the `open()` builtin.

            This exists because in Python 3, `open()` returns an
            `io.BufferedReader` by default.  This is bad, because
            `io.BufferedReader` doesn't support random access, which we need in
            some cases.  In the Python 3 case (implemented in the py3compat module)
            we must call open with buffering=0 to get a raw random-access file
            reader.
            """

>       return open(filename, mode, buffering=0)
E       FileNotFoundError: [Errno 2] No such file or directory: '/home/kingj/Software/gammapy-extra/test_datasets/spectrum/dummy_exclusion.fits'

/Users/deil/Library/Python/3.4/lib/python/site-packages/astropy-1.1.dev13977-py3.4-macosx-10.11-x86_64.egg/astropy/io/fits/util.py:400: FileNotFoundError

The issue is that you've hardcoded an absolute path here:
https://github.com/gammapy/gammapy-extra/blob/master/test_datasets/spectrum/spectrum_analysis_example_reflected.yaml#L28

excluded_regions:
  file: /home/kingj/Software/gammapy-extra/test_datasets/spectrum/dummy_exclusion.fits 

I have the GAMMAPY_EXTRA environment variable set on my computer, but putting $GAMMAPY_EXTRA/test_datasets/spectrum/dummy_exclusion.fits doesn't work, it's not automatically expanded.

Not sure how to make this work ... @kingj90 any thoughts?
Maybe we should do it by calling os.path.expandvars explicitly?
I guess putting relative filenames to the YAML file is not a good idea ... supporting environment variables in paths would be better?

@cdeil cdeil changed the title Reflected Regions Add gammapy.region and reflected region computation Dec 12, 2015
@cdeil
Copy link
Member

@cdeil cdeil commented Dec 12, 2015

@joleroi - I'd like to finally make the Gammapy 0.4 release next week and started some cleanup.

The PR description (i.e. the first comment on this issue) still has a task list with open points and the description is a bit cryptic. This happens to me also all the time ... the content and scope of PRs changes while I work on them. I think it helps a bit if the PR title and description and task list is cleaned up a bit once it's merged and the changelog entry is made. I don't mean spend a lot of time on it, just ~ 1 minute, because the fraction of people that will look at the PR in the future is small (mostly Gammapy maintainers that try to find out "why was this done this way?").

Can you please edit the description of this PR and your future PRs?

@joleroi
Copy link
Contributor Author

@joleroi joleroi commented Dec 14, 2015

There should be a changelog entry for this PR
Why don't we skip 0.4 an go for 0.5 directly? Is due son anyways, right?

@cdeil
Copy link
Member

@cdeil cdeil commented Dec 14, 2015

There should be a changelog entry for this PR

It's here: https://gammapy.readthedocs.org/en/latest/changelog.html#gammapy-0p4-release

Why don't we skip 0.4 an go for 0.5 directly? Is due son anyways, right?

I've adapted the milestone dates on GH.
I don't see a reason to skip 0.4

It's kind of bad how the schedule slips, but it's also normal for open-source projects, where the devs are scientists that work on this when they get a chance.
Having milestones with a planned release date is still better than no schedule at all, i.e. suddenly saying "I'm making a release now", no?

@joleroi
Copy link
Contributor Author

@joleroi joleroi commented Dec 14, 2015

I guess this happens to all projects (no only open source with scientist as devs). But I agree, that having milestones is a good thing!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

2 participants
You can’t perform that action at this time.