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
index 3d76a824acacdeee6c36bcb6537ce32dedcf10c5..61b1212a8bd9c49d3ed33e4bb77a91b8978ece62 100644
--- a/Paper/Examples/single_glitch.py
+++ b/Paper/Examples/single_glitch.py
@@ -80,8 +80,8 @@ theta_prior = {'F0': {'type': 'unif', 'lower': F0-DeltaF0/2.,
'Delta': Delta,
'tglitch': {'type': 'unif', 'lower': tstart+0.1*duration,
'upper': tend-0.1*duration},
- 'delta_F0': {'type': 'halfnorm', 'loc': 0, 'scale': 1e-3*F0},
- 'delta_F1': {'type': 'norm', 'loc': 0, 'scale': 1e-3*abs(F1)},
+ 'delta_F0': {'type': 'halfnorm', 'loc': 0, 'scale': DeltaF0},
+ 'delta_F1': {'type': 'norm', 'loc': 0, 'scale': DeltaF1},
}
ntemps = 3
log10temperature_min = -0.1
@@ -94,6 +94,6 @@ glitch_mcmc = pyfstat.MCMCGlitchSearch(
nwalkers=nwalkers, ntemps=ntemps,
log10temperature_min=log10temperature_min)
glitch_mcmc.run()
-glitch_mcmc.plot_corner(figsize=(3.2, 3.2))
+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
index c5c54c3c2fc318318e7b11c7b1c14025b6ded97f..8df9b12b4627b8b7a1367443ae3d3b5c7af37ca2 100644
--- a/Paper/Examples/transient_search_using_MCMC.py
+++ b/Paper/Examples/transient_search_using_MCMC.py
@@ -87,5 +87,5 @@ mcmc = pyfstat.MCMCTransientSearch(
nwalkers=nwalkers, ntemps=ntemps,
log10temperature_min=log10temperature_min)
mcmc.run()
-mcmc.plot_corner(add_prior=True)
+mcmc.plot_corner()
mcmc.print_summary()
diff --git a/Paper/bibliography.bib b/Paper/bibliography.bib
index d94502f196599b4c2949f72782063f7a438b09fc..8996364e5942dcdda3b9e46be021a9ddae43e934 100644
--- a/Paper/bibliography.bib
+++ b/Paper/bibliography.bib
@@ -428,3 +428,39 @@ year = {2015}
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/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 83a18b854cf46a18e1978c62492b74da814672d3..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}}},
@@ -1110,14 +1111,84 @@ 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:}
+\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_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