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

Improve ReflectedRegionsFinder class #2089

Merged
merged 41 commits into from Jun 26, 2019
Merged

Improve ReflectedRegionsFinder class #2089

merged 41 commits into from Jun 26, 2019

Conversation

@bkhelifi
Copy link
Contributor

@bkhelifi bkhelifi commented Mar 26, 2019

Hi all (after the closing of PR #2079),

As discussed during the last Coding Sprint, I have extended with @registerrier the features of this class in order to better analyse data (and the HESS data).

The class has been generalised to accept any astropy.regions, except the polygons. This shape could be a topic of upgrade. For a Circle shape, the new statistics for the ON region are identical to the previous implementation of the class; with this implementation I got slightly more OFF regions (+1). Note that a deepcopy of some astropy.regions classes does not work (problem of implementation of these classes? e.g. EllipseAnnulusSkyRegion): a trick should have been found that is hardly readable.

Together with the class, the test functions have been enlarged, and a new notebook associated, aiming to be a tutorial .

The possible improvements after this PR are:

  • Usage of PolygonSkyRegion,
  • Push the astropy.region developers to implement the panda shapes (ie the 'wedges', really needed for some PWN analysis, or the GC), and upgrade the RefelectedBackground class accordingly.

Good luck to the nice and cool reviewers;)
Bruno

@bkhelifi bkhelifi added this to the 0.12 milestone Mar 26, 2019
@bkhelifi bkhelifi self-assigned this Mar 26, 2019
@bkhelifi bkhelifi added this to To do in gammapy.makers via automation Mar 26, 2019
@bkhelifi bkhelifi requested review from adonath and registerrier Mar 26, 2019
Copy link
Contributor

@registerrier registerrier left a comment

Thanks @bkhelifi ! This adds an important functionality.

I have added some comments.

Is the event selection based on Map good enough or would we be better off with doing directly the selection on SkyRegion?

How about polygonalization of the regions, could it provide a more general solution? I assume this would require changes on the regions package itself.

Also, the change breaks most the tests using ReflectedRegionFinder, it the change is expected/acceptable the assert should be modified accordingly.

_test_pos = self._compute_xy(self._pix_center, self._offset, curr_angle)

# The function copy.deepcopy does not work for these regions classes
dict = { _ : getattr(self._pix_region, _) for _ in self._pix_region._repr_params}
Copy link
Contributor

@registerrier registerrier Mar 26, 2019

You should not call a dict dict this will mess with the built-in dict.

Copy link
Contributor Author

@bkhelifi bkhelifi Apr 1, 2019

ok... thanks

def _get_bounding_size(region, center):
"""Return the typical radius of the bounding size"""

