Commit b335e468 authored by Leo Pound Singer's avatar Leo Pound Singer
Browse files

Use robust estimators for location and eccentricity of ellipses

This produces better results on a wide variety of injections.
parent 5f18e1b6
Pipeline #13511 passed with stages
in 30 minutes and 48 seconds
......@@ -777,19 +777,28 @@ def posterior_max(prob, nest=False):
*hp.pix2ang(nside, i, nest=nest, lonlat=True), unit=u.deg)
# FIXME: use weighted percentile function.
# See https://github.com/numpy/numpy/pull/9211
def wmedian(x, weights):
i = np.argsort(x)
c = np.cumsum(weights[i])
return np.interp(0.5, c / c[-1], x[i])
def find_ellipse(prob, cl=90, projection='ARC', nest=False):
"""For a HEALPix map, find an ellipse that contains a given probability.
The center of the ellipse is the mean a posteriori sky position. The length
and orientation of the semi-major and semi-minor axes are measured as
follows:
The center of the ellipse is the median a posteriori sky position. The
length and orientation of the semi-major and semi-minor axes are measured
as follows:
1. The sky map is transformed to a WCS projection that may be specified by
the caller. The default projection is ``ARC`` (zenithal equidistant), in
which radial distances are proportional to the physical angular
separation from the center point.
2. The 1-sigma ellipse is found in that projection by principal component
analysis.
2. A 1-sigma ellipse is estimated by calculating the covariance matrix in
the projected image plane using three rounds of sigma clipping to reject
distant outlier points.
3. The 1-sigma ellipse is inflated until it encloses an integrated
probability of ``cl`` (default: 90%).
......@@ -863,7 +872,7 @@ def find_ellipse(prob, cl=90, projection='ARC', nest=False):
>>> ra, dec, a, b, pa, area = find_ellipse(prob)
>>> print(*np.around([ra, dec, a, b, pa, area], 5))
195.12382 -19.6692 7.58123 1.3988 64.66371 33.28661
195.03733 -19.29358 8.66547 1.1793 63.61698 32.07665
>>> s = 'fk5;ellipse({},{},{},{},{})'.format(ra, dec, a, b, pa)
>>> open('ds9.reg', 'w').write(s) # doctest: +SKIP
......@@ -920,7 +929,7 @@ def find_ellipse(prob, cl=90, projection='ARC', nest=False):
>>> skymap_moc = read_sky_map(healpix_hdu, moc=True)
>>> ellipse = find_ellipse(skymap_moc)
>>> print(*np.around(ellipse, 5))
195.12382 -19.66919 7.59245 1.40026 64.66363 33.03595
195.03715 -19.27587 8.67609 1.18167 63.60452 32.08015
Example 3
~~~~~~~~~
......@@ -948,12 +957,12 @@ def find_ellipse(prob, cl=90, projection='ARC', nest=False):
>>> prob = make_uniform_in_sin_theta(1)
>>> ra, dec, a, b, pa, area = np.around(find_ellipse(prob), 5) + 0
>>> print(dec, a, b, area)
90.0 0.82241 0.82241 2.36052
89.90863 0.87034 0.87034 2.37888
>>> prob = make_uniform_in_sin_theta(10)
>>> ra, dec, a, b, pa, area = np.around(find_ellipse(prob), 5) + 0
>>> print(dec, a, b, area)
90.0 9.05512 9.05512 255.10577
89.90828 9.02485 9.02484 255.11972
>>> prob = make_uniform_in_sin_theta(120)
>>> ra, dec, a, b, pa, area = np.around(find_ellipse(prob), 5) + 0
......@@ -992,7 +1001,7 @@ def find_ellipse(prob, cl=90, projection='ARC', nest=False):
...
>>> ra, dec, a, b, pa, area = np.around(find_ellipse(prob), 5) + 0
>>> print(ra, dec, a, b, area)
45.0 0.0 2.14209 2.14209 14.4677
44.99999 0.0 2.14242 2.14209 14.4677
This one is centered at RA=45°, Dec=0°, and is elongated in the north-south
direction.
......@@ -1003,7 +1012,7 @@ def find_ellipse(prob, cl=90, projection='ARC', nest=False):
...
>>> ra, dec, a, b, pa, area = np.around(find_ellipse(prob), 5) + 0
>>> print(ra, dec, a, b, pa, area)
45.0 0.0 13.44746 2.1082 90.0 88.55474
44.99998 0.0 13.58769 2.08298 90.0 88.57797
This one is centered at RA=0°, Dec=0°, and is elongated in the east-west
direction.
......@@ -1014,7 +1023,7 @@ def find_ellipse(prob, cl=90, projection='ARC', nest=False):
...
>>> ra, dec, a, b, pa, area = np.around(find_ellipse(prob), 5) + 0
>>> print(dec, a, b, pa, area)
0.0 13.4194 2.1038 0.0 88.61007
0.0 13.58392 2.08238 0.0 88.54623
This one is centered at RA=0°, Dec=0°, and has its long axis tilted about
10° to the west of north.
......@@ -1027,7 +1036,7 @@ def find_ellipse(prob, cl=90, projection='ARC', nest=False):
...
>>> ra, dec, a, b, pa, area = np.around(find_ellipse(prob), 5) + 0
>>> print(dec, a, b, pa, area)
0.0 63.82809 34.00824 80.78253 6379.84049
0.0 64.77133 33.50754 80.78231 6372.34466
This one is centered at RA=0°, Dec=0°, and has its long axis tilted about
10° to the east of north.
......@@ -1040,7 +1049,7 @@ def find_ellipse(prob, cl=90, projection='ARC', nest=False):
...
>>> ra, dec, a, b, pa, area = np.around(find_ellipse(prob), 5) + 0
>>> print(dec, a, b, pa, area)
0.0 63.82809 34.00824 99.21747 6379.84049
0.0 64.77133 33.50754 99.21769 6372.34466
This one is centered at RA=0°, Dec=0°, and has its long axis tilted about
80° to the east of north.
......@@ -1053,7 +1062,7 @@ def find_ellipse(prob, cl=90, projection='ARC', nest=False):
...
>>> ra, dec, a, b, pa, area = np.around(find_ellipse(prob), 5) + 0
>>> print(dec, a, b, pa, area)
0.0 63.82533 34.00677 170.78252 6379.81219
0.0 64.77564 33.50986 170.78252 6372.42573
This one is centered at RA=0°, Dec=0°, and has its long axis tilted about
80° to the west of north.
......@@ -1066,7 +1075,7 @@ def find_ellipse(prob, cl=90, projection='ARC', nest=False):
...
>>> ra, dec, a, b, pa, area = np.around(find_ellipse(prob), 5) + 0
>>> print(dec, a, b, pa, area)
0.0 63.82533 34.00677 9.21748 6379.81219
0.0 64.77564 33.50986 9.21748 6372.42573
"""
try:
prob['UNIQ']
......@@ -1084,9 +1093,9 @@ def find_ellipse(prob, cl=90, projection='ARC', nest=False):
area *= np.square(180 / np.pi)
nest = True
# Find mean right ascension and declination.
xyz0 = (hp.pix2vec(nside, ipix, nest=nest) * prob).sum(axis=1)
(ra,), (dec,) = hp.vec2ang(xyz0, lonlat=True)
# Find median a posteriori sky position.
xyz0 = [wmedian(x, prob) for x in hp.pix2vec(nside, ipix, nest=nest)]
(ra,), (dec,) = hp.vec2ang(np.asarray(xyz0), lonlat=True)
# Construct WCS with the specified projection
# and centered on mean direction.
......@@ -1106,8 +1115,13 @@ def find_ellipse(prob, cl=90, projection='ARC', nest=False):
prob = prob[keep]
area = area[keep]
# Find covariance matrix.
c = cov(xy, aweights=prob, rowvar=False)
# Find covariance matrix, performing three rounds of sigma-clipping
# to reject outliers.
keep = np.ones(len(xy), dtype=bool)
for _ in range(3):
c = cov(xy[keep], aweights=prob[keep], rowvar=False)
nsigmas = np.sqrt(np.sum(xy.T * np.linalg.solve(c, xy.T), axis=0))
keep &= (nsigmas < 3)
# If each point is n-sigma from the center, find n.
nsigmas = np.sqrt(np.sum(xy.T * np.linalg.solve(c, xy.T), axis=0))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment