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

Add option to rasterize to a custom order

parent 9f589d78
......@@ -12,6 +12,11 @@ Changelog
- The ``ligo-skymap-plot-airmass`` tool will now use the color map's full
dynamic range.
- Add ``order`` option to ``ligo.skymap.moc.rasterize`` and
``ligo.skymap.bayestar.rasterize`` and ``--nside`` option to
``ligo-skymap-flatten`` to support flattening multi-resolution HEALPix
datasets to specified resolutions.
0.1.6 (2019-03-26)
==================
......
......@@ -413,8 +413,38 @@ def localize(
return skymap
def rasterize(skymap):
skymap = Table(moc.rasterize(skymap), meta=skymap.meta)
def rasterize(skymap, order=None):
orig_order, _ = moc.uniq2nest(skymap['UNIQ'])
orig_order = orig_order.max()
if order is None or order < 0 or order >= orig_order \
or 'DISTMU' not in skymap.dtype.fields.keys():
skymap = Table(moc.rasterize(skymap, order=order), meta=skymap.meta)
else:
skymap = Table(skymap, copy=True, meta=skymap.meta)
probdensity = skymap['PROBDENSITY']
distmu = skymap.columns.pop('DISTMU')
distsigma = skymap.columns.pop('DISTSIGMA')
bad = ~(np.isfinite(distmu) & np.isfinite(distsigma))
distmean, diststd, _ = distance.parameters_to_moments(
distmu, distsigma)
distmean[bad] = np.nan
diststd[bad] = np.nan
skymap['DISTMEAN'] = probdensity * distmean
skymap['DISTVAR'] = probdensity * (
np.square(diststd) + np.square(distmean))
skymap = Table(moc.rasterize(skymap, order=order), meta=skymap.meta)
distmean = skymap.columns.pop('DISTMEAN') / skymap['PROBDENSITY']
diststd = np.sqrt(
skymap.columns.pop('DISTVAR') / skymap['PROBDENSITY']
- np.square(distmean))
skymap['DISTMU'], skymap['DISTSIGMA'], skymap['DISTNORM'] = \
distance.moments_to_parameters(
distmean, diststd)
skymap.rename_column('PROBDENSITY', 'PROB')
skymap['PROB'] *= 4 * np.pi / len(skymap)
skymap['PROB'].unit = u.pixel ** -1
......
......@@ -28,6 +28,8 @@ References
Recommendation <http://ivoa.net/documents/MOC/>.
"""
from astropy import table
import numpy as np
from .core import nest2uniq, uniq2nest, uniq2order, uniq2pixarea, uniq2ang
from .core import rasterize as _rasterize
......@@ -101,7 +103,7 @@ ipix : `numpy.ndarray`
""")
def rasterize(moc_data):
def rasterize(moc_data, order=None):
"""Convert a multi-order HEALPix dataset to fixed-order NESTED ordering.
Parameters
......@@ -111,6 +113,9 @@ def rasterize(moc_data):
first column is called UNIQ and contains the NUNIQ pixel index. Every
point on the unit sphere must be contained in exactly one pixel in the
dataset.
order : int, optional
The desired output resolution order, or :obj:`None` for the maximum
resolution present in the dataset.
Returns
-------
......@@ -118,7 +123,29 @@ def rasterize(moc_data):
A fixed-order, NESTED-ordering HEALPix dataset with all of the columns
that were in moc_data, with the exception of the UNIQ column.
"""
return _rasterize(moc_data)
if order is None or order < 0:
order = -1
else:
orig_order, orig_nest = uniq2nest(moc_data['UNIQ'])
to_downsample = order < orig_order
if np.any(to_downsample):
to_keep = table.Table(moc_data[~to_downsample])
orig_order = orig_order[to_downsample]
orig_nest = orig_nest[to_downsample]
to_downsample = table.Table(moc_data[to_downsample])
ratio = 1 << (2 * np.uint64(orig_order - order))
weights = 1.0 / ratio
for colname, column in to_downsample.columns.items():
if colname != 'UNIQ':
column *= weights
to_downsample['UNIQ'] = nest2uniq(order, orig_nest // ratio)
to_downsample = to_downsample.group_by(
'UNIQ').groups.aggregate(np.sum)
moc_data = table.vstack((to_keep, to_downsample))
return _rasterize(moc_data, order=order)
del add_newdoc_ufunc
from astropy import table
import healpy as hp
import numpy as np
import pytest
from .. import moc
def input_skymap(order1, d_order, fraction):
"""Construct a test multi-resolution sky map, with values that are
proportional to the NESTED pixel index.
To make the test more interesting by mixing together multiple resolutions,
part of the sky map is refined to a higher order.
Parameters
----------
order1 : int
The HEALPix resolution order.
d_order : int
The increase in orer for part of the sky map.
fraction : float
The fraction of the original pixels to refine.
"""
order2 = order1 + d_order
npix1 = hp.nside2npix(hp.order2nside(order1))
npix2 = hp.nside2npix(hp.order2nside(order2))
ipix1 = np.arange(npix1)
ipix2 = np.arange(npix2)
data1 = table.Table({
'UNIQ': moc.nest2uniq(order1, ipix1.astype(np.uint64)),
'VALUE': ipix1.astype(float),
'VALUE2': np.pi * ipix1.astype(float)
})
data2 = table.Table({
'UNIQ': moc.nest2uniq(order2, ipix2.astype(np.uint64)),
'VALUE': np.repeat(ipix1, npix2 // npix1).astype(float),
'VALUE2': np.pi * np.repeat(ipix1, npix2 // npix1).astype(float)
})
n = int(npix1 * (1 - fraction))
return table.vstack((data1[:n], data2[n * npix2 // npix1:]))
@pytest.mark.parametrize('order_in', [6])
@pytest.mark.parametrize('d_order_in', range(3))
@pytest.mark.parametrize('fraction_in', [0, 0.25, 0.5, 1])
@pytest.mark.parametrize('order_out', range(6))
def test_rasterize_downsample(order_in, d_order_in, fraction_in, order_out):
npix_in = hp.nside2npix(hp.order2nside(order_in))
npix_out = hp.nside2npix(hp.order2nside(order_out))
skymap_in = input_skymap(order_in, d_order_in, fraction_in)
skymap_out = moc.rasterize(skymap_in, order_out)
assert len(skymap_out) == npix_out
reps = npix_in // npix_out
expected = np.mean(np.arange(npix_in).reshape(-1, reps), axis=1)
np.testing.assert_array_equal(skymap_out['VALUE'], expected)
@pytest.mark.parametrize('order_in', [2])
@pytest.mark.parametrize('d_order_in', range(3))
@pytest.mark.parametrize('fraction_in', [0, 0.25, 0.5, 1])
@pytest.mark.parametrize('order_out', range(3, 9))
def test_rasterize_upsample(order_in, d_order_in, fraction_in, order_out):
npix_in = hp.nside2npix(hp.order2nside(order_in))
npix_out = hp.nside2npix(hp.order2nside(order_out))
skymap_in = input_skymap(order_in, d_order_in, fraction_in)
skymap_out = moc.rasterize(skymap_in, order_out)
assert len(skymap_out) == npix_out
ipix = np.arange(npix_in)
reps = npix_out // npix_in
for i in range(reps):
np.testing.assert_array_equal(skymap_out['VALUE'][i::reps], ipix)
@pytest.mark.parametrize('order', range(3))
def test_rasterize_default(order):
npix = hp.nside2npix(hp.order2nside(order))
skymap_in = input_skymap(order, 0, 0)
skymap_out = moc.rasterize(skymap_in)
assert len(skymap_out) == npix
......@@ -27,6 +27,7 @@ def parser():
type=FileType('rb'), help='Input FITS file')
parser.add_argument('output', metavar='OUTPUT.fits',
type=FileType('wb'), help='Output FITS file')
parser.add_argument('--nside', type=int, help='Output HEALPix resolution')
return parser
......@@ -35,9 +36,15 @@ def main(args=None):
import warnings
from astropy.io import fits
import healpy as hp
from ..io import read_sky_map, write_sky_map
from ..bayestar import rasterize
if args.nside is None:
order = None
else:
order = hp.nside2order(args.nside)
hdus = fits.open(args.input)
ordering = hdus[1].header['ORDERING']
expected_ordering = 'NUNIQ'
......@@ -45,4 +52,4 @@ def main(args=None):
msg = 'Expected the FITS file {} to have ordering {}, but it is {}'
warnings.warn(msg.format(args.input.name, expected_ordering, ordering))
table = read_sky_map(hdus, moc=True)
write_sky_map(args.output.name, rasterize(table), nest=True)
write_sky_map(args.output.name, rasterize(table, order=order), nest=True)
from astropy import table
import healpy as hp
import numpy as np
import pytest
from . import run_entry_point
from ...distance import moments_to_parameters, parameters_to_marginal_moments
from ... import moc
from ...io.fits import read_sky_map, write_sky_map
def input_skymap(order1, d_order, fraction):
"""Construct a test multi-resolution sky map, with values that are
proportional to the NESTED pixel index.
To make the test more interesting by mixing together multiple resolutions,
part of the sky map is refined to a higher order.
Parameters
----------
order1 : int
The HEALPix resolution order.
d_order : int
The increase in orer for part of the sky map.
fraction : float
The fraction of the original pixels to refine.
"""
order2 = order1 + d_order
npix1 = hp.nside2npix(hp.order2nside(order1))
npix2 = hp.nside2npix(hp.order2nside(order2))
ipix1 = np.arange(npix1)
ipix2 = np.arange(npix2)
# Create a random sky map.
area = hp.nside2pixarea(hp.order2nside(order1))
probdensity = np.random.uniform(0, 1, npix1)
prob = probdensity * area
normalization = prob.sum()
prob /= normalization
probdensity /= normalization
distmean = np.random.uniform(100, 110, npix1)
diststd = np.random.uniform(0, 1 / np.sqrt(3) - 0.1, npix1) * distmean
distmu, distsigma, distnorm = moments_to_parameters(distmean, diststd)
assert np.all(np.isfinite(distmu))
data1 = table.Table({
'UNIQ': moc.nest2uniq(order1, ipix1.astype(np.uint64)),
'PROBDENSITY': probdensity,
'DISTMU': distmu,
'DISTSIGMA': distsigma,
'DISTNORM': distnorm
})
# Add some upsampled pixels.
data2 = table.Table(np.repeat(data1, npix2 // npix1))
data2['UNIQ'] = moc.nest2uniq(order2, ipix2.astype(np.uint64))
n = int(npix1 * (1 - fraction))
result = table.vstack((data1[:n], data2[n * npix2 // npix1:]))
# Add marginal distance mean and standard deviation.
rbar = (prob * distmean).sum()
r2bar = (prob * (np.square(diststd) + np.square(distmean))).sum()
result.meta['distmean'] = rbar
result.meta['diststd'] = np.sqrt(r2bar - np.square(rbar))
return result
@pytest.mark.parametrize('order_in', [2])
@pytest.mark.parametrize('d_order_in', range(3))
@pytest.mark.parametrize('fraction_in', [0, 0.25, 0.5, 1])
@pytest.mark.parametrize('nside_out', [None, 1, 2, 4, 8, 512])
def test_flatten(tmpdir, order_in, d_order_in, fraction_in, nside_out):
"""Test ligo-skyamp-flatten."""
input_filename = str(tmpdir / 'bayestar.fits')
output_filename = str(tmpdir / 'bayestar.fits.gz')
skymap = input_skymap(order_in, d_order_in, fraction_in)
write_sky_map(input_filename, skymap, moc=True)
expected_distmean = skymap.meta['distmean']
expected_diststd = skymap.meta['diststd']
args = ['ligo-skymap-flatten', input_filename, output_filename]
if nside_out is not None:
args.extend(['--nside', str(nside_out)])
run_entry_point(*args)
(prob, distmu, distsigma, distnorm), _ = read_sky_map(
output_filename, distances=True)
distmean, diststd = parameters_to_marginal_moments(
prob, distmu, distsigma)
if nside_out is not None:
assert len(prob) == hp.nside2npix(nside_out)
assert prob.sum() == pytest.approx(1)
assert distmean == pytest.approx(expected_distmean)
assert diststd == pytest.approx(expected_diststd)
......@@ -76,15 +76,18 @@ void uniq2ang64(uint64_t uniq, double *theta, double *phi)
}
void *moc_rasterize64(const void *pixels, size_t offset, size_t itemsize, size_t len, size_t *npix)
void *moc_rasterize64(
const void *pixels, size_t offset, size_t itemsize, size_t len,
size_t *npix, int8_t order)
{
/* Calculate pixel size. */
const size_t pixelsize = offset + itemsize;
/* Find maximum order. Note: normally MOC datasets are stored in order
* of ascending MOC index, so the last pixel should have the highest order.
* However, our rasterization algorithm doesn't depend on this sorting,
* so let's just do a linear search for the maximum order. */
/* If the parameter order >= 0, then rasterize at that order.
* Otherwise, find maximum order. Note: normally MOC datasets are stored in
* order of ascending MOC index, so the last pixel should have the highest
* order. However, our rasterization algorithm doesn't depend on this
* sorting, so let's just do a linear search for the maximum order. */
int8_t max_order;
{
uint64_t max_uniq = 0;
......@@ -98,6 +101,14 @@ void *moc_rasterize64(const void *pixels, size_t offset, size_t itemsize, size_t
max_order = uniq2order64(max_uniq);
}
/* Don't handle downsampling here, because we don't know how to do
* reduction across pixels without more knowledge of the pixel datatype and
* contents. */
if (order >= max_order)
max_order = order;
else if (order >= 0)
GSL_ERROR_NULL("downsampling not implemented", GSL_EUNIMPL);
/* Allocate output. */
*npix = 12 * ((size_t) 1 << 2 * max_order);
void *ret = calloc(*npix, itemsize);
......@@ -109,7 +120,7 @@ void *moc_rasterize64(const void *pixels, size_t offset, size_t itemsize, size_t
{
const void *pixel = (const char *) pixels + i * pixelsize;
uint64_t nest;
const int8_t order = uniq2nest64(*(const uint64_t *) pixel, &nest);
order = uniq2nest64(*(const uint64_t *) pixel, &nest);
const size_t reps = (size_t) 1 << 2 * (max_order - order);
for (size_t j = 0; j < reps; j ++)
memcpy((char *) ret + (nest * reps + j) * itemsize,
......
......@@ -40,7 +40,7 @@ int8_t uniq2nest64(uint64_t uniq, uint64_t *nest);
void uniq2ang64(uint64_t uniq, double *theta, double *phi);
void *moc_rasterize64(const void *pixels, size_t offset, size_t itemsize, size_t len, size_t *npix);
void *moc_rasterize64(const void *pixels, size_t offset, size_t itemsize, size_t len, size_t *npix, int8_t order);
#endif /* __cplusplus */
......
......@@ -409,8 +409,24 @@ static void capsule_free(PyObject *self)
}
static PyObject *rasterize(PyObject *NPY_UNUSED(module), PyObject *arg)
static PyObject *rasterize(
PyObject *NPY_UNUSED(module), PyObject *args, PyObject *kwargs)
{
PyObject *arg;
int order = -1;
/* Names of arguments */
static const char *keywords[] = {"array", "order", NULL};
/* Parse arguments */
/* FIXME: PyArg_ParseTupleAndKeywords should expect keywords to be const */
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wincompatible-pointer-types"
if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O|i",
keywords, &arg, &order))
return NULL;
#pragma GCC diagnostic pop
PyArrayObject *arr = (PyArrayObject *) PyArray_FromAny(
arg, NULL, 1, 1, NPY_ARRAY_CARRAY_RO, NULL);
PyObject *uniq_key = PyUnicode_FromString("UNIQ");
......@@ -489,7 +505,7 @@ static PyObject *rasterize(PyObject *NPY_UNUSED(module), PyObject *arg)
void *out;
Py_BEGIN_ALLOW_THREADS
out = moc_rasterize64(pixels, offset, itemsize, len, &npix);
out = moc_rasterize64(pixels, offset, itemsize, len, &npix, order);
Py_END_ALLOW_THREADS
if (!out)
goto done;
......@@ -1022,7 +1038,8 @@ static void *const no_ufunc_data[] = {NULL};
static const char modulename[] = "core";
static PyMethodDef methods[] = {
{"rasterize", rasterize, METH_O, "fill me in"},
{"rasterize", (PyCFunction)rasterize,
METH_VARARGS | METH_KEYWORDS, "fill me in"},
{"toa_phoa_snr", (PyCFunction)sky_map_toa_phoa_snr,
METH_VARARGS | METH_KEYWORDS, "fill me in"},
{"test", (PyCFunction)test,
......
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