Complete run of MC results

Use N=1000 for the recoveries and N=10000 for the test in noise

Paper/AllSkyMC/repeat.sh → Paper/AllSkyMC/AllSkyMC_repeat.sh

@@ -4,9 +4,8 @@
 export PATH="/home/gregory.ashton/anaconda2/bin:$PATH"
 export MPLCONFIGDIR=/home/gregory.ashton/.config/matplotlib
 
-rm /local/user/gregory.ashton/MCResults*txt 
-for ((n=0;n<10;n++))
+for ((n=0;n<1;n++))
 do
-/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --quite --no-template-counting
+/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --no-template-counting --no-interactive
 done
-cp /local/user/gregory.ashton/MCResults*txt /home/gregory.ashton/PyFstat/Paper/AllSkyMC/CollectedOutput
+cp /local/user/gregory.ashton/MCResults_"$1".txt $(pwd)/CollectedOutput

Paper/AllSkyMC/generate_data.py

@@ -12,14 +12,14 @@ data_label = '{}_data'.format(label)
 results_file_name = '{}/MCResults_{}.txt'.format(outdir, ID)
 
 # Properties of the GW data
-sqrtSX = 2e-23
+sqrtSX = 1e-23
 tstart = 1000000000
 Tspan = 100*86400
 tend = tstart + Tspan
 
 # Fixed properties of the signal
 F0_center = 30
-F1_center = 1e-10
+F1_center = -1e-10
 F2 = 0
 tref = .5*(tstart+tend)
 
@@ -31,7 +31,7 @@ DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)
 DeltaAlpha = 0.02
 DeltaDelta = 0.02
 
-depths = np.linspace(100, 400, 7)
+depths = np.linspace(100, 400, 9)
 
 nsteps = 50
 run_setup = [((nsteps, 0), 20, False),
@@ -45,7 +45,7 @@ for depth in depths:
     h0 = sqrtSX / float(depth)
     F0 = F0_center + np.random.uniform(-0.5, 0.5)*DeltaF0
     F1 = F1_center + np.random.uniform(-0.5, 0.5)*DeltaF1
-    Alpha_center = np.random.uniform(0, 2*np.pi)
+    Alpha_center = np.random.uniform(DeltaAlpha, 2*np.pi-DeltaAlpha)
     Delta_center = np.arccos(2*np.random.uniform(0, 1)-1)-np.pi/2
     Alpha = Alpha_center + np.random.uniform(-0.5, 0.5)*DeltaAlpha
     Delta = Delta_center + np.random.uniform(-0.5, 0.5)*DeltaDelta
@@ -59,7 +59,6 @@ for depth in depths:
         Delta=Delta, h0=h0, sqrtSX=sqrtSX, psi=psi, phi=phi, cosi=cosi,
         detector='H1,L1')
     data.make_data()
-    predicted_twoF = data.predict_fstat()
 
     startTime = time.time()
     theta_prior = {'F0': {'type': 'unif',
@@ -95,6 +94,6 @@ for depth in depths:
     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, runTime))
+        f.write('{} {:1.8e} {:1.8e} {:1.8e} {:1.8e} {}\n'
+                .format(depth, h0, dF0, dF1, maxtwoF, runTime))
     os.system('rm {}/*{}*'.format(outdir, label))

Paper/AllSkyMC/plot_data.py

@@ -31,8 +31,7 @@ def binomialConfidenceInterval(N, K, confidence=0.95):
 df_list = []
 for fn in filenames:
     df = pd.read_csv(
-        fn, sep=' ', names=['depth', 'h0', 'dF0', 'dF1', 'twoF_predicted',
-                            'twoF', 'runTime'])
+        fn, sep=' ', names=['depth', 'h0', 'dF0', 'dF1', 'twoF', 'runTime'])
     df['CLUSTER_ID'] = fn.split('_')[1]
     df_list.append(df)
 df = pd.concat(df_list)
