diff --git a/Paper/VolumeConvergenceInvestigation.png b/Paper/VolumeConvergenceInvestigation.png new file mode 100644 index 0000000000000000000000000000000000000000..bb30ba721db0fa546b945861ce479d1cecb7854c Binary files /dev/null and b/Paper/VolumeConvergenceInvestigation.png differ diff --git a/Paper/VolumeConvergenceInvestigation/VolumeConvergence_repeater.sh b/Paper/VolumeConvergenceInvestigation/VolumeConvergence_repeater.sh new file mode 100755 index 0000000000000000000000000000000000000000..7a50c2490f56a23d0d9f2eaa37218d16e9ba8bcd --- /dev/null +++ b/Paper/VolumeConvergenceInvestigation/VolumeConvergence_repeater.sh @@ -0,0 +1,11 @@ +#!/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 diff --git a/Paper/VolumeConvergenceInvestigation/generate_data.py b/Paper/VolumeConvergenceInvestigation/generate_data.py new file mode 100644 index 0000000000000000000000000000000000000000..413027d569979609ece20d3707d2320cb7d1ccf2 --- /dev/null +++ b/Paper/VolumeConvergenceInvestigation/generate_data.py @@ -0,0 +1,99 @@ +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)) diff --git a/Paper/VolumeConvergenceInvestigation/plot_data.py b/Paper/VolumeConvergenceInvestigation/plot_data.py new file mode 100644 index 0000000000000000000000000000000000000000..a88534298f0081a3f7a91898858ebec4386f0735 --- /dev/null +++ b/Paper/VolumeConvergenceInvestigation/plot_data.py @@ -0,0 +1,68 @@ +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') diff --git a/Paper/VolumeConvergenceInvestigation/submitfile b/Paper/VolumeConvergenceInvestigation/submitfile new file mode 100644 index 0000000000000000000000000000000000000000..e656190d806b6217d70673dab20228e08f733316 --- /dev/null +++ b/Paper/VolumeConvergenceInvestigation/submitfile @@ -0,0 +1,12 @@ +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 diff --git a/Paper/definitions.tex b/Paper/definitions.tex index decc3ee8621b2f800498523b7047a43a4e6b3e19..6001da63333382190bce9c7a4c7baaa2486a027d 100644 --- a/Paper/definitions.tex +++ b/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})} diff --git a/Paper/paper_cw_mcmc.tex b/Paper/paper_cw_mcmc.tex index 8512a596d3108ffaf85e57f97b2b398c2d6bfe43..8b0e9a55fb901469be2c46155568b7b8200aa33f 100644 --- a/Paper/paper_cw_mcmc.tex +++ b/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}