diff --git a/Paper/Examples/follow_up.py b/Paper/Examples/follow_up.py
new file mode 100644
index 0000000000000000000000000000000000000000..3093daa0ca7e9219e695c87f74a8ee1249833b4b
--- /dev/null
+++ b/Paper/Examples/follow_up.py
@@ -0,0 +1,76 @@
+import pyfstat
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib
+
+F0 = 30.0
+F1 = -1e-10
+F2 = 0
+Alpha = 1.0
+Delta = 0.5
+
+# Properties of the GW data
+sqrtSX = 1e-23
+tstart = 1000000000
+duration = 100*86400
+tend = tstart+duration
+tref = .5*(tstart+tend)
+
+depth = 50
+data_label = 'follow_up'
+
+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()
+
+# The predicted twoF, given by lalapps_predictFstat can be accessed by
+twoF = data.predict_fstat()
+print 'Predicted twoF value: {}\n'.format(twoF)
+
+# Search
+VF0 = VF1 = 500
+DeltaF0 = VF0 * np.sqrt(3)/(np.pi*duration)
+DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*duration**2)
+DeltaAlpha = 1e-1
+DeltaDelta = 1e-1
+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': {'type': 'unif', 'lower': Alpha-DeltaAlpha,
+ 'upper': Alpha+DeltaAlpha},
+ 'Delta': {'type': 'unif', 'lower': Delta-DeltaDelta,
+ 'upper': Delta+DeltaDelta},
+ }
+
+ntemps = 3
+log10temperature_min = -0.5
+nwalkers = 100
+scatter_val = 1e-10
+nsteps = [200, 200]
+
+mcmc = pyfstat.MCMCFollowUpSearch(
+ label='follow_up', outdir='data',
+ sftfilepath='data/*'+data_label+'*sft', theta_prior=theta_prior, tref=tref,
+ minStartTime=tstart, maxStartTime=tend, nwalkers=nwalkers, nsteps=nsteps,
+ ntemps=ntemps, log10temperature_min=log10temperature_min,
+ scatter_val=scatter_val)
+
+fig, axes = plt.subplots(nrows=2, ncols=2)
+fig, axes = mcmc.run(
+ R=10, Nsegs0=100, subtractions=[F0, F1, Alpha, Delta], labelpad=0.01,
+ fig=fig, axes=axes, plot_det_stat=False, return_fig=True)
+axes[3].set_xlabel(r'$\textrm{Number of steps}$', labelpad=0.1)
+for ax in axes:
+ ax.set_xlim(0, axes[0].get_xlim()[-1])
+ ax.xaxis.set_major_locator(matplotlib.ticker.MaxNLocator(5))
+fig.tight_layout()
+fig.savefig('{}/{}_walkers.png'.format(mcmc.outdir, mcmc.label), dpi=400)
+
+mcmc.plot_corner(add_prior=True)
+mcmc.print_summary()
diff --git a/Paper/Examples/matplotlibrc b/Paper/Examples/matplotlibrc
new file mode 100644
index 0000000000000000000000000000000000000000..c2f237718bc52015e0edaf8a286e6f1371d99973
--- /dev/null
+++ b/Paper/Examples/matplotlibrc
@@ -0,0 +1,36 @@
+figure.figsize: 3.4, 2.5
+figure.facecolor: w
+
+axes.labelsize: 6
+axes.titlesize: 6
+xtick.labelsize: 4
+ytick.labelsize: 4
+legend.fontsize: 6
+font.size: 6
+
+grid.linewidth: 0.8
+lines.linewidth: 1.0
+patch.linewidth: 0.24
+lines.markersize: 3.6
+lines.markeredgewidth: 0
+
+xtick.major.size: 2.0
+ytick.major.size: 2.0
+xtick.minor.size: 1.5
+xtick.minor.visible: True
+ytick.minor.size: 1.5
+ytick.minor.visible: True
+
+xtick.major.pad: 2.5
+ytick.major.pad: 2.5
+
+pdf.compression: 9
+savefig.format: png
+savefig.dpi: 300
+
+font.family: serif
+font.serif: Computer Modern
+text.usetex: True
+axes.formatter.use_mathtext: True
+axes.formatter.useoffset: False
+axes.formatter.limits : -3, 4
diff --git a/Paper/follow_up_run_setup.tex b/Paper/follow_up_run_setup.tex
new file mode 100644
index 0000000000000000000000000000000000000000..711cb402182e5f06fc668f4b400c092986be34f3
--- /dev/null
+++ b/Paper/follow_up_run_setup.tex
@@ -0,0 +1,11 @@
+\begin{tabular}{c|cccccc}
+Stage & $\Nseg$ & $\Tcoh^{\rm days}$ &$\Nsteps$ & $\V$ & $\Vsky$ & $\Vpe$ \\ \hline
+0 & 100 & 1.0 & 200 & $1{\times}10^{2}$ & 7.0 & 10.0 \\
+1 & 56 & 1.8 & 200 & $1{\times}10^{3}$ & 20.0 & 40.0 \\
+2 & 31 & 3.2 & 200 & $1{\times}10^{4}$ & 70.0 & $1{\times}10^{2}$ \\
+3 & 17 & 5.9 & 200 & $1{\times}10^{5}$ & $2{\times}10^{2}$ & $5{\times}10^{2}$ \\
+4 & 9 & 11.1 & 200 & $1{\times}10^{6}$ & $8{\times}10^{2}$ & $2{\times}10^{3}$ \\
+5 & 5 & 20.0 & 200 & $1{\times}10^{7}$ & $2{\times}10^{3}$ & $5{\times}10^{3}$ \\
+6 & 2 & 50.0 & 200 & $1{\times}10^{8}$ & $3{\times}10^{3}$ & $3{\times}10^{4}$ \\
+7 & 1 & 100.0 & 200,200 & $3{\times}10^{8}$ & $4{\times}10^{3}$ & $6{\times}10^{4}$ \\
+\end{tabular}
diff --git a/Paper/follow_up_walkers.png b/Paper/follow_up_walkers.png
new file mode 100644
index 0000000000000000000000000000000000000000..9f8511b98ebfde60cae5ffaa875c3d909c54582d
Binary files /dev/null and b/Paper/follow_up_walkers.png differ
diff --git a/Paper/paper_cw_mcmc.tex b/Paper/paper_cw_mcmc.tex
index 3f947a548bb7336aab7fba35481efcb0a1fa1c4e..b336515a9d04338fbad17001ccb2524310e176e5 100644
--- a/Paper/paper_cw_mcmc.tex
+++ b/Paper/paper_cw_mcmc.tex
@@ -824,13 +824,13 @@ 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_weak_signal_follow_up}. In addition, we show the
+layed out in Table~\ref{tab_signal_follow_up}. 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_weak_signal_follow_up}
-\input{weak_signal_follow_up_run_setup}
+\label{tab_follow_up}
+\input{follow_up_run_setup}
\end{table}
The choice of $\mathcal{R}$ and $\V^{\rm min}$ is a compromise between the
@@ -842,7 +842,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_weak_signal_follow_up}, during
+the follow-up. As expected from Table~\ref{tab_follow_up}, 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
@@ -851,13 +851,13 @@ however this is due to the changing way that the Gaussian noise adds to the sign
Eventually, the walkers all converge to the true signal.
\begin{figure}[htb]
\centering
-\includegraphics[width=0.4\textwidth]{weak_signal_follow_up_walkers}
+\includegraphics[width=0.5\textwidth]{follow_up_walkers}
\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_weak_signal_follow_up}.}
+are listed in Table~\ref{tab_follow_up}.}
\label{fig_follow_up}
\end{figure}
@@ -877,8 +877,37 @@ 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]$.
+
+
+
\subsection{Follow-up candidates from a directed search}
+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
+
+
\begin{table}[htb]
\caption{Run-setup for the directed follow-up Monte-Carlo study, generated with
$\mathcal{R}=10$ and $\V^{\rm min}=100$}
@@ -894,6 +923,25 @@ come from random draws searches using the setup described in
Table~\ref{tab_directed_MC_follow_up}.}
\label{fig_directed_MC_follow_up}
\end{figure}
+
+\subsection{Follow-up candidates from an all-sky search}
+
+\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$}
+\label{tab_allsky_MC_follow_up}
+\input{AllSky_run_setup}
+\end{table}
+
+\begin{figure}[htb]
+\centering
+\includegraphics[width=0.45\textwidth]{allsky_recovery}
+\caption{Recovery fraction for the all-sky follow-up. The Monte-Carlo results
+come from random draws searches using the setup described in
+Table~\ref{tab_directed_MC_follow_up}.}
+\label{fig_allsky_MC_follow_up}
+\end{figure}
+
\section{Alternative waveform models}
In a grided search, the template bank is constructed to ensure that a canonical
diff --git a/Paper/weak_signal_follow_up_run_setup.tex b/Paper/weak_signal_follow_up_run_setup.tex
deleted file mode 100644
index a43fb4a45d09ee85d526acbd63a7e59264aa5d92..0000000000000000000000000000000000000000
--- a/Paper/weak_signal_follow_up_run_setup.tex
+++ /dev/null
@@ -1,9 +0,0 @@
-\begin{tabular}{c|cccccc}
-Stage & $\Nseg$ & $\Tcoh^{\rm days}$ &$\Nsteps$ & $\V$ & $\Vsky$ & $\Vpe$ \\ \hline
-0 & 93 & 1.1 & 100 & 20.0 & 2.0 & 10.0 \\
-1 & 43 & 2.3 & 100 & $2{\times}10^{2}$ & 10.0 & 20.0 \\
-2 & 20 & 5.0 & 100 & $2{\times}10^{3}$ & 50.0 & 50.0 \\
-3 & 9 & 11.1 & 100 & $2{\times}10^{4}$ & $2{\times}10^{2}$ & $1{\times}10^{2}$ \\
-4 & 4 & 25.0 & 100 & $2{\times}10^{5}$ & $1{\times}10^{3}$ & $2{\times}10^{2}$ \\
-5 & 1 & 100.0 & 100,100 & $2{\times}10^{6}$ & $3{\times}10^{3}$ & $9{\times}10^{2}$ \\
-\end{tabular}
diff --git a/Paper/weak_signal_follow_up_walkers.png b/Paper/weak_signal_follow_up_walkers.png
deleted file mode 100644
index d901ea5dce2cc0b572741a0c426c12d435ebe792..0000000000000000000000000000000000000000
Binary files a/Paper/weak_signal_follow_up_walkers.png and /dev/null differ
diff --git a/examples/weak_signal_follow_up.py b/examples/weak_signal_follow_up.py
index 426981b850fa0fce7922f4d1093c6e06b6280fb3..5aa1e1d11f909c59a9e7c1de4d7eb44e86b98f1c 100644
--- a/examples/weak_signal_follow_up.py
+++ b/examples/weak_signal_follow_up.py
@@ -1,4 +1,6 @@
import pyfstat
+import numpy as np
+import matplotlib.pyplot as plt
F0 = 30.0
F1 = -1e-10
@@ -29,15 +31,20 @@ twoF = data.predict_fstat()
print 'Predicted twoF value: {}\n'.format(twoF)
# Search
-theta_prior = {'F0': {'type': 'unif', 'lower': F0*(1-1e-6),
- 'upper': F0*(1+1e-6)},
- 'F1': {'type': 'unif', 'lower': F1*(1+1e-2),
- 'upper': F1*(1-1e-2)},
+VF0 = VF1 = 500
+DeltaF0 = VF0 * np.sqrt(3)/(np.pi*duration)
+DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*duration**2)
+DeltaAlpha = 1e-1
+DeltaDelta = 1e-1
+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': {'type': 'unif', 'lower': Alpha-1e-2,
- 'upper': Alpha+1e-2},
- 'Delta': {'type': 'unif', 'lower': Delta-1e-2,
- 'upper': Delta+1e-2},
+ 'Alpha': {'type': 'unif', 'lower': Alpha-DeltaAlpha,
+ 'upper': Alpha+DeltaAlpha},
+ 'Delta': {'type': 'unif', 'lower': Delta-DeltaDelta,
+ 'upper': Delta+DeltaDelta},
}
ntemps = 3
@@ -52,6 +59,11 @@ mcmc = pyfstat.MCMCFollowUpSearch(
minStartTime=tstart, maxStartTime=tend, nwalkers=nwalkers, nsteps=nsteps,
ntemps=ntemps, log10temperature_min=log10temperature_min,
scatter_val=scatter_val)
-mcmc.run(R=10, Nsegs0=50, subtractions=[F0, Alpha, Delta], context='paper')
+
+fig, axes = plt.subplots(nrows=2, ncols=2)
+mcmc.run(
+ R=10, Nsegs0=100, subtractions=[F0, F1, Alpha, Delta], context='paper',
+ fig=fig, axes=axes, plot_det_stat=False, return_fig=True)
+
mcmc.plot_corner(add_prior=True)
mcmc.print_summary()
diff --git a/pyfstat.py b/pyfstat.py
index 1ddb98f9ae455e7b021c9f410ee7c4fbf341d6d3..df2e112610a91d7a1b65e1ed4579348e6abf4b64 100755
--- a/pyfstat.py
+++ b/pyfstat.py
@@ -1386,9 +1386,12 @@ class MCMCSearch(BaseSearchClass):
def plot_walkers(self, sampler, symbols=None, alpha=0.4, color="k", temp=0,
lw=0.1, burnin_idx=None, add_det_stat_burnin=False,
fig=None, axes=None, xoffset=0, plot_det_stat=True,
- context='classic', subtractions=None):
+ context='classic', subtractions=None, labelpad=0.05):
""" Plot all the chains from a sampler """
+ if np.ndim(axes) > 1:
+ axes = axes.flatten()
+
shape = sampler.chain.shape
if len(shape) == 3:
nwalkers, nsteps, ndim = shape
@@ -1419,8 +1422,6 @@ class MCMCSearch(BaseSearchClass):
if ndim > 1:
for i in range(ndim):
axes[i].ticklabel_format(useOffset=False, axis='y')
- if i < ndim-1:
- axes[i].set_xticklabels([])
cs = chain[:, :, i].T
if burnin_idx:
axes[i].plot(xoffset+idxs[:burnin_idx],
@@ -1432,10 +1433,11 @@ class MCMCSearch(BaseSearchClass):
color="k", alpha=alpha, lw=lw)
if symbols:
if subtractions[i] == 0:
- axes[i].set_ylabel(symbols[i])
+ axes[i].set_ylabel(symbols[i], labelpad=labelpad)
else:
axes[i].set_ylabel(
- symbols[i]+'$-$'+symbols[i]+'$_0$')
+ symbols[i]+'$-$'+symbols[i]+'$_0$',
+ labelpad=labelpad)
else:
axes[0].ticklabel_format(useOffset=False, axis='y')
@@ -1446,12 +1448,13 @@ class MCMCSearch(BaseSearchClass):
axes[0].plot(idxs[burnin_idx:], cs[burnin_idx:], color="k",
alpha=alpha, lw=lw)
if symbols:
- axes[0].set_ylabel(symbols[0])
+ axes[0].set_ylabel(symbols[0], labelpad=labelpad)
- if len(axes) == ndim:
- axes.append(fig.add_subplot(ndim+1, 1, ndim+1))
if plot_det_stat:
+ if len(axes) == ndim:
+ axes.append(fig.add_subplot(ndim+1, 1, ndim+1))
+
lnl = sampler.lnlikelihood[temp, :, :]
if burnin_idx and add_det_stat_burnin:
burn_in_vals = lnl[:, :burnin_idx].flatten()
@@ -1478,7 +1481,7 @@ class MCMCSearch(BaseSearchClass):
xfmt.set_powerlimits((-4, 4))
axes[-1].xaxis.set_major_formatter(xfmt)
- axes[-2].set_xlabel(r'$\textrm{Number of steps}$', labelpad=0.1)
+ axes[-2].set_xlabel(r'$\textrm{Number of steps}$', labelpad=0.2)
return fig, axes
def apply_corrections_to_p0(self, p0):
@@ -2374,8 +2377,8 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch):
return run_setup
def run(self, run_setup=None, proposal_scale_factor=2, R=10, Nsegs0=None,
- create_plots=True, log_table=True, gen_tex_table=True,
- **kwargs):
+ create_plots=True, log_table=True, gen_tex_table=True, fig=None,
+ axes=None, return_fig=False, **kwargs):
""" Run the follow-up with the given run_setup
Parameters
@@ -2403,8 +2406,6 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch):
self.nsegs = run_setup[-1][1]
return
- fig = None
- axes = None
nsteps_total = 0
for j, ((nburn, nprod), nseg, reset_p0) in enumerate(run_setup):
if j == 0:
@@ -2445,17 +2446,10 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch):
fig, axes = self.plot_walkers(
sampler, symbols=self.theta_symbols, fig=fig, axes=axes,
burnin_idx=nburn, xoffset=nsteps_total, **kwargs)
- for ax in axes[:-1]:
- ax.axvline(nsteps_total, color='k', ls='--')
- nsteps_total += nburn+nprod
+ for ax in axes[:self.ndim]:
+ ax.axvline(nsteps_total, color='k', ls='--', lw=0.25)
- if create_plots:
- try:
- fig.tight_layout()
- except ValueError as e:
- logging.warning('Tight layout encountered {}'.format(e))
- fig.savefig('{}/{}_walkers.png'.format(
- self.outdir, self.label), dpi=200)
+ nsteps_total += nburn+nprod
samples = sampler.chain[0, :, nburn:, :].reshape((-1, self.ndim))
lnprobs = sampler.lnprobability[0, :, nburn:].reshape((-1))
@@ -2466,6 +2460,17 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch):
self.lnlikes = lnlikes
self.save_data(sampler, samples, lnprobs, lnlikes)
+ if create_plots:
+ try:
+ fig.tight_layout()
+ except (ValueError, RuntimeError) as e:
+ logging.warning('Tight layout encountered {}'.format(e))
+ if return_fig:
+ return fig, axes
+ else:
+ fig.savefig('{}/{}_walkers.png'.format(
+ self.outdir, self.label), dpi=200)
+
class MCMCTransientSearch(MCMCSearch):
""" MCMC search for a transient signal using the ComputeFstat """