@@ -52,9 +51,9 @@ for d in depths:
 
 yerr = np.abs(recovery_fraction - np.array(recovery_fraction_CI).T)
 fig, ax = plt.subplots()
-ax.errorbar(depths, recovery_fraction, yerr=yerr, fmt='sk', marker='s', ms=2,
+ax.errorbar(depths, recovery_fraction, yerr=yerr, fmt='sr', marker='s', ms=2,
             capsize=1, capthick=0.5, elinewidth=0.5,
-            label='Monte-Carlo result')
+            label='Monte-Carlo result', zorder=10)
 
 fname = 'analytic_data.txt'
 if os.path.isfile(fname):

Paper/AllSkyMC/submitfile

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

Paper/AllSkyMCNoiseOnly/generate_data.py

@@ -9,17 +9,17 @@ outdir = sys.argv[2]
 
 label = 'run_{}'.format(ID)
 data_label = '{}_data'.format(label)
-results_file_name = '{}/MCResults_{}.txt'.format(outdir, ID)
+results_file_name = '{}/NoiseOnlyMCResults_{}.txt'.format(outdir, ID)
 
 # Properties of the GW data
-sqrtSX = 2e-23
+sqrtSX = 1e-23
 tstart = 1000000000
 Tspan = 100*86400
 tend = tstart + Tspan
 
 # Fixed properties of the signal
 F0_center = 30
-F1_center = 1e-10
+F1_center = -1e-10
 F2 = 0
 tref = .5*(tstart+tend)
 
@@ -28,6 +28,9 @@ VF0 = VF1 = 100
 DeltaF0 = VF0 * np.sqrt(3)/(np.pi*Tspan)
 DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)
 
+DeltaAlpha = 0.02
+DeltaDelta = 0.02
+
 nsteps = 50
 run_setup = [((nsteps, 0), 20, False),
              ((nsteps, 0), 11, False),
@@ -35,21 +38,14 @@ run_setup = [((nsteps, 0), 20, False),
              ((nsteps, 0), 3, False),
              ((nsteps, nsteps), 1, False)]
 
-DeltaAlpha = 0.05
-DeltaDelta = 0.05
 
 h0 = 0
 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)
-Alpha_min = Alpha - DeltaAlpha*(1-fAlpha)
-Alpha_max = Alpha + DeltaAlpha*fAlpha
-fDelta = np.random.uniform(0, 1)
-Delta_min = Delta - DeltaDelta*(1-fDelta)
-Delta_max = Delta + DeltaDelta*fDelta
-
+Alpha_center = np.random.uniform(DeltaAlpha, 2*np.pi-DeltaAlpha)
+Delta_center = np.arccos(2*np.random.uniform(0, 1)-1)-np.pi/2
+Alpha = Alpha_center + np.random.uniform(-0.5, 0.5)*DeltaAlpha
+Delta = Delta_center + np.random.uniform(-0.5, 0.5)*DeltaDelta
 psi = np.random.uniform(-np.pi/4, np.pi/4)
 phi = np.random.uniform(0, 2*np.pi)
 cosi = np.random.uniform(-1, 1)
@@ -60,25 +56,24 @@ data = pyfstat.Writer(
     Delta=Delta, h0=h0, sqrtSX=sqrtSX, psi=psi, phi=phi, cosi=cosi,
     detector='H1,L1')
 data.make_data()
