From ffdeeb226c7b4686b3bf99b501c0476b91390ae5 Mon Sep 17 00:00:00 2001
From: Sylvia Biscoveanu <sylvia.biscoveanu@ligo.org>
Date: Wed, 24 Apr 2019 01:02:51 -0500
Subject: [PATCH] Calibration plotting
---
bilby/core/utils.py | 25 ++++++++++++++
bilby/gw/result.py | 81 +++++++++++++++++++++++++++++++++++++++++++++
bilby/gw/utils.py | 71 ++++++++++++++++++++++++++++++++++++++-
3 files changed, 176 insertions(+), 1 deletion(-)
diff --git a/bilby/core/utils.py b/bilby/core/utils.py
index 5e97a641b..5636a594a 100644
--- a/bilby/core/utils.py
+++ b/bilby/core/utils.py
@@ -703,6 +703,31 @@ def logtrapzexp(lnf, dx):
return np.log(dx / 2.) + logsumexp([logsumexp(lnf[:-1]), logsumexp(lnf[1:])])
+def credible_interval(samples, confidence_level=.9, lower=True):
+ """
+ Return location of lower or upper confidence levels
+ Based on lalinference.bayespputils.cred_interval
+
+ Parameters
+ ----------
+ x: List of samples.
+ cl: Confidence level to return the bound of.
+ lower: If ``True``, return the lower bound, otherwise return the upper bound.
+
+ Returns
+ -------
+ float: the upper or lower confidence level
+
+ """
+ def cred_level(cl, x):
+ return np.sort(x, axis=0)[int(cl * len(x))]
+
+ if lower:
+ return cred_level((1. - confidence_level) / 2, samples)
+ else:
+ return cred_level((1. + confidence_level) / 2, samples)
+
+
def run_commandline(cl, log_level=20, raise_error=True, return_output=True):
"""Run a string cmd as a subprocess, check for errors and return output.
diff --git a/bilby/gw/result.py b/bilby/gw/result.py
index 8f38de813..ed9e0ada1 100644
--- a/bilby/gw/result.py
+++ b/bilby/gw/result.py
@@ -1,7 +1,13 @@
from __future__ import division
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+
from ..core.result import Result as CoreResult
from ..core.utils import logger
+from .utils import (plot_spline_pos,
+ spline_angle_xform)
class CompactBinaryCoalesenceResult(CoreResult):
@@ -80,5 +86,80 @@ class CompactBinaryCoalesenceResult(CoreResult):
logger.info("No injection for detector {}".format(detector))
return None
+ def plot_calibration_posterior(self, level=.9):
+ """ Plots the calibration amplitude and phase uncertainty.
+ Adapted from the LALInference version in bayespputils
+
+ Parameters
+ ----------
+ level: float, percentage for confidence levels
+
+ Returns
+ -------
+ saves a plot to outdir+label+calibration.png
+
+ """
+ fig, [ax1, ax2] = plt.subplots(2, 1, figsize=(15, 15), dpi=500)
+ posterior = self.posterior
+
+ font_size = 32
+ outdir = self.outdir
+
+ parameters = posterior.keys()
+ ifos = np.unique([param.split('_')[1] for param in parameters if 'recalib_' in param])
+ if ifos.size == 0:
+ logger.info("No calibration parameters found. Aborting calibration plot.")
+ return
+
+ for ifo in ifos:
+ if ifo == 'H1':
+ color = 'r'
+ elif ifo == 'L1':
+ color = 'g'
+ elif ifo == 'V1':
+ color = 'm'
+ else:
+ color = 'c'
+
+ # Assume spline control frequencies are constant
+ freq_params = np.sort([param for param in parameters if
+ 'recalib_{0}_frequency_'.format(ifo) in param])
+
+ logfreqs = np.log([posterior[param].iloc[0] for param in freq_params])
+
+ # Amplitude calibration model
+ plt.sca(ax1)
+ amp_params = np.sort([param for param in parameters if
+ 'recalib_{0}_amplitude_'.format(ifo) in param])
+ if len(amp_params) > 0:
+ amplitude = 100 * np.column_stack([posterior[param] for param in amp_params])
+ plot_spline_pos(logfreqs, amplitude, color=color, level=level,
+ label="{0} (mean, {1}%)".format(ifo.upper(), int(level * 100)))
+
+ # Phase calibration model
+ plt.sca(ax2)
+ phase_params = np.sort([param for param in parameters if
+ 'recalib_{0}_phase_'.format(ifo) in param])
+ if len(phase_params) > 0:
+ phase = np.column_stack([posterior[param] for param in phase_params])
+ plot_spline_pos(logfreqs, phase, color=color, level=level,
+ label="{0} (mean, {1}%)".format(ifo.upper(), int(level * 100)),
+ xform=spline_angle_xform)
+
+ ax1.tick_params(labelsize=.75 * font_size)
+ ax2.tick_params(labelsize=.75 * font_size)
+ plt.legend(loc='upper right', prop={'size': .75 * font_size}, framealpha=0.1)
+ ax1.set_xscale('log')
+ ax2.set_xscale('log')
+
+ ax2.set_xlabel('Frequency (Hz)', fontsize=font_size)
+ ax1.set_ylabel('Amplitude (%)', fontsize=font_size)
+ ax2.set_ylabel('Phase (deg)', fontsize=font_size)
+
+ filename = os.path.join(outdir, self.label + '_calibration.png')
+ fig.tight_layout()
+ fig.savefig(filename, bbox_inches='tight')
+ plt.close(fig)
+
CBCResult = CompactBinaryCoalesenceResult
diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py
index 3bb38eb7a..934bd3914 100644
--- a/bilby/gw/utils.py
+++ b/bilby/gw/utils.py
@@ -3,10 +3,13 @@ import os
import json
import numpy as np
+from scipy import interpolate
+import matplotlib.pyplot as plt
from ..core.utils import (gps_time_to_gmst, ra_dec_to_theta_phi,
speed_of_light, logger, run_commandline,
- check_directory_exists_and_if_not_mkdir)
+ check_directory_exists_and_if_not_mkdir,
+ credible_interval)
try:
from gwpy.timeseries import TimeSeries
@@ -856,3 +859,69 @@ def lalsim_SimInspiralWaveformParamsInsertTidalLambda2(
return lalsim.SimInspiralWaveformParamsInsertTidalLambda2(
waveform_dictionary, lambda_2)
+
+
+def spline_angle_xform(delta_psi):
+ """
+ Returns the angle in degrees corresponding to the spline
+ calibration parameters delta_psi.
+ Based on the same function in lalinference.bayespputils
+
+ Parameters
+ ----------
+ delta_psi: calibration phase uncertainity
+
+ Returns
+ -------
+ float: delta_psi in degrees
+
+ """
+ rotation = (2.0 + 1.0j * delta_psi) / (2.0 - 1.0j * delta_psi)
+
+ return 180.0 / np.pi * np.arctan2(np.imag(rotation), np.real(rotation))
+
+
+def plot_spline_pos(log_freqs, samples, nfreqs=100, level=0.9, color='k', label=None, xform=None):
+ """
+ Plot calibration posterior estimates for a spline model in log space.
+ Adapted from the same function in lalinference.bayespputils
+
+ Parameters
+ ----------
+ log_freqs: The (log) location of spline control points.
+ samples: List of posterior draws of function at control points ``log_freqs``
+ nfreqs: Number of points to evaluate spline at for plotting.
+ level: Credible level to fill in.
+ color: Color to plot with.
+ label: Label for plot.
+ xform: Function to transform the spline into plotted values.
+
+ """
+ freq_points = np.exp(log_freqs)
+ freqs = np.logspace(min(log_freqs), max(log_freqs), nfreqs, base=np.exp(1))
+
+ data = np.zeros((samples.shape[0], nfreqs))
+
+ if xform is None:
+ scaled_samples = samples
+ else:
+ scaled_samples = xform(samples)
+
+ mu = np.mean(scaled_samples, axis=0)
+ lower_confidence_level = mu - credible_interval(scaled_samples, level, lower=True)
+ upper_confidence_level = credible_interval(scaled_samples, level, lower=False) - mu
+ plt.errorbar(freq_points, mu, yerr=[lower_confidence_level, upper_confidence_level],
+ fmt='.', color=color, lw=4, alpha=0.5, capsize=0)
+
+ for i, sample in enumerate(samples):
+ temp = interpolate.spline(log_freqs, sample, np.log(freqs))
+ if xform is None:
+ data[i] = temp
+ else:
+ data[i] = xform(temp)
+
+ line, = plt.plot(freqs, np.mean(data, axis=0), color=color, label=label)
+ color = line.get_color()
+ plt.fill_between(freqs, credible_interval(data, level), credible_interval(data, level, lower=False),
+ color=color, alpha=.1, linewidth=0.1)
+ plt.xlim(freq_points.min() - .5, freq_points.max() + 50)
--
GitLab