Commit a3d23060 authored by Leo Pound Singer's avatar Leo Pound Singer Committed by Adam Mercer
Browse files

Port the BAYESTAR signal amplitude model unit test to Python

This unit test constituted most of the calls to other LALSuite
functions in C.
parent 25afcc41
......@@ -23,6 +23,7 @@
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wcast-qual"
#include <numpy/arrayobject.h>
#include <numpy/ufuncobject.h>
#pragma GCC diagnostic pop
#include <chealpix.h>
#include <gsl/gsl_errno.h>
......@@ -355,6 +356,27 @@ fail: /* Cleanup */
};
static void signal_amplitude_model_loop(
char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(data))
{
const npy_intp n = dimensions[0];
for (npy_intp i = 0; i < n; i ++)
{
/* FIXME: args must be void ** to avoid alignment warnings */
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wcast-align"
*(double complex *) &args[4][i * steps[4]] =
bayestar_signal_amplitude_model(
*(double complex *) &args[0][i * steps[0]],
*(double complex *) &args[1][i * steps[1]],
*(double *) &args[2][i * steps[2]],
*(double *) &args[3][i * steps[3]]);
#pragma GCC diagnostic pop
}
}
static PyObject *test(
PyObject *NPY_UNUSED(module), PyObject *NPY_UNUSED(arg))
{
......@@ -386,6 +408,15 @@ static PyModuleDef moduledef = {
};
static const PyUFuncGenericFunction
signal_amplitude_model_loops[] = {signal_amplitude_model_loop};
static const char signal_amplitude_model_ufunc_types[] = {
NPY_CDOUBLE, NPY_CDOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_CDOUBLE};
static void *const no_ufunc_data[] = {NULL};
PyMODINIT_FUNC PyInit__sky_map(void); /* Silence -Wmissing-prototypes */
PyMODINIT_FUNC PyInit__sky_map(void)
{
......@@ -393,6 +424,7 @@ PyMODINIT_FUNC PyInit__sky_map(void)
gsl_set_error_handler_off();
import_array();
import_umath();
sky_map_descr = sky_map_create_descr();
if (!sky_map_descr)
......@@ -402,6 +434,22 @@ PyMODINIT_FUNC PyInit__sky_map(void)
if (!module)
goto done;
/* Ignore warnings in Numpy API */
#pragma GCC diagnostic push
#ifdef __clang__
#pragma GCC diagnostic ignored "-Wincompatible-pointer-types-discards-qualifiers"
#else
#pragma GCC diagnostic ignored "-Wdiscarded-qualifiers"
#endif
PyModule_AddObject(
module, "signal_amplitude_model", PyUFunc_FromFuncAndData(
signal_amplitude_model_loops, no_ufunc_data,
signal_amplitude_model_ufunc_types, 1, 4, 1, PyUFunc_None,
"signal_amplitude_model", NULL, 0));
#pragma GCC diagnostic pop
done:
return module;
}
......
......@@ -81,14 +81,11 @@
#include <lal/cubic_interp.h>
#include <lal/DetResponse.h>
#include <lal/InspiralInjectionParams.h>
#include <lal/FrequencySeries.h>
#include <lal/LALSimInspiral.h>
#include <lal/LALSimulation.h>
#include <chealpix.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_math.h>
......@@ -534,8 +531,8 @@ static double complex complex_antenna_factor(
/* Expression for complex amplitude on arrival (without 1/distance factor) */
static double complex signal_amplitude_model(
double complex F, /* Antenna factor */
double complex bayestar_signal_amplitude_model(
double complex F, /* Complex antenna factor */
double complex exp_i_twopsi, /* e^(i*2*psi), for polarization angle psi */
double u, /* cos(inclination) */
double u2 /* cos^2(inclination */
......@@ -768,7 +765,7 @@ static void bayestar_sky_map_toa_phoa_snr_pixel(
for (unsigned int iifo = 0; iifo < nifos; iifo ++)
{
p2 += cabs2(
z_times_r[iifo] = signal_amplitude_model(
z_times_r[iifo] = bayestar_signal_amplitude_model(
F[iifo], exp_i_twopsi, u, u2));
}
p2 *= 0.5;
......@@ -1076,7 +1073,7 @@ double bayestar_log_likelihood_toa_phoa_snr(
responses[iifo], ra, dec, gmst) * horizons[iifo];
const double complex z_times_r =
signal_amplitude_model(F, exp_i_twopsi, u, u2);
bayestar_signal_amplitude_model(F, exp_i_twopsi, u, u2);
i0arg_complex_times_r += conj(z_times_r)
* eval_snr(snrs[iifo], nsamples, (t - dt[iifo]) * sample_rate - 0.5 * (nsamples - 1));
......@@ -1183,154 +1180,6 @@ static void test_eval_snr(void)
}
static void test_signal_amplitude_model(
double ra,
double dec,
double inclination,
double polarization,
unsigned long epoch_nanos,
const char *instrument
) {
const LALDetector *detector = XLALDetectorPrefixToLALDetector(
instrument);
LIGOTimeGPS epoch;
double gmst;
XLALINT8NSToGPS(&epoch, epoch_nanos);
gmst = XLALGreenwichMeanSiderealTime(&epoch);
double complex result;
{
double complex F;
const double complex exp_i_twopsi = exp_i(2 * polarization);
const double u = cos(inclination);
const double u2 = gsl_pow_2(u);
XLALComputeDetAMResponse(
(double *)&F, /* Type-punned real part */
1 + (double *)&F, /* Type-punned imag part */
detector->response, ra, dec, 0, gmst);
result = signal_amplitude_model(F, exp_i_twopsi, u, u2);
}
float abs_expected;
{
SimInspiralTable sim_inspiral;
memset(&sim_inspiral, 0, sizeof(sim_inspiral));
sim_inspiral.distance = 1;
sim_inspiral.longitude = ra;
sim_inspiral.latitude = dec;
sim_inspiral.inclination = inclination;
sim_inspiral.polarization = polarization;
sim_inspiral.geocent_end_time = epoch;
sim_inspiral.end_time_gmst = gmst;
XLALInspiralSiteTimeAndDist(
&sim_inspiral, detector, &epoch, &abs_expected);
abs_expected = 1 / abs_expected;
}
/* This is the *really* slow way of working out the signal amplitude:
* generate a frequency-domain waveform and inject it. */
double complex expected = 0;
{
COMPLEX16FrequencySeries
*Htemplate = NULL, *Hsignal = NULL, *Hcross = NULL;
double h = 0, Fplus, Fcross;
int ret;
LALDict *params = XLALCreateDict();
assert(params);
XLALSimInspiralWaveformParamsInsertPNPhaseOrder(params, LAL_PNORDER_NEWTONIAN);
XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(params, LAL_PNORDER_NEWTONIAN);
/* Calculate antenna factors */
XLALComputeDetAMResponse(
&Fplus, &Fcross, detector->response, ra, dec, polarization, gmst);
/* Check that the way I compute the antenna factors matches */
{
double complex F;
XLALComputeDetAMResponse(
(double *)&F, 1 + (double *)&F, detector->response, ra, dec, 0, gmst);
F *= exp_i(-2 * polarization);
gsl_test_abs(creal(F), Fplus, 4 * GSL_DBL_EPSILON,
"testing Fplus(ra=%0.03f, dec=%0.03f, "
"inc=%0.03f, pol=%0.03f, epoch_nanos=%lu, detector=%s)",
ra, dec, inclination, polarization, epoch_nanos, instrument);
gsl_test_abs(cimag(F), Fcross, 4 * GSL_DBL_EPSILON,
"testing Fcross(ra=%0.03f, dec=%0.03f, "
"inc=%0.03f, pol=%0.03f, epoch_nanos=%lu, detector=%s)",
ra, dec, inclination, polarization, epoch_nanos, instrument);
}
/* "Template" waveform with inclination angle of zero */
ret = XLALSimInspiralFD(
&Htemplate, &Hcross, 1.4 * LAL_MSUN_SI, 1.4 * LAL_MSUN_SI,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 100, 101, 100, params,
TaylorF2);
assert(ret == XLAL_SUCCESS);
/* Discard any non-quadrature phase component of "template" */
for (size_t i = 0; i < Htemplate->data->length; i ++)
Htemplate->data->data[i] += I * Hcross->data->data[i];
/* Free Hcross */
XLALDestroyCOMPLEX16FrequencySeries(Hcross);
Hcross = NULL;
/* Normalize "template" */
for (size_t i = 0; i < Htemplate->data->length; i ++)
h += cabs2(Htemplate->data->data[i]);
h = 2 / h;
for (size_t i = 0; i < Htemplate->data->length; i ++)
Htemplate->data->data[i] *= h;
/* "Signal" waveform with requested inclination angle */
ret = XLALSimInspiralFD(
&Hsignal, &Hcross, 1.4 * LAL_MSUN_SI, 1.4 * LAL_MSUN_SI,
0, 0, 0, 0, 0, 0, 1, inclination, 0, 0, 0, 0, 1, 100, 101, 100,
params, TaylorF2);
assert(ret == XLAL_SUCCESS);
/* Project "signal" using antenna factors */
for (size_t i = 0; i < Hsignal->data->length; i ++)
Hsignal->data->data[i] = Fplus * Hsignal->data->data[i]
+ Fcross * Hcross->data->data[i];
/* Free Hcross */
XLALDestroyCOMPLEX16FrequencySeries(Hcross);
Hcross = NULL;
/* Work out complex amplitude by comparing with "template" waveform */
for (size_t i = 0; i < Htemplate->data->length; i ++)
expected += conj(Htemplate->data->data[i]) * Hsignal->data->data[i];
XLALDestroyDict(params);
XLALDestroyCOMPLEX16FrequencySeries(Hsignal);
XLALDestroyCOMPLEX16FrequencySeries(Htemplate);
Hsignal = Htemplate = NULL;
}
/* Test to nearly float (32-bit) precision because
* XLALInspiralSiteTimeAndDist returns result as float. */
gsl_test_abs(cabs(result), abs_expected, 1.5 * GSL_FLT_EPSILON,
"testing abs(signal_amplitude_model(ra=%0.03f, dec=%0.03f, "
"inc=%0.03f, pol=%0.03f, epoch_nanos=%lu, detector=%s))",
ra, dec, inclination, polarization, epoch_nanos, instrument);
gsl_test_abs(creal(result), creal(expected), 4 * GSL_DBL_EPSILON,
"testing re(signal_amplitude_model(ra=%0.03f, dec=%0.03f, "
"inc=%0.03f, pol=%0.03f, epoch_nanos=%lu, detector=%s))",
ra, dec, inclination, polarization, epoch_nanos, instrument);
gsl_test_abs(cimag(result), cimag(expected), 4 * GSL_DBL_EPSILON,
"testing im(signal_amplitude_model(ra=%0.03f, dec=%0.03f, "
"inc=%0.03f, pol=%0.03f, epoch_nanos=%lu, detector=%s))",
ra, dec, inclination, polarization, epoch_nanos, instrument);
}
static void test_log_radial_integral(
double expected, double tol, double r1, double r2, double p2, double b, int k)
{
......@@ -1437,15 +1286,6 @@ int bayestar_test(void)
test_complex_catrom();
test_eval_snr();
for (double ra = -M_PI; ra <= M_PI; ra += 0.4 * M_PI)
for (double dec = -M_PI_2; dec <= M_PI_2; dec += 0.4 * M_PI)
for (double inc = -M_PI; inc <= M_PI; inc += 0.4 * M_PI)
for (double pol = -M_PI; pol <= M_PI; pol += 0.4 * M_PI)
for (unsigned long epoch = 1000000000000000000;
epoch <= 1000086400000000000;
epoch += 7200000000000)
test_signal_amplitude_model(ra, dec, inc, pol, epoch, "H1");
/* Tests of radial integrand with p2=0, b=0. */
test_log_radial_integral(0, 0, 0, 1, 0, 0, 0);
test_log_radial_integral(0, 0, exp(1), exp(2), 0, 0, -1);
......
......@@ -125,6 +125,17 @@ double bayestar_log_likelihood_toa_phoa_snr(
const double *horizons /* SNR=1 horizon distances for each detector */
);
/* Expression for complex amplitude on arrival (without 1/distance factor).
* This is more of an internal function, but it's *really* important that
* it agrees with LAL conventions, so we expose it in the interface in order
* to be to validate it in Python against the LALSimluation SWIG bindings. */
double complex bayestar_signal_amplitude_model(
double complex F, /* Complex antenna factor */
double complex exp_i_twopsi, /* e^(i*2*psi), for polarization angle psi */
double u, /* cos(inclination) */
double u2 /* cos^2(inclination */
);
/* Unit test suite. Return EXIT_SUCCESS if tests passed,
* or otherwise EXIT_FAILURE. */
int bayestar_test(void);
......
......@@ -18,6 +18,7 @@ test_programs += test_cubic_interp
# Disable test_multiband.sh for now
#test_scripts = test_multiband.sh
if SWIG_BUILD_PYTHON
test_scripts += test_bayestar_signal_amplitude_model.py
test_scripts += test_plot.py
test_scripts += test_detframe.py
test_scripts += test_io_events.py
......
from __future__ import division
from __future__ import print_function
import os
import sys
try:
from astropy.tests.helper import pytest
except ImportError:
print('these tests require pytest', file=sys.stderr)
raise SystemExit(77)
from pytest import approx, mark
import lal
import lalinspiral
import lalmetaio
import lalsimulation
from lalinference.bayestar._sky_map import signal_amplitude_model
import numpy as np
def get_eff_dist(detector, ra, dec, inclination, polarization, epoch, gmst):
sim_inspiral = lalmetaio.SimInspiralTable()
sim_inspiral.distance = 1
sim_inspiral.longitude = ra
sim_inspiral.latitude = dec
sim_inspiral.inclination = inclination
sim_inspiral.polarization = polarization
sim_inspiral.geocent_end_time = epoch
sim_inspiral.end_time_gmst = gmst
_ = lal.LIGOTimeGPS()
_, eff_dist = lalinspiral.InspiralSiteTimeAndDist(
sim_inspiral, detector, _)
return eff_dist
def get_complex_antenna(response, ra, dec, gmst):
Fplus, Fcross = lal.ComputeDetAMResponse(response, ra, dec, 0, gmst)
return Fplus + 1j * Fcross
@mark.parametrize('ra', np.arange(-np.pi, 1.1 * np.pi, 0.8 * np.pi))
@mark.parametrize('dec', np.arange(-np.pi, 1.1 * np.pi, 0.8 * np.pi))
@mark.parametrize('inclination', np.arange(-np.pi, 1.1 * np.pi, 0.8 * np.pi))
@mark.parametrize('polarization', np.arange(-np.pi, 1.1 * np.pi, 0.8 * np.pi))
@mark.parametrize('epoch', np.arange(1000000000.0, 1000086400.0, 14400.0))
@mark.parametrize('instrument', ['H1'])
def test_bayestar_signal_amplitude_model(ra, dec, inclination, polarization,
epoch, instrument):
"""Test BAYESTAR signal amplitude model against LAL injection code."""
detector = lalsimulation.DetectorPrefixToLALDetector(instrument)
epoch = lal.LIGOTimeGPS(epoch)
gmst = lal.GreenwichMeanSiderealTime(epoch)
exp_i_twopsi = np.exp(2j * polarization)
u = np.cos(inclination)
u2 = np.square(u)
F = get_complex_antenna(detector.response, ra, dec, gmst)
result = signal_amplitude_model(F, exp_i_twopsi, u, u2)
abs_expected = 1 / get_eff_dist(
detector, ra, dec, inclination, polarization, epoch, gmst)
# This is the *really* slow way of working out the signal amplitude:
# generate a frequency-domain waveform and inject it.
params = lal.CreateDict()
lalsimulation.SimInspiralWaveformParamsInsertPNPhaseOrder(
params, lalsimulation.PNORDER_NEWTONIAN)
lalsimulation.SimInspiralWaveformParamsInsertPNAmplitudeOrder(
params, lalsimulation.PNORDER_NEWTONIAN)
# Calculate antenna factors
Fplus, Fcross = lal.ComputeDetAMResponse(
detector.response, ra, dec, polarization, gmst)
# Check that the way I compute the antenna factors matches
F = get_complex_antenna(detector.response, ra, dec, gmst)
F *= np.exp(-2j * polarization)
assert F.real == approx(Fplus, abs=4 * np.finfo(np.float64).eps)
assert F.imag == approx(Fcross, abs=4 * np.finfo(np.float64).eps)
# "Template" waveform with inclination angle of zero
Htemplate, Hcross = lalsimulation.SimInspiralFD(
1.4 * lal.MSUN_SI, 1.4 * lal.MSUN_SI,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 100, 101, 100, params,
lalsimulation.TaylorF2)
# Discard any non-quadrature phase component of "template"
Htemplate.data.data += 1j * Hcross.data.data
# Normalize "template"
h = np.sum(np.square(Htemplate.data.data.real) +
np.square(Htemplate.data.data.imag))
h = 2 / h
Htemplate.data.data *= h
# "Signal" waveform with requested inclination angle
Hsignal, Hcross = lalsimulation.SimInspiralFD(
1.4 * lal.MSUN_SI, 1.4 * lal.MSUN_SI,
0, 0, 0, 0, 0, 0, 1, inclination, 0, 0, 0, 0, 1, 100, 101, 100,
params, lalsimulation.TaylorF2)
# Project "signal" using antenna factors
Hsignal.data.data = Fplus * Hsignal.data.data + Fcross * Hcross.data.data
# Work out complex amplitude by comparing with "template" waveform
expected = np.sum(Htemplate.data.data.conj() * Hsignal.data.data)
# Test to nearly float (32-bit) precision because
# lalinspiral.InspiralSiteTimeAndDist returns result as float.
assert abs(expected) == approx(abs_expected,
abs=1.5 * np.finfo(np.float32).eps)
assert abs(result) == approx(abs_expected,
abs=1.5 * np.finfo(np.float32).eps)
assert result.real == approx(expected.real,
abs=4 * np.finfo(np.float64).eps)
assert result.imag == approx(expected.imag,
abs=4 * np.finfo(np.float64).eps)
if __name__ == '__main__':
raise SystemExit(pytest.main(['-x', __file__]))
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