HEALPix Trees (`ligo.skymap.healpix_tree`)¶

Multiresolution HEALPix trees

class ligo.skymap.healpix_tree.HEALPixTree(samples, max_samples_per_pixel, max_order, order=0, needs_sort=True)[source]¶

Bases: object

Data structure used internally by the function adaptive_healpix_histogram().

property flat_bitmap¶: Return flattened HEALPix representation.

static key_for_order(order)[source]¶: Create a function that downsamples full-resolution pixel indices.

property order¶: Return the maximum HEALPix order required to represent this tree, which is the same as the tree depth.

visit(order='depthfirst', extra=True)[source]¶

Traverse the leaves of the HEALPix tree.

Parameters:

orderstring, optional: Traversal order: ‘depthfirst’ (the default) or ‘breadthfirst’.
extrabool: Whether to output extra information about the pixel (default is True).

Yields:

nsideint: The HEALPix resolution of the node.
full_nsideint, present if extra=True: The HEALPix resolution of the deepest node in the tree.
ipixint: The nested HEALPix index of the node.
ipix0int, present if extra=True: The start index of the range of pixels spanned by the node at the resolution full_nside.
ipix1int, present if extra=True: The end index of the range of pixels spanned by the node at the resolution full_nside.
sampleslist, present if extra=True: The list of samples contained in the node.

Examples

