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] [edit on github]

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] [edit on github]

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] [edit on github]

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] [edit on github]

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] [edit on github]

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] [edit on github]

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)
>>> [tuple(_) for _ in result]
[(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)
>>> [tuple(_) for _ in result]
[(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)
>>> [tuple(_) for _ in result]
[(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)
>>> [tuple(_) for _ in result]
[(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)
>>> [tuple(_) for _ in result]
[(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)
>>> [tuple(_) for _ in result]
[(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)]