Multiple detectors, O1 data#
Here we compare lalpulsar_parameter_estimation_nested
with cwinpy
in the case of real
gravitational-wave data from two detectors (H1 and L1). This data (you can downloaded the
H1
and L1
data) comes from the O1 run and has been heterodyned using
the parameters for the Crab pulsar given in this file: J0534+2200.par
. The parameters being estimated are \(h_0\), \(\phi_0\),
\(\psi\) and \(\cos{\iota}\), all with uniform priors.
The script for this comparison, using the dynesty nested sampling algorithm, is shown at the bottom of the page. It produces the following comparison data:
Method |
\(\ln{(Z)}\) |
\(\ln{(Z)}\) noise |
\(\ln{}\) Odds |
---|---|---|---|
|
22220423.101 |
22220428.879 |
-5.779±0.078 |
|
22121194.510 |
22121200.131 |
-5.621±0.071 |
|
22121194.427 |
-5.704 |
Method |
\(h_0\) |
\(\phi_0\) (rad) |
\(\psi\) (rad) |
\(\cos{\iota}\) |
---|---|---|---|---|
|
1.98±1.48×10-26 |
1.90±1.09 |
0.84±0.42 |
-0.00±0.49 |
90% credible intervals |
[0.19, 4.87]×10-26 |
[0.07, 3.07] |
[0.10, 1.46] |
[-0.83, 0.86] |
|
1.97±1.48×10-26 |
1.98±1.07 |
0.87±0.42 |
0.04±0.46 |
90% credible intervals |
[0.16, 4.98]×10-26 |
[0.09, 3.08] |
[0.12, 1.46] |
[-0.78, 0.81] |
Method |
\(h_0\) |
\(\phi_0\) (rad) |
\(\psi\) (rad) |
\(\cos{\iota}\) |
\(\ln{(L)}\) max |
---|---|---|---|---|---|
|
4.07×10-26 |
2.91 |
1.07 |
0.00 |
22220431.18 |
|
4.20×10-26 |
2.92 |
1.13 |
0.09 |
22121202.43 |
#!/usr/bin/env python
"""
Compare cwinpy with lalpulsar_parameter_estimation_nested for O1 data
from two detectors (H1 and L1). This uses data for the Crab pulsar.
"""
import os
import subprocess as sp
import h5py
import matplotlib
import numpy as np
from astropy.utils.data import download_file
from bilby.core.prior import Uniform
from comparitors import comparisons
from lalinference import LALInferenceHDF5PosteriorSamplesDatasetName
from lalinference.io import read_samples
from matplotlib import pyplot as plt
from cwinpy.pe import pe
from cwinpy.plot import Plot
matplotlib.use("Agg")
# URL for ephemeris files
DOWNLOAD_URL = "https://git.ligo.org/lscsoft/lalsuite/raw/master/lalpulsar/lib/{}"
label = "multi_detector_O1_data"
outdir = "outputs"
if not os.path.isdir(outdir):
os.makedirs(outdir)
# set the par file
parfile = os.path.join("data", "J0534+2200.par")
# set the data files
hetfiles = [
os.path.join("data", "O1_Crab_H1.txt.gz"),
os.path.join("data", "O1_Crab_L1.txt.gz"),
]
# set the detector names
detectors = ["H1", "L1"]
# create priors
phi0range = [0.0, np.pi]
psirange = [0.0, np.pi / 2.0]
cosiotarange = [-1.0, 1.0]
h0range = [0.0, 1e-23]
# set prior for lalpulsar_parameter_estimation_nested
priorfile = os.path.join(outdir, "{}_prior.txt".format(label))
priorcontent = """H0 uniform {} {}
PHI0 uniform {} {}
PSI uniform {} {}
COSIOTA uniform {} {}
"""
with open(priorfile, "w") as fp:
fp.write(priorcontent.format(*(h0range + phi0range + psirange + cosiotarange)))
# set prior for bilby
priors = {}
priors["h0"] = Uniform(h0range[0], h0range[1], "h0", latex_label=r"$h_0$")
priors["phi0"] = Uniform(
phi0range[0], phi0range[1], "phi0", latex_label=r"$\phi_0$", unit="rad"
)
priors["psi"] = Uniform(
psirange[0], psirange[1], "psi", latex_label=r"$\psi$", unit="rad"
)
priors["cosiota"] = Uniform(
cosiotarange[0], cosiotarange[1], "cosiota", latex_label=r"$\cos{\iota}$"
)
# run lalpulsar_parameter_estimation_nested
try:
execpath = os.environ["CONDA_PREFIX"]
except KeyError:
raise KeyError(
"Please work in a conda environment with lalsuite and cwinpy installed"
)
execpath = os.path.join(execpath, "bin")
lppen = os.path.join(execpath, "lalpulsar_parameter_estimation_nested")
n2p = os.path.join(execpath, "lalinference_nest2pos")
Nlive = 1000 # number of nested sampling live points
Nmcmcinitial = 0 # set to 0 so that prior samples are not resampled
outfile = os.path.join(outdir, "{}_nest.hdf".format(label))
# set ephemeris files
efile = download_file(DOWNLOAD_URL.format("earth00-40-DE405.dat.gz"), cache=True)
sfile = download_file(DOWNLOAD_URL.format("sun00-40-DE405.dat.gz"), cache=True)
tfile = download_file(DOWNLOAD_URL.format("te405_2000-2040.dat.gz"), cache=True)
# set the command line arguments
runcmd = " ".join(
[
lppen,
"--verbose",
"--input-files",
",".join(hetfiles),
"--detectors",
",".join(detectors),
"--par-file",
parfile,
"--prior-file",
priorfile,
"--Nlive",
"{}".format(Nlive),
"--Nmcmcinitial",
"{}".format(Nmcmcinitial),
"--outfile",
outfile,
"--ephem-earth",
efile,
"--ephem-sun",
sfile,
"--ephem-timecorr",
tfile,
]
)
with sp.Popen(
runcmd,
stdout=sp.PIPE,
stderr=sp.PIPE,
shell=True,
bufsize=1,
universal_newlines=True,
) as p:
for line in p.stderr:
print(line, end="")
# convert nested samples to posterior samples
outpost = os.path.join(outdir, "{}_post.hdf".format(label))
runcmd = " ".join([n2p, "-p", outpost, outfile])
with sp.Popen(
runcmd,
stdout=sp.PIPE,
stderr=sp.PIPE,
shell=True,
bufsize=1,
universal_newlines=True,
) as p:
for line in p.stdout:
print(line, end="")
# get posterior samples
post = read_samples(outpost, tablename=LALInferenceHDF5PosteriorSamplesDatasetName)
lp = len(post["H0"])
postsamples = np.zeros((lp, len(priors)))
for i, p in enumerate(priors.keys()):
postsamples[:, i] = post[p.upper()]
# get evidence
hdf = h5py.File(outpost, "r")
a = hdf["lalinference"]["lalinference_nest"]
evsig = a.attrs["log_evidence"]
evnoise = a.attrs["log_noise_evidence"]
hdf.close()
# run bilby via the pe interface
runner = pe(
data_file=hetfiles,
par_file=parfile,
prior=priors,
detector=detectors,
outdir=outdir,
label=label,
)
result = runner.result
# evaluate the likelihood on a grid
gridpoints = 30
grid_size = dict()
for p in priors.keys():
grid_size[p] = np.linspace(
np.min(result.posterior[p]), np.max(result.posterior[p]), gridpoints
)
grunner = pe(
data_file=hetfiles,
par_file=parfile,
prior=priors,
detector=detectors,
outdir=outdir,
label=label,
grid=True,
grid_kwargs={"grid_size": grid_size},
numba=True,
)
grid = grunner.grid
# output comparisons
comparisons(label, outdir, grid, priors, cred=0.9)
# create results plot
allresults = {
"lalpulsar_parameter_estimation_nested": outpost,
"cwinpy_pe": result,
"cwinpy_pe (grid)": grid,
}
colors = {
key: plt.rcParams["axes.prop_cycle"].by_key()["color"][i]
for i, key in enumerate(allresults.keys())
}
plot = Plot(
results=allresults,
parameters=list(priors.keys()),
plottype="corner",
)
plot.plot(
bins=50,
smooth=0.9,
quantiles=[0.16, 0.84],
levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.0)),
fill_contours=True,
colors=colors,
)
plot.savefig(os.path.join(outdir, "{}_corner.png".format(label)), dpi=150)