From 03b42dc70de207e66b53327f4e8e125f86191e6c Mon Sep 17 00:00:00 2001 From: Melissa DeLucchi Date: Mon, 18 Nov 2024 10:50:54 -0500 Subject: [PATCH 1/2] Remove margin fine filtering. --- .github/workflows/pre-commit-ci.yml | 2 +- .github/workflows/testing-and-coverage.yml | 2 +- src/hats/pixel_math/__init__.py | 1 - src/hats/pixel_math/healpix_shim.py | 20 --- src/hats/pixel_math/margin_bounding.py | 139 ------------------ tests/hats/pixel_math/test_margin_bounding.py | 109 -------------- tests/hats/pixel_math/test_spatial_index.py | 19 ++- 7 files changed, 11 insertions(+), 281 deletions(-) delete mode 100644 src/hats/pixel_math/margin_bounding.py delete mode 100644 tests/hats/pixel_math/test_margin_bounding.py diff --git a/.github/workflows/pre-commit-ci.yml b/.github/workflows/pre-commit-ci.yml index f31b3ee6..642e7b9d 100644 --- a/.github/workflows/pre-commit-ci.yml +++ b/.github/workflows/pre-commit-ci.yml @@ -7,7 +7,7 @@ on: push: branches: [ main ] pull_request: - branches: [ main ] + branches: [ main, margin ] jobs: pre-commit-ci: diff --git a/.github/workflows/testing-and-coverage.yml b/.github/workflows/testing-and-coverage.yml index 3d3b867f..8438b944 100644 --- a/.github/workflows/testing-and-coverage.yml +++ b/.github/workflows/testing-and-coverage.yml @@ -7,7 +7,7 @@ on: push: branches: [ main ] pull_request: - branches: [ main ] + branches: [ main, margin ] jobs: build: diff --git a/src/hats/pixel_math/__init__.py b/src/hats/pixel_math/__init__.py index 132eaa84..327d83f9 100644 --- a/src/hats/pixel_math/__init__.py +++ b/src/hats/pixel_math/__init__.py @@ -2,7 +2,6 @@ from .healpix_pixel import HealpixPixel from .healpix_pixel_convertor import HealpixInputTypes, get_healpix_pixel -from .margin_bounding import check_margin_bounds from .partition_stats import empty_histogram, generate_alignment, generate_histogram from .pixel_margins import get_margin from .spatial_index import compute_spatial_index, spatial_index_to_healpix diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index 766bb5d7..44aed38b 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -20,18 +20,10 @@ def npix2nside(param): return hp.npix2nside(param) -def nside2npix(param): - return hp.nside2npix(param) - - def order2npix(param): return hp.order2npix(param) -def npix2order(param): - return hp.npix2order(param) - - def nside2resol(*args, **kwargs): return hp.nside2resol(*args, **kwargs) @@ -63,10 +55,6 @@ def query_polygon(*args, **kwargs): return hp.query_polygon(*args, **kwargs) -def boundaries(*args, **kwargs): - return hp.boundaries(*args, **kwargs) - - def get_all_neighbours(*args, **kwargs): return hp.get_all_neighbours(*args, **kwargs) @@ -78,14 +66,6 @@ def ang2vec(*args, **kwargs): return hp.ang2vec(*args, **kwargs) -def pix2ang(*args, **kwargs): - return hp.pix2ang(*args, **kwargs) - - -def vec2dir(*args, **kwargs): - return hp.vec2dir(*args, **kwargs) - - def pix2xyf(*args, **kwargs): return hp.pix2xyf(*args, **kwargs) diff --git a/src/hats/pixel_math/margin_bounding.py b/src/hats/pixel_math/margin_bounding.py deleted file mode 100644 index 8234f740..00000000 --- a/src/hats/pixel_math/margin_bounding.py +++ /dev/null @@ -1,139 +0,0 @@ -"""Utilities to build bounding boxes around healpixels that include a neighor margin.""" - -import numpy as np -from astropy.coordinates import SkyCoord -from numba import njit - -import hats.pixel_math.healpix_shim as hp - - -def check_margin_bounds(r_asc, dec, pixel_order, pixel, margin_threshold, step=100, chunk_size=1000): - # pylint: disable=too-many-locals - """Check a set of points in ra and dec space to see if they fall within the - `margin_threshold` of a given healpixel. - - Args: - r_asc (:obj:`np.array`): one dimmensional array representing the ra of a - given set of points. - dec (:obj:`np.array`): one dimmensional array representing the dec of a - given set of points. - pixel_order (int): the order of the healpixel to be checked against. - pixel (int): the healpixel to be checked against. - margin_threshold (double): the size of the margin cache in arcseconds. - step (int): the amount of samples we'll get from the healpixel boundary - to test the datapoints against. the higher the step the more - granular the point checking and therefore the more accurate, however - using a smaller step value might be helpful for super large datasets - to save compute time if accuracy is less important. - chunk_size (int): the size of the chunk of data points we process on a single - thread at any given time. We only process a certain amount of the data - at once as trying to check all data points at once can lead to - memory issues. - Returns: - :obj:`numpy.array` of boolean values corresponding to the ra and dec - coordinates checked against whether a given point is within any of the - provided polygons. - """ - num_points = len(r_asc) - if num_points != len(dec): - raise ValueError( - f"length of r_asc ({num_points}) is different \ - from length of dec ({len(dec)})" - ) - - if num_points == 0: - return np.array([]) - - bounds = hp.vec2dir(hp.boundaries(2**pixel_order, pixel, step=step, nest=True), lonlat=True) - margin_check = np.full((num_points,), False) - - distances = _get_boundary_point_distances(pixel_order, pixel, step) - margin_threshold_deg = margin_threshold / 3600.0 - - num_bounds = step * 4 - ra_chunks = np.array_split(r_asc, chunk_size) - dec_chunks = np.array_split(dec, chunk_size) - index = 0 - for chunk_index, ra_chunk in enumerate(ra_chunks): - dec_chunk = dec_chunks[chunk_index] - - chunk_len = len(ra_chunk) - - if chunk_len > 0: - ra_matrix = np.repeat(np.reshape([ra_chunk], (-1, 1)), num_bounds, axis=1) - dec_matrix = np.repeat(np.reshape([dec_chunk], (-1, 1)), num_bounds, axis=1) - - bounds_ra_matrix = np.repeat(np.array([bounds[0]]), len(ra_chunk), axis=0) - bounds_dec_matrix = np.repeat(np.array([bounds[1]]), len(dec_chunk), axis=0) - - points_coords = SkyCoord(ra=ra_matrix, dec=dec_matrix, unit="deg") - bounds_coords = SkyCoord(ra=bounds_ra_matrix, dec=bounds_dec_matrix, unit="deg") - - separations = points_coords.separation(bounds_coords).deg - - bisectors = np.apply_along_axis( - _find_minimum_distance, - 1, - separations, - distances=distances, - margin_threshold=margin_threshold_deg, - ) - margin_check[index : index + chunk_len] = bisectors <= margin_threshold_deg - index += chunk_len - - del ra_matrix, dec_matrix, bounds_ra_matrix, bounds_dec_matrix - del points_coords, bounds_coords, separations - - del distances - - return margin_check - - -# numba jit compiler doesn't count for coverage tests, so we'll set no cover. -@njit -def _find_minimum_distance(separations, distances, margin_threshold): # pragma: no cover - """Find the minimum distance between a given datapoint and a healpixel""" - minimum_index = np.argmin(separations) - - if separations[minimum_index] <= margin_threshold: - return separations[minimum_index] - - plus_one, minus_one = minimum_index + 1, minimum_index - 1 - num_seps = len(separations) - if minimum_index == 0: - minus_one = num_seps - 1 - if minimum_index == num_seps - 1: - plus_one = 0 - - if separations[minus_one] <= separations[plus_one]: - other_index = minus_one - else: - other_index = plus_one - - side_x = np.radians(separations[minimum_index]) - side_y = np.radians(separations[other_index]) - if (0 in (minimum_index, other_index)) and (num_seps - 1 in (minimum_index, other_index)): - side_z = distances[-1] - else: - side_z = distances[min(minimum_index, other_index)] - side_z *= np.pi / 180.0 - - ang = np.arccos((np.cos(side_y) - (np.cos(side_x) * np.cos(side_z))) / (np.sin(side_x) * np.sin(side_z))) - bisector = np.degrees(np.arcsin(np.sin(ang) * np.sin(side_x))) - - return bisector - - -def _get_boundary_point_distances(order, pixel, step): - """Get the distance of the segments between healpixel boundary points""" - boundary_points = hp.vec2dir(hp.boundaries(2**order, pixel, step=step, nest=True), lonlat=True) - - # shift forward all the coordinates by one - shift_ra = np.roll(boundary_points[0], 1) - shift_dec = np.roll(boundary_points[1], 1) - - coord_set_1 = SkyCoord(ra=boundary_points[0], dec=boundary_points[1], unit="deg") - coord_set_2 = SkyCoord(ra=shift_ra, dec=shift_dec, unit="deg") - - separations = coord_set_1.separation(coord_set_2).degree - return separations diff --git a/tests/hats/pixel_math/test_margin_bounding.py b/tests/hats/pixel_math/test_margin_bounding.py deleted file mode 100644 index 98945f2d..00000000 --- a/tests/hats/pixel_math/test_margin_bounding.py +++ /dev/null @@ -1,109 +0,0 @@ -"""Tests of margin bounding utility functions""" - -import numpy as np -import numpy.testing as npt -import pytest - -import hats.pixel_math as pm - - -def test_check_margin_bounds(): - """Make sure check_margin_bounds works properly""" - - ras = np.array([42.4704538, 56.25, 56.25, 50.61225197]) - decs = np.array([1.4593925, 9.6, 10.0, 14.4767556]) - - checks = pm.check_margin_bounds(ras, decs, 3, 4, 360.0) - - expected = np.array([False, True, False, True]) - - npt.assert_array_equal(checks, expected) - - -def test_check_margin_bounds_low_order(): - """Make sure check_margin_bounds works when using a low order healpixel""" - - ras = np.array( - [ - 142.4704538, - 45.09, - 45.2, - 37.31343956517391, - 42.649354753311535, - 32.62796809350278, - 39.89468227954832, - 27.718121934039974, - ] - ) - - decs = np.array( - [ - 1.4593925, - 0.0, - 0.0, - 6.566326903165274, - 2.005185097251452, - 10.597884275167646, - 4.465967883812584, - 14.959672304191956, - ] - ) - - checks = pm.check_margin_bounds(ras, decs, 0, 4, 360.0) - - expected = np.array([False, True, False, True, True, True, True, True]) - - npt.assert_array_equal(checks, expected) - - -def test_check_polar_margin_bounds_north(): - """Make sure check_polar_margin_bounds works at the north pole""" - order = 0 - pix = 1 - - ras = np.array([89, -179, -45]) - decs = np.array([89.9, 89.9, 85.0]) - - vals = pm.check_margin_bounds(ras, decs, order, pix, 360.0) - - expected = np.array([True, True, False]) - - npt.assert_array_equal(vals, expected) - - -def test_check_polar_margin_bounds_south(): - """Make sure check_polar_margin_bounds works at the south pole""" - order = 0 - pix = 9 - - ras = np.array([89, -179, -45]) - decs = np.array([-89.9, -89.9, -85.0]) - - vals = pm.check_margin_bounds(ras, decs, order, pix, 360.0) - - expected = np.array([True, True, False]) - - npt.assert_array_equal(vals, expected) - - -def test_check_margin_bounds_uneven(): - """Ensure check_margin_bounds fails when ra and dec arrays are unbalanced.""" - - ras = np.array([42.4704538, 56.25, 56.25, 50.61225197]) - decs = np.array([1.4593925, 9.6, 10.0]) - - with pytest.raises(ValueError, match="length of r_asc"): - pm.check_margin_bounds(ras, decs, 3, 4, 360.0) - - -def test_check_margin_bounds_empty(): - """Ensure check_margin_bounds works when passed empty coordinate arrays.""" - - ras = np.array([]) - decs = np.array([]) - - checks = pm.check_margin_bounds(ras, decs, 3, 4, 360.0) - - expected = np.array([]) - - npt.assert_array_equal(checks, expected) diff --git a/tests/hats/pixel_math/test_spatial_index.py b/tests/hats/pixel_math/test_spatial_index.py index 587482f0..ca9b04eb 100644 --- a/tests/hats/pixel_math/test_spatial_index.py +++ b/tests/hats/pixel_math/test_spatial_index.py @@ -132,11 +132,11 @@ def test_spatial_index_to_healpix_low_order(): def test_healpix_to_spatial_index_single(): orders = [3, 3, 4, 1] pixels = [0, 12, 1231, 11] - pixels_at_high_order = [p * (4 ** (SPATIAL_INDEX_ORDER - o)) for o, p in zip(orders, pixels)] - lon, lat = hp.pix2ang( - [2**SPATIAL_INDEX_ORDER] * len(orders), pixels_at_high_order, nest=True, lonlat=True - ) - actual_spatial_indices = compute_spatial_index(lon, lat) + + ra = [45.0, 45.0, 0.0, 225.0] + dec = [7.11477952e-08, 1.94712207e01, 1.44775123e01, 4.18103150e01] + + actual_spatial_indices = compute_spatial_index(ra, dec) test_spatial_indices = [healpix_to_spatial_index(o, p) for o, p in zip(orders, pixels)] assert np.all(test_spatial_indices == actual_spatial_indices) @@ -144,10 +144,9 @@ def test_healpix_to_spatial_index_single(): def test_healpix_to_spatial_index_array(): orders = [3, 3, 4, 1] pixels = [0, 12, 1231, 11] - pixels_at_high_order = [p * (4 ** (SPATIAL_INDEX_ORDER - o)) for o, p in zip(orders, pixels)] - lon, lat = hp.pix2ang( - [2**SPATIAL_INDEX_ORDER] * len(orders), pixels_at_high_order, nest=True, lonlat=True - ) - actual_spatial_indices = compute_spatial_index(lon, lat) + + ra = [45.0, 45.0, 0.0, 225.0] + dec = [7.11477952e-08, 1.94712207e01, 1.44775123e01, 4.18103150e01] + actual_spatial_indices = compute_spatial_index(ra, dec) test_spatial_indices = healpix_to_spatial_index(orders, pixels) assert np.all(test_spatial_indices == actual_spatial_indices) From 2dca8ec18b1d867abd29b1ca340fd09a19d30301 Mon Sep 17 00:00:00 2001 From: Melissa DeLucchi Date: Mon, 18 Nov 2024 11:08:24 -0500 Subject: [PATCH 2/2] Remove unused unseen. --- src/hats/pixel_math/healpix_shim.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index 44aed38b..bee05e47 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -40,10 +40,6 @@ def ring2nest(*args, **kwargs): return hp.ring2nest(*args, **kwargs) -def unseen_pixel(): - return hp.pixelfunc.UNSEEN - - ## Query