Adds examples for fully-coherent searches on glitching data

README.md

@@ -10,7 +10,8 @@ All examples can be run from their source scripts in [examples](examples), or
 for each example there is descriptive documentation:
 
 * [Making fake data with and without glitches](docs/make_fake_data.md)
-* [Fully coherent searches MCMC](docs/fully_coherent_search.md)
+* [Fully coherent MCMC search](docs/fully_coherent_search.md)
+* [Fully coherent MCMC search on data containing glitching signals](docs/fully_coherent_search_on_glitching_data.md)
 
 ## Installation
 

docs/fully_coherent_search_on_glitching_data.md

+# Fully coherent search on glitching data using MCMC
+
+This example applies the basic [fully coherent
+search](fully_coherent_search.md), to the glitching signal data set created in
+[make fake data](make_fake_data.md]). The aim here is to illustrate the effect
+such a signal can have on a fully-coherent search. The complete script for this
+example canbe found
+[here](../example/fully_cohrent_search_on_glitching_data.py).
+
+
+We use the same prior as in the basic fully-coherent search, except that we
+will modify the prior on `F0` to a flat uniform prior. The reason for this is
+to highlight the multimodal nature of the posterior which results from the
+glitch (a normal prior centered on one of the modes will bias one mode over
+the other). So our initial set up is
+
+```
+from pyfstat import MCMCSearch
+
+F0 = 30.0
+F1 = -1e-10
+F2 = 0
+Alpha = 5e-3
+Delta = 6e-2
+tref = 362750407.0
+
+tstart = 1000000000
+duration = 100*86400
+tend = tstart = duration
+
+theta_prior = {'F0': {'type': 'unif', 'lower': F0-5e-5,
+                      'upper': F0+5e-5},
+               'F1': {'type': 'norm', 'loc': F1, 'scale': abs(1e-6*F1)},
+               'F2': F2,
+               'Alpha': Alpha,
+               'Delta': Delta
+               }
+```
+
+Next, we will use 10 temperatures and a larger number of walkers - these have
+been tuned to illustrate the bimodal nature of the posterior
+
+```
+ntemps = 10
+betas = np.logspace(0, -30, ntemps)
+nwalkers = 500
+nsteps = [100, 100, 100]
+
+mcmc = MCMCSearch('fully_coherent_on_glitching_data', 'data',
+                  sftlabel='glitch', sftdir='data',
+                  theta_prior=theta_prior, tref=tref, tstart=tstart, tend=tend,
+                  nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps, betas=betas,
+                  scatter_val=1e-6)
+
+mcmc.run()
+mcmc.plot_corner(add_prior=True)
+```
+
+Running this takes slightly longer than the basic example (a few minutes) and
+produces a multimodal posterior:
+![](img/fully_coherent_on_glitching_data_corner.png)
+
+Clearly one central peak pertains to the original frequency `F0=30`. At
+`30+0.4e-5`, we find the second largest peak - this is the mode corresponding
+to the second half of the data. In reality, we would expect both peaks to be
+of equal size since the glitch occurs exactly half way through (see [how the
+data was made](make_fake_data.md)). We will confirm this later on by performing
+a grid search.
+
+Finally, the maximum twoF value found is
+
+```
+>>> mcmc.print_summary()
+Max twoF: 411.595
+```
+That is, compared to the basic search (on a smooth signal) which had a twoF of
+`1756.44177246` (in agreement with the predicted twoF), we have lost a large
+fraction of the SNR due to the glitch.
+

docs/make_fake_data.md

@@ -74,6 +74,7 @@ data.run_makefakedata()
 In fact, the previous two commands are wrapped together by a single call to
 `data.make_data()` which we will use from now on.
 
+
 ## Glitching signal
 
 We now want to generate a set of data which contains a *glitching signal*. We
@@ -83,8 +84,8 @@ another `Writer` instance called `glitch_data`, and then run `make_data()`
 
 ```
 dtglitch = duration/2.0
-delta_F0 = 1e-6 * F0
-delta_F1 = 1e-5 * F1
+delta_F0 = 0.4e-5
+delta_F1 = 0
 
 glitch_data = Writer(
     label='glitch', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,
@@ -117,9 +118,9 @@ Delta = 5.999999999999999778e-02
 h0 = 9.999999999999999604e-24
 cosi = 0.000000000000000000e+00
 psi = 0.000000000000000000e+00
-phi0 = -1.222261350197196007e+05
-Freq = 3.000003064156959098e+01
-f1dot = -1.000009999999999993e-10
+phi0 = -1.612440256772935390e+04
+Freq = 3.000000400000000056e+01
+f1dot = -1.000000000000000036e-10
 f2dot = 0.000000000000000000e+00
 refTime = 362750407.000000
 transientWindowType=rect
@@ -130,4 +131,14 @@ transientTauDays=50.000
 The glitch config file uses transient windows to create two non-overlapping,
 but continuous signals.
 
+## Expected twoF
+
+Finally, the `Writer` class also provides a wrapper of `lalapps_PredictFstat`.
+So calling
+
+```
+>>> print data.predict_fstat()
+1721.1
+```
 
+Notice that the predicted value will be the same for both sets of data.

examples/fully_coherent_search_on_glitching_data.py

+from pyfstat import MCMCSearch
+import numpy as np
+
+F0 = 30.0
+F1 = -1e-10
+F2 = 0
+Alpha = 5e-3
+Delta = 6e-2
+tref = 362750407.0
+
+tstart = 1000000000
+duration = 100*86400
+tend = tstart = duration
+
+theta_prior = {'F0': {'type': 'unif', 'lower': F0-5e-5,
+                      'upper': F0+5e-5},
+               'F1': {'type': 'norm', 'loc': F1, 'scale': abs(1e-6*F1)},
+               'F2': F2,
+               'Alpha': Alpha,
+               'Delta': Delta
+               }
+
+ntemps = 10
+betas = np.logspace(0, -30, ntemps)
+nwalkers = 500
+nsteps = [100, 100, 100]
+
+mcmc = MCMCSearch('fully_coherent_on_glitching_data', 'data',
+                  sftlabel='glitch', sftdir='data',
+                  theta_prior=theta_prior, tref=tref, tstart=tstart, tend=tend,
+                  nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps, betas=betas,
+                  scatter_val=1e-6)
+mcmc.run()
+mcmc.plot_corner(add_prior=True)
+mcmc.print_summary()

examples/make_fake_data.py

@@ -27,8 +27,8 @@ data.make_data()
 
 # Next, taking the same signal parameters, we include a glitch half way through
 dtglitch = duration/2.0
-delta_F0 = 1e-6 * F0
-delta_F1 = 1e-5 * F1
+delta_F0 = 0.4e-5
+delta_F1 = 0
 
 glitch_data = Writer(
     label='glitch', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,
@@ -36,3 +36,7 @@ glitch_data = Writer(
     dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1)
 glitch_data.make_data()
 
+
+# The predicted twoF, given by lalapps_predictFstat can be accsessed by
+
+print data.predict_fstat()