From b354669abadbbb5c0a708ccb96af6519ccf7d9a0 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Mon, 23 Dec 2024 17:26:16 +0100 Subject: [PATCH 1/4] modify triggers get par --- bin/all_sky_search/pycbc_fit_sngls_binned | 17 +++++++---------- bin/all_sky_search/pycbc_fit_sngls_by_template | 8 +++----- .../pycbc_fit_sngls_over_multiparam | 3 +-- pycbc/events/triggers.py | 14 +++++++++----- 4 files changed, 20 insertions(+), 22 deletions(-) diff --git a/bin/all_sky_search/pycbc_fit_sngls_binned b/bin/all_sky_search/pycbc_fit_sngls_binned index 126ed493fc3..00d1ca53b53 100644 --- a/bin/all_sky_search/pycbc_fit_sngls_binned +++ b/bin/all_sky_search/pycbc_fit_sngls_binned @@ -184,21 +184,20 @@ for veto_file, veto_segment_name in zip(args.veto_file, args.veto_segment_name): ### Functions for doing the pruning (removal of trigs at loudest times) -def get_pars(args, tag, m1, m2, s1z, s2z): +def get_pars(args, tag, bank, tid): # here used for both pruning and binning params paramarg = getattr(args, tag+'_param') try: # will fail if m1 is a float rather than a sequence - logging.info('Getting %s values for %i triggers' % (paramarg, len(m1))) + logging.info('Getting %s values for %i triggers' % (paramarg, len(bank['mass1']))) except: pass - return triggers.get_param(paramarg, args, m1, m2, s1z, s2z) + return triggers.get_param(paramarg, args, bank, tid) if args.prune_param: logging.info('Getting min and max param values') - prpars = get_pars(args, 'prune', - templatef['mass1'][:], templatef['mass2'][:], - templatef['spin1z'][:], templatef['spin2z'][:]) + templateid = np.arange(len(templatef['template_hash'][:])) + prpars = get_pars(args, 'prune', templatef, templateid) minprpar = min(prpars) maxprpar = max(prpars) del prpars @@ -232,8 +231,7 @@ if args.prune_param: lstat = statpruneall[loudest] ltid = tidpruneall[loudest] ltime = timepruneall[loudest] - m1, m2, s1z, s2z = triggers.get_mass_spin(templatef, ltid) - lbin = trstats.which_bin(get_pars(args, 'prune', m1, m2, s1z, s2z), + lbin = trstats.which_bin(get_pars(args, 'prune', templatef, ltid), minprpar, maxprpar, args.prune_bins, log=args.log_prune_param) # is the bin where the loudest trigger lives full already? @@ -275,8 +273,7 @@ if args.prune_param: if trig_dur: binpars = tdur + args.min_duration else: - m1, m2, s1z, s2z = triggers.get_mass_spin(templatef, tid) - binpars = get_pars(args, 'bin', m1, m2, s1z, s2z) + binpars = get_pars(args, 'bin', templatef, tid) logging.info("Parameter range of triggers: %f - %f" % (min(binpars), max(binpars))) diff --git a/bin/all_sky_search/pycbc_fit_sngls_by_template b/bin/all_sky_search/pycbc_fit_sngls_by_template index c7816bb7e23..d666725efd1 100755 --- a/bin/all_sky_search/pycbc_fit_sngls_by_template +++ b/bin/all_sky_search/pycbc_fit_sngls_by_template @@ -261,9 +261,8 @@ total_time = abs(all_segments) # do pruning (removal of trigs at N loudest times defined over param bins) if args.prune_param: logging.info('Getting min and max param values') - pars = triggers.get_param(args.prune_param, args, - templatef['mass1'][:], templatef['mass2'][:], - templatef['spin1z'][:], templatef['spin2z'][:]) + templateid = np.arange(len(templatef['template_hash'][:])) + pars = triggers.get_param(args.prune_param, args, templatef, templateid) minpar = min(pars) maxpar = max(pars) del pars @@ -297,9 +296,8 @@ if args.prune_param: lstat = statpruneall[loudest] ltid = tidpruneall[loudest] ltime = timepruneall[loudest] - m1, m2, s1z, s2z = triggers.get_mass_spin(templatef, ltid) lbin = trstats.which_bin(triggers.get_param(args.prune_param, args, - m1, m2, s1z, s2z), + templatef, ltid), minpar, maxpar, args.prune_bins, log=args.log_prune_param) # is the bin where the loudest trigger lives full already? diff --git a/bin/all_sky_search/pycbc_fit_sngls_over_multiparam b/bin/all_sky_search/pycbc_fit_sngls_over_multiparam index 169dcf44d9e..876db56861d 100755 --- a/bin/all_sky_search/pycbc_fit_sngls_over_multiparam +++ b/bin/all_sky_search/pycbc_fit_sngls_over_multiparam @@ -388,13 +388,12 @@ if args.smoothing_method == 'n_closest' and n_required > nabove.sum(): logging.info('Calculating template parameter values') bank = HFile(args.bank_file, 'r') -m1, m2, s1z, s2z = triggers.get_mass_spin(bank, tid) parvals = [] parnames = [] for param, slog in zip(args.fit_param, args.log_param): - data = triggers.get_param(param, args, m1, m2, s1z, s2z) + data = triggers.get_param(param, args, bank, tid) if slog in ['false', 'False', 'FALSE']: logging.info('Using param: %s', param) parvals.append(data) diff --git a/pycbc/events/triggers.py b/pycbc/events/triggers.py index c21cda49c5a..6b4a04423e0 100644 --- a/pycbc/events/triggers.py +++ b/pycbc/events/triggers.py @@ -104,7 +104,7 @@ def get_mass_spin(bank, tid): return m1, m2, s1z, s2z -def get_param(par, args, m1, m2, s1z, s2z): +def get_param(par, args, bank, tid): """ Helper function @@ -114,14 +114,16 @@ def get_param(par, args, m1, m2, s1z, s2z): Name of parameter to calculate args : Namespace object returned from ArgumentParser instance Calling code command line options, used for f_lower value - m1 : float or array of floats - First binary component mass (etc.) + bank: + + tid: Returns ------- parvals : float or array of floats Calculated parameter values """ + m1, m2, s1z, s2z = get_mass_spin(bank, tid) if par == 'mchirp': parvals = conversions.mchirp_from_mass1_mass2(m1, m2) elif par == 'mtotal': @@ -131,10 +133,12 @@ def get_param(par, args, m1, m2, s1z, s2z): elif par in ['chi_eff', 'effective_spin']: parvals = conversions.chi_eff(m1, m2, s1z, s2z) elif par == 'template_duration': + if 'template_duration' in bank: + parvals = bank['template_duration'][:][tid] # default to SEOBNRv4 duration function - if not hasattr(args, 'approximant') or args.approximant is None: + elif not hasattr(args, 'approximant') or args.approximant is None: args.approximant = "SEOBNRv4" - parvals = pnutils.get_imr_duration(m1, m2, s1z, s2z, args.f_lower, + parvals = pnutils.get_imr_duration(m1, m2, s1z, s2z, args.f_lower, args.approximant) if args.min_duration: parvals += args.min_duration From abd1ae867579c985cee6936d381dbb246dc38267 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Sat, 15 Feb 2025 00:22:05 +0100 Subject: [PATCH 2/4] add eccentricity to smooth over --- pycbc/events/triggers.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pycbc/events/triggers.py b/pycbc/events/triggers.py index 6b4a04423e0..ac1f17a37cd 100644 --- a/pycbc/events/triggers.py +++ b/pycbc/events/triggers.py @@ -130,8 +130,12 @@ def get_param(par, args, bank, tid): parvals = m1 + m2 elif par == 'eta': parvals = conversions.eta_from_mass1_mass2(m1, m2) + elif par == 'q': + parvals = conversions.q_from_mass1_mass2(m1, m2) elif par in ['chi_eff', 'effective_spin']: parvals = conversions.chi_eff(m1, m2, s1z, s2z) + elif par == 'eccentricity': + parvals = bank['eccentricity'][:][tid] elif par == 'template_duration': if 'template_duration' in bank: parvals = bank['template_duration'][:][tid] From 7f1a64c440f85220e09ccea560164aef69914758 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Fri, 7 Mar 2025 19:08:34 +0100 Subject: [PATCH 3/4] check fix over parameters not having 0 --- bin/all_sky_search/pycbc_fit_sngls_over_multiparam | 3 +++ 1 file changed, 3 insertions(+) diff --git a/bin/all_sky_search/pycbc_fit_sngls_over_multiparam b/bin/all_sky_search/pycbc_fit_sngls_over_multiparam index 876db56861d..1a41dcb25f9 100755 --- a/bin/all_sky_search/pycbc_fit_sngls_over_multiparam +++ b/bin/all_sky_search/pycbc_fit_sngls_over_multiparam @@ -501,6 +501,9 @@ else: **kwarg_dict ) +if not isinstance(nabove_smoothed, numpy.ndarray): + nabove_smoothed = numpy.array(nabove_smoothed) +assert all(nabove_smoothed > 0), "Some nabove_smoothed values are not above zero" logging.info("Writing output") outfile = HFile(args.output, 'w') outfile['template_id'] = tid From c411d750cdb4f186d03fe2fe6037afe015b6ec19 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Tue, 14 Apr 2026 13:24:44 +0200 Subject: [PATCH 4/4] clean up and polish --- bin/all_sky_search/pycbc_fit_sngls_binned | 4 ++-- bin/all_sky_search/pycbc_fit_sngls_by_template | 2 +- .../pycbc_fit_sngls_over_multiparam | 12 +++++++++--- pycbc/events/triggers.py | 16 ++++++++-------- 4 files changed, 20 insertions(+), 14 deletions(-) diff --git a/bin/all_sky_search/pycbc_fit_sngls_binned b/bin/all_sky_search/pycbc_fit_sngls_binned index 00d1ca53b53..e75adc8579a 100644 --- a/bin/all_sky_search/pycbc_fit_sngls_binned +++ b/bin/all_sky_search/pycbc_fit_sngls_binned @@ -188,7 +188,7 @@ def get_pars(args, tag, bank, tid): # here used for both pruning and binning params paramarg = getattr(args, tag+'_param') try: - # will fail if m1 is a float rather than a sequence + # will fail if mass1 is not in the bank logging.info('Getting %s values for %i triggers' % (paramarg, len(bank['mass1']))) except: pass @@ -196,7 +196,7 @@ def get_pars(args, tag, bank, tid): if args.prune_param: logging.info('Getting min and max param values') - templateid = np.arange(len(templatef['template_hash'][:])) + templateid = np.arange(len(templatef['mass1'])) # get the full range prpars = get_pars(args, 'prune', templatef, templateid) minprpar = min(prpars) maxprpar = max(prpars) diff --git a/bin/all_sky_search/pycbc_fit_sngls_by_template b/bin/all_sky_search/pycbc_fit_sngls_by_template index d666725efd1..edd7c47cfd7 100755 --- a/bin/all_sky_search/pycbc_fit_sngls_by_template +++ b/bin/all_sky_search/pycbc_fit_sngls_by_template @@ -261,7 +261,7 @@ total_time = abs(all_segments) # do pruning (removal of trigs at N loudest times defined over param bins) if args.prune_param: logging.info('Getting min and max param values') - templateid = np.arange(len(templatef['template_hash'][:])) + templateid = np.arange(len(templatef['mass1'])) # get the full range pars = triggers.get_param(args.prune_param, args, templatef, templateid) minpar = min(pars) maxpar = max(pars) diff --git a/bin/all_sky_search/pycbc_fit_sngls_over_multiparam b/bin/all_sky_search/pycbc_fit_sngls_over_multiparam index 1a41dcb25f9..3a8611eca6b 100755 --- a/bin/all_sky_search/pycbc_fit_sngls_over_multiparam +++ b/bin/all_sky_search/pycbc_fit_sngls_over_multiparam @@ -193,6 +193,7 @@ parser.add_argument("--fit-param", nargs='+', help="Parameter(s) over which to regress the background " "fit coefficients. Required. Either read from " "template fit file or choose from mchirp, mtotal, " + "q, eccentricity, " "chi_eff, eta, tau_0, tau_3, template_duration, " "a frequency cutoff in pnutils or a frequency function" "in LALSimulation. To regress the background over " @@ -501,9 +502,14 @@ else: **kwarg_dict ) -if not isinstance(nabove_smoothed, numpy.ndarray): - nabove_smoothed = numpy.array(nabove_smoothed) -assert all(nabove_smoothed > 0), "Some nabove_smoothed values are not above zero" +non_positive_mask = smoothed_vals[:, 0] <= 0 +if numpy.any(non_positive_mask): + n_bad = numpy.sum(non_positive_mask) + raise ValueError( + f"{n_bad} template(s) have zero triggers above threshold. " + "This will cause problems for the fit coefficient values, and is likely " + "due to too small smoothing width or too large fitting threshold. " + ) logging.info("Writing output") outfile = HFile(args.output, 'w') outfile['template_id'] = tid diff --git a/pycbc/events/triggers.py b/pycbc/events/triggers.py index ac1f17a37cd..d3a7779dd3f 100644 --- a/pycbc/events/triggers.py +++ b/pycbc/events/triggers.py @@ -106,7 +106,8 @@ def get_mass_spin(bank, tid): def get_param(par, args, bank, tid): """ - Helper function + Helper function to extract parameters from bank and calculate + derived parameters Parameters ---------- @@ -114,9 +115,10 @@ def get_param(par, args, bank, tid): Name of parameter to calculate args : Namespace object returned from ArgumentParser instance Calling code command line options, used for f_lower value - bank: - - tid: + bank : h5py File object + Bank parameter file + tid : integer or array of int + Indices of the entries to be returned Returns ------- @@ -137,12 +139,10 @@ def get_param(par, args, bank, tid): elif par == 'eccentricity': parvals = bank['eccentricity'][:][tid] elif par == 'template_duration': - if 'template_duration' in bank: - parvals = bank['template_duration'][:][tid] # default to SEOBNRv4 duration function - elif not hasattr(args, 'approximant') or args.approximant is None: + if not hasattr(args, 'approximant') or args.approximant is None: args.approximant = "SEOBNRv4" - parvals = pnutils.get_imr_duration(m1, m2, s1z, s2z, args.f_lower, + parvals = pnutils.get_imr_duration(m1, m2, s1z, s2z, args.f_lower, args.approximant) if args.min_duration: parvals += args.min_duration