Adds volume convergence investigation

Paper/VolumeConvergenceInvestigation/VolumeConvergence_repeater.sh

+#!/bin/bash
+
+. /home/gregory.ashton/lalsuite-install/etc/lalapps-user-env.sh
+export PATH="/home/gregory.ashton/anaconda2/bin:$PATH"
+export MPLCONFIGDIR=/home/gregory.ashton/.config/matplotlib
+
+for ((n=0;n<1;n++))
+do
+/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --no-template-counting --no-interactive
+done
+mv /local/user/gregory.ashton/VolumeConvergenceResults_"$1".txt $(pwd)/CollectedOutput

Paper/VolumeConvergenceInvestigation/generate_data.py

+import pyfstat
+import numpy as np
+import os
+import sys
+import time
+
+ID = sys.argv[1]
+outdir = sys.argv[2]
+
+label = 'VCrun_{}'.format(ID)
+data_label = 'VCrunData_{}_'.format(ID)
+results_file_name = '{}/VolumeConvergenceResults_{}.txt'.format(outdir, ID)
+
+# Properties of the GW data
+sqrtSX = 1e-23
+tstart = 1000000000
+Tspan = 100*86400
+tend = tstart + Tspan
+
+# Fixed properties of the signal
+F0 = 30
+F1 = -1e-10
+F2 = 0
+Alpha = np.radians(83.6292)
+Delta = np.radians(22.0144)
+tref = .5*(tstart+tend)
+
+depth = 100
+
+h0 = sqrtSX / float(depth)
+
+psi = np.random.uniform(-np.pi/4, np.pi/4)
+phi = np.random.uniform(0, 2*np.pi)
+cosi = np.random.uniform(-1, 1)
+
+data = pyfstat.Writer(
+    label=data_label, outdir=outdir, tref=tref,
+    tstart=tstart, F0=F0, F1=F1, F2=F2, duration=Tspan, Alpha=Alpha,
+    Delta=Delta, h0=h0, sqrtSX=sqrtSX, psi=psi, phi=phi, cosi=cosi,
+    detectors='L1')
+data.make_data()
+twoF_PM = data.predict_fstat()
+
+nsteps = [500, 500]
+
+Vs = np.arange(50, 550, 50)
+Vs = [50, 150, 200, 250, 350, 400, 450]
+
+for V in Vs:
+
+    DeltaF0 = V * np.sqrt(3)/(np.pi*Tspan)
+    DeltaF1 = V * np.sqrt(45/4.)/(np.pi*Tspan**2)
+
+    startTime = time.time()
+    theta_prior = {'F0': {'type': 'unif',
+                          'lower': F0-DeltaF0/2.0,
+                          'upper': F0+DeltaF0/2.0},
+                   'F1': {'type': 'unif',
+                          'lower': F1-DeltaF1/2.0,
+                          'upper': F1+DeltaF1/2.0},
+                   'F2': F2,
+                   'Alpha': Alpha,
+                   'Delta': Delta
+                   }
+
+    ntemps = 3
+    log10temperature_min = -0.5
+    nwalkers = 100
+
+    mcmc = pyfstat.MCMCSearch(
+        label=label, outdir=outdir, nsteps=nsteps,
+        sftfilepath='{}/*{}*sft'.format(outdir, data_label),
+        theta_prior=theta_prior,
+        tref=tref, minStartTime=tstart, maxStartTime=tend,
+        nwalkers=nwalkers, ntemps=ntemps,
+        log10temperature_min=log10temperature_min)
+    mcmc.setup_convergence_testing(
+        convergence_period=10, convergence_length=10,
+        convergence_burnin_fraction=0.1, convergence_threshold_number=5,
+        convergence_threshold=1.1, convergence_early_stopping=False)
+    mcmc.run(create_plots=False,
+             log_table=False, gen_tex_table=False
+             )
+    mcmc.print_summary()
+    cdF0, cdF1 = mcmc.convergence_diagnostic[-1]
+    d, maxtwoF = mcmc.get_max_twoF()
+    d_med_std = mcmc.get_median_stds()
+    dF0 = F0 - d['F0']
+    dF1 = F1 - d['F1']
+    F0_std = d_med_std['F0_std']
+    F1_std = d_med_std['F1_std']
+    runTime = time.time() - startTime
+    with open(results_file_name, 'a') as f:
+        f.write('{} {:1.8e} {:1.8e} {:1.8e} {:1.8e} {:1.8e} {:1.8e} {:1.8e} {:1.8e} {:1.5e} {:1.5e} {}\n'
+                .format(V, dF0, dF1, F0_std, F1_std, DeltaF0, DeltaF1, maxtwoF, twoF_PM, cdF0, cdF1,
+                        runTime))
+    os.system('rm {}/*{}*'.format(outdir, label))
+
+os.system('rm {}/*{}*'.format(outdir, data_label))

