Update All Sky MC table and plots

Paper/AllSkyMC/generate_data.py

@@ -2,6 +2,7 @@ import pyfstat
 import numpy as np
 import os
 import sys
+import time
 
 ID = sys.argv[1]
 outdir = sys.argv[2]
@@ -29,21 +30,20 @@ DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)
 
 depths = np.linspace(100, 400, 7)
 
-run_setup = [((100, 0), 27, False),
-             ((100, 0), 15, False),
-             ((100, 0), 8, False),
-             ((100, 0), 4, False),
-             ((50, 50), 1, False)]
+nsteps = 50
+run_setup = [((nsteps, 0), 20, False),
+             ((nsteps, 0), 11, False),
+             ((nsteps, 0), 6, False),
+             ((nsteps, 0), 3, False),
+             ((nsteps, nsteps), 1, False)]
 
 DeltaAlpha = 0.05
 DeltaDelta = 0.05
 
 for depth in depths:
     h0 = sqrtSX / float(depth)
-    r = np.random.uniform(0, 1)
-    theta = np.random.uniform(0, 2*np.pi)
-    F0 = F0_center + 3*np.sqrt(r)*np.cos(theta)/(np.pi**2 * Tspan**2)
-    F1 = F1_center + 45*np.sqrt(r)*np.sin(theta)/(4*np.pi**2 * Tspan**4)
+    F0 = F0_center + np.random.uniform(-0.5, 0.5)*DeltaF0
+    F1 = F1_center + np.random.uniform(-0.5, 0.5)*DeltaF1
     Alpha = np.random.uniform(0, 2*np.pi)
     Delta = np.arccos(2*np.random.uniform(0, 1)-1)-np.pi/2
     fAlpha = np.random.uniform(0, 1)
@@ -65,6 +65,7 @@ for depth in depths:
     data.make_data()
     predicted_twoF = data.predict_fstat()
 
+    startTime = time.time()
     theta_prior = {'F0': {'type': 'unif',
                           'lower': F0-DeltaF0/2.,
                           'upper': F0+DeltaF0/2.},
@@ -96,7 +97,8 @@ for depth in depths:
     d, maxtwoF = mcmc.get_max_twoF()
     dF0 = F0 - d['F0']
     dF1 = F1 - d['F1']
+    runTime = time.time() - startTime
     with open(results_file_name, 'a') as f:
-        f.write('{} {:1.8e} {:1.8e} {:1.8e} {:1.8e} {:1.8e}\n'
-                .format(depth, h0, dF0, dF1, predicted_twoF, maxtwoF))
+        f.write('{} {:1.8e} {:1.8e} {:1.8e} {:1.8e} {:1.8e} {}\n'
+                .format(depth, h0, dF0, dF1, predicted_twoF, maxtwoF, runTime))
     os.system('rm {}/*{}*'.format(outdir, label))

Paper/AllSkyMC/generate_table.py

@@ -3,7 +3,7 @@ import numpy as np
 
 outdir = 'data'
 
-label = 'AllSky'
+label = 'allsky_setup'
 data_label = '{}_data'.format(label)
 
 # Properties of the GW data
@@ -23,32 +23,30 @@ VF0 = VF1 = 100
 DeltaF0 = VF0 * np.sqrt(3)/(np.pi*Tspan)
 DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)
 
-depths = np.linspace(100, 400, 7)
-
-run_setup = [((100, 0), 27, False),
-             ((100, 0), 15, False),
-             ((100, 0), 8, False),
-             ((100, 0), 4, False),
-             ((50, 50), 1, False)]
-
 DeltaAlpha = 0.05
 DeltaDelta = 0.05
 
 depth = 100
 
+nsteps = 50
+run_setup = [((nsteps, 0), 20, False),
+             ((nsteps, 0), 11, False),
+             ((nsteps, 0), 6, False),
+             ((nsteps, 0), 3, False),
+             ((nsteps, nsteps), 1, False)]
+
 h0 = sqrtSX / float(depth)
-F0 = F0_center
-F1 = F1_center
+r = np.random.uniform(0, 1)
+theta = np.random.uniform(0, 2*np.pi)
+F0 = F0_center + 3*np.sqrt(r)*np.cos(theta)/(np.pi**2 * Tspan**2)
+F1 = F1_center + 45*np.sqrt(r)*np.sin(theta)/(4*np.pi**2 * Tspan**4)
+
 Alpha = 0
 Delta = 0