-predicted_twoF = data.predict_fstat()
 
 startTime = time.time()
 theta_prior = {'F0': {'type': 'unif',
-                      'lower': F0-DeltaF0/2.,
-                      'upper': F0+DeltaF0/2.},
+                      'lower': F0_center-DeltaF0,
+                      'upper': F0_center+DeltaF0},
                'F1': {'type': 'unif',
-                      'lower': F1-DeltaF1/2.,
-                      'upper': F1+DeltaF1/2.},
+                      'lower': F1_center-DeltaF1,
+                      'upper': F1_center+DeltaF1},
                'F2': F2,
                'Alpha': {'type': 'unif',
-                         'lower': Alpha_min,
-                         'upper': Alpha_max},
+                         'lower': Alpha_center-DeltaAlpha,
+                         'upper': Alpha_center+DeltaAlpha},
                'Delta': {'type': 'unif',
-                         'lower': Delta_min,
-                         'upper': Delta_max},
+                         'lower': Delta_center-DeltaDelta,
+                         'upper': Delta_center+DeltaDelta},
                }
 
-ntemps = 1
+ntemps = 2
 log10temperature_min = -1
 nwalkers = 100
 
@@ -96,6 +91,6 @@ 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} {}\n'
-            .format(dF0, dF1, predicted_twoF, maxtwoF, runTime))
+    f.write('{:1.8e} {:1.8e} {:1.8e} {}\n'
+            .format(dF0, dF1, maxtwoF, runTime))
 os.system('rm {}/*{}*'.format(outdir, label))

Paper/AllSkyMCNoiseOnly/plot_data.py

@@ -24,11 +24,11 @@ def maxtwoFinNoise(twoF, Ntrials):
 df_list = []
 for fn in filenames:
     df = pd.read_csv(
-        fn, sep=' ', names=['dF0', 'dF1', 'twoF_predicted',
-                            'twoF', 'runTime'])
+        fn, sep=' ', names=['dF0', 'dF1', 'twoF', 'runTime'])
     df['CLUSTER_ID'] = fn.split('_')[1]
     df_list.append(df)
 df = pd.concat(df_list)
+print 'Number of samples = ', len(df)
 
 fig, ax = plt.subplots()
 ax.hist(df.twoF, bins=50, histtype='step', color='k', normed=True, linewidth=1)

Paper/AllSkyMCNoiseOnly/submitfile

-Executable= repeat.sh
+Executable= AllSkyMCNoiseOnly_repeat.sh
 Arguments= $(Cluster)_$(Process)
 Universe=vanilla
 Input=/dev/null
 accounting_group = ligo.dev.o2.cw.explore.test
-Output=CollectedOutput/out.$(Process)
-Error=CollectedOutput/err.$(Process)
-Log=CollectedOutput/log.$(Process)
+Output=CollectedOutput/out.$(Cluster).$(Process)
+Error=CollectedOutput/err.$(Cluster).$(Process)
+Log=CollectedOutput/log.$(Cluster).$(Process)
 request_cpus = 1
 request_memory = 16 GB
 

Paper/DirectedMC/repeat.sh → Paper/DirectedMC/DirectedMC_repeat.sh

@@ -4,9 +4,8 @@
 export PATH="/home/gregory.ashton/anaconda2/bin:$PATH"
 export MPLCONFIGDIR=/home/gregory.ashton/.config/matplotlib
 
-rm /local/user/gregory.ashton/MCResults*txt 
-for ((n=0;n<10;n++))
+for ((n=0;n<1;n++))
 do
-/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --quite --no-template-counting
+/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --no-template-counting --no-interactive
 done
-cp /local/user/gregory.ashton/MCResults*txt /home/gregory.ashton/PyFstat/Paper/DirectedMC/CollectedOutput
+cp /local/user/gregory.ashton/MCResults_"$1".txt $(pwd)/CollectedOutput

Paper/DirectedMC/generate_data.py

@@ -12,17 +12,17 @@ data_label = '{}_data'.format(label)
 results_file_name = '{}/MCResults_{}.txt'.format(outdir, ID)
 
 # Properties of the GW data
-sqrtSX = 2e-23
+sqrtSX = 1e-23
 tstart = 1000000000
 Tspan = 100*86400
 tend = tstart + Tspan
 
 # Fixed properties of the signal
 F0_center = 30
-F1_center = 1e-10
+F1_center = -1e-10
 F2 = 0
