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.pywill 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