From 2a3b8e7c6927b0d39a552bb3f1f4d615505bf6e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marlin=20Sch=C3=A4fer?= Date: Thu, 26 May 2022 11:14:37 +0200 Subject: [PATCH] Fixed bug where injections could be counted multiple times. Changes: -Fixed bug that would cause injections to be counted multiple times if there were multiple true-positive foreground events for these injections. --- evaluate.py | 46 +++++++++++++++++++++++++++++++++++++--------- 1 file changed, 37 insertions(+), 9 deletions(-) diff --git a/evaluate.py b/evaluate.py index d3cc9fb..4a4fc0b 100755 --- a/evaluate.py +++ b/evaluate.py @@ -7,6 +7,7 @@ import h5py import os import logging +from tqdm import tqdm def find_injection_times(fgfiles, injfile, padding_start=0, padding_end=0): @@ -153,11 +154,14 @@ def get_stats(fgevents, bgevents, injparams, duration=None, logging.info('Finding injection times closest to event times') idxs = find_closest_index(injtimes, fgevents[0]) + # print(idxs) diff = np.abs(injtimes[idxs] - fgevents[0]) logging.info('Finding true- and false-positives') - tpidxs = np.arange(len(fgevents[0]))[diff <= fgevents[2]] - fpidxs = np.arange(len(fgevents[0]))[diff > fgevents[2]] + tpbidxs = diff <= fgevents[2] + tpidxs = np.arange(len(fgevents[0]))[tpbidxs] + fpbidxs = diff > fgevents[2] + fpidxs = np.arange(len(fgevents[0]))[fpbidxs] tpevents = fgevents.T[tpidxs].T fpevents = fgevents.T[fpidxs].T @@ -190,14 +194,37 @@ def get_stats(fgevents, bgevents, injparams, duration=None, far = far / duration ret['far'] = far + # Find best true-positive for each injection + verbose = logging.root.level is logging.INFO + found_injections = [] + tmpsidxs = idxs.argsort() + sorted_idxs = idxs[tmpsidxs] + iidxs = np.full(len(idxs), False) + for i in tqdm(range(len(injtimes)), ascii=True, disable=not verbose, + desc='Determining found injections'): + L = np.searchsorted(sorted_idxs, i, side='left') + if L >= len(idxs) or sorted_idxs[L] != i: + continue + R = np.searchsorted(sorted_idxs, i, side='right') + # All indices that point to the same injection + iidxs[tmpsidxs[L:R]] = True + # Indices of the true-positives that belong to the same injection + eidxs = np.logical_and(iidxs[tmpsidxs[L:R]], + tpbidxs[tmpsidxs[L:R]]) + if eidxs.any(): + found_injections.append([i, + np.max(fgevents[1][tmpsidxs[L:R]][eidxs])]) + iidxs[tmpsidxs[L:R]] = False + + found_injections = np.array(found_injections).T + # Calculate sensitivity # CARE! THIS APPLIES ONLY WHEN THE DISTRIBUTION IS CHOSEN CORRECTLY logging.info('Calculating sensitivity') - sidxs = tpevents[1].argsort() - tp_sort = tpevents[1][sidxs] + sidxs = found_injections[1].argsort() + found_injections = found_injections.T[sidxs].T # Sort found injections if chirp_distance: - found_mchirp_total = massc[idxs[tpidxs]] - found_mchirp_total = found_mchirp_total[sidxs] + found_mchirp_total = massc[found_injections[0].astype(int)] mchirp_max = massc.max() max_distance = dist.max() vtot = (4. / 3.) * np.pi * max_distance**3. @@ -208,11 +235,12 @@ def get_stats(fgevents, bgevents, injparams, duration=None, mc_norm = Ninj prefactor = vtot / mc_norm - nfound = len(tp_sort) - np.searchsorted(tp_sort, noise_stats, - side='right') + nfound = len(found_injections[1]) - np.searchsorted(found_injections[1], + noise_stats, + side='right') if chirp_distance: # Get found chirp-mass indices for given threshold - fidxs = np.searchsorted(tp_sort, noise_stats, side='right') + fidxs = np.searchsorted(found_injections[1], noise_stats, side='right') found_mchirp_total = np.flip(found_mchirp_total) # Calculate sum(found_mchirp ** (5/2))