def from_cli(opt, dyn_range_fac=1, precision='single',
             inj_filter_rejector=None):
    """Parses the CLI options related to strain data reading and conditioning.

    Parameters
    ----------
    opt : object
        Result of parsing the CLI with OptionParser, or any object with the
        required attributes  (gps-start-time, gps-end-time, strain-high-pass,
        pad-data, sample-rate, (frame-cache or frame-files), channel-name,
        fake-strain, fake-strain-seed, fake-strain-from-file, gating_file).
    dyn_range_fac : {float, 1}, optional
        A large constant to reduce the dynamic range of the strain.
    precision : string
        Precision of the returned strain ('single' or 'double').
    inj_filter_rejector : InjFilterRejector instance; optional, default=None
        If given send the InjFilterRejector instance to the inject module so
        that it can store a reduced representation of injections if
        necessary.

    Returns
    -------
    strain : TimeSeries
        The time series containing the conditioned strain data.
    """
    gating_info = {}

    if opt.frame_cache or opt.frame_files or opt.frame_type:
        if opt.frame_cache:
            frame_source = opt.frame_cache
        if opt.frame_files:
            frame_source = opt.frame_files

        logging.info("Reading Frames")

        if hasattr(opt, 'frame_sieve') and opt.frame_sieve:
            sieve = opt.frame_sieve
        else:
            sieve = None

        if opt.frame_type:
            strain = pycbc.frame.query_and_read_frame(
                    opt.frame_type, opt.channel_name,
                    start_time=opt.gps_start_time-opt.pad_data,
                    end_time=opt.gps_end_time+opt.pad_data,
                    sieve=sieve)
        else:
            strain = pycbc.frame.read_frame(
                    frame_source, opt.channel_name,
                    start_time=opt.gps_start_time-opt.pad_data,
                    end_time=opt.gps_end_time+opt.pad_data,
                    sieve=sieve)

        if opt.zpk_z and opt.zpk_p and opt.zpk_k:
            logging.info("Highpass Filtering")
            strain = highpass(strain, frequency=opt.strain_high_pass)

            logging.info("Applying zpk filter")
            z = numpy.array(opt.zpk_z)
            p = numpy.array(opt.zpk_p)
            k = float(opt.zpk_k)
            strain = filter_zpk(strain.astype(numpy.float64), z, p, k)

        if opt.normalize_strain:
            logging.info("Dividing strain by constant")
            l = opt.normalize_strain
            strain = strain / l

        if opt.injection_file:
            logging.info("Applying injections")
            injector = InjectionSet(opt.injection_file)
            injections = \
                injector.apply(strain, opt.channel_name[0:2],
                               distance_scale=opt.injection_scale_factor,
                               inj_filter_rejector=inj_filter_rejector)

        if opt.sgburst_injection_file:
            logging.info("Applying sine-Gaussian burst injections")
            injector = SGBurstInjectionSet(opt.sgburst_injection_file)
            injector.apply(strain, opt.channel_name[0:2],
                             distance_scale=opt.injection_scale_factor)

        logging.info("Highpass Filtering")
        strain = highpass(strain, frequency=opt.strain_high_pass)

        if precision == 'single':
            logging.info("Converting to float32")
            strain = (strain * dyn_range_fac).astype(pycbc.types.float32)
        elif precision == "double":
            logging.info("Converting to float64")
            strain = (strain * dyn_range_fac).astype(pycbc.types.float64)
        else:
            raise ValueError("Unrecognized precision {}".format(precision))

        if opt.gating_file is not None:
            logging.info("Gating times contained in gating file")
            gate_params = numpy.loadtxt(opt.gating_file)
            if len(gate_params.shape) == 1:
                gate_params = [gate_params]
            strain = gate_data(strain, gate_params)
            gating_info['file'] = \
                    [gp for gp in gate_params \
                     if (gp[0] + gp[1] + gp[2] >= strain.start_time) \
                     and (gp[0] - gp[1] - gp[2] <= strain.end_time)]

        if opt.autogating_threshold is not None:
            # the + 0 is for making a copy
            glitch_times = detect_loud_glitches(
                    strain + 0., threshold=opt.autogating_threshold,
                    cluster_window=opt.autogating_cluster,
                    low_freq_cutoff=opt.strain_high_pass,
                    high_freq_cutoff=opt.sample_rate/2,
                    corrupt_time=opt.pad_data+opt.autogating_pad)
            gate_params = [[gt, opt.autogating_width, opt.autogating_taper] \
                           for gt in glitch_times]
            if len(glitch_times) > 0:
                logging.info('Autogating at %s',
                             ', '.join(['%.3f' % gt for gt in glitch_times]))
            strain = gate_data(strain, gate_params)
            gating_info['auto'] = gate_params

        logging.info("Resampling data")
        strain = resample_to_delta_t(strain, 1.0/opt.sample_rate, method='ldas')

        logging.info("Highpass Filtering")
        strain = highpass(strain, frequency=opt.strain_high_pass)

        if hasattr(opt, 'witness_frame_type') and opt.witness_frame_type:
            stilde = strain.to_frequencyseries()
            import h5py
            tf_file = h5py.File(opt.witness_tf_file)
            for key in tf_file:
                witness = pycbc.frame.query_and_read_frame(opt.witness_frame_type, str(key),
                       start_time=strain.start_time, end_time=strain.end_time)
                witness = (witness * dyn_range_fac).astype(strain.dtype)
                tf = pycbc.types.load_frequencyseries(opt.witness_tf_file, group=key)
                tf = tf.astype(stilde.dtype)

                flen = int(opt.witness_filter_length * strain.sample_rate)
                tf = pycbc.psd.interpolate(tf, stilde.delta_f)

                tf_time = tf.to_timeseries()
                window = Array(numpy.hanning(flen*2), dtype=strain.dtype)
                tf_time[0:flen] *= window[flen:]
                tf_time[len(tf_time)-flen:] *= window[0:flen]
                tf = tf_time.to_frequencyseries()

                kmax = min(len(tf), len(stilde)-1)
                stilde[:kmax] -= tf[:kmax] * witness.to_frequencyseries()[:kmax]

            strain = stilde.to_timeseries()


        logging.info("Remove Padding")
        start = opt.pad_data*opt.sample_rate
        end = len(strain)-opt.sample_rate*opt.pad_data
        strain = strain[start:end]

    if opt.fake_strain or opt.fake_strain_from_file:
        logging.info("Generating Fake Strain")
        if not opt.low_frequency_cutoff:
            raise ValueError('Please provide low frequency cutoff to '
                             'generate a fake strain')
        duration = opt.gps_end_time - opt.gps_start_time
        tlen = duration * opt.sample_rate
        pdf = 1.0/128
        plen = int(opt.sample_rate / pdf) / 2 + 1

        if opt.fake_strain_from_file:
            logging.info("Reading ASD from file")
            strain_psd = pycbc.psd.from_txt(opt.fake_strain_from_file, plen, pdf,
                                            opt.low_frequency_cutoff, is_asd_file=True)
        elif opt.fake_strain != 'zeroNoise':
            logging.info("Making PSD for strain")
            strain_psd = pycbc.psd.from_string(opt.fake_strain, plen, pdf,
                                               opt.low_frequency_cutoff)

        if opt.fake_strain == 'zeroNoise':
            logging.info("Making zero-noise time series")
            strain = TimeSeries(pycbc.types.zeros(tlen),
                                delta_t=1.0/opt.sample_rate,
                                epoch=opt.gps_start_time)
        else:
            logging.info("Making colored noise")
            from pycbc.noise.reproduceable import colored_noise
            lowfreq = opt.low_frequency_cutoff / 2.
            strain = colored_noise(strain_psd, opt.gps_start_time,
                                          opt.gps_end_time,
                                          seed=opt.fake_strain_seed,
                                          low_frequency_cutoff=lowfreq)
            strain = resample_to_delta_t(strain, 1.0/opt.sample_rate)

        if not opt.channel_name and (opt.injection_file \
                                     or opt.sgburst_injection_file):
            raise ValueError('Please provide channel names with the format '
                             'ifo:channel (e.g. H1:CALIB-STRAIN) to inject '
                             'simulated signals into fake strain')

        if opt.injection_file:
            logging.info("Applying injections")
            injector = InjectionSet(opt.injection_file)
            injections = \
                injector.apply(strain, opt.channel_name[0:2],
                               distance_scale=opt.injection_scale_factor,
                               inj_filter_rejector=inj_filter_rejector)

        if opt.sgburst_injection_file:
            logging.info("Applying sine-Gaussian burst injections")
            injector =  SGBurstInjectionSet(opt.sgburst_injection_file)
            injector.apply(strain, opt.channel_name[0:2],
                             distance_scale=opt.injection_scale_factor)

        if precision == 'single':
            logging.info("Converting to float32")
            strain = (dyn_range_fac * strain).astype(pycbc.types.float32)
        elif precision == 'double':
            logging.info("Converting to float64")
            strain = (dyn_range_fac * strain).astype(pycbc.types.float64)
        else:
            raise ValueError("Unrecognized precision {}".format(precision))

    if opt.taper_data:
        logging.info("Tapering data")
        # Use auto-gating stuff for this, a one-sided gate is a taper
        pd_taper_window = opt.taper_data
        gate_params = [(strain.start_time, 0., pd_taper_window)]
        gate_params.append( (strain.end_time, 0.,
                             pd_taper_window) )
        gate_data(strain, gate_params)

    if opt.injection_file:
        strain.injections = injections
    strain.gating_info = gating_info

    return strain