>>> ipix = np.arange(12) * HEALPIX_MACHINE_NSIDE**2
>>> tree = HEALPixTree(ipix, max_samples_per_pixel=1, max_order=1)
>>> [tuple(_) for _ in tree.visit(extra=False)]
[(1, 0), (1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (1, 7), (1, 8), (1, 9), (1, 10), (1, 11)]

ligo.skymap.healpix_tree.adaptive_healpix_histogram(theta, phi, max_samples_per_pixel, nside=-1, max_nside=-1, nest=False)[source]¶: Adaptively histogram the posterior samples represented by the (theta, phi) points using a recursively subdivided HEALPix tree. Nodes are subdivided until each leaf contains no more than max_samples_per_pixel samples. Finally, the tree is flattened to a fixed-resolution HEALPix image with a resolution appropriate for the depth of the tree. If nside is specified, the result is resampled to another desired HEALPix resolution.

ligo.skymap.healpix_tree.interpolate_nested(m, nest=False)[source]¶

Apply bilinear interpolation to a multiresolution HEALPix map, assuming that runs of pixels containing identical values are nodes of the tree. This smooths out the stair-step effect that may be noticeable in contour plots.

Here is how it works. Consider a coarse tile surrounded by base tiles, like this:

        +---+---+
        |   |   |
        +-------+
        |   |   |
+---+---+---+---+---+---+
|   |   |       |   |   |
+-------+       +-------+
|   |   |       |   |   |
+---+---+---+---+---+---+
        |   |   |
        +-------+
        |   |   |
        +---+---+

The value within the central coarse tile is computed by downsampling the sky map (averaging the fine tiles), upsampling again (with bilinear interpolation), and then finally copying the interpolated values within the coarse tile back to the full-resolution sky map. This process is applied recursively at all successive HEALPix resolutions.

Note that this method suffers from a minor discontinuity artifact at the edges of regions of coarse tiles, because it temporarily treats the bordering fine tiles as constant. However, this artifact seems to have only a minor effect on generating contour plots.

Parameters:

m: `~numpy.ndarray`: a HEALPix array
nest: bool, default: False: Whether the input array is stored in the NESTED indexing scheme (True) or the RING indexing scheme (False).

ligo.skymap.healpix_tree.reconstruct_nested(m, order='depthfirst', extra=True)[source]¶

Reconstruct the leaves of a multiresolution tree.

Parameters:

mndarray: A HEALPix array in the NESTED ordering scheme.
order{‘depthfirst’, ‘breadthfirst’}, optional: Traversal order: ‘depthfirst’ (the default) or ‘breadthfirst’.
extrabool: Whether to output extra information about the pixel (default is True).

Yields:

nsideint: The HEALPix resolution of the node.
full_nsideint, present if extra=True: The HEALPix resolution of the deepest node in the tree.
ipixint: The nested HEALPix index of the node.
ipix0int, present if extra=True: The start index of the range of pixels spanned by the node at the resolution full_nside.
ipix1int, present if extra=True: The end index of the range of pixels spanned by the node at the resolution full_nside.
valuelist, present if extra=True: The value of the map at the node.

Examples

An nside=1 array of all zeros:

>>> m = np.zeros(12)
>>> result = reconstruct_nested(m, order='breadthfirst', extra=False)
>>> np.asarray(list(result))
array([[ 1,  0],
       [ 1,  1],
       [ 1,  2],
       [ 1,  3],
       [ 1,  4],
       [ 1,  5],
       [ 1,  6],
       [ 1,  7],
       [ 1,  8],
       [ 1,  9],
       [ 1, 10],
       [ 1, 11]])

An nside=1 array of distinct values:

>>> m = range(12)
>>> result = reconstruct_nested(m, order='breadthfirst', extra=False)
>>> np.asarray(list(result))
array([[ 1,  0],
       [ 1,  1],
       [ 1,  2],
       [ 1,  3],
       [ 1,  4],
       [ 1,  5],
       [ 1,  6],
       [ 1,  7],
       [ 1,  8],
       [ 1,  9],
       [ 1, 10],
       [ 1, 11]])

An nside=8 array of zeros:

>>> m = np.zeros(768)
>>> result = reconstruct_nested(m, order='breadthfirst', extra=False)
>>> np.asarray(list(result))
array([[ 1,  0],
       [ 1,  1],
       [ 1,  2],
       [ 1,  3],
       [ 1,  4],
       [ 1,  5],
       [ 1,  6],
       [ 1,  7],
       [ 1,  8],
       [ 1,  9],
       [ 1, 10],
       [ 1, 11]])

An nside=2 array, all zeros except for four consecutive distinct elements:

>>> m = np.zeros(48); m[:4] = range(4)
>>> result = reconstruct_nested(m, order='breadthfirst', extra=False)
>>> np.asarray(list(result))
array([[ 1,  1],
       [ 1,  2],
       [ 1,  3],
       [ 1,  4],
       [ 1,  5],
       [ 1,  6],
       [ 1,  7],
       [ 1,  8],
       [ 1,  9],
       [ 1, 10],
       [ 1, 11],
       [ 2,  0],
       [ 2,  1],
       [ 2,  2],
       [ 2,  3]])

Same, but in depthfirst order:

>>> result = reconstruct_nested(m, order='depthfirst', extra=False)
>>> np.asarray(list(result))
array([[ 2,  0],
       [ 2,  1],
       [ 2,  2],
       [ 2,  3],
       [ 1,  1],
       [ 1,  2],
       [ 1,  3],
       [ 1,  4],
       [ 1,  5],
       [ 1,  6],
       [ 1,  7],
       [ 1,  8],
       [ 1,  9],
       [ 1, 10],
       [ 1, 11]])

An nside=2 array, all elements distinct except for four consecutive zeros:

>>> m = np.arange(48); m[:4] = 0
>>> result = reconstruct_nested(m, order='breadthfirst', extra=False)
>>> np.asarray(list(result))
array([[ 1,  0],
       [ 2,  4],
       [ 2,  5],
       [ 2,  6],
       [ 2,  7],
       [ 2,  8],
       [ 2,  9],
       [ 2, 10],
       [ 2, 11],
       [ 2, 12],
       [ 2, 13],
       [ 2, 14],
       [ 2, 15],
       [ 2, 16],
       [ 2, 17],
       [ 2, 18],
       [ 2, 19],
       [ 2, 20],
       [ 2, 21],
       [ 2, 22],
       [ 2, 23],
       [ 2, 24],
       [ 2, 25],
       [ 2, 26],
       [ 2, 27],
       [ 2, 28],
       [ 2, 29],
       [ 2, 30],
       [ 2, 31],
       [ 2, 32],
       [ 2, 33],
       [ 2, 34],
       [ 2, 35],
       [ 2, 36],
       [ 2, 37],
       [ 2, 38],
       [ 2, 39],
       [ 2, 40],
       [ 2, 41],
       [ 2, 42],
       [ 2, 43],
       [ 2, 44],
       [ 2, 45],
       [ 2, 46],
       [ 2, 47]])

Navigation

HEALPix Trees (ligo.skymap.healpix_tree)¶

HEALPix Trees (`ligo.skymap.healpix_tree`)¶