WIP: Added binary script to generate mock data

bin/pycbc_generate_mock_data

@@ -0,0 +1,191 @@
+#! /usr/bin/env python
+
+# Copyright (C) 2025 Sumit K., Shichao Wu
+#
+# 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.
+__description__ = \
+"""Generate the mock data frame files for a given network of detectors.
+It contains Gaussian noise, population injections, and glitches.
+"""
+
+import os
+import argparse
+import numpy
+
+from pycbc.types import MultiDetOptionAction
+import logging
+
+
+parser = argparse.ArgumentParser()
+
+parser.add_argument("--ifo-list", nargs="+", type=str,
+                        help="List of detector names, e.g. H1 L1 V1",
+                        required=True)
+
+parser.add_argument('--low-frequency-cutoff', nargs='+',
+                        action=MultiDetOptionAction,
+                        metavar='DETECTOR:COLUMN', type=float,
+                        help='For each detector, provide minimum frequency for'
+                        'the data generation.')
+
+
+parser.add_argument('--channel-name', nargs='+', action=MultiDetOptionAction,
+                        metavar='DETECTOR:COLUMN', type=str,
+                        help='For each detector, provide channel names for the'
+                        'frame files.')
+
+parser.add_argument('--psd-model', nargs='+', action=MultiDetOptionAction,
+                        metavar='DETECTOR:COLUMN', type=str,
+                        help='For each detector, provide PSD option for'
+                        'the noise generation.')
+
+parser.add_argument('--fake-strain-seed', nargs='+',
+                        action=MultiDetOptionAction,
+                        metavar='DETECTOR:COLUMN', type=int,
+                        help='For each detector, provide a random number seed')
+
+parser.add_argument('--gps-start-time',
+                        help='GPS start time for frame files',
+                        required=True, type=int)
+
+parser.add_argument('--gps-end-time',
+                        help='GPS end time for frame files',
+                        required=True, type=int)
+
+parser.add_argument('--sample-rate', type=float,
+		        help='Sample rate of the frames to be generated')
+
+parser.add_argument('--injection-file', type=str,
+                        help='Injection file containing parameters of'
+                        'injection population')
+
+parser.add_argument('--output-path',
+                        help="Path for the output frame files",
+                        default='.', type=str)
+
+parser.add_argument('--tag', help='Provide your tag for naming frame file',
+                        default=None)
+
+parser.add_argument('--len-arm', type=float, default=None,
+                        help="The arm length of LISA, in the unit of 'm'")
+
+parser.add_argument('--acc-noise-level', type=float,
+                        help="The level of acceleration noise")
+
+parser.add_argument('--oms-noise-level', type=float,
+                        help="The level of OMS noise")
+
+parser.add_argument('--tdi', type=str,
+                        help="The version of TDI. Choose from '1.5' or '2.0'")
+
+parser.add_argument('--duration', type=float,
+                        help="The duration of observation, between 0 and 10,"
+                             "in the unit of years")
+
+# parse the command line
+opts = parser.parse_args()
+
+
+# Function to generate command line script
+def create_command_opts(ifo, rseed, output_file):
+    if 'LISA' in ifo:
+        detector_tag = 'space'
+    else:
+        detector_tag = 'ground'
+
+    if detector_tag=='ground':
+        args_opt = [
+             'pycbc_condition_strain',
+             '--fake-strain %s'%(opts.psd_model[ifo]),
+             '--fake-strain-seed %d'%rseed,
+             '--sample-rate %d'%(opts.sample_rate),
+             '--gps-start-time %d'%(opts.gps_start_time),
+             '--gps-end-time %d'%(opts.gps_end_time),
+             '--channel-name %s:%s'%(ifo, opts.channel_name[ifo]),
+             '--output-strain-file %s'%(output_file)]
+    elif detector_tag=='space':
+        # Sanity Check
+        if opts.len_arm is None:
+            raise ValueError("Space-borne missions require len-arm")
+        # Build extra args for LISA
+        extra_args = []
+        if opts.len_arm is not None:
+            extra_args.append('len_arm:%s' % opts.len_arm)
+        if opts.acc_noise_level is not None:
+            extra_args.append('acc_noise_level:%s' % opts.acc_noise_level)
+        if opts.oms_noise_level is not None:
+            extra_args.append('oms_noise_level:%s' % opts.oms_noise_level)
+        if opts.tdi is not None:
+            extra_args.append('tdi:%s' % opts.tdi)
+        if opts.duration is not None:
+            extra_args.append('duration:%s' % opts.duration)
+
+        args_opt = [
+            'pycbc_condition_strain',
+            '--fake-strain %s'%(opts.psd_model[ifo]),
+            '--fake-strain-seed %d'%rseed,
+            '--sample-rate %f'%(opts.sample_rate),
+            '--gps-start-time %d'%(opts.gps_start_time),
+            '--gps-end-time %d'%(opts.gps_end_time),
+            '--channel-name %s:%s'%(ifo, opts.channel_name[ifo]),
+            '--output-strain-file %s'%(output_file)
+        ]
+
+        # Add extra args if any
+        if extra_args:
+            args_opt.append('--fake-strain-extra-args %s' % ' '.join(extra_args))
+    return args_opt
+
+
+# Sanity Checks
+if set(opts.ifo_list) != set(opts.low_frequency_cutoff.keys()):
+    raise ValueError("Low frequency cutoff for each detector is required")
+if set(opts.ifo_list) != set(opts.psd_model.keys()):
+    raise ValueError("PSD model for each detector is required")
+
+
+for ifo in opts.ifo_list:
+    # Sanity Checks
+    if opts.fake_strain_seed[ifo] is not None:
+        if set(opts.ifo_list) != set(opts.fake_strain_seed.keys()):
+            raise ValueError("Fake strain seed for each detector is required")
+
+    logging.info("Generating strain for the Detector:%s"%ifo)
+    if opts.fake_strain_seed[ifo] is None:
+        rseed = numpy.random.randint(1e8)
+    else:
+        rseed = opts.fake_strain_seed[ifo]
+
+    time_duration = opts.gps_end_time - opts.gps_start_time
+    if opts.tag is None:
+        output_file = '%s/%s-SIMULATED_STRAIN-%d-%d.gwf'%(
+                    opts.output_path, ifo, opts.gps_start_time, time_duration)
+    else:
+        output_file = '%s/%s-SIMULATED_STRAIN-%s-%d-%d.gwf'%(
+                    opts.output_path, ifo, opts.tag, opts.gps_start_time,
+                    time_duration)
+
+    # Get the base command options
+    args_opt = create_command_opts(ifo, rseed, output_file)
+    if opts.injection_file is not None:
+        args_opt.append('--injection-file %s'%(opts.injection_file))
+    if opts.low_frequency_cutoff is not None:
+        args_opt.append('--fake-strain-flow %f'%(
+                        opts.low_frequency_cutoff[ifo]))
+
+
+    cmd = ' '.join(args_opt)
+    print(cmd)
+    os.system(cmd)

