Glitch robust search

This is an example of running a glitch-robust search: a search which allows the signal to undergo one or more jumps in the frequency and frequency derivatives.

Generating the data

In order to demonstrate the method working, we will simulate a data set containing a glitching signal. This can be done be running the example script glitch_robust_search_make_simulated_data.py. This script contains some setup of times and signal properties, then uses the Writer to generate the simulated data:

import numpy as np
import pyfstat

outdir = 'data'
label = 'simulated_glitching_signal'

# Properties of the GW data
tstart = 1000000000
Tspan = 60 * 86400

tref = tstart + .5 * Tspan

# Fixed properties of the signal
F0s = 30
F1s = -1e-8
F2s = 0
Alpha = np.radians(83.6292)
Delta = np.radians(22.0144)
h0 = 1e-25
sqrtSX = 1e-24
psi = -0.1
phi = 0
cosi = 0.5

# Glitch properties
dtglitch = 0.45 * Tspan  # time (in secs) after minStartTime
dF0 = 5e-6
dF1 = 1e-12


detectors = 'H1'

glitch_data = pyfstat.Writer(
    label=label, outdir=outdir, tref=tref, tstart=tstart,
    F0=F0s, F1=F1s, F2=F2s, duration=Tspan, Alpha=Alpha,
    Delta=Delta, sqrtSX=sqrtSX, dtglitch=dtglitch,
    h0=h0, cosi=cosi, phi=phi, psi=psi,
    delta_F0=dF0, delta_F1=dF1, add_noise=True)

glitch_data.make_data()

Running this script, e.g.

$ python glitch_robust_search_make_simulated_data.py

will produce data/simulated_glitching_signal.cff (the config file passed to lalapps_Makefakedata_v5) and an sft file (also in the data directory).

Running the search

The script examples/glitch_robust_search.py contains the basic setup and commands to run the glitch-robust search for a single glitch:

import numpy as np
import pyfstat

outdir = 'data'
label = 'glitch_robust_search'

# Properties of the GW data
tstart = 1000000000
Tspan = 60 * 86400

# Fixed properties of the signal
F0s = 30
F1s = -1e-8
F2s = 0
Alpha = np.radians(83.6292)
Delta = np.radians(22.0144)

tref = tstart + .5 * Tspan

sftfilepath = 'data/*glitching_signal*sft'

F0_width = np.sqrt(3)/(np.pi*Tspan)
F1_width = np.sqrt(45/4.)/(np.pi*Tspan**2)
DeltaF0 = 50 * F0_width
DeltaF1 = 50 * F1_width

theta_prior = {'F0': {'type': 'unif',
                      'lower': F0s-DeltaF0,
                      'upper': F0s+DeltaF0},
               'F1': {'type': 'unif',
                      'lower': F1s-DeltaF1,
                      'upper': F1s+DeltaF1},
               'F2': F2s,
               'delta_F0': {'type': 'unif',
                            'lower': 0,
                            'upper': 1e-5},
               'delta_F1': {'type': 'unif',
                            'lower': -1e-11,
                            'upper': 1e-11},
               'tglitch': {'type': 'unif',
                           'lower': tstart+0.1*Tspan,
                           'upper': tstart+0.9*Tspan},
               'Alpha': Alpha,
               'Delta': Delta,
               }

search = pyfstat.MCMCGlitchSearch(
    label=label, outdir=outdir, sftfilepath=sftfilepath,
    theta_prior=theta_prior, nglitch=1, tref=tref, nsteps=[500, 500],
    ntemps=3, log10temperature_min=-0.5, minStartTime=tstart,
    maxStartTime=tstart+Tspan)
search.run()
search.plot_corner()
search.print_summary()

In this example, we of course know the properties of the simulated signal. In a real search, the prior will be based on the candidate one is following up. Here for example, we are performing a directed follow up, one could search over Alpha and Delta by specifying a prior distribution (rather than a fixed value).

The search.run() command runs the sampler (which may take a few minutes) and will save an image data/glitch_robust_search_walkers.png:

then we plot a corner plot

and print a summary to the terminal:

10:59 INFO    : Summary:
10:59 INFO    : theta0 index: 0
10:59 INFO    : Max twoF: 5155.82543945 with parameters:
  F0         = 3.000000001e+01
  F1         = -9.999994683e-09
  delta_F0   = 4.993886517e-06
  delta_F1   = 9.876309376e-13
  tglitch    = 1.002318834e+09
10:59 INFO    : Median +/- std for production values
10:59 INFO    :   F0         = 3.000000001e+01 +/- 1.698754802e-08
10:59 INFO    :   F1         = -9.999993595e-09 +/- 1.144523267e-14
10:59 INFO    :   delta_F0   = 4.996919647e-06 +/- 2.021047780e-08
10:59 INFO    :   delta_F1   = 9.869144054e-13 +/- 1.365274015e-14
10:59 INFO    :   tglitch    = 1.002319104e+09 +/- 9.289483961e+03