_map = WcsNDMap.create(
Copy link
Contributor

@registerrier registerrier Mar 26, 2019

Why create a new map here? You only need a MapGeom, not the full array of data that comes with the Map. You could also pass it to the method, since you already use one.

Copy link
Contributor Author

@bkhelifi bkhelifi Apr 17, 2019

ok. done


if 'ra' in reg_center.representation_component_names:
_plotmap = WcsNDMap.create(
skydir=reg_center, binsz=self.binsz, width=10., coordsys="CEL", proj="TAN"
Copy link
Contributor

@registerrier registerrier Mar 26, 2019

Why does the Map have a default size of 10 degrees? Shouldn't it adapt either to user input or to the actual distribution of reflected regions?

Copy link
Contributor Author

@bkhelifi bkhelifi Apr 1, 2019

For the SCT, one might have a 10deg FoV. The user might not known which telescope type is used, so I used 10deg as the worst case.

dx = offset * np.sin(angle)
dy = offset * np.cos(angle)
"""Compute x, y position for a given position angle and offset."""
dx = offset * np.cos(angle)
Copy link
Contributor

@registerrier registerrier Mar 26, 2019

Why the change cos -> sin.

Copy link
Contributor Author

@bkhelifi bkhelifi Apr 1, 2019

Because it was a bug!

_mask = geom.region_mask(self.reflected_regions, inside=True)
self.off_reference_map = WcsNDMap(geom=geom, data=_mask)
# TODO: remove this lign after the correction of the Issue #2074
self.off_reference_map.data = self.off_reference_map.data.astype(float)
Copy link
Contributor

@registerrier registerrier Mar 26, 2019

Would a SkyRegion based event selection make things easier regarding bin size?

Copy link
Contributor Author

@bkhelifi bkhelifi Apr 1, 2019

This functionality was not present at the time of the PR. If it exists, one could use also the proposed solution

log.debug("Placing reflected region\n{}".format(refl_region))
reflected_regions.append(refl_region)
_test_reg = self._create_rotated_reg(curr_angle)
if _pixel_excluded is False or not np.any(_test_reg.contains(_excluded_pixcoords)):
Copy link
Contributor

@registerrier registerrier Mar 26, 2019

Is this approach rapid enough? Does it impact the performances for exclusion masks with many excluded pixels?

How about testing if pixels inside the rotated region are excluded or not?

Copy link
Contributor Author

@bkhelifi bkhelifi Apr 1, 2019

Can you precise me how to quantify this and what is the threshold of 'rapid enough'?
I do not understand the second question

@@ -62,29 +63,39 @@ def test_find_reflected_regions(exclusion_mask, on_region):
)
fregions.run()
regions = fregions.reflected_regions
assert len(regions) == 14
assert_quantity_allclose(regions[3].center.icrs.ra, Angle("83.674 deg"), rtol=1e-2)
assert len(regions) == 15
Copy link
Contributor

@registerrier registerrier Mar 26, 2019

Where does the change come from? Bin size?

Copy link
Contributor Author

@bkhelifi bkhelifi Apr 17, 2019

My guess is that it comes from the fact that the logics to determine the starting angle of OFF regions has been changed.


# Test without exclusion
fregions.exclusion_mask = None
fregions.run()
regions = fregions.reflected_regions
assert len(regions) == 16
assert len(regions) == 17
Copy link
Contributor

@registerrier registerrier Mar 26, 2019

A change without mask is bit surprising no? Where does it come from?

Copy link
Contributor Author

@bkhelifi bkhelifi Apr 17, 2019

See above

assert len(regions) == 16
assert_quantity_allclose(regions[3].center.icrs.ra, Angle("83.674 deg"), rtol=1e-2)
assert len(regions) == 17
assert_quantity_allclose(regions[14].center.icrs.ra, Angle("83.868 deg"), rtol=1e-2)
Copy link
Contributor

@registerrier registerrier Mar 26, 2019

This is due to the change of rotation angle in the region definition?

Copy link
Contributor Author

@bkhelifi bkhelifi Apr 17, 2019

yes, I think so

fregions.exclusion_mask = small_mask
fregions.run()
regions = fregions.reflected_regions
assert len(regions) == 16
assert_quantity_allclose(regions[3].center.icrs.ra, Angle("83.674 deg"), rtol=1e-2)
assert len(regions) == 17
Copy link
Contributor

@registerrier registerrier Mar 26, 2019

Why 16 - > 17?

Copy link
Contributor Author

@bkhelifi bkhelifi Apr 17, 2019

See above

Copy link
Member

@cdeil cdeil left a comment

@bkhelifi - Apologies for not responding here sooner.

This PR was briefly discussed in the Gammapy dev call today. As I mentioned, I'll open a PIG for region handling in Gammapy today that's related. There's a few design decisions concerning whether pixel or sky regions are used, and specifically for reflected regions how the rotation is computed. Probably this will be controversial and the final region solution for Gammapy 1.0 will require quite some discussion in the next weeks, and then ~ 1 month of work to implement.

@registerrier mentioned that he'd like to finish this PR and get it merged. In this case, please see the inline comments with suggestions. There's already a merge conflict, please do it soon.

from .background_estimate import BackgroundEstimate

__all__ = ["ReflectedRegionsFinder", "ReflectedRegionsBackgroundEstimator"]

log = logging.getLogger(__name__)


def _compute_distance_image(mask_map):
Copy link
Member

@cdeil cdeil May 3, 2019

This PR removes the distance image implementation. It is used to make reflected region finding fast for the common case of circular regions. I think this speed-up really is needed: one of my first contributions to the HESS software was to add this, after running spectral analyses for HGPS sources. What I found was that reflected region computation was dominating analysis time overall completely - even event reconstruction, gamma-hadron separation, spectral likelihood fitting, everything else was much faster.

Can you please put it back and keep reflected region analysis fast for the common case of having circles (like for AGN)?

Copy link
Contributor Author

@bkhelifi bkhelifi May 3, 2019

I will make comparison of speed tests, put the results here, which will justify the retained solution.

Copy link
Contributor Author

@bkhelifi bkhelifi May 6, 2019

Your method using the optimised function '_compute_distance_image' is based on the assumption that a region has a radius (cf. _is_inside_exclusion). If one used an outer radius, more pixels than in reality will be considered as excluded, which is not optimal in this regards. In this PR, one uses the function 'contains' of a region, that garanties to have exact pixels.
Do you have an idea to adapt the function '_is_inside_exclusion' to be more accurate?

Copy link
Member

@cdeil cdeil May 6, 2019

The distance image algorithm is only for the case of a circular region.

Copy link
Contributor Author

@bkhelifi bkhelifi May 6, 2019

Ok. So, should I keep these algo only for the case of circular regions?
The pros: most of the case, we will have circular region, and so we gain some computing time.
The cons: more lines for just one case.
Any preference?

Copy link
Member

@cdeil cdeil May 6, 2019

Yes, please keep the fast case for circular regions.

Copy link
Contributor Author

@bkhelifi bkhelifi May 6, 2019

I have re-implemented this algorithm. The timing measurement are on the notebook intro_reflected_background using a circular region:
Christoph algo: total CPUTime=25.4s, Wall time=7.75s
Current algo using PixRegion.contains: total CPUTime=25.2s, Wall time=8.07s
The time different is tiny!
Should we thus keep an implementation using a different method depending on the On region type? My personal feeling is thus no.

Copy link
Member

@cdeil cdeil May 7, 2019

@bkhelifi - I see is https://github.com/gammapy/gammapy/blob/22eda59cf6feff6a29b421afc2a30a0158ec86f2/tutorials/intro_reflected_background.ipynb but it has no timings. Please share a link to the notebook with your benchmark.

A good way to share a notebook or some other file that contains something that's not supposed to be merged with a PR is on the side, e.g. in a repo in your account, or I often post a gist: https://gist.github.com/ .

Not sure what test case you measured. Here's one that could be used: https://docs.gammapy.org/0.11/background/reflected.html
The time will depend on the number of off regions tested, so a smaller on region (e.g. 0.1 deg) and larger FOV offset (e.g. 2 deg) so that there are ~ 100 off regions would make a good test case.
With the current algorithm in Gammapy master it will be fast, with that removed I would expect it to be slow.

Copy link
Contributor Author

@bkhelifi bkhelifi May 7, 2019

By using the code given in the web page that you mentioned, I made a fast notebook. Here are the results, first using your method, then the new one.

For each of the finder given in the notebook, the results for your method are:
CPU times: user 12.6 s, sys: 400 ms, total: 13 s; Wall time: 4.37 s
CPU times: user 12.1 s, sys: 291 ms, total: 12.4 s; Wall time: 4.04 s
CPU times: user 12.5 s, sys: 307 ms, total: 12.8 s; Wall time: 4.13 s

With the new method:
CPU times: user 12 s, sys: 413 ms, total: 12.4 s; Wall time: 4.32 s
CPU times: user 11.8 s, sys: 385 ms, total: 12.2 s; Wall time: 4.15 s
CPU times: user 11.9 s, sys: 356 ms, total: 12.3 s; Wall time: 4.18 s

When using an offset of 2deg, 23 reflected regions are built (r=0.3deg). The execution times are the the first finder:
Your method: CPU times: user 16 s, sys: 941 ms, total: 16.9 s; Wall time: 6.95 s
The new method: CPU times: user 15.5 s, sys: 837 ms, total: 16.3 s; Wall time: 7.3 s

The conclusions are the identical. The current algorithm is almost as fast as the previous one

Copy link
Member

@cdeil cdeil May 7, 2019

@bkhelifi - Can you please share the notebook in a gist so that I can try it out as well?

I'll be at a meeting the next days, will only be able to look at this again Monday.

@registerrier - If you have time to review this, please go ahead. If you're happy with the PR, feel free to merge, no need to wait for me.

# Make the ON reference map
_mask = geom.region_mask([self.region], inside=True)
self.on_reference_map = WcsNDMap(geom=geom, data=_mask)
# TODO: remove this lign after the correction of the Issue #2074
Copy link
Member

@cdeil cdeil May 3, 2019

Resolve this TODO?

Copy link
Contributor Author

@bkhelifi bkhelifi May 3, 2019

This ToDo is a reference towards an on-going work, that is not absolutely necessary for this class. It is just a fair mention of possible future code clean-up.
I would not like that this PR is delayed just because of that... This class should be in the next gammapy version.

Copy link
Member

@cdeil cdeil May 3, 2019

The TODO references issue #2074, but this issue was resolved and closed by @adonath via #2078 on March 28. So these TODOs should be resolved in this PR before merging.

@registerrier - usually what I do is to push commits to a PR branch and then hit the merge button. This is better than merging first and then follow-up commits in master or new PRs, because it's cleaner in the changelog and triggers CI builds in a better way (don't break master).

