diff --git a/Paper/AllSky_run_setup.tex b/Paper/AllSky_run_setup.tex
deleted file mode 100644
index b54e483bd7b1183179784ac25cff08a97a682026..0000000000000000000000000000000000000000
--- a/Paper/AllSky_run_setup.tex
+++ /dev/null
@@ -1,8 +0,0 @@
-\begin{tabular}{c|cccccc}
-Stage & $\Nseg$ & $\Tcoh^{\rm days}$ &$\Nsteps$ & $\V$ & $\Vsky$ & $\Vpe$ \\ \hline
-0 & 27 & 3.7 & 100 & 1.0 & 1.0 & 1.0 \\
-1 & 15 & 6.7 & 100 & 1.0 & 1.0 & 1.0 \\
-2 & 8 & 12.5 & 100 & 1.0 & 1.0 & 1.0 \\
-3 & 4 & 25.0 & 100 & 1.0 & 1.0 & 1.0 \\
-4 & 1 & 100.0 & 50,50 & 1.0 & 1.0 & 1.0 \\
-\end{tabular}
diff --git a/Paper/DirectedMC/generate_table.py b/Paper/DirectedMC/generate_table.py
index c1e8274ffde31c10f65ce7de994ed94819f0af34..3877cf5c14401f297f58f08319d8fe2b5234bf75 100644
--- a/Paper/DirectedMC/generate_table.py
+++ b/Paper/DirectedMC/generate_table.py
@@ -27,7 +27,7 @@ DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)
depth = 100
-nsteps = 50
+nsteps = 20
run_setup = [((nsteps, 0), 20, False),
((nsteps, 0), 7, False),
((nsteps, 0), 2, False),
diff --git a/Paper/Examples/Makefile b/Paper/Examples/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..3557f777a60519dacd872c5ed0c83adf44d031da
--- /dev/null
+++ b/Paper/Examples/Makefile
@@ -0,0 +1,10 @@
+frequency_grid_files = grided_frequency_search_1D.png
+fully_coherent_files = fully_coherent_search_using_MCMC_walkers.png
+follow_up_files = follow_up_run_setup.tex follow_up_walkers.png
+transient_files = transient_search_initial_stage_twoFcumulative.png transient_search_corner.png
+glitch_files = single_glitch_F0F1_grid_2D.png single_glitch_glitchSearch_corner.png
+
+all_files = $(frequency_grid_files) $(fully_coherent_files) $(follow_up_files) $(transient_files) $(glitch_files)
+
+copyfiles:
+ cd data; cp $(all_files) ../../
diff --git a/Paper/Examples/fully_coherent_search_using_MCMC.py b/Paper/Examples/fully_coherent_search_using_MCMC.py
new file mode 100644
index 0000000000000000000000000000000000000000..fb4464ec9de22ab5f73f53bf11d5366597c85e54
--- /dev/null
+++ b/Paper/Examples/fully_coherent_search_using_MCMC.py
@@ -0,0 +1,61 @@
+import pyfstat
+import numpy as np
+
+# Properties of the GW data
+sqrtSX = 1e-23
+tstart = 1000000000
+duration = 100*86400
+tend = tstart + duration
+
+# Properties of the signal
+F0 = 30.0
+F1 = -1e-10
+F2 = 0
+Alpha = 5e-3
+Delta = 6e-2
+tref = .5*(tstart+tend)
+
+depth = 10
+h0 = sqrtSX / depth
+data_label = 'fully_coherent_search_using_MCMC'
+
+data = pyfstat.Writer(
+ label=data_label, outdir='data', tref=tref,
+ tstart=tstart, F0=F0, F1=F1, F2=F2, duration=duration, Alpha=Alpha,
+ Delta=Delta, h0=h0, sqrtSX=sqrtSX)
+data.make_data()
+
+# The predicted twoF, given by lalapps_predictFstat can be accessed by
+twoF = data.predict_fstat()
+print 'Predicted twoF value: {}\n'.format(twoF)
+
+DeltaF0 = 1e-7
+DeltaF1 = 1e-13
+VF0 = (np.pi * duration * DeltaF0)**2 / 3.0
+VF1 = (np.pi * duration**2 * DeltaF1)**2 * 4/45.
+print '\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1)
+
+theta_prior = {'F0': {'type': 'unif',
+ 'lower': F0-DeltaF0/2.,
+ 'upper': F0+DeltaF0/2.},
+ 'F1': {'type': 'unif',
+ 'lower': F1-DeltaF1/2.,
+ 'upper': F1+DeltaF1/2.},
+ 'F2': F2,
+ 'Alpha': Alpha,
+ 'Delta': Delta
+ }
+
+ntemps = 1
+log10temperature_min = -1
+nwalkers = 1000
+nsteps = [50, 50]
+
+mcmc = pyfstat.MCMCSearch(
+ label='fully_coherent_search_using_MCMC', outdir='data',
+ sftfilepath='data/*'+data_label+'*sft', theta_prior=theta_prior, tref=tref,
+ minStartTime=tstart, maxStartTime=tend, nsteps=nsteps, nwalkers=nwalkers,
+ ntemps=ntemps, log10temperature_min=log10temperature_min)
+mcmc.run(context='paper', subtractions=[30, -1e-10])
+mcmc.plot_corner(add_prior=True)
+mcmc.print_summary()
diff --git a/Paper/Examples/grided_frequency_search.py b/Paper/Examples/grided_frequency_search.py
new file mode 100644
index 0000000000000000000000000000000000000000..5e4a423ce8a080081de364bf058287d184df9b01
--- /dev/null
+++ b/Paper/Examples/grided_frequency_search.py
@@ -0,0 +1,64 @@
+import pyfstat
+import numpy as np
+import matplotlib.pyplot as plt
+
+plt.style.use('paper')
+
+F0 = 30.0
+F1 = 0
+F2 = 0
+Alpha = 1.0
+Delta = 1.5
+
+# Properties of the GW data
+sqrtSX = 1e-23
+tstart = 1000000000
+duration = 100*86400
+tend = tstart+duration
+tref = .5*(tstart+tend)
+
+depth = 70
+data_label = 'grided_frequency_depth_{:1.0f}'.format(depth)
+
+h0 = sqrtSX / depth
+
+data = pyfstat.Writer(
+ label=data_label, outdir='data', tref=tref,
+ tstart=tstart, F0=F0, F1=F1, F2=F2, duration=duration, Alpha=Alpha,
+ Delta=Delta, h0=h0, sqrtSX=sqrtSX)
+data.make_data()
+
+m = 0.001
+dF0 = np.sqrt(12*m)/(np.pi*duration)
+DeltaF0 = 800*dF0
+F0s = [F0-DeltaF0/2., F0+DeltaF0/2., dF0]
+F1s = [F1]
+F2s = [F2]
+Alphas = [Alpha]
+Deltas = [Delta]
+search = pyfstat.GridSearch(
+ 'grided_frequency_search', 'data', 'data/*'+data_label+'*sft', F0s, F1s,
+ F2s, Alphas, Deltas, tref, tstart, tend)
+search.run()
+
+fig, ax = plt.subplots()
+xidx = search.keys.index('F0')
+frequencies = np.unique(search.data[:, xidx])
+twoF = search.data[:, -1]
+
+#mismatch = np.sign(x-F0)*(duration * np.pi * (x - F0))**2 / 12.0
+ax.plot(frequencies, twoF, 'k', lw=1)
+DeltaF = frequencies - F0
+sinc = np.sin(np.pi*DeltaF*duration)/(np.pi*DeltaF*duration)
+A = np.abs((np.max(twoF)-4)*sinc**2 + 4)
+ax.plot(frequencies, A, '-r', lw=1)
+ax.set_ylabel('$\widetilde{2\mathcal{F}}$')
+ax.set_xlabel('Frequency')
+ax.set_xlim(F0s[0], F0s[1])
+dF0 = np.sqrt(12*1)/(np.pi*duration)
+xticks = [F0-10*dF0, F0, F0+10*dF0]
+ax.set_xticks(xticks)
+xticklabels = ['$f_0 {-} 10\Delta f$', '$f_0$', '$f_0 {+} 10\Delta f$']
+ax.set_xticklabels(xticklabels)
+plt.tight_layout()
+fig.savefig('{}/{}_1D.png'.format(search.outdir, search.label), dpi=300)
diff --git a/Paper/Examples/single_glitch.py b/Paper/Examples/single_glitch.py
new file mode 100644
index 0000000000000000000000000000000000000000..61b1212a8bd9c49d3ed33e4bb77a91b8978ece62
--- /dev/null
+++ b/Paper/Examples/single_glitch.py
@@ -0,0 +1,99 @@
+import pyfstat
+import numpy as np
+import matplotlib.pyplot as plt
+
+sqrtSX = 1e-22
+tstart = 1000000000
+duration = 100*86400
+tend = tstart+duration
+
+# Define parameters of the Crab pulsar as an example
+tref = .5*(tstart + tend)
+F0 = 30.0
+F1 = -1e-10
+F2 = 0
+Alpha = 5e-3
+Delta = 6e-2
+
+# Signal strength
+depth = 10
+h0 = sqrtSX / depth
+
+VF0 = VF1 = 200
+dF0 = np.sqrt(3)/(np.pi*duration)
+dF1 = np.sqrt(45/4.)/(np.pi*duration**2)
+DeltaF0 = VF0 * dF0
+DeltaF1 = VF1 * dF1
+
+# Next, taking the same signal parameters, we include a glitch half way through
+dtglitch = duration/2.0
+delta_F0 = 0.25*DeltaF0
+delta_F1 = -0.1*DeltaF1
+
+glitch_data = pyfstat.Writer(
+ label='single_glitch', outdir='data', tref=tref, tstart=tstart, F0=F0,
+ F1=F1, F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0,
+ sqrtSX=sqrtSX, dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1)
+glitch_data.make_data()
+
+F0s = [F0-DeltaF0/2., F0+DeltaF0/2., 1*dF0]
+F1s = [F1-DeltaF1/2., F1+DeltaF1/2., 1*dF1]
+F2s = [F2]
+Alphas = [Alpha]
+Deltas = [Delta]
+search = pyfstat.GridSearch(
+ 'single_glitch_F0F1_grid', 'data', 'data/*single_glitch*sft', F0s, F1s,
+ F2s, Alphas, Deltas, tref, tstart, tend)
+search.run()
+search.plot_2D('F0', 'F1')
+
+theta_prior = {'F0': {'type': 'unif', 'lower': F0-DeltaF0/2.,
+ 'upper': F0+DeltaF0/2},
+ 'F1': {'type': 'unif', 'lower': F1-DeltaF1/2.,
+ 'upper': F1+DeltaF1/2},
+ 'F2': F2,
+ 'Alpha': Alpha,
+ 'Delta': Delta
+ }
+ntemps = 3
+log10temperature_min = -0.05
+nwalkers = 100
+nsteps = [500, 500]
+
+mcmc = pyfstat.MCMCSearch(
+ 'single_glitch', 'data', sftfilepath='data/*_single_glitch*.sft',
+ theta_prior=theta_prior, tref=tref, minStartTime=tstart, maxStartTime=tend,
+ nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
+ log10temperature_min=log10temperature_min)
+
+mcmc.run()
+mcmc.plot_corner(figsize=(3.2, 3.2))
+mcmc.print_summary()
+
+
+theta_prior = {'F0': {'type': 'unif', 'lower': F0-DeltaF0/2.,
+ 'upper': F0+DeltaF0/2},
+ 'F1': {'type': 'unif', 'lower': F1-DeltaF1/2.,
+ 'upper': F1+DeltaF1/2},
+ 'F2': F2,
+ 'Alpha': Alpha,
+ 'Delta': Delta,
+ 'tglitch': {'type': 'unif', 'lower': tstart+0.1*duration,
+ 'upper': tend-0.1*duration},
+ 'delta_F0': {'type': 'halfnorm', 'loc': 0, 'scale': DeltaF0},
+ 'delta_F1': {'type': 'norm', 'loc': 0, 'scale': DeltaF1},
+ }
+ntemps = 3
+log10temperature_min = -0.1
+nwalkers = 100
+nsteps = [1000, 1000]
+glitch_mcmc = pyfstat.MCMCGlitchSearch(
+ 'single_glitch_glitchSearch', 'data',
+ sftfilepath='data/*_single_glitch*.sft', theta_prior=theta_prior,
+ tref=tref, minStartTime=tstart, maxStartTime=tend, nsteps=nsteps,
+ nwalkers=nwalkers, ntemps=ntemps,
+ log10temperature_min=log10temperature_min)
+glitch_mcmc.run()
+glitch_mcmc.plot_corner(figsize=(6, 6))
+glitch_mcmc.print_summary()
+
diff --git a/Paper/Examples/transient_search_using_MCMC.py b/Paper/Examples/transient_search_using_MCMC.py
new file mode 100644
index 0000000000000000000000000000000000000000..8df9b12b4627b8b7a1367443ae3d3b5c7af37ca2
--- /dev/null
+++ b/Paper/Examples/transient_search_using_MCMC.py
@@ -0,0 +1,91 @@
+import pyfstat
+import numpy as np
+import matplotlib.pyplot as plt
+
+plt.style.use('thesis')
+
+F0 = 30.0
+F1 = -1e-10
+F2 = 0
+Alpha = 5e-3
+Delta = 6e-2
+
+tstart = 1000000000
+duration = 100*86400
+data_tstart = tstart - duration
+data_tend = data_tstart + 3*duration
+tref = .5*(data_tstart+data_tend)
+
+h0 = 1e-23
+sqrtSX = 1e-22
+
+transient = pyfstat.Writer(
+ label='transient', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,
+ F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX,
+ minStartTime=data_tstart, maxStartTime=data_tend)
+transient.make_data()
+print transient.predict_fstat()
+
+
+
+DeltaF0 = 6e-7
+DeltaF1 = 1e-13
+VF0 = (np.pi * duration * DeltaF0)**2 / 3.0
+VF1 = (np.pi * duration**2 * DeltaF1)**2 * 4/45.
+print '\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1)
+
+theta_prior = {'F0': {'type': 'unif',
+ 'lower': F0-DeltaF0/2.,
+ 'upper': F0+DeltaF0/2.},
+ 'F1': {'type': 'unif',
+ 'lower': F1-DeltaF1/2.,
+ 'upper': F1+DeltaF1/2.},
+ 'F2': F2,
+ 'Alpha': Alpha,
+ 'Delta': Delta
+ }
+
+ntemps = 3
+log10temperature_min = -1
+nwalkers = 100
+nsteps = [750, 250]
+
+mcmc = pyfstat.MCMCSearch(
+ label='transient_search_initial_stage', outdir='data',
+ sftfilepath='data/*transient*sft', theta_prior=theta_prior, tref=tref,
+ minStartTime=data_tstart, maxStartTime=data_tend, nsteps=nsteps,
+ nwalkers=nwalkers, ntemps=ntemps,
+ log10temperature_min=log10temperature_min)
+mcmc.run()
+mcmc.plot_cumulative_max()
+mcmc.print_summary()
+
+theta_prior = {'F0': {'type': 'unif',
+ 'lower': F0-DeltaF0/2.,
+ 'upper': F0+DeltaF0/2.},
+ 'F1': {'type': 'unif',
+ 'lower': F1-DeltaF1/2.,
+ 'upper': F1+DeltaF1/2.},
+ 'F2': F2,
+ 'Alpha': Alpha,
+ 'Delta': Delta,
+ 'transient_tstart': {'type': 'unif',
+ 'lower': data_tstart,
+ 'upper': data_tend},
+ 'transient_duration': {'type': 'halfnorm',
+ 'loc': 0,
+ 'scale': 0.5*duration}
+ }
+
+nwalkers = 500
+nsteps = [200, 200]
+
+mcmc = pyfstat.MCMCTransientSearch(
+ label='transient_search', outdir='data',
+ sftfilepath='data/*transient*sft', theta_prior=theta_prior, tref=tref,
+ minStartTime=data_tstart, maxStartTime=data_tend, nsteps=nsteps,
+ nwalkers=nwalkers, ntemps=ntemps,
+ log10temperature_min=log10temperature_min)
+mcmc.run()
+mcmc.plot_corner()
+mcmc.print_summary()
diff --git a/Paper/allsky_setup_run_setup.tex b/Paper/allsky_setup_run_setup.tex
index 16033d738786dd62c55554defc247c9ed9c38794..1d765771ee4de15b1a45e08e3ca6ab86da3dc83b 100644
--- a/Paper/allsky_setup_run_setup.tex
+++ b/Paper/allsky_setup_run_setup.tex
@@ -1,8 +1,8 @@
\begin{tabular}{c|cccccc}
Stage & $\Nseg$ & $\Tcoh^{\rm days}$ &$\Nsteps$ & $\V$ & $\Vsky$ & $\Vpe$ \\ \hline
-0 & 20 & 5.0 & 100 & $2{\times}10^{2}$ & 10.0 & 10.0 \\
-1 & 11 & 9.1 & 100 & $2{\times}10^{3}$ & 40.0 & 50.0 \\
-2 & 6 & 16.7 & 100 & $2{\times}10^{4}$ & $1{\times}10^{2}$ & $2{\times}10^{2}$ \\
-3 & 3 & 33.3 & 100 & $1{\times}10^{5}$ & $2{\times}10^{2}$ & $6{\times}10^{2}$ \\
-4 & 1 & 100.0 & 100,100 & $8{\times}10^{5}$ & $3{\times}10^{2}$ & $2{\times}10^{3}$ \\
+0 & 20 & 5.0 & 50 & $2{\times}10^{2}$ & 10.0 & 10.0 \\
+1 & 11 & 9.1 & 50 & $2{\times}10^{3}$ & 40.0 & 50.0 \\
+2 & 6 & 16.7 & 50 & $2{\times}10^{4}$ & $1{\times}10^{2}$ & $2{\times}10^{2}$ \\
+3 & 3 & 33.3 & 50 & $1{\times}10^{5}$ & $2{\times}10^{2}$ & $6{\times}10^{2}$ \\
+4 & 1 & 100.0 & 50,50 & $8{\times}10^{5}$ & $3{\times}10^{2}$ & $2{\times}10^{3}$ \\
\end{tabular}
diff --git a/Paper/bibliography.bib b/Paper/bibliography.bib
index 770b86af8acaba6c4a87971141cc69193a138755..8996364e5942dcdda3b9e46be021a9ddae43e934 100644
--- a/Paper/bibliography.bib
+++ b/Paper/bibliography.bib
@@ -413,3 +413,54 @@ volume = {91},
year = {2015}
}
+@article{wette2012,
+ arxivId = {1111.5650},
+ author = {Wette, Karl},
+ doi = {10.1103/PhysRevD.85.042003},
+ eprint = {1111.5650},
+ file = {:home/greg/Dropbox/Papers/Wette{\_}2012.pdf:pdf},
+ isbn = {0556-2821},
+ issn = {15507998},
+ journal = {Physical Review D},
+ number = {4},
+ title = {{Estimating the sensitivity of wide-parameter-space searches for gravitational-wave pulsars}},
+ volume = {85},
+ year = {2012}
+}
+
+@article{melatos2008,
+ archivePrefix = {arXiv},
+ arxivId = {0710.1021},
+ author = {Melatos, A. and Peralta, C. and Wyithe, J. S. B.},
+ doi = {10.1086/523349},
+ eprint = {0710.1021},
+ file = {:home/greg/Dropbox/Papers/Mealtos{\_}2007.pdf:pdf},
+ issn = {0004-637X},
+ journal = {The Astrophysical Journal},
+ keywords = {dense matter — pulsars: general — stars: interiors},
+ number = {Jensen 1998},
+ pages = {1103},
+ title = {{Avalanche dynamics of radio pulsar glitches}},
+ url = {http://arxiv.org/abs/0710.1021},
+ volume = {672},
+ year = {2008}
+}
+
+@article{espinoza2011,
+ archivePrefix = {arXiv},
+ arxivId = {1102.1743},
+ author = {Espinoza, Crist{\'{o}}bal and Lyne, Andrew and Stappers, Ben and Kramer, Michael},
+ doi = {10.1063/1.3615093},
+ eprint = {1102.1743},
+ file = {:home/greg/Dropbox/Papers/Espinoza{\_}2011.pdf:pdf},
+ isbn = {9780735409156},
+ issn = {0094243X},
+ journal = {AIP Conference Proceedings},
+ keywords = {general,glitches,pulsars,stars: neutron},
+ number = {March},
+ pages = {117--120},
+ title = {{Glitches in the rotation of pulsars}},
+ volume = {1357},
+ year = {2011}
+}
+
diff --git a/Paper/definitions.tex b/Paper/definitions.tex
index 4eb752ef1eb7d4ed0e5e82b6f14c6e8cc61f8e6a..233a66e55fe61797cf43f0482788a9ac978ce1b1 100644
--- a/Paper/definitions.tex
+++ b/Paper/definitions.tex
@@ -3,6 +3,10 @@
\newcommand{\A}{\boldsymbol{\mathcal{A}}}
\newcommand{\blambda}{\boldsymbol{\mathbf{\lambda}}}
\newcommand{\blambdaSignal}{\boldsymbol{\mathbf{\lambda}}^{\rm s}}
+\newcommand{\tglitch}{t_{\rm glitch}}
+\newcommand{\tstart}{t_{\rm start}}
+\newcommand{\tend}{t_{\rm end}}
+\newcommand{\Nglitches}{N_{\rm glitches}}
\newcommand{\Tspan}{T_{\rm span}}
\newcommand{\Tcoh}{T_{\rm coh}}
\newcommand{\tref}{t_{\rm ref}}
diff --git a/Paper/directed_setup_run_setup.tex b/Paper/directed_setup_run_setup.tex
new file mode 100644
index 0000000000000000000000000000000000000000..bf81a4a6451a8a16f1de11fb11edb47edf81b0f2
--- /dev/null
+++ b/Paper/directed_setup_run_setup.tex
@@ -0,0 +1,7 @@
+\begin{tabular}{c|cccc}
+Stage & $\Nseg$ & $\Tcoh^{\rm days}$ &$\Nsteps$ & $\Vpe$ \\ \hline
+0 & 20 & 5.0 & 20 & 10.0 \\
+1 & 7 & 14.3 & 20 & $1{\times}10^{2}$ \\
+2 & 2 & 50.0 & 20 & $1{\times}10^{3}$ \\
+3 & 1 & 100.0 & 20,20 & $2{\times}10^{3}$ \\
+\end{tabular}
diff --git a/Paper/fully_coherent_search_using_MCMC_walkers.png b/Paper/fully_coherent_search_using_MCMC_walkers.png
index d338b9fcbc9d77a39c38c7ccab75710570e59759..d7f82cc39211b04ecaf3a421e705103e1ef1b97d 100644
Binary files a/Paper/fully_coherent_search_using_MCMC_walkers.png and b/Paper/fully_coherent_search_using_MCMC_walkers.png differ
diff --git a/Paper/grided_frequency_search_1D.png b/Paper/grided_frequency_search_1D.png
index 7df348058ea5b974603b629a9680e50ebc8893fe..d63b0c346443932d7c848c5edea6cdd1b201e524 100644
Binary files a/Paper/grided_frequency_search_1D.png and b/Paper/grided_frequency_search_1D.png differ
diff --git a/Paper/paper_cw_mcmc.tex b/Paper/paper_cw_mcmc.tex
index b336515a9d04338fbad17001ccb2524310e176e5..e1648c22abcd80293274a59a9f9d82afeb2c5684 100644
--- a/Paper/paper_cw_mcmc.tex
+++ b/Paper/paper_cw_mcmc.tex
@@ -296,6 +296,7 @@ computed analytically as
\int
\mathcal{L}(x ;\A, \blambda)
P(\A| \Hs, I) d\A
+\label{eqn_B_S_N}
\\
& = \frac{C (2\pi)^{2} e^{\F(x| \blambda)}}
{\sqrt{\textrm{det} \mathcal{M}}},
@@ -790,6 +791,13 @@ fraction of the volume at adjacent stages. That is we define
\end{equation}
where $\mathcal{R} \ge 1$ as $\Tcoh^{i+1} > \Tcoh^{i}$.
+For the MCMC simulations to be succesful, this initial bounding
+box, given the segment setup which produced the candidate, must be small
+compared to the width of the signal (at that segment setup). If we start out
+follow-up with the same search setup used in the source search, i.e. we set
+$\Nseg^0$ equal to the number of segments used in the input search, then for
+the MCMC simulation to work, we require that $\V(\Nseg^{0}) \sim \mathcal{O}(100)$.
+
Given a fixed prior on the Doppler parameters and a fixed span of data, the
metric volume $\V_{\Nstages}$ for $\Tcoh^{\Nstages} = \Tspan$ is fixed, or in
terms of the number of segments, $\Nseg^{\Nstages} = 1$.
@@ -824,12 +832,12 @@ in the noise.
First, we must define the setup for the run. Using $\mathcal{R}=10$ and
$\V^{\rm min}=100$ our optimisation procedure is run and proposes the setup
-layed out in Table~\ref{tab_signal_follow_up}. In addition, we show the
+layed out in Table~\ref{tab_follow_up_run_setup}. In addition, we show the
number of steps taken at each stage.
\begin{table}[htb]
\caption{The search setup used in Figure~\ref{fig_follow_up}, generated with
$\mathcal{R}=10$ and $\V^{\rm min}=100$.}
-\label{tab_follow_up}
+\label{tab_follow_up_run_setup}
\input{follow_up_run_setup}
\end{table}
@@ -842,7 +850,7 @@ found to result in the MCMC simulations `loosing' the peaks between stages, we
conservatively opt for 10 here, but values as large as 100 where also successful.
In Figure~\ref{fig_follow_up} we show the progress of the MCMC sampler during
-the follow-up. As expected from Table~\ref{tab_follow_up}, during
+the follow-up. As expected from Table~\ref{tab_follow_up_run_setup}, during
the initial stage the signal peak is broad with respect to the size of the
prior volume, therefore the MCMC simulation quickly converges to it. Subsequently,
each time the number of segments is reduced, the peak narrows and the samplers
@@ -855,9 +863,8 @@ Eventually, the walkers all converge to the true signal.
\caption{In the top three panels we show the progress of the 500 parallel
walkers (see Figure~\ref{fig_MCMC_simple_example} for a description) during the
MCMC simulation for each of the search parameters, frequency $f$,
-right-ascension $\alpha$ and declination $\delta$. Each vertical dashed line
-indicates the start of a new stage of the search, the parameters for all stages
-are listed in Table~\ref{tab_follow_up}.}
+right-ascension $\alpha$ and declination $\delta$. Each vertical dashed line indicates the start of a new stage of the search, the parameters for all stages
+are listed in Table~\ref{tab_follow_up_run_setup}.}
\label{fig_follow_up}
\end{figure}
@@ -869,52 +876,89 @@ chains explore the other `noise peaks' in the data.
\section{Monte Carlo studies}
-In order to understand how well the MCMC follow-up method works, we will now
-study the recovery fraction as a function of signal depth. This will be done in
-a Monte Carlo study, with independent random realisations of the Guassian
-noise, amplitude, and Doppler parameters in suitable ranges. Such a method is
-analagous to the studies performed in \citet{shaltev2013}, except that we
-present results as a function of the fixed signal depth, rather than the
-squared SNR.
-
-In particular we will generate \comment{N} realisations of Gaussian noise data
-lasting for 100 days, each with a simulated CW signal. We choose the parameters
-of the signal in such a way to model the candidates generated from directed and
-all-sky searches by drawing the signal parameters from appropriate
-distributions. However, we do not draw $h_0$ randomly, but instead run the MC
-study at a number of selected values chosen such that given the fixed
-$\sqrt{S_n}=2\times10^{3}$, the signals are injected with a depth $\mathcal{D}
-\in [100, 400]$.
-
-To simulate an isotropic distribution of sources, we draw the remaining
-amplitude parameters for each signal uniformly from $\phi \in [0, 2\pi]$, $\psi
-\in [-\pi/4, \pi/4]$, and $\cos\iota \in [-1, 1]$.
-
+In order to understand how well the MCMC follow-up method works, we will test
+its ability to succesfully identify simulated signals in Gaussian. This will be
+done in a Monte Carlo study, with independent random realisations of the
+Guassian noise, amplitude, and Doppler parameters in suitable ranges. Such a
+method is analagous to the studies performed in \citet{shaltev2013}, except
+that we present results as a function of the fixed injected signal depth,
+rather than the squared SNR.
+
+In particular we will generate a number of 100-day data sets containing
+independent realisations of Gaussian noise and a simulated CW signal. We choose
+the parameters of the signal in such a way to model the candidates generated
+from directed and all-sky searches by drawing the signal parameters from
+appropriate distributions. However, we do not draw $h_0$ randomly, but instead
+run the MC study at a number of selected values chosen such that given the
+fixed $\sqrt{S_n}=2\times10^{3}$, the signals are injected with a depth
+$\mathcal{D} \in [100, 400]$. To simulate an isotropic distribution of
+sources, we draw the remaining amplitude parameters for each signal uniformly
+from $\phi \in [0, 2\pi]$, $\psi \in [-\pi/4, \pi/4]$, and $\cos\iota \in [-1,
+1]$.
+
+To provide a reference, we will compare the MC follow-up study against the
+expected maximum theoretical detection probability for an infinitly dense
+fully-coherent search of data containing isotropically-distributed signals as
+calculated by \citet{wette2012}. Note however that we will parameterise with
+respect to the signal depth (i.e. using Equation~(3.8) of \citet{wette2012} to
+relate the averaged-SNR to the depth). The probability is maximal in the
+sense that signals are `lost' during the follow-up due simply to the fact that
+they are not sufficiently strong to rise above the noise.
\subsection{Follow-up candidates from a directed search}
+\label{sec_directed_follow_up}
In a directed search, the sky location parameters $\alpha$ and $\delta$ are
fixed - in our study we fix them on the location of the Crab pulsar, but this
choice is arbitrary and holds no particular significance. The ouput of a
-gridded directed search would provide the grid-point with the highest detection
-statistic, and some estimate of the iso-mismatch contours in which the
-candidate is expected to exist. To simulate this, we will define a frequency
-and spindown of $f_0=30$~Hz and $\dot{f}_0=10^{-10}$Hz/s and a surrounding box
-$\Delta f$ and $\Delta \dot{f}$ which correspoinds to a fully-coherent
-$\V=10^{4}$ with $\V_{f}=\V_{\dot{f}}$. Then, we pick a candidate we first pick
-a point randomly in the unit circle, then using the PE-phase-metric we convert
-this into a random point in an isomismatch contour. In addition, we also select
-the amplitude parameters We then select a set of particular values of
-
+semi-coherent gridded directed search would provide the grid-point with the
+highest detection statistic, and some box bounding the candidate given some
+uncertainty. We will assume in this section that given the search setup
+$\Nseg^{0}$ of the input search, the bounding box is sufficiently small to
+ensure that $\V(\Nseg^0)\sim \mathcal{O}(100)$. This is not a limiting
+assumption as any search can (quite cheaply) increase the density of grid
+points around any interesting candidates in order to better constrain their
+uncertainty.
+
+The behaviour of the follow-up is independent of the exact frequency and
+spin-down values used (this is not true for an all-sky follow-up as discussed
+in Section~\ref{sec_all_sky_follow_up}). As such, we can, without loss of
+generality define our Monte-Carlo follow-up in the following way. First, we
+select an arbitrary frequency and spindown of $f_0=30$~Hz and
+$\dot{f}_0=10^{-10}$Hz/s and a surrounding uncertainty box $\Delta f$ and
+$\Delta \dot{f}$ chosen such that the uncertainty in frequency and spindown are
+roughly equivalent in terms of mismatch. Then, we pick a candidate uniformly
+from within this uncertainty box; this choice reflects the fact that the grid
+is chosen such that the probability distribution of candidate signals is
+uniform.
+
+
+
+Having generated the data given the prescription above, we proceed to perform a
+hierarchical MCMC follow-up. Given the data span and initial bounding box, we
+compute the optimal heirarchical setup, the details of which are given in
+Table~\ref{tab_directed_MC_follow_up}, this setup was used for all MC
+simulations. This table also lists the number of steps, which in this case we
+chose to be 10 - quite a small number of steps for a typical MCMC simulation.
\begin{table}[htb]
\caption{Run-setup for the directed follow-up Monte-Carlo study, generated with
-$\mathcal{R}=10$ and $\V^{\rm min}=100$}
+$\mathcal{R}=10$ and $\Nseg^0=20$.}
\label{tab_directed_MC_follow_up}
\input{directed_setup_run_setup}
\end{table}
+This process yeilds a maximum detection statistic $\widetilde{2\F}^{\rm max}$.
+The signal is considered `detected' if $\widetilde{2\F}^{\rm max} >
+\widetilde{2\F}^{\rm th}$, where we set the threshold at $2\F^{\rm th}=60$,
+corresponding to a p-value of \comment{Finish section}. In
+Figure~\ref{fig_directed_MC_follow_up} we plot the the fraction of the MC
+simulations which where recovered and compare this against the theoretical
+maximum, given the threshold. The figure demonstrates that the recovery power
+of the MCMC follow-up shows only a small margin of loss compared to the
+theoretical maximum. \comment{Probably best to run with more steps and see if
+this doesn't remove the loss}
\begin{figure}[htb]
\centering
\includegraphics[width=0.45\textwidth]{directed_recovery}
@@ -924,13 +968,39 @@ Table~\ref{tab_directed_MC_follow_up}.}
\label{fig_directed_MC_follow_up}
\end{figure}
+
\subsection{Follow-up candidates from an all-sky search}
+\label{sec_all_sky_follow_up}
+
+We now test the follow-up method when applied to candidates from an all-sky
+search which, by definition, have uncertain sky-position parameters $\alpha$
+and $\delta$. Searching over these parameters, in addition to the frequency
+and spin-down not only increases the parameter space volume that needs to be
+searched, but also adds difficulty due to correlations between the sky-position
+and spin-down \comment{Does some of Karl's work discuss this?}.
+
+To replicate the condition of candidates from an all-sky search, we draw the
+candidate positions isotropically from the unit sphere (i.e. $\alpha$ uniform
+on $[0, 2\pi]$ and $\delta = \cos^{-1}(2u{-}1){- }\pi/2$ where $u$ is uniform
+on $[0, 1]$). We then place an uncertainty box containing the candidates with a
+width $\Delta\alpha=\Delta\delta=0.05$; this box is not centered on the
+candidate, but chosen in such a way that the location of the candidate has a
+uniform probability distrubution within the box. This neglects the non-uniform
+variation in $\delta$ over the sky-pacth, but given the small area should not
+cause any significant bias. The frequency, spin-down, and amplitude parameters
+are chosen in the same way as for the directed search
+(Section~\ref{sec_directed_follow_up}).
+
+Producing \CHECK{1000} indepedant MC simulations we the perform a follow-up on
+each using the setup given in Table~\ref{tab_allsky_MC_follow_up}. The
+resulting recovery fraction as a function of the injected signal depth is given
+in Figure~\ref{fig_allsky_MC_follow_up}.
\begin{table}[htb]
\caption{Run-setup for the all-sky follow-up Monte-Carlo study, generated with
-$\mathcal{R}=10$ and $\V^{\rm min}=100$}
+$\mathcal{R}=10$ and $\Nseg^0=20$.}
\label{tab_allsky_MC_follow_up}
-\input{AllSky_run_setup}
+\input{allsky_setup_run_setup}
\end{table}
\begin{figure}[htb]
@@ -1043,6 +1113,82 @@ a simulated transient signal and Gaussian noise.}
\subsection{Glitches}
\label{sec_glitches}
+Observations of radio pulsars show that occasionally neutron stars undergo
+sudden glitch events during which the pulsation frequency suddenly
+\emph{increases} quite distinct from the usual spin-down due to EM or GW
+torques (see \citet{espinoza2011} for a review). These events are typically
+Poisson distributed \citet{melatos2008}, but the exact mechanism is not yet
+fully understood. However, it seems plausible that the glitch mechanism will
+similarly effect any CW emission, i.e. CW signals will also undergo glitches.
+In \citet{ashton2016} we demonstrated that if this is the case, CW sources in
+directed and all-sky searches have a reasonable probability of undergoing one
+or more glitches during typical observations spans and that these glitches may
+be sufficiently large to disrupt the detection; this is particuarly true for
+areas of parameter space with large spindowns. However, the effect is mitigated
+by the use of incoherent searches. As such, we may expect all-sky and directed
+searches to correctly identify glitching signals as candidates, but
+subseuqnetly `loose' them during the follow-up process. As such, we will now
+discuss how the $\F$-statistic optimisation routine discussed in this work can
+be used to search for glitching signals.
+
+We will model a glitch as a sudden discontinuous change in the Doppler
+parameters $\delta \blambda$ at some time $\tglitch$: that the change is only
+in the frequency and spin-down Doppler parameters, i.e. $\delta\blambda = [0,
+0, \delta f, \delta \fdot, \ldots]$. This is to be interpretted, as is done for
+radio-pulsar glitches, as the instantaneous change in frequency and spindown at
+the time of the glitch. Moreover, multiple glitches can be included (we will
+define $\Nglitches$ as the number of glitches) by adding
+an index to the glitch magnitude and glitch time.
+
+An ideal glitch search would include these new parameters into the phase model
+itself and compute the $\F$-statistic fully-coherently i.e. rewriting
+Equation~\eqref{eqn_B_S_N} as
+\begin{multline}
+\Bsn(x| \tglitch, \delta\blambda, \Pic, \blambda, I) \equiv \\
+\int
+\mathcal{L}(x ;\A, \blambda, \tglitch, \delta\blambda)
+P(\A| \Hs, I) d\A
+\end{multline}
+where $\tglitch, \delta\blambda$ modify the phase model. However, we propose the
+following pragmatic approach: let us instead use an semi-coherent detection
+statistic with the epoch between \emph{segments} as $\tglitch$, i.e.
+\begin{multline}
+\Bsn(x| \tglitch, \delta\blambda, \Pic, \blambda, I) \equiv \\
+\int_{\tstart}^{\tglitch}
+\mathcal{L}(x ;\A, \blambda)
+P(\A| \Hs, I) d\A \\ +
+\int_{\tglitch}^{\tend}
+\mathcal{L}(x ;\A, \blambda{+}\delta\blambda)
+P(\A| \Hs, I) d\A
+\end{multline}
+
+This simplistic method leverages readily available and tested code with a loss
+of sensitivity over the fully-coherent search by a factor
+$\sqrt{\Nglitches+1}$. Such a loss is acceptable, given that the signals must
+be sufficiently strong to have been identified by the initial all-sky or directed
+search.
+
+As an illustration of the use of this method, we simulate a signal in Gaussian
+noise which undergoes a glitch. Firstly in Figure~\ref{fig_glitch_grid} we show
+the fully-coherent $2\F$ value computed for a grid of points: this distinctly
+shows two peaks, a feature indicative of a glitching signal, with a maximum
+value of $\sim\CHECK{400}$. Using our incoherent glitch statistic, we plot
+the resulting posterior distributions in Figure~\ref{fig_glitch_posteriors}.
+
+\begin{figure}[htb]
+\centering
+\includegraphics[width=0.5\textwidth]{single_glitch_F0F1_grid_2D}
+\caption{}
+\label{fig_glitch_grid}
+\end{figure}
+
+\begin{figure}[htb]
+\centering
+\includegraphics[width=0.5\textwidth]{single_glitch_glitchSearch_corner}
+\caption{}
+\label{fig_glitch_posteriors}
+\end{figure}
+
\section{Conclusion}
\label{sec_conclusion}
diff --git a/Paper/single_glitch_F0F1_grid_2D.png b/Paper/single_glitch_F0F1_grid_2D.png
new file mode 100644
index 0000000000000000000000000000000000000000..e13b851e3e86640cf2eef82b7ae93939f233211b
Binary files /dev/null and b/Paper/single_glitch_F0F1_grid_2D.png differ
diff --git a/Paper/single_glitch_glitchSearch_corner.png b/Paper/single_glitch_glitchSearch_corner.png
new file mode 100644
index 0000000000000000000000000000000000000000..11d87222720c04248ebd3c6f08187ec2bedc2b74
Binary files /dev/null and b/Paper/single_glitch_glitchSearch_corner.png differ
diff --git a/Paper/transient_search_corner.png b/Paper/transient_search_corner.png
index 7507dd26a823a3a5776dcd2db0d23e5458fb719a..cc829501625520d77e5917b6f6e6be7460f962e8 100644
Binary files a/Paper/transient_search_corner.png and b/Paper/transient_search_corner.png differ
diff --git a/Paper/transient_search_initial_stage_twoFcumulative.png b/Paper/transient_search_initial_stage_twoFcumulative.png
index 991babe08f271c1793806d365385cdd912eb4d2b..987c54ae52727f6cf61c2136a1545ed95fb00fc9 100644
Binary files a/Paper/transient_search_initial_stage_twoFcumulative.png and b/Paper/transient_search_initial_stage_twoFcumulative.png differ
diff --git a/pyfstat.py b/pyfstat.py
index df2e112610a91d7a1b65e1ed4579348e6abf4b64..2017133c0e0ae86e0ddd9a7489c0fc410ba3f481 100755
--- a/pyfstat.py
+++ b/pyfstat.py
@@ -2365,8 +2365,8 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch):
nsteps = rs[0][0]
else:
nsteps = '{},{}'.format(*rs[0])
- line = line.format(i, rs[1], Tcoh, nsteps,
- texify_float(Vpe))
+ line = line.format(i, rs[1], '{:1.1f}'.format(Tcoh),
+ nsteps, texify_float(Vpe))
f.write(line)
f.write(r'\end{tabular}' + '\n')