Paper/VolumeConvergenceInvestigation/plot_data.py

+import matplotlib.pyplot as plt
+import pandas as pd
+import numpy as np
+import os
+from tqdm import tqdm
+from oct2py import octave
+import glob
+import scipy.stats
+
+filenames = glob.glob("CollectedOutput/*.txt")
+
+plt.style.use('paper')
+
+Tspan = 100 * 86400
+
+df_list = []
+for fn in filenames:
+    if '3356003' in fn:
+        continue
+    else:
+        df = pd.read_csv(
+            fn, sep=' ', names=['V', 'dF0', 'dF1', 'F0_std', 'F1_std', 'DeltaF0',
+                                'DeltaF1', 'MaxtwoF', 'twoF_PM', 'cdF0', 'cdF1', 'runTime'])
+        df['mismatch'] = (df.twoF_PM - df.MaxtwoF)/df.twoF_PM
+        df['F0_fraction'] = df.F0_std / df.DeltaF0
+        df['F1_fraction'] = df.F1_std / df.DeltaF1
+        df['CLUSTER_ID'] = fn.split('_')[1]
+        df['RUN_ID'] = fn.split('_')[2]
+        df_list.append(df)
+
+df = pd.concat(df_list)
+
+df.sort_values('V', inplace=True)
+
+df['cd_ave'] = df[['cdF0', 'cdF1']].mean(axis=1)
+df['fraction_ave'] = df[['F0_fraction', 'F1_fraction']].mean(axis=1)
+df = df[df.V <= 500]
+
+fig, (ax2, ax3) = plt.subplots(nrows=2, figsize=(3.2, 4), sharex=True)
+
+print df.groupby('V')['dF0'].count()
+
+
+Vs = df['V'].unique()
+mismatch_ave = df.groupby('V')['mismatch'].mean()
+ave_cd_ave = df.groupby('V')['cd_ave'].apply(scipy.stats.mstats.gmean)
+cdF0_ave = df.groupby('V')['cdF0'].mean()
+cdF1_ave = df.groupby('V')['cdF1'].mean()
+ave_fraction_ave = df.groupby('V')['fraction_ave'].median()
+F0_fraction_ave = df.groupby('V')['F0_fraction'].mean()
+F1_fraction_ave = df.groupby('V')['F1_fraction'].mean()
+
+ax2.plot(df.V, df.cd_ave, 'o', color='k', alpha=0.5)
+#ax2.plot(Vs, ave_cd_ave, color='r')
+ax2.set_ylim(0, 10)
+ax2.set_ylabel(r'$\langle \textrm{PSRF} \rangle$')
+
+ax3.plot(df.V, df.fraction_ave, 'o', color='k', alpha=0.2)
+#ax3.plot(Vs, ave_fraction_ave, color='r')
+ax3.set_ylabel(r'$\langle \textrm{posterior std. / prior width} \rangle$')
+ax3.set_xlabel(
+    r'$\mathcal{V}_\mathrm{PE}^{(0)}=\mathcal{V}_\mathrm{PE}^{(1)}=\sqrt{\mathcal{V}}$')
+ax3.set_xlim(0, 525)
+
+ax3.set_xticks(Vs)
+
+fig.tight_layout()
+fig.savefig('VolumeConvergenceInvestigation.png')

Paper/VolumeConvergenceInvestigation/submitfile

