# Path: redmapper/utilities.py
# def get_healsparse_subpix_indices(subpix_nside, subpix_hpix, subpix_border, coverage_nside):
#     """
#     Retrieve the coverage pixels that intersect the region, with a border.
# 
#     Parameters
#     ----------
#     subpix_nside: `int`
#        Nside for the subregion
#     subpix_hpix: `list`
#        Pixel numbers for the subregion (ring format).
#     subpix_border: `float`
#        Border radius to cover outside subpix_hpix
#     coverage_nside: `int`
#        Nside of the healsparse coverage map
#     """
# 
#     if len(subpix_hpix) == 0:
#         raise RuntimeError("subpix_hpix cannot be an empty list.")
#     if subpix_border > 0.0 and len(subpix_hpix) > 1:
#         raise NotImplementedError("Cannot read multiple subpixels with a border.")
# 
#     # First, we need to know which pixel(s) from nside_coverage are covered by
#     # subpix_hpix
# 
#     if subpix_nside == 0:
#         # Special case for the full sky
#         return None
#     elif subpix_nside == coverage_nside:
#         covpix = hp.ring2nest(subpix_nside, subpix_hpix)
#     elif subpix_nside > coverage_nside:
#         # what pixels are these contained in?
#         theta, phi = hp.pix2ang(subpix_nside, subpix_hpix, nest=False)
#         covpix = hp.ang2pix(coverage_nside, theta, phi, nest=True)
#     else:
#         # This is subpix_nside < coverage_nside
#         # what coverage pixels are contained in subpix_hpix?
#         subpix_hpix_nest = hp.ring2nest(subpix_nside, subpix_hpix)
#         bit_shift = 2 * int(np.round(np.log(coverage_nside / subpix_nside) / np.log(2)))
#         n_pix_per_pix = 2**bit_shift
#         covpix = np.zeros(len(subpix_hpix_nest) * n_pix_per_pix, dtype=np.int64)
#         for i, hpix in enumerate(subpix_hpix_nest):
#             covpix[i * n_pix_per_pix: (i + 1) * n_pix_per_pix] = np.left_shift(hpix, bit_shift) + np.arange(n_pix_per_pix)
# 
#     # And now if we have a border...(only if the length is 1)
#     if subpix_border > 0.0:
#         nside_testing = max([coverage_nside * 4, subpix_nside * 4])
#         boundaries = hp.boundaries(subpix_nside, subpix_hpix[0], step=nside_testing/subpix_nside)
# 
#         extrapix = np.zeros(0, dtype=np.int64)
# 
#         # These are pixels that touch the boundary
#         for i in range(boundaries.shape[1]):
#             pixint = hp.query_disc(nside_testing, boundaries[:, i],
#                                    np.radians(subpix_border), inclusive=True, fact=8)
#             extrapix = np.append(extrapix, pixint)
# 
#         extrapix = np.unique(extrapix)
#         theta, phi = hp.pix2ang(nside_testing, extrapix)
#         covpix = np.unique(np.append(covpix, hp.ang2pix(coverage_nside, theta, phi, nest=True)))
# 
#     return covpix


import unittest
import numpy.testing as testing
import numpy as np
import healpy as hp
import os
import esutil
from redmapper.utilities import get_healsparse_subpix_indices

        covpix = get_healsparse_subpix_indices(subpix_nside, subpix_hpix, subpix_border, coverage_nside)

        # Test that we have the right number
        testing.assert_equal(covpix.size, int(hp.nside2npix(coverage_nside) / hp.nside2npix(subpix_nside)))

        # Test that they are all in the correct pixel
        theta, phi = hp.pix2ang(coverage_nside, covpix, nest=True)
        ipring = hp.ang2pix(subpix_nside, theta, phi)

        testing.assert_equal(ipring, subpix_hpix[0])

        covpix_temp = covpix

        # Test with two pixels
        subpix_hpix = [16, 17]

        covpix = get_healsparse_subpix_indices(subpix_nside, subpix_hpix, subpix_border, coverage_nside)

        # Test that they are all in the correct pixel
        theta, phi = hp.pix2ang(coverage_nside, covpix, nest=True)
        ipring = hp.ang2pix(subpix_nside, theta, phi)

        testing.assert_equal(ipring[0: 256], subpix_hpix[0])
        testing.assert_equal(ipring[256: ], subpix_hpix[1])

        # Test with border (most typical value)...
        subpix_hpix = [16]
        subpix_border = 2.5


