Skip to content

Commit

Permalink
Improvements to single fit plots (gwastro#4509)
Browse files Browse the repository at this point in the history 
* Improvements to single fit plots

* Apply suggestions from Gareth

Co-authored-by: Gareth S Cabourn Davies <[email protected]>

---------

Co-authored-by: Gareth S Cabourn Davies <[email protected]>
  • Loading branch information
@titodalcanton @GarethCabournDavies
titodalcanton and GarethCabournDavies committed Feb 2, 2024
1 parent 6c2fb27 commit 4081be4
Show file tree
Hide file tree
Showing 2 changed files with 62 additions and 69 deletions.
45 changes: 20 additions & 25 deletions bin/live/pycbc_live_plot_combined_single_fits
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,18 +12,23 @@
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.

import h5py, numpy as np, argparse
"""Plot the time evolution of fit parameters of PyCBC Live triggers.
"""

import argparse
import logging
import numpy as np
import h5py
import matplotlib
matplotlib.use('agg')
from matplotlib import pyplot as plt
import logging

from lal import gpstime
import pycbc

import pycbc


parser = argparse.ArgumentParser(usage="",
description="Combine fitting parameters from several different files")
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--verbose", action="store_true",
help="Print extra debugging information", default=False)
parser.add_argument("--combined-fits-file", required=True,
Expand All @@ -45,7 +50,7 @@ parser.add_argument("--colormap", default="rainbow_r", choices=plt.colormaps(),
parser.add_argument("--log-colormap", action='store_true',
help="Use log spacing for choosing colormap values "
"based on duration bins.")
args=parser.parse_args()
args = parser.parse_args()

if '{ifo}' not in args.output_plot_name_format or \
'{type}' not in args.output_plot_name_format:
Expand Down Expand Up @@ -94,21 +99,6 @@ with h5py.File(args.combined_fits_file, 'r') as cff:
bin_starts = bins_edges[:-1]
bin_ends = bins_edges[1:]

bin_max = max(bin_ends)
bin_min = min(bin_starts)

def bin_proportion(upper, lower, log_spacing=False):
if log_spacing:
ll = np.log(lower)
ul = np.log(lower)
centl = (ll + ul) / 2.
minl = np.log(bin_min)
maxl = np.log(bin_max)
return (centl - minl) / (maxl - minl)

else:
return ((lower + upper) / 2. - bin_min) / (bin_max - bin_min)

# Set up the x ticks - note that these are rounded to the nearest
# midnight, so may not line up exactly with the data
min_start = min([separate_starts[ifo].min() for ifo in ifos])
Expand Down Expand Up @@ -157,8 +147,7 @@ for ifo in ifos:
mr = mean_count[ifo][i] / live_total[ifo]
cr = cons_count[ifo][i] / live_total[ifo]

bin_prop = bin_proportion(bu, bl,
log_spacing=args.log_colormap)
bin_prop = i / len(bin_starts)
bin_colour = plt.get_cmap(args.colormap)(bin_prop)
bin_label = f"duration {bl:.2f}-{bu:.2f}"
alpha_lines += ax_alpha.plot(separate_starts[ifo], alphas, c=bin_colour,
Expand Down Expand Up @@ -199,6 +188,12 @@ for ifo in ifos:
ax.grid(zorder=-30)

fig_count.tight_layout()
fig_count.savefig(args.output_plot_name_format.format(ifo=ifo, type='counts'))
fig_count.savefig(
args.output_plot_name_format.format(ifo=ifo, type='counts')
)
fig_alpha.tight_layout()
fig_alpha.savefig(args.output_plot_name_format.format(ifo=ifo, type='fit_coeffs'))
fig_alpha.savefig(
args.output_plot_name_format.format(ifo=ifo, type='fit_coeffs')
)

logging.info("Done")
86 changes: 42 additions & 44 deletions bin/live/pycbc_live_plot_single_trigger_fits
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,17 +12,24 @@
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.

"""Plot histograms of PyCBC Live triggers split over various parameters, and
the corresponding fits.
"""

import argparse
import logging
import numpy as np
import pycbc
from pycbc import bin_utils
from pycbc.events import trigger_fits as trstats
import h5py
import matplotlib
matplotlib.use('agg')
from matplotlib import pyplot as plt
import argparse, logging, h5py

parser = argparse.ArgumentParser(usage="",
description="Plot histograms of triggers split over various parameters")
import pycbc
from pycbc.bin_utils import IrregularBins
from pycbc.events.trigger_fits import cum_fit as eval_cum_fit


parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--verbose", action="store_true",
help="Print extra debugging information", default=False)
parser.add_argument("--trigger-fits-file", required=True,
Expand All @@ -39,11 +46,14 @@ parser.add_argument("--colormap", default="rainbow_r", choices=plt.colormaps(),
parser.add_argument("--log-colormap", action='store_true',
help="Use log spacing for choosing colormap values "
"based on duration bins.")
parser.add_argument("--x-lim-lower", type=float,
help="Add a lower limit to the x-axis of the plot")
parser.add_argument("--x-lim-upper", type=float,
help="Add an upper limit to the x-axis of the plot")
parser.add_argument("--y-lim-lower", type=float,
help="Add a lower limit to the y-axis of the plot")
parser.add_argument("--y-lim-upper", type=float,
help="Add an upper limit to the y-axis of the plot")
#parser.add_argument("--", default="", help="")

#Add some input sanitisation
args = parser.parse_args()
Expand Down Expand Up @@ -80,33 +90,15 @@ with h5py.File(args.trigger_fits_file, 'r') as trfit_f:
bu = trfit_f['bins_upper'][:]
bl = trfit_f['bins_lower'][:]

bin_max = max(bu)
bin_min = min(bl)

def bin_proportion(upper, lower, log_spacing=False):
if log_spacing:
ll = np.log(lower)
ul = np.log(lower)
centl = (ll + ul) / 2.
minl = np.log(bin_min)
maxl = np.log(bin_max)
return (centl - minl) / (maxl - minl)

else:
return ((lower + upper) / 2. - bin_min) / (bin_max - bin_min)


duration_bin_edges = list(bl) + [bu[-1]]

tbins = bin_utils.IrregularBins(duration_bin_edges)
tbins = IrregularBins(duration_bin_edges)

logger = logging.getLogger()
init_level = logger.level

logging.info("Plotting fits")

for ifo in all_ifos:
# Skip if no triggers in this IFO
fig, ax = plt.subplots(1)
oput_plot = args.output_plot_name_format.format(ifo=ifo)

Expand All @@ -133,19 +125,19 @@ for ifo in all_ifos:
maxstat = stats[ifo].max()
max_rate = 0

plotbins = np.linspace(fit_threshold, 1.05 * maxstat)
plotbins = np.linspace(fit_threshold, 1.05 * maxstat, 400)

logging.info("Putting events into bins")
event_bin = np.array([tbins[d] for d in durations[ifo]])

for bin_num, lower_upper in enumerate(zip(duration_bin_edges[:-1],
duration_bin_edges[1:])):
for bin_num, lower_upper in enumerate(zip(tbins.lower(), tbins.upper())):
lower, upper = lower_upper
binlabel = f"{lower:.3g} - {upper:.3g}"

inbin = event_bin == bin_num
# Skip if there are no triggers in this bin in this IFO
if not any(inbin): continue
if not any(inbin):
continue
binned_sngl_stats = stats[ifo][event_bin == bin_num]

# Histogram the triggers
Expand All @@ -155,33 +147,39 @@ for ifo in all_ifos:

max_rate = max(max_rate, cum_rate[0])

cum_fit = counts[ifo][bin_num] / live_time[ifo] * \
trstats.cum_fit(fit_function, plotbins,
alphas[ifo][bin_num], fit_threshold)
ecf = eval_cum_fit(
fit_function,
plotbins,
alphas[ifo][bin_num],
fit_threshold
)
cum_fit = counts[ifo][bin_num] / live_time[ifo] * ecf

# Get the colour from the centre of the bin vs the full range
# of all bins
bin_prop = bin_proportion(upper, lower,
log_spacing=args.log_colormap)
bin_prop = bin_num / len(tbins)
bin_colour = plt.get_cmap(args.colormap)(bin_prop)
ax.plot(edges[:-1], cum_rate, linewidth=2,
color=bin_colour, label=binlabel, alpha=0.6)
ax.plot(plotbins, cum_fit, "--", color=bin_colour,
label=r"$\alpha = $%.2f" % alphas[ifo][bin_num])
ax.semilogy()
ax.grid()
x_upper = args.x_lim_upper or 1.05 * maxstat
y_upper = args.y_lim_upper or 1.5 * max_rate
ax.set_xlim(fit_threshold, x_upper)
ax.set_ylim(0.5 / live_time[ifo], y_upper)
ax.set_xlim(
fit_threshold if args.x_lim_lower is None else args.x_lim_lower,
1.05 * maxstat if args.x_lim_upper is None else args.x_lim_upper
)
ax.set_ylim(
0.5 / live_time[ifo] if args.y_lim_lower is None else args.y_lim_lower,
1.5 * max_rate if args.y_lim_upper is None else args.y_lim_upper
)
ax.set_xlabel(sngl_ranking)
ax.set_ylabel("Number of louder triggers per live time")
title = f"{ifo} {analysis_date} trigger fits"
ax.set_title(title)
ax.legend(loc='upper right')
ax.legend(loc='center left', bbox_to_anchor=(1.01, 0.5))
logging.info(f"Saving {oput_plot}")
# Save initial logging level
logger.setLevel(logging.WARNING)
fig.savefig(oput_plot)
fig.savefig(oput_plot, bbox_inches="tight")
logger.setLevel(init_level)


logging.info("Done")

0 comments on commit 4081be4

Please sign in to comment.