Nr mode power and opt

bin/numrel/pycbc_measure_nr_mode_power

@@ -0,0 +1,219 @@
+#!/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 distutils
+
+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 <[email protected]>"
+__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)
+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)

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

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',