# Path: redmapper/utilities.py
# def get_healsparse_subpix_indices(subpix_nside, subpix_hpix, subpix_border, coverage_nside):
#     """
#     Retrieve the coverage pixels that intersect the region, with a border.
# 
#     Parameters
#     ----------
#     subpix_nside: `int`
#        Nside for the subregion
#     subpix_hpix: `list`
#        Pixel numbers for the subregion (ring format).
#     subpix_border: `float`
#        Border radius to cover outside subpix_hpix
#     coverage_nside: `int`
#        Nside of the healsparse coverage map
#     """
# 
#     if len(subpix_hpix) == 0:
#         raise RuntimeError("subpix_hpix cannot be an empty list.")
#     if subpix_border > 0.0 and len(subpix_hpix) > 1:
#         raise NotImplementedError("Cannot read multiple subpixels with a border.")
# 
#     # First, we need to know which pixel(s) from nside_coverage are covered by
#     # subpix_hpix
# 
#     if subpix_nside == 0:
#         # Special case for the full sky
#         return None
#     elif subpix_nside == coverage_nside:
#         covpix = hp.ring2nest(subpix_nside, subpix_hpix)
#     elif subpix_nside > coverage_nside:
#         # what pixels are these contained in?
#         theta, phi = hp.pix2ang(subpix_nside, subpix_hpix, nest=False)
#         covpix = hp.ang2pix(coverage_nside, theta, phi, nest=True)
#     else:
#         # This is subpix_nside < coverage_nside
#         # what coverage pixels are contained in subpix_hpix?
#         subpix_hpix_nest = hp.ring2nest(subpix_nside, subpix_hpix)
#         bit_shift = 2 * int(np.round(np.log(coverage_nside / subpix_nside) / np.log(2)))
#         n_pix_per_pix = 2**bit_shift
#         covpix = np.zeros(len(subpix_hpix_nest) * n_pix_per_pix, dtype=np.int64)
#         for i, hpix in enumerate(subpix_hpix_nest):
#             covpix[i * n_pix_per_pix: (i + 1) * n_pix_per_pix] = np.left_shift(hpix, bit_shift) + np.arange(n_pix_per_pix)
# 
#     # And now if we have a border...(only if the length is 1)
#     if subpix_border > 0.0:
#         nside_testing = max([coverage_nside * 4, subpix_nside * 4])
#         boundaries = hp.boundaries(subpix_nside, subpix_hpix[0], step=nside_testing/subpix_nside)
# 
#         extrapix = np.zeros(0, dtype=np.int64)
# 
#         # These are pixels that touch the boundary
#         for i in range(boundaries.shape[1]):
#             pixint = hp.query_disc(nside_testing, boundaries[:, i],
#                                    np.radians(subpix_border), inclusive=True, fact=8)
#             extrapix = np.append(extrapix, pixint)
# 
#         extrapix = np.unique(extrapix)
#         theta, phi = hp.pix2ang(nside_testing, extrapix)
#         covpix = np.unique(np.append(covpix, hp.ang2pix(coverage_nside, theta, phi, nest=True)))
# 
#     return covpix


# Path: tests/test_subpix_indices.py
import unittest
import numpy.testing as testing
import numpy as np
import healpy as hp
import os
import esutil
from redmapper.utilities import get_healsparse_subpix_indices


        covpix = get_healsparse_subpix_indices(subpix_nside, subpix_hpix, subpix_border, coverage_nside)

        # Test that we have the right number
        testing.assert_equal(covpix.size, int(hp.nside2npix(coverage_nside) / hp.nside2npix(subpix_nside)))

        # Test that they are all in the correct pixel
        theta, phi = hp.pix2ang(coverage_nside, covpix, nest=True)
        ipring = hp.ang2pix(subpix_nside, theta, phi)

        testing.assert_equal(ipring, subpix_hpix[0])

        covpix_temp = covpix

        # Test with two pixels
        subpix_hpix = [16, 17]

        covpix = get_healsparse_subpix_indices(subpix_nside, subpix_hpix, subpix_border, coverage_nside)

        # Test that they are all in the correct pixel
        theta, phi = hp.pix2ang(coverage_nside, covpix, nest=True)
        ipring = hp.ang2pix(subpix_nside, theta, phi)

        testing.assert_equal(ipring[0: 256], subpix_hpix[0])
        testing.assert_equal(ipring[256: ], subpix_hpix[1])

        # Test with border (most typical value)...
        subpix_hpix = [16]
        subpix_border = 2.5


        covpix2 = get_healsparse_subpix_indices(subpix_nside, subpix_hpix, subpix_border, coverage_nside)