Merge branch 'master' of gitlab.aei.uni-hannover.de:GregAshton/PyFstat

Paper/Examples/follow_up.py

+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()

Paper/Examples/matplotlibrc

+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

Paper/follow_up_run_setup.tex

+\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}

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

Paper/weak_signal_follow_up_run_setup.tex

-\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}

examples/weak_signal_follow_up.py

 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()

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 """