-Alpha = 5e-3
-Delta = 6e-2
+Alpha = np.radians(83.6292)
+Delta = np.radians(22.0144)
 tref = .5*(tstart+tend)
 
 
@@ -30,9 +30,9 @@ 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)
+depths = np.linspace(100, 400, 9)
 
-nsteps = 50
+nsteps = 25
 run_setup = [((nsteps, 0), 20, False),
              ((nsteps, 0), 7, False),
              ((nsteps, 0), 2, False),
@@ -53,7 +53,6 @@ for depth in depths:
         Delta=Delta, h0=h0, sqrtSX=sqrtSX, psi=psi, phi=phi, cosi=cosi,
         detector='H1,L1')
     data.make_data()
-    predicted_twoF = data.predict_fstat()
 
     startTime = time.time()
     theta_prior = {'F0': {'type': 'unif',
@@ -80,11 +79,12 @@ for depth in depths:
         log10temperature_min=log10temperature_min)
     mcmc.run(run_setup=run_setup, create_plots=False, log_table=False,
              gen_tex_table=False)
+    mcmc.print_summary()
     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, runTime))
+        f.write('{} {:1.8e} {:1.8e} {:1.8e} {:1.8e} {}\n'
+                .format(depth, h0, dF0, dF1, maxtwoF, runTime))
     os.system('rm {}/*{}*'.format(outdir, label))

Paper/DirectedMC/plot_data.py

@@ -31,10 +31,12 @@ def binomialConfidenceInterval(N, K, confidence=0.95):
 df_list = []
 for fn in filenames:
     df = pd.read_csv(
-        fn, sep=' ', names=['depth', 'h0', 'dF0', 'dF1', 'twoF_predicted',
-                            'twoF', 'runTime'])
+        fn, sep=' ', names=['depth', 'h0', 'dF0', 'dF1', 'twoF', 'runTime'])
     df['CLUSTER_ID'] = fn.split('_')[1]
-    df_list.append(df)
+    if len(df) != 9:
+        print len(df), fn
+    else:
+        df_list.append(df)
 df = pd.concat(df_list)
 
 twoFstar = 60
@@ -52,9 +54,9 @@ for d in depths:
 
 yerr = np.abs(recovery_fraction - np.array(recovery_fraction_CI).T)
 fig, ax = plt.subplots()
-ax.errorbar(depths, recovery_fraction, yerr=yerr, fmt='sk', marker='s', ms=2,
+ax.errorbar(depths, recovery_fraction, yerr=yerr, fmt='sr', marker='s', ms=2,
             capsize=1, capthick=0.5, elinewidth=0.5,
-            label='Monte-Carlo result')
+            label='Monte-Carlo result', zorder=10)
 
 fname = 'analytic_data.txt'
 if os.path.isfile(fname):

Paper/DirectedMC/submitfile

-Executable= repeat.sh
-Arguments= $(Cluster)_$(Process)
+Executable=DirectedMC_repeat.sh
+Arguments=$(Cluster)_$(Process)
 Universe=vanilla
 Input=/dev/null
 accounting_group = ligo.dev.o2.cw.explore.test
-Output=CollectedOutput/out.$(Process)
-Error=CollectedOutput/err.$(Process)
-Log=CollectedOutput/log.$(Process)
+Output=CollectedOutput/out.$(Cluster).$(Process)
+Error=CollectedOutput/err.$(Cluster).$(Process)
+Log=CollectedOutput/log.$(Cluster).$(Process)
 request_cpus = 1
-request_memory = 8 GB
+request_memory = 16 GB
 
-Queue 100
+Queue 6

Paper/AllSkyMCNoiseOnly/repeat.sh → Paper/DirectedMCNoiseOnly/DirectedMCNoiseOnly_repeat.sh

@@ -4,9 +4,8 @@
 export PATH="/home/gregory.ashton/anaconda2/bin:$PATH"
 export MPLCONFIGDIR=/home/gregory.ashton/.config/matplotlib
 