def from_cli ( opt , dyn_range_fac = 1 , precision = 'single' , inj_filter_rejector = None ) : gating_info = { } if opt . frame_cache or opt . frame_files or opt . frame_type : if opt . frame_cache : frame_source = opt . frame_cache if opt . frame_files : frame_source = opt . frame_files logging . info ( "Reading Frames" ) if hasattr ( opt , 'frame_sieve' ) and opt . frame_sieve : sieve = opt . frame_sieve else : sieve = None if opt . frame_type : strain = pycbc . frame . query_and_read_frame ( opt . frame_type , opt . channel_name , start_time = opt . gps_start_time - opt . pad_data , end_time = opt . gps_end_time + opt . pad_data , sieve = sieve ) else : strain = pycbc . frame . read_frame ( frame_source , opt . channel_name , start_time = opt . gps_start_time - opt . pad_data , end_time = opt . gps_end_time + opt . pad_data , sieve = sieve ) if opt . zpk_z and opt . zpk_p and opt . zpk_k : logging . info ( "Highpass Filtering" ) strain = highpass ( strain , frequency = opt . strain_high_pass ) logging . info ( "Applying zpk filter" ) z = numpy . array ( opt . zpk_z ) p = numpy . array ( opt . zpk_p ) k = float ( opt . zpk_k ) strain = filter_zpk ( strain . astype ( numpy . float64 ) , z , p , k ) if opt . normalize_strain : logging . info ( "Dividing strain by constant" ) l = opt . normalize_strain strain = strain / l if opt . injection_file : logging . info ( "Applying injections" ) injector = InjectionSet ( opt . injection_file ) injections = injector . apply ( strain , opt . channel_name [ 0 : 2 ] , distance_scale = opt . injection_scale_factor , inj_filter_rejector = inj_filter_rejector ) if opt . sgburst_injection_file : logging . info ( "Applying sine-Gaussian burst injections" ) injector = SGBurstInjectionSet ( opt . sgburst_injection_file ) injector . apply ( strain , opt . channel_name [ 0 : 2 ] , distance_scale = opt . injection_scale_factor ) logging . info ( "Highpass Filtering" ) strain = highpass ( strain , frequency = opt . strain_high_pass ) if precision == 'single' : logging . info ( "Converting to float32" ) strain = ( strain * dyn_range_fac ) . astype ( pycbc . types . float32 ) elif precision == "double" : logging . info ( "Converting to float64" ) strain = ( strain * dyn_range_fac ) . astype ( pycbc . types . float64 ) else : raise ValueError ( "Unrecognized precision {}" . format ( precision ) ) if opt . gating_file is not None : logging . info ( "Gating times contained in gating file" ) gate_params = numpy . loadtxt ( opt . gating_file ) if len ( gate_params . shape ) == 1 : gate_params = [ gate_params ] strain = gate_data ( strain , gate_params ) gating_info [ 'file' ] = [ gp for gp in gate_params if ( gp [ 0 ] + gp [ 1 ] + gp [ 2 ] >= strain . start_time ) and ( gp [ 0 ] - gp [ 1 ] - gp [ 2 ] <= strain . end_time ) ] if opt . autogating_threshold is not None : glitch_times = detect_loud_glitches ( strain + 0. , threshold = opt . autogating_threshold , cluster_window = opt . autogating_cluster , low_freq_cutoff = opt . strain_high_pass , high_freq_cutoff = opt . sample_rate / 2 , corrupt_time = opt . pad_data + opt . autogating_pad ) gate_params = [ [ gt , opt . autogating_width , opt . autogating_taper ] for gt in glitch_times ] if len ( glitch_times ) > 0 : logging . info ( 'Autogating at %s' , ', ' . join ( [ '%.3f' % gt for gt in glitch_times ] ) ) strain = gate_data ( strain , gate_params ) gating_info [ 'auto' ] = gate_params logging . info ( "Resampling data" ) strain = resample_to_delta_t ( strain , 1.0 / opt . sample_rate , method = 'ldas' ) logging . info ( "Highpass Filtering" ) strain = highpass ( strain , frequency = opt . strain_high_pass ) if hasattr ( opt , 'witness_frame_type' ) and opt . witness_frame_type : stilde = strain . to_frequencyseries ( ) import h5py tf_file = h5py . File ( opt . witness_tf_file ) for key in tf_file : witness = pycbc . frame . query_and_read_frame ( opt . witness_frame_type , str ( key ) , start_time = strain . start_time , end_time = strain . end_time ) witness = ( witness * dyn_range_fac ) . astype ( strain . dtype ) tf = pycbc . types . load_frequencyseries ( opt . witness_tf_file , group = key ) tf = tf . astype ( stilde . dtype ) flen = int ( opt . witness_filter_length * strain . sample_rate ) tf = pycbc . psd . interpolate ( tf , stilde . delta_f ) tf_time = tf . to_timeseries ( ) window = Array ( numpy . hanning ( flen * 2 ) , dtype = strain . dtype ) tf_time [ 0 : flen ] *= window [ flen : ] tf_time [ len ( tf_time ) - flen : ] *= window [ 0 : flen ] tf = tf_time . to_frequencyseries ( ) kmax = min ( len ( tf ) , len ( stilde ) - 1 ) stilde [ : kmax ] -= tf [ : kmax ] * witness . to_frequencyseries ( ) [ : kmax ] strain = stilde . to_timeseries ( ) logging . info ( "Remove Padding" ) start = opt . pad_data * opt . sample_rate end = len ( strain ) - opt . sample_rate * opt . pad_data strain = strain [ start : end ] if opt . fake_strain or opt . fake_strain_from_file : logging . info ( "Generating Fake Strain" ) if not opt . low_frequency_cutoff : raise ValueError ( 'Please provide low frequency cutoff to ' 'generate a fake strain' ) duration = opt . gps_end_time - opt . gps_start_time tlen = duration * opt . sample_rate pdf = 1.0 / 128 plen = int ( opt . sample_rate / pdf ) / 2 + 1 if opt . fake_strain_from_file : logging . info ( "Reading ASD from file" ) strain_psd = pycbc . psd . from_txt ( opt . fake_strain_from_file , plen , pdf , opt . low_frequency_cutoff , is_asd_file = True ) elif opt . fake_strain != 'zeroNoise' : logging . info ( "Making PSD for strain" ) strain_psd = pycbc . psd . from_string ( opt . fake_strain , plen , pdf , opt . low_frequency_cutoff ) if opt . fake_strain == 'zeroNoise' : logging . info ( "Making zero-noise time series" ) strain = TimeSeries ( pycbc . types . zeros ( tlen ) , delta_t = 1.0 / opt . sample_rate , epoch = opt . gps_start_time ) else : logging . info ( "Making colored noise" ) from pycbc . noise . reproduceable import colored_noise lowfreq = opt . low_frequency_cutoff / 2. strain = colored_noise ( strain_psd , opt . gps_start_time , opt . gps_end_time , seed = opt . fake_strain_seed , low_frequency_cutoff = lowfreq ) strain = resample_to_delta_t ( strain , 1.0 / opt . sample_rate ) if not opt . channel_name and ( opt . injection_file or opt . sgburst_injection_file ) : raise ValueError ( 'Please provide channel names with the format ' 'ifo:channel (e.g. H1:CALIB-STRAIN) to inject ' 'simulated signals into fake strain' ) if opt . injection_file : logging . info ( "Applying injections" ) injector = InjectionSet ( opt . injection_file ) injections = injector . apply ( strain , opt . channel_name [ 0 : 2 ] , distance_scale = opt . injection_scale_factor , inj_filter_rejector = inj_filter_rejector ) if opt . sgburst_injection_file : logging . info ( "Applying sine-Gaussian burst injections" ) injector = SGBurstInjectionSet ( opt . sgburst_injection_file ) injector . apply ( strain , opt . channel_name [ 0 : 2 ] , distance_scale = opt . injection_scale_factor ) if precision == 'single' : logging . info ( "Converting to float32" ) strain = ( dyn_range_fac * strain ) . astype ( pycbc . types . float32 ) elif precision == 'double' : logging . info ( "Converting to float64" ) strain = ( dyn_range_fac * strain ) . astype ( pycbc . types . float64 ) else : raise ValueError ( "Unrecognized precision {}" . format ( precision ) ) if opt . taper_data : logging . info ( "Tapering data" ) pd_taper_window = opt . taper_data gate_params = [ ( strain . start_time , 0. , pd_taper_window ) ] gate_params . append ( ( strain . end_time , 0. , pd_taper_window ) ) gate_data ( strain , gate_params ) if opt . injection_file : strain . injections = injections strain . gating_info = gating_info return strain
Parses the CLI options related to strain data reading and conditioning.

    Parameters
    ----------
    opt : object
        Result of parsing the CLI with OptionParser, or any object with the
        required attributes  (gps-start-time, gps-end-time, strain-high-pass,
        pad-data, sample-rate, (frame-cache or frame-files), channel-name,
        fake-strain, fake-strain-seed, fake-strain-from-file, gating_file).
    dyn_range_fac : {float, 1}, optional
        A large constant to reduce the dynamic range of the strain.
    precision : string
        Precision of the returned strain ('single' or 'double').
    inj_filter_rejector : InjFilterRejector instance; optional, default=None
        If given send the InjFilterRejector instance to the inject module so
        that it can store a reduced representation of injections if
        necessary.

    Returns
    -------
    strain : TimeSeries
        The time series containing the conditioned strain data.