|
|
# Glitch robust search
|
|
|
|
|
|
TBD |
|
|
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:
|
|
|
|
|
|
```python
|
|
|
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:
|
|
|
|
|
|
```python
|
|
|
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`:
|
|
|
![](img/glitch_robust_search_walkers.png)
|
|
|
|
|
|
then we plot a corner plot
|
|
|
![](img/glitch_robust_search_corner.png)
|
|
|
|
|
|
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
|
|
|
``` |