docs/apps.rst

@@ -21,3 +21,4 @@ template banks) should read the documentation at:
    faithsim
    upload_to_gracedb
    pycbc_condition_strain
+   mock_data_generation

docs/mock_data_generation.rst

@@ -0,0 +1,40 @@
+##################################################
+``pycbc_generate_mock_data``: Generating mock data
+##################################################
+
+``pycbc_generate_mock_data`` is a command-line tool which can be used for
+generating mock data for a network of detectors. It also has a functionality
+to choose from a population model, merger rate, and glitch distribution.
+The basic idea of the tool is to do the following operations:
+
+1. Generate noise realizations for a given network of detectors.
+2. Inject signal population in the given noise realizations.
+3. Provide zero noise realizations of signal and glitch population.
+
+The following sections show concrete examples of common use cases.  For a
+detailed list of all options, and option-specific help, see the output of::
+
+    pycbc_generate_mock_data --help
+
+=====================================================
+Generate noise realization for a network of detectors
+=====================================================
+
+The following command generate 16 second of data for a given noise realiztion::
+
+    pycbc_generate_mock_data \
+        --ifo-list H1 L1 V1 \
+        --low-frequency-cutoff H1:20 L1:20 V1:20 \
+        --psd-model H1:aLIGOZeroDetHighPower \
+                    L1:aLIGOZeroDetHighPower \
+                    V1:aLIGOZeroDetHighPower \
+        --fake-strain-seed H1:1234 L1:2345 V1:3456 \
+        --gps-start-time 1257294808 \
+        --gps-end-time 1257294908 \
+        --sample-rate 2048 \
+        --channel-name H1:SIMULATED_STRAIN \
+                       L1:SIMULATED_STRAIN \
+                       V1:SIMULATED_STRAIN \
+        --tag HLV_NOISE
+
+