-rm /local/user/gregory.ashton/MCResults*txt 
 for ((n=0;n<100;n++))
 do
-/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --quite --no-template-counting
+/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --no-template-counting --no-interactive
 done
-cp /local/user/gregory.ashton/MCResults*txt /home/gregory.ashton/PyFstat/Paper/AllSkyMCNoiseOnly/CollectedOutput
+cp /local/user/gregory.ashton/NoiseOnlyMCResults_"$1".txt $(pwd)/CollectedOutput

Paper/DirectedMCNoiseOnly/generate_data.py

@@ -10,27 +10,27 @@ outdir = sys.argv[2]
 
 label = 'run_{}'.format(ID)
 data_label = '{}_data'.format(label)
-results_file_name = '{}/MCResults_{}.txt'.format(outdir, ID)
+results_file_name = '{}/NoiseOnlyMCResults_{}.txt'.format(outdir, ID)
 
 # Properties of the GW data
-sqrtSX = 2e-23
+sqrtSX = 1e-23
 tstart = 1000000000
 Tspan = 100*86400
 tend = tstart + Tspan
 
 # Fixed properties of the signal
 F0_center = 30
-F1_center = 1e-10
+F1_center = -1e-10
 F2 = 0
-Alpha = 5e-3
-Delta = 6e-2
+Alpha = np.radians(83.6292)
+Delta = np.radians(22.0144)
 tref = .5*(tstart+tend)
 
 VF0 = VF1 = 100
 DeltaF0 = VF0 * np.sqrt(3)/(np.pi*Tspan)
 DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)
 
-nsteps = 20
+nsteps = 25
 run_setup = [((nsteps, 0), 20, False),
              ((nsteps, 0), 7, False),
              ((nsteps, 0), 2, False),
@@ -50,21 +50,20 @@ data = pyfstat.Writer(
     Delta=Delta, h0=h0, sqrtSX=sqrtSX, psi=psi, phi=phi, cosi=cosi,
     detector='H1,L1')
 data.make_data()
-predicted_twoF = data.predict_fstat()
 
 startTime = time.time()
 theta_prior = {'F0': {'type': 'unif',
-                      'lower': F0-DeltaF0/2.,
-                      'upper': F0+DeltaF0/2.},
+                      'lower': F0_center-DeltaF0,
+                      'upper': F0_center+DeltaF0},
                'F1': {'type': 'unif',
-                      'lower': F1-DeltaF1/2.,
-                      'upper': F1+DeltaF1/2.},
+                      'lower': F1_center-DeltaF1,
+                      'upper': F1_center+DeltaF1},
                'F2': F2,
                'Alpha': Alpha,
                'Delta': Delta
                }
 
-ntemps = 1
+ntemps = 2
 log10temperature_min = -1
 nwalkers = 100
 
@@ -77,11 +76,12 @@ mcmc = pyfstat.MCMCFollowUpSearch(
     log10temperature_min=log10temperature_min)
 mcmc.run(run_setup=run_setup, create_plots=False, log_table=False,
          gen_tex_table=False)
+mcmc.print_summary()
 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} {}\n'
-            .format(dF0, dF1, predicted_twoF, maxtwoF, runTime))
+    f.write('{:1.8e} {:1.8e} {:1.8e} {}\n'
+            .format(dF0, dF1, maxtwoF, runTime))
 os.system('rm {}/*{}*'.format(outdir, label))

Paper/DirectedMCNoiseOnly/repeat.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
-
-rm /local/user/gregory.ashton/MCResults*txt 
-for ((n=0;n<100;n++))
-do
-/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --quite --no-template-counting
-done
-cp /local/user/gregory.ashton/MCResults*txt /home/gregory.ashton/PyFstat/Paper/DirectedMCNoiseOnly/CollectedOutput

Paper/DirectedMCNoiseOnly/submitfile

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