min_ang = np.max(angles)-np.min(angles)
if min_ang.value > np.pi:
new_angles = np.zeros_like(angles)
new_angles[angles.value > 0] = angles[angles.value > 0] - np.pi*u.rad
Copy link
Member

@cdeil cdeil May 3, 2019

Can this be simplified? astropy.coordinates.Angle has methods that implement such wrapping logic?

Copy link
Contributor Author

@bkhelifi bkhelifi May 6, 2019

Except using the class Longitude, I do not see how to to... Should I use this class?

Copy link
Member

@cdeil cdeil May 6, 2019

Try to use the Angle.wrap_at method.

>>> from astropy.coordinates import Angle
>>> a = Angle([-200, -300, 0, 100, 200, 300, 400], unit='deg')
>>> a.wrap_at('180 deg')
<Angle [ 160.,   60.,    0.,  100., -160.,  -60.,   40.] deg>

Copy link
Contributor Author

@bkhelifi bkhelifi May 6, 2019

done

curr_angle = curr_angle + self._min_ang
if self.max_region_number <= len(reflected_regions):
break
else:
curr_angle = curr_angle + self.angle_increment

print("Found {0} reflected regions for the Obs #{1}".format(len(reflected_regions), self.obs_id))
Copy link
Member

@cdeil cdeil May 3, 2019