examples/mdc_generation/injections.ini

@@ -0,0 +1,44 @@
+[variable_params]
+; waveform parameters that will vary in MCMC
+;distance =
+;comoving_volume =
+;inclination =
+polarization =
+ra =
+dec =
+distance = 
+inclination = 
+tc =
+
+[static_params]
+; waveform parameters that will not change in MCMC
+approximant = IMRPhenomXPHM
+f_lower = 4
+f_ref = 4
+mass1 = 39
+mass2 = 32
+
+[prior-ra+dec]
+; sky position prior
+name = uniform_sky
+
+[prior-polarization]
+; polarization prior
+name = uniform_angle
+
+[prior-inclination]
+;inclination prior
+name = sin_angle
+
+[prior-distance]
+;prior distance
+name = uniform_radius
+min-distance = 10
+max-distance = 2000
+
+[prior-tc]
+;prior tc
+name = uniform
+min-tc = 1257294808 
+max-tc = 1257294908
+

examples/mdc_generation/injections_lisa.ini

@@ -0,0 +1,28 @@
+[variable_params]
+
+[static_params]
+; This assumes all those values are in LISA frame.
+; You can set "ref_frame = SSB", but then you should also add it to
+; "static_params" section in PE .ini file.
+ref_frame = LISA
+approximant = BBHX_PhenomD
+; You can use "1.5" or "2.0" for TDI.
+; Please use the same TDI version for PSD and static_params in the PE .ini file.
+tdi = 2.0
+mass1 = 1015522.4376
+mass2 = 796849.1091
+spin1z = 0.597755394865021
+spin2z = 0.36905807298613247
+distance = 17758.367941273442
+inclination = 1.5970175301911231
+coa_phase = 4.275929308696054
+eclipticlongitude = 5.4431083771985165
+eclipticlatitude = -1.2734504596198182
+polarization = 0.22558110042980073
+tc = 4799624.274911478
+t_obs_start = 31536000
+; Put LISA behind the Earth by ~20 degrees.
+t_offset = 7365189.431698299
+f_lower = 1e-4
+f_ref = 1e-4
+f_final = 0.1

examples/mdc_generation/make_injections.sh

@@ -0,0 +1,9 @@
+
+N_INJ=10
+
+pycbc_create_injections \
+       --ninjections $N_INJ \
+       --output-file injection_10inj.hdf\
+       --config-files injections.ini\
+       --force \
+       --seed 233

examples/mdc_generation/make_injections_lisa.sh

@@ -0,0 +1,10 @@
+#!/bin/sh
+pycbc_create_injections --verbose \
+        --config-files injections_lisa.ini \
+        --ninjections 1 \
+        --seed 10 \
+        --output-file injections_lisa.hdf \
+        --variable-params-section variable_params \
+        --static-params-section static_params \
+        --dist-section prior \
+        --force

examples/mdc_generation/run_generate_mdc.sh

@@ -0,0 +1,14 @@
+#!/bin/bash
+
+# Generate injection file
+./make_injections.sh
+
+# Generate Noise without signals frame files
+pycbc_generate_mock_data --ifo-list H1 L1 V1 --low-frequency-cutoff H1:20 L1:20 V1:20 --psd-model H1:aLIGOZeroDetHighPower L1:aLIGOZeroDetHighPower V1:aLIGOZeroDetHighPower --fake-strain-seed H1:1234 L1:2345 V1:3456 --gps-start-time 1257294808 --gps-end-time 1257294908 --sample-rate 2048 --channel-name H1:SIMULATED_STRAIN L1:SIMULATED_STRAIN V1:SIMULATED_STRAIN --tag HLV_NOISE 
+
+# Generate Noise + Signals frame files
+pycbc_generate_mock_data --ifo-list H1 L1 V1 --low-frequency-cutoff H1:20 L1:20 V1:20 --psd-model H1:aLIGOZeroDetHighPower L1:aLIGOZeroDetHighPower V1:aLIGOZeroDetHighPower --fake-strain-seed H1:1234 L1:2345 V1:3456 --gps-start-time 1257294808 --gps-end-time 1257294908 --sample-rate 2048 --channel-name H1:SIMULATED_STRAIN L1:SIMULATED_STRAIN V1:SIMULATED_STRAIN --tag HLV_NOISE_PLUS_SIGNAL --injection-file injection_10inj.hdf
+
+# Generate Signal (with Zeronoise) frame files 
+pycbc_generate_mock_data --ifo-list H1 L1 V1 --low-frequency-cutoff H1:20 L1:20 V1:20 --psd-model H1:zeroNoise L1:zeroNoise V1:zeroNoise --fake-strain-seed H1:1234 L1:2345 V1:3456 --gps-start-time 1257294808 --gps-end-time 1257294908 --sample-rate 2048 --channel-name H1:SIMULATED_STRAIN L1:SIMULATED_STRAIN V1:SIMULATED_STRAIN --tag HLV_ZERONOISE_SIGNAL --injection-file injection_10inj.hdf
+