-Alpha_min = Alpha - DeltaAlpha/2
-Alpha_max = Alpha + DeltaAlpha/2
-Delta_min = Delta - DeltaDelta/2
-Delta_max = Delta + DeltaDelta/2
 
-psi = 0
-phi = 0
-cosi = 0
+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,
@@ -66,11 +64,11 @@ theta_prior = {'F0': {'type': 'unif',
                       'upper': F1+DeltaF1/2.},
                'F2': F2,
                'Alpha': {'type': 'unif',
-                         'lower': Alpha_min,
-                         'upper': Alpha_max},
+                         'lower': Alpha-DeltaAlpha/2.,
+                         'upper': Alpha+DeltaAlpha/2.},
                'Delta': {'type': 'unif',
-                         'lower': Delta_min,
-                         'upper': Delta_max},
+                         'lower': Delta-DeltaDelta/2.,
+                         'upper': Delta+DeltaDelta/2.},
                }
 
 ntemps = 1
@@ -84,4 +82,4 @@ mcmc = pyfstat.MCMCFollowUpSearch(
     tref=tref, minStartTime=tstart, maxStartTime=tend,
     nwalkers=nwalkers, ntemps=ntemps,
     log10temperature_min=log10temperature_min)
-mcmc.run(run_setup=run_setup)
+mcmc.run(run_setup)

Paper/AllSkyMC/plot_data.py

@@ -4,6 +4,9 @@ import numpy as np
 import os
 from tqdm import tqdm
 from oct2py import octave
+import glob
+
+filenames = glob.glob("CollectedOutput/*.txt")
 
 plt.style.use('paper')
 
@@ -25,11 +28,14 @@ def binomialConfidenceInterval(N, K, confidence=0.95):
     [l, u] =  octave.eval(cmd, verbose=False, return_both=True)[0].split('\n')
     return float(l.split('=')[1]), float(u.split('=')[1])
 
-results_file_name = 'MCResults.txt'
-
-df = pd.read_csv(
-    results_file_name, sep=' ', names=['depth', 'h0', 'dF0', 'dF1',
-                                       'twoF_predicted', 'twoF'])
+df_list = []
+for fn in filenames:
+    df = pd.read_csv(
+        fn, sep=' ', names=['depth', 'h0', 'dF0', 'dF1', 'twoF_predicted',
+                            'twoF', 'runTime'])
+    df['CLUSTER_ID'] = fn.split('_')[1]
+    df_list.append(df)
+df = pd.concat(df_list)
 
 twoFstar = 60
 depths = np.unique(df.depth.values)
@@ -71,3 +77,10 @@ ax.legend(loc=1, frameon=False)
 
 fig.tight_layout()
 fig.savefig('allsky_recovery.png')
+
+
+fig, ax = plt.subplots()
+ax.hist(df.runTime, bins=20)
+ax.set_xlabel('runTime per follow-up [s]')
+fig.savefig('runTimeHist.png')
+

Paper/AllSkyMC/submitfile

@@ -9,4 +9,4 @@ Log=CollectedOutput/log.$(Process)
 request_cpus = 1
 request_memory = 16 GB
 
-Queue 10
+Queue 90

Paper/allsky_setup_run_setup.tex

+\begin{tabular}{c|cccccc}
+Stage & $\Nseg$ & $\Tcoh^{\rm days}$ &$\Nsteps$ & $\V$ & $\Vsky$ & $\Vpe$ \\ \hline
+0 & 20 & 5.0 & 100 & $2{\times}10^{2}$ & 10.0 & 10.0 \\
+1 & 11 & 9.1 & 100 & $2{\times}10^{3}$ & 40.0 & 50.0 \\
+2 & 6 & 16.7 & 100 & $2{\times}10^{4}$ & $1{\times}10^{2}$ & $2{\times}10^{2}$ \\
+3 & 3 & 33.3 & 100 & $1{\times}10^{5}$ & $2{\times}10^{2}$ & $6{\times}10^{2}$ \\
+4 & 1 & 100.0 & 100,100 & $8{\times}10^{5}$ & $3{\times}10^{2}$ & $2{\times}10^{3}$ \\
+\end{tabular}