Commit a6eb6943 authored by Michal Was's avatar Michal Was
Browse files

Merge branch 'master' of ligo-vcs.phys.uwm.edu:/usr/local/git/lalsuite

Original: d7e676db363ec7b59f0df49c786c2f61c566e9c2
parents 38c49d81 b80708e6
......@@ -1452,20 +1452,29 @@ int XLALSetFlatLatticeEllipticalBounds(
// Check input
XLAL_CHECK(tiling != NULL, XLAL_EFAULT);
XLAL_CHECK(x_semi > 0.0, XLAL_EINVAL);
XLAL_CHECK(y_semi > 0.0, XLAL_EINVAL);
XLAL_CHECK(x_semi >= 0.0, XLAL_EINVAL);
XLAL_CHECK(y_semi >= 0.0, XLAL_EINVAL);
// Allocate and set bounds data
double* bounds = XLALCalloc(4, sizeof(double));
XLAL_CHECK(bounds != NULL, XLAL_ENOMEM);
bounds[0] = x_centre;
bounds[1] = y_centre;
bounds[2] = x_semi;
bounds[3] = y_semi;
// Set parameter space bound
// Set parameter space X bound
XLAL_CHECK(XLALSetFlatLatticeConstantBound(tiling, x_dimension, x_centre - x_semi, x_centre + x_semi) == XLAL_SUCCESS, XLAL_EFAILED);
XLAL_CHECK(XLALSetFlatLatticeBound(tiling, x_dimension + 1, false, EllipticalYBound, (void*)bounds) == XLAL_SUCCESS, XLAL_EFAILED);
// Set parameter space Y bound
if (x_semi == 0.0 || y_semi == 0.0) {
XLAL_CHECK(XLALSetFlatLatticeConstantBound(tiling, x_dimension + 1, y_centre - y_semi, y_centre + y_semi) == XLAL_SUCCESS, XLAL_EFAILED);
} else {
// Allocate and set bounds data
double* bounds = XLALCalloc(4, sizeof(double));
XLAL_CHECK(bounds != NULL, XLAL_ENOMEM);
bounds[0] = x_centre;
bounds[1] = y_centre;
bounds[2] = x_semi;
bounds[3] = y_semi;
// Set parameter space bound
XLAL_CHECK(XLALSetFlatLatticeBound(tiling, x_dimension + 1, false, EllipticalYBound, (void*)bounds) == XLAL_SUCCESS, XLAL_EFAILED);
}
return XLAL_SUCCESS;
......
......@@ -321,10 +321,13 @@ double XLALGetLocalRate(LALCosmologicalRateParameters *rate)
* Implements the fit to the SFR in Eq.7 of Coward, Burman 2005 ( http://arxiv.org/abs/astro-ph/0505181 )
* See also Porciani & Madau ( http://arxiv.org/pdf/astro-ph/0008294v2.pdf ) and references therein
**/
double XLALStarFormationDensity(double z, void *rate)
double XLALStarFormationDensity(double z, void *params)
{
LALCosmologicalRateParameters *p = (LALCosmologicalRateParameters *)rate;
return p->r0*(1.0+p->W)*exp(p->Q*z)/(exp(p->R*z)+p->W);
LALCosmologicalParametersAndRate *p = (LALCosmologicalParametersAndRate *)params;
double hz = XLALHubbleParameter(z,p->omega);
double x = 1.0/sqrt(1.+z);
double hz0 = x*x*x;
return hz0/hz*p->rate->r0*(1.0+p->rate->W)*exp(p->rate->Q*z)/(exp(p->rate->R*z)+p->rate->W);
}
/**
* Returns the Rate weighted uniform comoving volume density
......@@ -333,7 +336,7 @@ double XLALRateWeightedUniformComovingVolumeDensity(double z, void *params)
{
LALCosmologicalParametersAndRate *p = (LALCosmologicalParametersAndRate *)params;
double dvdz = XLALUniformComovingVolumeDensity(z,p->omega);
double ez = XLALStarFormationDensity(z,p->rate);
double ez = XLALStarFormationDensity(z,p);
return ez*dvdz;
}
......
......@@ -49,7 +49,10 @@ AC_CONFIG_FILES([\
src/online/Makefile \
src/hwinjection/Makefile \
])
AM_INIT_AUTOMAKE([1.11 foreign color-tests parallel-tests])
# FIXME: -Wno-unsupported is a hack which silences warnings regarding
# subdir-objects in automake-1.14+. this will be the default, and only,
# behaviour in automake-2.0+ so this is only a temporary measure
AM_INIT_AUTOMAKE([1.11 foreign color-tests parallel-tests -Wno-unsupported])
AH_TOP([
#ifndef CONFIG_H
#define CONFIG_H])
......
......@@ -28,6 +28,7 @@ summary values for each sky map:
* searched posterior probability
* angle between true sky location and maximum a posteriori estimate
* runtime in seconds
* (optional) areas of specified probability contours
The filenames of the sky maps may be provided as positional command line
arguments, and may also be provided as globs (such as '*.fits.gz').
......@@ -48,7 +49,11 @@ if __name__ == '__main__':
Option("-o", "--output", default="/dev/stdout",
help="Name of output file [default: %default]"),
Option("-j", "--jobs", default=1, type=int,
help="Number of threads [default: %default]")
help="Number of threads [default: %default]"),
Option("-p", "--contour", default=[], action="append",
type=float, metavar="PERCENT",
help="Report the area of the smallest contour containing this "
+ "much probability. Can be repeated mulitple times")
]
)
opts, args = parser.parse_args()
......@@ -70,9 +75,10 @@ from lalinference.bayestar import fits
from lalinference.bayestar import postprocess
def startup(dbfilename):
global db
def startup(dbfilename, opts_contour):
global db, contours
db = sqlite3.connect(dbfilename)
contours = opts_contour
def process(fitsfilename):
......@@ -94,10 +100,10 @@ def process(fitsfilename):
WHERE cem1.table_name = 'sim_inspiral'
AND cem2.table_name = 'coinc_event' AND cem2.event_id = ?""",
(coinc_event_id,)).fetchone()
searched_area, searched_prob, offset = postprocess.find_injection(
sky_map, true_ra, true_dec)
searched_area, searched_prob, offset, contour_areas = postprocess.find_injection(
sky_map, true_ra, true_dec, contours=[0.01 * p for p in contours])
return coinc_event_id, simulation_id, far, searched_area, searched_prob, offset, runtime
return [coinc_event_id, simulation_id, far, searched_area, searched_prob, offset, runtime] + contour_areas
if __name__ == '__main__':
......@@ -109,17 +115,19 @@ if __name__ == '__main__':
progress.update(-1, 'spawning {0} workers'.format(opts.jobs))
if opts.jobs == 1:
from itertools import imap
startup(dbfilename)
else:
import multiprocessing
imap = multiprocessing.Pool(opts.jobs, startup, (dbfilename,)).imap_unordered
imap = multiprocessing.Pool(opts.jobs, startup, (dbfilename, opts.contour)).imap_unordered
startup(dbfilename, opts.contour)
progress.update(-1, 'obtaining filenames of sky maps')
fitsfilenames = tuple(itertools.chain.from_iterable(glob.iglob(fitsfileglob)
for fitsfileglob in fitsfileglobs))
print('coinc_event_id', 'simulation_id', 'far', 'searched_area', 'searched_prob', 'offset', 'runtime',
sep="\t", file=outfile)
colnames = ['coinc_event_id', 'simulation_id', 'far', 'searched_area',
'searched_prob', 'offset', 'runtime'] + ["area({0:g})".format(p)
for p in contours]
print(*colnames, sep="\t", file=outfile)
count_records = 0
progress.max = len(fitsfilenames)
......
......@@ -83,7 +83,10 @@ log = logging.getLogger('BAYESTAR')
from lalinference.bayestar import fits
from lalinference.bayestar.ligolw_sky_map import gracedb_sky_map
sky_map, epoch, elapsed_time, instruments = gracedb_sky_map(open(infilename, "rb"), open(opts.psd, "rb"), opts.waveform, opts.f_low, opts.min_distance, opts.max_distance, opts.prior_distance_power, reference_frequency=opts.reference_frequency, nside=opts.nside)
sky_map, epoch, elapsed_time, instruments = gracedb_sky_map(
open(infilename, "rb"), open(opts.psd, "rb"), opts.waveform, opts.f_low,
opts.min_distance, opts.max_distance, opts.prior_distance_power,
reference_frequency=opts.reference_frequency, nside=opts.nside)
# Write sky map
fits.write_sky_map(opts.output, sky_map, gps_time=float(epoch),
......
......@@ -72,18 +72,17 @@ def ascii_trigger(line, columns=KLEINEWELLE_COLUMNS):
@returns a `SnglBurst` built from the `ASCII` data
"""
if isinstance(line, str):
dat = map(float, _delim.split(line.rstrip()))
else:
dat = map(float, line)
t = lsctables.SnglBurst()
t.search = u"kleinewelle"
dat = line.rstrip().split()
if len(dat) == 9:
channel = re.sub("_", ":", dat.pop(-1), 1)
if re.search("_\d+_\d+\Z", c):
if re.search("_\d+_\d+\Z", channel):
channel = channel.rsplit("_", 2)[0]
if 'channel' in columns:
t.channel = channel
elif len(dat) == 8:
if len(dat) == 8:
(start, stop, peak, freq, energy,
amplitude, n_pix, sig) = list(map(float, dat))
start = LIGOTimeGPS(start)
......@@ -96,8 +95,6 @@ def ascii_trigger(line, columns=KLEINEWELLE_COLUMNS):
raise ValueError("Wrong number of columns in ASCII line. "
"Cannot read.")
t = lsctables.SnglBurst()
t.search = u"kleinewelle"
# set times
if 'start_time' in columns:
t.start_time = start.gpsSeconds
......
......@@ -34,9 +34,6 @@ def ligolw_sky_map(sngl_inspirals, approximant, amplitude_order, phase_order, f_
"""Convenience function to produce a sky map from LIGO-LW rows. Note that
min_distance and max_distance should be in Mpc."""
if method == "toa_snr" and prior_distance_power is None:
raise ValueError("For method='toa_snr', the argument prior_distance_power is required.")
ifos = [sngl_inspiral.ifo for sngl_inspiral in sngl_inspirals]
# Extract masses from the table.
......
......@@ -28,7 +28,8 @@ import healpy as hp
def angle_distance(theta0, phi0, theta1, phi1):
"""Angular separation in radians between two points on the unit sphere."""
cos_angle_distance = np.cos(phi1 - phi0) * np.sin(theta0) * np.sin(theta1) + np.cos(theta0) * np.cos(theta1)
cos_angle_distance = (np.cos(phi1 - phi0) * np.sin(theta0) * np.sin(theta1)
+ np.cos(theta0) * np.cos(theta1))
if cos_angle_distance > 1:
return 0.
elif cos_angle_distance < -1:
......@@ -37,17 +38,20 @@ def angle_distance(theta0, phi0, theta1, phi1):
return np.arccos(cos_angle_distance)
def find_injection(sky_map, true_ra, true_dec):
def find_injection(sky_map, true_ra, true_dec, contours=()):
"""
Given a sky map and the true right ascension and declination (in radians),
find the smallest area in deg^2 that would have to be searched to find the
source, the smallest posterior mass, and the angular offset in degrees from
the true location to the maximum (mode) of the posterior.
the true location to the maximum (mode) of the posterior. Optionally, also
compute the areas of the smallest contours containing a given total
probability.
"""
# Compute the HEALPix lateral resolution parameter for this sky map.
npix = len(sky_map)
nside = hp.npix2nside(npix)
deg2perpix = hp.nside2pixarea(nside, degrees=True)
# Convert from ra, dec to conventional spherical polar coordinates.
true_theta = 0.5 * np.pi - true_dec
......@@ -69,21 +73,29 @@ def find_injection(sky_map, true_ra, true_dec):
# Find the index of the true location in the cumulative distribution.
idx = (i for i, pix in enumerate(indices) if pix == true_pix).next()
# Find the smallest area that would have to be searched to find the true
# location.
searched_area = (idx + 1) * hp.nside2pixarea(nside, degrees=True)
# Find the smallest area that would have to be searched to find
# the true location. Note that 1 is added to the index because we want
# the **length** of the array up to and including the idx'th element,
# not the index itself.
searched_area = (idx + 1) * deg2perpix
# Find the smallest posterior mass that would have to be searched to find
# the true location.
searched_prob = cum_sky_map[idx]
# Permute the cumulative distribution so that it is indexed the same way
# as the original sky map.
cum_sky_map[indices] = cum_sky_map
# Reverse pixel order so that we can use np.searchsorted() to locate
# the index where a given confidence level ends.
cum_sky_map = cum_sky_map[::-1]
# For each of the given confidence levels, compute the area of the
# smallest region containing that probability.
contour_areas = [(np.searchsorted(cum_sky_map, p, side='right') + 1)
* deg2perpix for p in contours]
# Find the angular offset between the mode and true locations.
offset = np.rad2deg(angle_distance(true_theta, true_phi, mode_theta, mode_phi))
offset = np.rad2deg(angle_distance(true_theta, true_phi,
mode_theta, mode_phi))
# Done.
return searched_area, searched_prob, offset
return searched_area, searched_prob, offset, contour_areas
......@@ -168,7 +168,7 @@ int XLALSetFlatLatticeFnDotConstantBound(
info->bounds[1] = GSL_MAX(bound1, bound2);
// Set parameter space bound
XLAL_CHECK(XLALSetFlatLatticeBound(tiling, dimension, false, FnDotConstantBound, (void*)info) == XLAL_SUCCESS, XLAL_EFAILED);
XLAL_CHECK(XLALSetFlatLatticeBound(tiling, dimension, bound1 == bound2, FnDotConstantBound, (void*)info) == XLAL_SUCCESS, XLAL_EFAILED);
return XLAL_SUCCESS;
......
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