#!/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')