Never print from a library. If you want to communicate with users, use log.

Copy link
Contributor Author

@bkhelifi bkhelifi May 3, 2019

ok. corrected

# Make the OFF reference map
_mask = geom.region_mask(self.reflected_regions, inside=True)
self.off_reference_map = WcsNDMap(geom=geom, data=_mask)
# TODO: remove this lign after the correction of the Issue #2074
Copy link
Member

@cdeil cdeil May 3, 2019

Resolve this TODO?

Copy link
Contributor Author

@bkhelifi bkhelifi May 6, 2019

done

dict = { _ : getattr(self._pix_region, _) for _ in self._pix_region._repr_params}
_test_reg = self._pix_region.__class__(_test_pos, **dict)

if hasattr(_test_reg, 'angle'):
Copy link
Member

@cdeil cdeil May 3, 2019

Yes, an error should be raised if this doesn't work for some regions (and a test added to make sure that it does).

$ python -c 'import this' | grep Error
Errors should never pass silently.

binsz = kwargs.get("binsz")
else:
binsz = self._get_bounding_size(on_region, observations[0].pointing_radec)/10.
if binsz > 0.01*u.deg:
Copy link
Member

@cdeil cdeil May 3, 2019

Do results really depend on bin size?
Probably because of the distance image?
Otherwise it shouldn't, no?