+Executable=VolumeConvergence_repeater.sh
+Arguments=$(Cluster)_$(Process)
+Universe=vanilla
+Input=/dev/null
+accounting_group=ligo.dev.o2.cw.directedisolatedother.semicoherent 
+Output=CollectedOutput/out.$(Cluster).$(Process)
+Error=CollectedOutput/err.$(Cluster).$(Process)
+Log=CollectedOutput/log.$(Cluster).$(Process)
+request_cpus = 1
+request_memory = 16 GB
+
+Queue 100

Paper/definitions.tex

@@ -32,4 +32,5 @@
 \newcommand{\smax}{s_{\textrm{max}}}
 \newcommand{\fmax}{f_{\textrm{max}}}
 \newcommand{\rhohatmax}{\hat{\rho}_{\mathrm{max}}}
+\newcommand{\logmintemp}{\log_{10}(T_\textrm{min})}
 

Paper/paper_cw_mcmc.tex

@@ -888,6 +888,61 @@ chains remain isolated from the main bulk can be easily countered by the use
 of parallel tempering (see Sec.~\ref{sec_parallel_tempering}). This was intenionally
 not used here to highlight the utility of the convergence statistic.
 
+\section{Choice of the prior metric volume}
+
+The prior volume estimates the approximate size of the parameter space relative
+to the size of a signal; any optimization routine, given finite resources,
+will fail if this volume is too large. We will now investigate precisely what
+this means for an MCMC search given a fixed search setup.
+
+In a data set of 100 days, we simulate Gaussian noise and a signal with a
+sensitivity depth of 100. Then, we perform a directed search (over $f$ and
+$\fdot$) with 500 burn-in and 500 production steps, 100 walkers, and using
+parallel tempering with $\Ntemps=3, \logmintemp=-3$. The prior for the search
+is uniform in frequency and spin-down with the width defined by an equal and
+fixed metric volume (i.e. $\Vpe^{(0)}=\Vpe^{(1)}$). This process is repeated
+500 times for a range of prior volumes so that we can infer the statistical
+behaviour. In the upper panel of Fig.~\ref{fig_volume_convergence} we plot the
+PSRF (calculated for a subset of the production samples) averaged over the two
+search dimensions fo each simulation; this demonstrates that for $\sqrt{\V}
+\gtrsim 200$, the majority of simulations where not converged (i.e. having
+PSRF < 1.2). For $\sqrt{\V}=100$, only one simulation of the 500 failed, while
+for $\sqrt{V}=50$ all simulations reached convergence.
+
+\begin{figure}[htb]
+\centering
+\includegraphics[width=0.5\textwidth]{VolumeConvergenceInvestigation}
+\caption{Results of a study into the convergence of MCMC simulations as a
+function of the prior volume}
+\label{fig_volume_convergence}
+\end{figure}
+
+It is reasonable to ask ``what happens for these cases that don't converge?''.
+For all the simulations considered here, the MCMC algorithim did find the peak
+in the likelhood: this was demonstrated by checking that the minimum mismatch
+for each simulation was $\lesssim 0.1$. However, for those simulations which
+where not converged, this was true only for a small subset of the walkers, the
+other walkers where sampling other points in the parameter space (away from the
+signal) and hence had not converged. This is demonstrated in the lower panel of
+Fig.~\ref{fig_volume_convergence} where we plot the ratio of the posterior
+standard deviation to the prior width (averaged over the two search
+dimensions). For the converged solutions, this is small indicating that the
+posterior is narrower than the prior. However, for the solution which did not
+converge, this ratio is larger; if none of the walkers had converged and hence
+the posterior samples where equivalent to the prior (i.e. uniform over some
+range), we would expect this ratio to be $\sim 0.29$. For the largest prior
+volumes considered here the ratio approaches this value.
+This equivalence of the posterior to the prior is expected behaviour since the
+walkers are initialised from the uniform prior (which is much wider than the
+signal); if we had extended the burn-in period, we would expect eventually for
+these walkers to converge: the performance at a fixed prior volume will depend
+on the setup.
+
+For a given setup, one can choose the optimal prior volume (i.e. balancing the
+risk of loosing the signal due to non-convergence against computing cost) by
+stipulating a failure rate; for this example one may choose $\sqrt{\V}=100$ based
+on an estimated failure rate of 1 in 500.
+
 
 \section{Zoom follow-up}
 \label{sec_follow_up}