examples/mdc_generation/run_generate_mdc_lisa.sh

@@ -0,0 +1,36 @@
+#!/bin/bash
+
+# Generate injection file
+./make_injections_lisa.sh
+
+# Generate noise without signals frame files
+pycbc_generate_mock_data --ifo-list LISA_A LISA_E LISA_T \
+                         --low-frequency-cutoff LISA_A:1e-4 LISA_E:1e-4 LISA_T:1e-4 \
+                         --psd-model LISA_A:analytical_psd_lisa_tdi_AE LISA_E:analytical_psd_lisa_tdi_AE LISA_T:analytical_psd_lisa_tdi_T \
+                         --fake-strain-seed LISA_A:100 LISA_E:150 LISA_T:200 \
+                         --gps-start-time 2121224.274911478 --gps-end-time 4799924.274911478 --sample-rate 0.2 \
+                         --len-arm 2.5e9 --acc-noise-level 2.4e-15 --oms-noise-level 7.9e-12 --tdi '2.0' \
+                         --channel-name LISA_A:SIMULATED_STRAIN LISA_E:SIMULATED_STRAIN LISA_T:SIMULATED_STRAIN \
+                         --tag LISA_NOISE
+
+# Generate noise + signals frame files
+pycbc_generate_mock_data --ifo-list LISA_A LISA_E LISA_T \
+                         --low-frequency-cutoff LISA_A:1e-4 LISA_E:1e-4 LISA_T:1e-4 \
+                         --psd-model LISA_A:analytical_psd_lisa_tdi_AE LISA_E:analytical_psd_lisa_tdi_AE LISA_T:analytical_psd_lisa_tdi_T \
+                         --fake-strain-seed LISA_A:100 LISA_E:150 LISA_T:200 \
+                         --gps-start-time 2121224.274911478 --gps-end-time 4799924.274911478 --sample-rate 0.2 \
+                         --len-arm 2.5e9 --acc-noise-level 2.4e-15 --oms-noise-level 7.9e-12 --tdi '2.0' \
+                         --channel-name LISA_A:SIMULATED_STRAIN LISA_E:SIMULATED_STRAIN LISA_T:SIMULATED_STRAIN \
+                         --injection-file injections_lisa.hdf \
+                         --tag LISA_NOISE_PLUS_SIGNAL
+
+# Generate Signal (with Zeronoise) frame files 
+pycbc_generate_mock_data --ifo-list LISA_A LISA_E LISA_T \
+                         --low-frequency-cutoff LISA_A:1e-4 LISA_E:1e-4 LISA_T:1e-4 \
+                         --psd-model LISA_A:zeroNoise LISA_E:zeroNoise LISA_T:zeroNoise \
+                         --fake-strain-seed LISA_A:100 LISA_E:150 LISA_T:200 \
+                         --gps-start-time 2121224.274911478 --gps-end-time 4799924.274911478 --sample-rate 0.2 \
+                         --len-arm 2.5e9 --acc-noise-level 2.4e-15 --oms-noise-level 7.9e-12 --tdi '2.0' \
+                         --channel-name LISA_A:SIMULATED_STRAIN LISA_E:SIMULATED_STRAIN LISA_T:SIMULATED_STRAIN \
+                         --injection-file injections_lisa.hdf \
+                         --tag LISA_ZERONOISE_SIGNAL