From a85ebf67c085b9206cb920b48042702275b7691c Mon Sep 17 00:00:00 2001 From: Ian Harry Date: Wed, 30 Mar 2016 13:52:13 +0200 Subject: [PATCH 1/3] Add ability to select specific modes + optimization --- pycbc/waveform/nr_waveform.py | 83 ++++++++++++++++++++++++----------- 1 file changed, 58 insertions(+), 25 deletions(-) diff --git a/pycbc/waveform/nr_waveform.py b/pycbc/waveform/nr_waveform.py index 30ea8baa46b..ce4fad82a4c 100644 --- a/pycbc/waveform/nr_waveform.py +++ b/pycbc/waveform/nr_waveform.py @@ -137,10 +137,31 @@ def get_rotation_angles_from_h5_file(filepointer, inclination, phi_ref): return theta, psi, calpha, salpha, source_vecs -def get_hplus_hcross_from_directory(hd5_file_name, template_params, delta_t): +def get_hplus_hcross_from_directory(hd5_file_name, template_params, delta_t, + limit_modes=None, cplx_mode_dict=None): """ Generate hplus, hcross from a NR hd5 file, for a given total mass and f_lower. + + Parameters + ----------- + hd5_file_name : string + Filename of the HDF5 file containing compressed NR data. + template_params : dict or namespace + Object containing values corresponding to the parameters of the + waveform to be generated. + delta_t : float + Timestep to generate waveform at. + limit_modes : list of tuples, default=None + If given (ie. not None) this should be a list of tuples ((2,2), (2,1)) + containing a list of (l,m) modes to generate the waveform with. + cplx_mode_dict : dictionary + **WARNING: DO NOT USE THIS UNLESS YOU REALLY UNDERSTAND WHAT IT DOES!** + If this option is given the Ylm time series will be stored in this + dictionary if they do not already exist. If a mode is in the dictionary + it will be used directly. This can be used if generating a NR waveform + numerous times with the same total mass but different inclination and + phi values. """ def get_param(value): try: @@ -248,35 +269,47 @@ def get_param(value): hp = numpy.zeros(len(time_series), dtype=float) hc = numpy.zeros(len(time_series), dtype=float) + # Identify which modes are needed. + if limit_modes is None: + mode_list = [] + # NOTE: Currently will only search up to l=8 by default. + for l in (2,3,4,5,6,7,8): + for m in range(-l,l+1): + mode_list.append((l,m)) + else: + mode_list = limit_modes + # Generate the waveform - # FIXME: should parse list of (l,m)-pairs - # IWH: Code currently checks for existence of all modes, and includes a - # mode if it is present is this not the right behaviour? If not, - # what is? - for l in (2,3,4,5,6,7,8): - for m in range(-l,l+1): - amp_key = 'amp_l%d_m%d' %(l,m) - phase_key = 'phase_l%d_m%d' %(l,m) - if amp_key not in fp.keys() or phase_key not in fp.keys(): - continue + for l,m in mode_list: + amp_key = 'amp_l%d_m%d' %(l,m) + phase_key = 'phase_l%d_m%d' %(l,m) + if amp_key not in fp.keys() or phase_key not in fp.keys(): + continue + if cplx_mode_dict is not None and \ + amp_key in cplx_mode_dict: + curr_h_real, curr_h_imag = cplx_mode_dict[amp_key] + else: curr_amp = get_data_from_h5_file(fp, time_series_M, amp_key) curr_phase = get_data_from_h5_file(fp, time_series_M, phase_key) curr_h_real = curr_amp * numpy.cos(curr_phase) curr_h_imag = curr_amp * numpy.sin(curr_phase) - curr_ylm = lal.SpinWeightedSphericalHarmonic(theta, psi, -2, l, m) - # Here is what 0709.0093 defines h_+ and h_x as. This defines the - # NR wave frame - curr_hp = curr_h_real * curr_ylm.real - curr_h_imag * curr_ylm.imag - curr_hc = -curr_h_real*curr_ylm.imag - curr_h_imag * curr_ylm.real - - # Correct for the "alpha" angle as given in T1600045 to translate - # from the NR wave frame to LAL wave-frame - hp_corr = (calpha*calpha - salpha*salpha) * curr_hp - hp_corr += 2 * calpha * salpha * curr_hc - hc_corr = - 2 * calpha * salpha * curr_hp - hc_corr += (calpha*calpha - salpha*salpha) * curr_hc - hp += hp_corr - hc += hc_corr + if cplx_mode_dict is not None: + cplx_mode_dict[amp_key] = (curr_h_real, curr_h_imag) + + curr_ylm = lal.SpinWeightedSphericalHarmonic(theta, psi, -2, l, m) + # Here is what 0709.0093 defines h_+ and h_x as. This defines the + # NR wave frame + curr_hp = curr_h_real * curr_ylm.real - curr_h_imag * curr_ylm.imag + curr_hc = -curr_h_real*curr_ylm.imag - curr_h_imag * curr_ylm.real + + # Correct for the "alpha" angle as given in T1600045 to translate + # from the NR wave frame to LAL wave-frame + hp_corr = (calpha*calpha - salpha*salpha) * curr_hp + hp_corr += 2 * calpha * salpha * curr_hc + hc_corr = - 2 * calpha * salpha * curr_hp + hc_corr += (calpha*calpha - salpha*salpha) * curr_hc + hp += hp_corr + hc += hc_corr # Scale by distance # The original NR scaling is 1M. The steps below scale the distance From c6be2c5ffe8ddccecca7f06168eff007115780df Mon Sep 17 00:00:00 2001 From: Ian Harry Date: Wed, 30 Mar 2016 13:52:58 +0200 Subject: [PATCH 2/3] Add a new code to compute power in modes --- bin/numrel/pycbc_measure_nr_mode_power | 215 +++++++++++++++++++++++++ setup.py | 1 + 2 files changed, 216 insertions(+) create mode 100644 bin/numrel/pycbc_measure_nr_mode_power diff --git a/bin/numrel/pycbc_measure_nr_mode_power b/bin/numrel/pycbc_measure_nr_mode_power new file mode 100644 index 00000000000..049a7d8473a --- /dev/null +++ b/bin/numrel/pycbc_measure_nr_mode_power @@ -0,0 +1,215 @@ +#!/usr/bin/env python + +# Copyright (C) 2016 Ian W. Harry, Patricia Schmidt +# +# 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 3 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. + +""" +Measure power in a set of NR modes. +""" +from __future__ import division + +import logging +import h5py +import argparse +import numpy +from numpy import cos, sin + +import pycbc.version +import pycbc.psd +import pycbc.strain +from pycbc import pnutils +from pycbc.waveform import get_hplus_hcross_from_directory +from pycbc.filter import sigma +from pycbc.types import TimeSeries, zeros, real_same_precision_as + + +__author__ = "Ian Harry " +__version__ = pycbc.version.git_verbose_msg +__date__ = pycbc.version.date +__program__ = "pycbc_measure_nr_mode_power" + +def measure_nr_snr(nr_file, total_mass, inclination, phi, polarizations, psds, + f_lower, f_upper, delta_t, modes_list=None, + cplx_mode_dict=None): + """ + Measure the SNR of the NR waveform given the input parameters. + """ + wav_params = {} + fp = h5py.File(nr_file, 'r') + eta = fp.attrs['eta'] + mass1, mass2 = pnutils.mtotal_eta_to_mass1_mass2(total_mass, eta) + wav_params['mass1'] = mass1 + wav_params['mass2'] = mass2 + # These are specified in the LAL frame. The NR code does frame rotations. + wav_params['inclination'] = inclination + wav_params['coa_phase'] = phi + wav_params['f_lower'] = f_lower + # These take arbitrary values + wav_params['end_time'] = 900000000 + wav_params['distance'] = 100. + fp.close() + + try: + hplus, hcross = get_hplus_hcross_from_directory(nr_file, wav_params, + delta_t, limit_modes=modes_list, + cplx_mode_dict=cplx_mode_dict) + except ValueError: + return [-1 for i in range(len(polarizations))] + + # We don't need to worry about sky-location here, so assume polarization + # is a simple function (this is equivalent to source overhead detector) + + # Needs to be binary length + new_length = 2**numpy.ceil(numpy.log2(len(hplus))) + psd_length = int(new_length / 2) + 1 + hoft_padded = TimeSeries(zeros(new_length), delta_t=hplus.delta_t, + dtype=real_same_precision_as(hplus)) + psd = psds[psd_length] + if len(hplus) > len(hoft_padded): + raise ValueError("PSD too short for waveform. You shouldn't see this.") + + snrs = [] + for polarization in polarizations: + hoft = hplus * cos(polarization) + hcross * sin(polarization) + hoft_padded[0:len(hoft)] = hoft + hoft_padded = hoft_padded * pycbc.DYN_RANGE_FAC + + snr = sigma(hoft_padded, psd=psd, low_frequency_cutoff=f_lower, + high_frequency_cutoff=f_upper) + snrs.append(snr) + return snrs + +parser = argparse.ArgumentParser(description=__doc__[1:]) + +# Begin with code specific options +parser.add_argument("--version", action="version", version=__version__) +parser.add_argument("--verbose", action="store_true", default=False, + help="verbose output") +parser.add_argument("--output-file", action="store", required=True, + help="Output file to dump data to.") +parser.add_argument("--nr-file", action="store", required=True, + help="Path to the NR file to use.") +parser.add_argument("--f-lower", action="store", required=True, type=float, + help="Set lower frequency cutoff to this value (in Hz). ") +parser.add_argument("--f-upper", action="store", default=2048., type=float, + help="Set upper frequency cutoff to this value (in Hz). " + "Default=2048Hz.") + +# Insert the PSD options +pycbc.psd.insert_psd_option_group(parser) + +# Insert the data reading options (needed if reading PSD from data) +pycbc.strain.insert_strain_option_group(parser) + +opts = parser.parse_args() + +if opts.verbose: + log_level = logging.DEBUG +else: + log_level = logging.WARN +log_format='%(asctime)s %(message)s' +logging.basicConfig(format=log_format, level=log_level) + +pycbc.psd.verify_psd_options(opts, parser) +if opts.psd_estimation: + # If reading data, check the data reading options + pycbc.strain.verify_strain_options(opts, parser) + +# If we are going to use h(t) to estimate a PSD we need to get h(t) +if opts.psd_estimation: + logging.info("Obtaining h(t) for PSD generation") + # NOTE: I am using DYN_RANGE_FAC. This will apeear in *every* number. + # Be careful with this!! + strain = pycbc.strain.from_cli(opts, pycbc.DYN_RANGE_FAC) +else: + strain = None + +# Get the PSD using the pycbc interface +logging.info("Obtaining PSD") +# Want the number of samples to be a binary number and Nyquist must be above +# opts.f_upper. All this assumes that 1 / deltaF is a binary number +nyquistFreq = 2**numpy.ceil(numpy.log2(opts.f_upper)) +delta_t = 1.0 / (2*nyquistFreq) +psds = {} +for psd_len in [1,2,4,8,16,32,64,128,256]: + delta_f = 1./psd_len + numSamples = int(round(nyquistFreq / delta_f)) + 1 + psd = pycbc.psd.from_cli(opts, length=numSamples, delta_f=delta_f, + low_frequency_cutoff=opts.f_lower, strain=strain, + dyn_range_factor=pycbc.DYN_RANGE_FAC) + psds[numSamples] = psd + + +phi_range = numpy.linspace(0, 2*numpy.pi, 20) +inc_range = numpy.linspace(0, numpy.pi, 20) +polarization_range = numpy.linspace(0, 2*numpy.pi, 20) +mass_range = numpy.arange(30,100.1, 5) + +# Define some mode combinations +l2_modes = [(2,-2),(2,-1),(2,0),(2,1),(2,2)] +# This is same as above but using python list magic +l3_modes = [(3,x) for x in range(-3,4)] +l4_modes = [(4,x) for x in range(-4,5)] +l5_modes = [(5,x) for x in range(-5,6)] + +output_data = [] + +for total_mass in mass_range: + print total_mass + # This *must* be reset this if the total mass changes + cplx_mode_dict = {} + for inclination in inc_range: + for phi in phi_range: + snrs_all = measure_nr_snr(opts.nr_file, total_mass, inclination, + phi, polarization_range, psds, opts.f_lower, opts.f_upper, + delta_t, modes_list=None, cplx_mode_dict=cplx_mode_dict) + snrs_l2 = measure_nr_snr(opts.nr_file, total_mass, inclination, + phi, polarization_range, psds, opts.f_lower, opts.f_upper, + delta_t, modes_list=l2_modes, cplx_mode_dict=cplx_mode_dict) + snrs_l3 = measure_nr_snr(opts.nr_file, total_mass, inclination, + phi, polarization_range, psds, opts.f_lower, opts.f_upper, + delta_t, modes_list=l3_modes, cplx_mode_dict=cplx_mode_dict) + snrs_l4 = measure_nr_snr(opts.nr_file, total_mass, inclination, + phi, polarization_range, psds, opts.f_lower, opts.f_upper, + delta_t, modes_list=l4_modes, cplx_mode_dict=cplx_mode_dict) + snrs_l5 = measure_nr_snr(opts.nr_file, total_mass, inclination, + phi, polarization_range, psds, opts.f_lower, opts.f_upper, + delta_t, modes_list=l5_modes, cplx_mode_dict=cplx_mode_dict) + snrs_22 = measure_nr_snr(opts.nr_file, total_mass, inclination, + phi, polarization_range, psds, opts.f_lower, opts.f_upper, + delta_t, modes_list=[(2,2)], cplx_mode_dict=cplx_mode_dict) + snrs_21 = measure_nr_snr(opts.nr_file, total_mass, inclination, + phi, polarization_range, psds, opts.f_lower, opts.f_upper, + delta_t, modes_list=[(2,1)], cplx_mode_dict=cplx_mode_dict) + snrs_33 = measure_nr_snr(opts.nr_file, total_mass, inclination, + phi, polarization_range, psds, opts.f_lower, opts.f_upper, + delta_t, modes_list=[(3,3)], cplx_mode_dict=cplx_mode_dict) + snrs_44 = measure_nr_snr(opts.nr_file, total_mass, inclination, + phi, polarization_range, psds, opts.f_lower, opts.f_upper, + delta_t, modes_list=[(4,4)], cplx_mode_dict=cplx_mode_dict) + snrs_55 = measure_nr_snr(opts.nr_file, total_mass, inclination, + phi, polarization_range, psds, opts.f_lower, opts.f_upper, + delta_t, modes_list=[(5,5)], cplx_mode_dict=cplx_mode_dict) + for pol_idx in range(len(polarization_range)): + data = [total_mass, inclination, phi, + polarization_range[pol_idx], snrs_all[pol_idx], + snrs_l2[pol_idx], snrs_l3[pol_idx], snrs_l4[pol_idx], + snrs_l5[pol_idx], snrs_22[pol_idx], snrs_21[pol_idx], + snrs_33[pol_idx], snrs_44[pol_idx], snrs_55[pol_idx]] + output_data.append(data) + +output_data = numpy.array(output_data) +numpy.savetxt(opts.output_file, output_data, header="Total mass, Inclination, Phi, Polarization, SNR all modes, SNR l=2 modes, SNR l=3 modes, SNR l=4 modes, SNR l=5 modes, SNR 22 mode, SNR 21 mode, SNR 33 mode, SNR 44 mode, SNR 55 mode") diff --git a/setup.py b/setup.py index 1cbe9ed07f8..87ec34a662f 100644 --- a/setup.py +++ b/setup.py @@ -447,6 +447,7 @@ def run(self): 'bin/hdfcoinc/pycbc_plot_bank_bins', 'bin/pygrb/pygrb_make_offline_workflow', 'bin/pygrb/pygrb_make_summary_page', + 'bin/numrel/pycbc_measure_nr_mode_power', ], packages = [ 'pycbc', From 16f8d4ce6aacec68f68c25305f17f00da3b3c533 Mon Sep 17 00:00:00 2001 From: Ian Harry Date: Wed, 30 Mar 2016 14:24:15 +0200 Subject: [PATCH 3/3] header jwarg not available in older numpy --- bin/numrel/pycbc_measure_nr_mode_power | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/bin/numrel/pycbc_measure_nr_mode_power b/bin/numrel/pycbc_measure_nr_mode_power index 049a7d8473a..68752958659 100644 --- a/bin/numrel/pycbc_measure_nr_mode_power +++ b/bin/numrel/pycbc_measure_nr_mode_power @@ -26,6 +26,7 @@ import h5py import argparse import numpy from numpy import cos, sin +import distutils import pycbc.version import pycbc.psd @@ -212,4 +213,7 @@ for total_mass in mass_range: output_data.append(data) output_data = numpy.array(output_data) -numpy.savetxt(opts.output_file, output_data, header="Total mass, Inclination, Phi, Polarization, SNR all modes, SNR l=2 modes, SNR l=3 modes, SNR l=4 modes, SNR l=5 modes, SNR 22 mode, SNR 21 mode, SNR 33 mode, SNR 44 mode, SNR 55 mode") +if distutils.version.LooseVersion(numpy.version.version) >= distutils.version.LooseVersion('1.7.0'): + numpy.savetxt(opts.output_file, output_data, header="Total mass, Inclination, Phi, Polarization, SNR all modes, SNR l=2 modes, SNR l=3 modes, SNR l=4 modes, SNR l=5 modes, SNR 22 mode, SNR 21 mode, SNR 33 mode, SNR 44 mode, SNR 55 mode") +else: + numpy.savetxt(opts.output_file, output_data)