If you really need to keep this parameter "0.01 deg" - make it a constant with a name and use it everywhere, don't duplicate the hard-coded number several times.

@cdeil cdeil mentioned this pull request May 3, 2019
"""Return the typical radius of the bounding size"""

geom = WcsGeom.create(skydir=center, width=10.*u.deg, binsz=Angle("0.01 deg"), coordsys="GAL", proj="TAN")

Copy link
Contributor

@registerrier registerrier May 21, 2019

After testing, this fails for regions and/or positions that located close to the poles because the ixmaxand ixmin
values are not defined.
For safety a test should be made to determine the most relevant coordinate system to use.

A possible solution is to determine if the bin size is sufficiently large in the ReflectedRegionFinder itself during the run execution, when the reference_map for the current observation is created.

@registerrier registerrier removed this from the 0.12 milestone May 27, 2019
@registerrier registerrier added this to the 0.13 milestone May 27, 2019
@cdeil
Copy link
Member

@cdeil cdeil commented Jun 17, 2019

I've added copy and rotate methods for all PixelRegion classes, and just now released regions v0.4 with this feature: https://astropy-regions.readthedocs.io/en/v0.4/changelog.html

My suggestion would be bump to regions>=0.4 for Gammapy now, and remove home-grown region copy and rotate code here.

@registerrier - Thoughts?

@cdeil
Copy link
Member

@cdeil cdeil commented Jun 17, 2019

Discussed with @registerrier and @adonath just now in a call: plan here is to bump the regions requirement to v0.4 and use region.rotate from there now and clean this up as much as possible and get it in.

@registerrier - Could you please pu regions>=0.4 in setup.py

'regions>=0.3',

and then we see if we can get CI green for this PR?

@registerrier
Copy link
Contributor

@registerrier registerrier commented Jun 26, 2019

This should be ready now.

  • I have made tests with heavy use of exclusion regions to test the time difference with and without the distance_image approach. The latter is a bit faster. Typically, 150 ms while the current implementation takes about 250-300 ms. This seems a rather acceptable difference. If this proves a problem, we might reintroduce it for CircleSkyRegion.
  • The code now uses region v0.4 and setup.py has been modified to ensure this version is used. The Region.rotate is used. The rotation is performed towards positive angles. This introduces a change compared to the previous implementation. I kept rotation along positive angles since this is more natural. I had to update some numerical results in the tests to adapt to the new. In one or two cases this resulted in a difference of one exclusion region found.
  • The content of the notebook has been moved into documentation. Several example files have been added to show how tu use regions and to show the result for some Rectangular regions.
  • I have added tests for several circular regions as well as non circular regions. I also added a test for a bad region where the center is the same as the pointing to ensure that no error is encountered and that the list of off regions is empty.
  • The ON and OFF selection is done using the EventList.select_region recently introduced in PR #2254 . Rather than passing a list of OFF region a union is created. This might be simplified if we introduce a helper function for that.

@cdeil cdeil removed the request for review from adonath Jun 26, 2019
@cdeil cdeil assigned registerrier and unassigned bkhelifi Jun 26, 2019
cdeil
cdeil approved these changes Jun 26, 2019
Copy link
Member

@cdeil cdeil left a comment

@bkhelifi and @registerrier - Thanks!

I'm merging this now.

Let's continue the work on regions-related code / tests / docs in future issues and pull requests.

@cdeil cdeil merged commit 44dd2ba into gammapy:master Jun 26, 2019
7 of 9 checks passed
gammapy.makers automation moved this from To do to Done Jun 26, 2019
@adonath adonath changed the title Improved Reflected Background class Improve ReflectedRegionsFinder class Jul 25, 2019
@bkhelifi bkhelifi deleted the reflected_bis branch Oct 27, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
Linked issues

Successfully merging this pull request may close these issues.

None yet

3 participants