#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Will Farr, Johnathon Merrit, and Derek Davis
# https://git.ligo.org/will-farr/lvkratesummaryplot/
# Copyright 2019
# JD Merritt jonathan.merritt@ligo.org
# Copyright 2020
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
# MA 02110-1301, USA.
import matplotlib
import matplotlib.pyplot as plt
import astropy.time as at
import h5py
import numpy as np
import matplotlib.ticker as ticker
from scipy.interpolate import interp1d
from scipy.stats import gamma, poisson
import seaborn as sns
matplotlib.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
#sns.set_context('talk')
sns.set_style('ticks')
sns.set_palette('colorblind')
#03a paper format convention:
matplotlib.use('agg')
matplotlib.rcParams['text.usetex'] = True
matplotlib.rcParams['font.size'] = 9
matplotlib.rcParams['savefig.dpi'] = 300
matplotlib.rcParams['legend.fontsize'] = 9
#Hardcoded event times for O1/O2/O3a
O1evts = [at.Time('2015-09-14'),
at.Time('2015-10-12'),
at.Time('2015-12-26')]
O2evts = [at.Time('2017-01-04'),
at.Time('2017-06-08'),
at.Time('2017-07-29'),
at.Time('2017-08-09'),
at.Time('2017-08-14'),
at.Time('2017-08-17'),
at.Time('2017-08-18'),
at.Time('2017-08-23')]
O3aevts = [at.Time('2019-09-30'),
at.Time('2019-09-29'),
at.Time('2019-09-24'),
at.Time('2019-09-15'),
at.Time('2019-09-10'),
at.Time('2019-09-09'),
at.Time('2019-08-28'),
at.Time('2019-08-28'),
at.Time('2019-08-14'),
at.Time('2019-08-03'),
at.Time('2019-07-31'),
at.Time('2019-07-28'),
at.Time('2019-07-27'),
at.Time('2019-07-20'),
at.Time('2019-07-19'),
at.Time('2019-07-08'),
at.Time('2019-07-07'),
at.Time('2019-07-06'),
at.Time('2019-07-01'),
at.Time('2019-06-30'),
at.Time('2019-06-20'),
at.Time('2019-06-02'),
at.Time('2019-05-27'),
at.Time('2019-05-21'),
at.Time('2019-05-21'),
at.Time('2019-05-19'),
at.Time('2019-05-17'),
at.Time('2019-05-14'),
at.Time('2019-05-13'),
at.Time('2019-05-12'),
at.Time('2019-05-03'),
at.Time('2019-04-26'),
at.Time('2019-04-25'),
at.Time('2019-04-24'),
at.Time('2019-04-21'),
at.Time('2019-04-13'),
at.Time('2019-04-13'),
at.Time('2019-04-12'),
at.Time('2019-04-08')]
allevts = np.concatenate((O1evts, O2evts, O3aevts))
allevts_jyear = [e.jyear for e in allevts]
O1VTs = np.genfromtxt('NETWORK-CUMULATIVE_TIME_VOLUME-O1.csv', delimiter=',')
O1ts = at.Time(O1VTs[:,0], format='gps')
O1VTs = O1VTs[:,1]
O2VTs = np.genfromtxt('NETWORK-CUMULATIVE_TIME_VOLUME-O2.csv', delimiter=',')
O2ts = at.Time(O2VTs[:,0], format='gps')
O2VTs = O2VTs[:,1]
O3aVTs = np.genfromtxt('NETWORK-CUMULATIVE_TIME_VOLUME-O3a.csv', delimiter=',')
O3ats = at.Time(O3aVTs[:,0], format='gps')
O3aVTs = O3aVTs[:,1]
allts_jyear = np.concatenate((O1ts.jyear, O2ts.jyear, O3ats.jyear))
allVTs = np.concatenate((O1VTs, O1VTs[-1]+O2VTs, O1VTs[-1]+O2VTs[-1]+O3aVTs))
jyear_of_VT = interp1d(allVTs, allts_jyear, bounds_error=False, fill_value=(allts_jyear[0], allts_jyear[-1]))
VT_of_jyear = interp1d(allts_jyear, allVTs, bounds_error=False, fill_value=(0.0, allVTs[-1]))
end_O1 = O1ts[-1].jyear
end_O2 = O2ts[-1].jyear
@ticker.FuncFormatter
def VT_to_year_formatter(VT, pos=None):
return '{:.1f}'.format(jyear_of_VT(VT))
vts = np.linspace(0, allVTs[-1], 16384)
ts = jyear_of_VT(vts)
#Generate synthetic detection histories based on a poisson likelihood.
rstate = np.random.get_state()
np.random.seed(971059855)
try:
N = 10000
p_L = gamma(len(allevts)+1)
Ns = poisson(p_L.rvs(N)).rvs()
synthetic_events = []
for N in Ns:
vs = np.sort(np.random.uniform(low=0, high=allVTs[-1], size=N))
synthetic_events.append(vs)
synthetic_events = np.array(synthetic_events)
nsyn = []
for syn in synthetic_events:
ns, _ = np.histogram(syn, bins=vts)
ns = np.concatenate(((0,), ns))
nsyn.append(np.cumsum(ns))
nsyn = np.array(nsyn)
finally:
np.random.set_state(rstate)
#plotting
fig, ax = plt.subplots(figsize=(3.375,3))
#Rsyn = nsyn / vts[np.newaxis,:]
nallts = np.array([np.count_nonzero(allevts_jyear < t) for t in ts])
#Rallts = nallts / vts
l, = ax.plot(vts, np.median(nsyn, axis=0))
ax.fill_between(vts, np.percentile(nsyn, 75, axis=0), np.percentile(nsyn, 25, axis=0), color=l.get_color(), alpha=0.25)
ax.fill_between(vts, np.percentile(nsyn, 95, axis=0), np.percentile(nsyn, 5, axis=0), color=l.get_color(), alpha=0.25)
ax.plot(vts, nallts, '-k')
ax.set_ylabel(r"$\mathrm{Cumulative \ Detections}$")
ax.set_xlim(0, np.max(vts))
ax.set_xlabel(r'$\mathrm{Effective \ BNS \ VT \ } [\mathrm{Mpc}^3 \ \mathrm{kyr}]$')
ymin,ymax = ax.get_ylim()
ymin = 0
ax.set_ylim(ymin, ymax)
ax.fill_betweenx([ymin, ymax], VT_of_jyear(end_O1), 0, color=sns.color_palette()[1], alpha=0.2)
ax.fill_betweenx([ymin, ymax], VT_of_jyear(end_O2), VT_of_jyear(end_O1), color=sns.color_palette()[2], alpha=0.2)
ax.fill_betweenx([ymin, ymax], vts[-1], VT_of_jyear(end_O2), color=sns.color_palette()[3], alpha=0.2)
ax.text(20, 50, r'$\mathrm{O1}$', fontsize=9)
ax.text(250, 50, r'$\mathrm{O2}$', fontsize=9)
ax.text(750, 50, r'$\mathrm{O3a}$', fontsize=9)
plt.tight_layout()
plt.savefig('BNS_VT.pdf')