Remove old paper directory

Paper/AllSkyMC/AllSkyMC_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
-
-for ((n=0;n<10;n++))
-do
-/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_"$1".txt $(pwd)/CollectedOutput

Paper/AllSkyMC/generate_data.py

-import pyfstat
-import numpy as np
-import os
-import sys
-import time
-
-ID = sys.argv[1]
-outdir = sys.argv[2]
-
-label = 'run_{}'.format(ID)
-data_label = '{}_data'.format(label)
-results_file_name = '{}/MCResults_{}.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_center = 30
-F1_center = -1e-10
-F2 = 0
-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)
-
-DeltaAlpha = 0.02
-DeltaDelta = 0.02
-
-depths = np.linspace(100, 400, 9)
-depths = [118.75, 156.25]
-
-nsteps = 50
-run_setup = [((nsteps, 0), 20, False),
-             ((nsteps, 0), 11, False),
-             ((nsteps, 0), 6, False),
-             ((nsteps, 0), 3, False),
-             ((nsteps, nsteps), 1, False)]
-
-
-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(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)
-
-    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,
-        detector='H1,L1')
-    data.make_data()
-
-    startTime = time.time()
-    theta_prior = {'F0': {'type': 'unif',
-                          'lower': F0_center-DeltaF0,
-                          'upper': F0_center+DeltaF0},
-                   'F1': {'type': 'unif',
-                          'lower': F1_center-DeltaF1,
-                          'upper': F1_center+DeltaF1},
-                   'F2': F2,
-                   'Alpha': {'type': 'unif',
-                             'lower': Alpha_center-DeltaAlpha,
-                             'upper': Alpha_center+DeltaAlpha},
-                   'Delta': {'type': 'unif',
-                             'lower': Delta_center-DeltaDelta,
-                             'upper': Delta_center+DeltaDelta},
-                   }
-
-    ntemps = 2
-    log10temperature_min = -1
-    nwalkers = 100
-
-    mcmc = pyfstat.MCMCFollowUpSearch(
-        label=label, outdir=outdir,
-        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.run(run_setup=run_setup, create_plots=False, log_table=False,
-             gen_tex_table=False)
-    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(depth, h0, dF0, dF1, maxtwoF, runTime))
-    os.system('rm {}/*{}*'.format(outdir, label))

Paper/AllSkyMC/generate_failures.py

-import pyfstat
-import numpy as np
-import os
-import time
-
-outdir = 'data'
-
-label = 'run_failures'
-data_label = '{}_data'.format(label)
-results_file_name = '{}/MCResults_failures.txt'.format(outdir)
-
-# Properties of the GW data
-sqrtSX = 2e-23
-tstart = 1000000000
-Tspan = 100*86400
-tend = tstart + Tspan
-
-# Fixed properties of the signal
-F0_center = 30
-F1_center = 1e-10
-F2 = 0
-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)
-
-DeltaAlpha = 0.02
-DeltaDelta = 0.02
-
-depths = [140]
-
-nsteps = 50
-run_setup = [((nsteps, 0), 20, False),
-             ((nsteps, 0), 11, False),
-             ((nsteps, 0), 6, False),
-             ((nsteps, 0), 3, False),
-             ((nsteps, nsteps), 1, False)]
-
-
-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)
-    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)
-
-    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,
-        detector='H1,L1')
-    data.make_data()
-    predicted_twoF = data.predict_fstat()
-
-    startTime = time.time()
-    theta_prior = {'F0': {'type': 'unif',
-                          'lower': F0_center-DeltaF0,
-                          'upper': F0_center+DeltaF0},
-                   'F1': {'type': 'unif',
-                          'lower': F1_center-DeltaF1,
-                          'upper': F1_center+DeltaF1},
-                   'F2': F2,
-                   'Alpha': {'type': 'unif',
-                             'lower': Alpha_center-DeltaAlpha,
-                             'upper': Alpha_center+DeltaAlpha},
-                   'Delta': {'type': 'unif',
-                             'lower': Delta_center-DeltaDelta,
-                             'upper': Delta_center+DeltaDelta},
-                   }
-
-    ntemps = 2
-    log10temperature_min = -1
-    nwalkers = 100
-
-    mcmc = pyfstat.MCMCFollowUpSearch(
-        label=label, outdir=outdir,
-        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.run(run_setup=run_setup, create_plots=True, log_table=False,
-             gen_tex_table=False)
-    d, maxtwoF = mcmc.get_max_twoF()
-    print 'MaxtwoF = {}'.format(maxtwoF)

Paper/AllSkyMC/generate_table.py

-import pyfstat
-import numpy as np
-
-outdir = 'data'
-
-label = 'allsky_setup'
-data_label = '{}_data'.format(label)
-
-# Properties of the GW data
-sqrtSX = 2e-23
-tstart = 1000000000
-Tspan = 100*86400
-tend = tstart + Tspan
-
-# Fixed properties of the signal
-F0_center = 30
-F1_center = 1e-10
-F2 = 0
-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)
-
-DeltaAlpha = 0.02
-DeltaDelta = 0.02
-
-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)
-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_center = 0
-Delta_center = 0
-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)
-
-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,
-    detector='H1,L1')
-data.make_data()
-predicted_twoF = data.predict_fstat()
-
-theta_prior = {'F0': {'type': 'unif',
-                      'lower': F0_center-DeltaF0,
-                      'upper': F0_center+DeltaF0},
-               'F1': {'type': 'unif',
-                      'lower': F1_center-DeltaF1,
-                      'upper': F1_center+DeltaF1},
-               'F2': F2,
-               'Alpha': {'type': 'unif',
-                         'lower': Alpha_center-DeltaAlpha,
-                         'upper': Alpha_center+DeltaAlpha},
-               'Delta': {'type': 'unif',
-                         'lower': Delta_center-DeltaDelta,
-                         'upper': Delta_center+DeltaDelta},
-               }
-
-ntemps = 2
-log10temperature_min = -1
-nwalkers = 100
-
-mcmc = pyfstat.MCMCFollowUpSearch(
-    label=label, outdir=outdir,
-    sftfilepath='{}/*{}*sft'.format(outdir, data_label),
-    theta_prior=theta_prior,
-    tref=tref, minStartTime=tstart, maxStartTime=tend,
-    nwalkers=nwalkers, ntemps=ntemps, nsteps=[nsteps, nsteps],
-    log10temperature_min=log10temperature_min)
-mcmc.run(Nsegs0=20, R=10)
-#mcmc.run(run_setup)

Paper/AllSkyMC/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
-
-filenames = glob.glob("CollectedOutput/*.txt")
-
-plt.style.use('paper')
-
-Tspan = 100 * 86400
-
-
-def Recovery(Tspan, Depth, twoFstar=60, detectors='H1,L1'):
-    numDetectors = len(detectors.split(','))
-    cmd = ("DetectionProbabilityStackSlide('Nseg', 1, 'Tdata', {},"
-           "'misHist', createDeltaHist(0), 'avg2Fth', {}, 'detectors', '{}',"
-           "'Depth', {})"
-           ).format(numDetectors*Tspan, twoFstar, detectors, Depth)
-    return octave.eval(cmd, verbose=False)
-
-
-def binomialConfidenceInterval(N, K, confidence=0.95):
-    cmd = '[fLow, fUpper] = binomialConfidenceInterval({}, {}, {})'.format(
-        N, K, confidence)
-    [l, u] =  octave.eval(cmd, verbose=False, return_both=True)[0].split('\n')
-    return float(l.split('=')[1]), float(u.split('=')[1])
-
-df_list = []
-for fn in filenames:
-    df = pd.read_csv(
-        fn, sep=' ', names=['depth', 'h0', 'dF0', 'dF1', 'twoF', 'runTime'])
-    df['CLUSTER_ID'] = fn.split('_')[1]
-    df_list.append(df)
-df = pd.concat(df_list)
-
-twoFstar = 70
-depths = np.unique(df.depth.values)
-recovery_fraction = []
-recovery_fraction_CI = []
-for d in depths:
-    twoFs = df[df.depth == d].twoF.values
-    N = len(twoFs)
-    K = np.sum(twoFs > twoFstar)
-    print d, N, K
-    recovery_fraction.append(K/float(N))
-    [fLower, fUpper] = binomialConfidenceInterval(N, K)
-    recovery_fraction_CI.append([fLower, fUpper])
-
-yerr = np.abs(recovery_fraction - np.array(recovery_fraction_CI).T)
-fig, ax = plt.subplots()
-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', zorder=10)
-
-fname = 'analytic_data_{}.txt'.format(twoFstar)
-if os.path.isfile(fname):
-    depths_smooth, recovery_analytic = np.loadtxt(fname)
-else:
-    depths_smooth = np.linspace(10, 550, 100)
-    recovery_analytic = []
-    for d in tqdm(depths_smooth):
-        recovery_analytic.append(Recovery(Tspan, d, twoFstar))
-    np.savetxt(fname, np.array([depths_smooth, recovery_analytic]))
-depths_smooth = np.concatenate(([0], depths_smooth))
-recovery_analytic = np.concatenate(([1], recovery_analytic))
-ax.plot(depths_smooth, recovery_analytic, '-k', label='Theoretical maximum')
-
-
-ax.set_ylim(0, 1.05)
-ax.set_xlabel(r'Sensitivity depth', size=10)
-ax.set_ylabel(r'Recovered fraction', size=10)
-ax.legend(loc=1, frameon=False)
-
-fig.tight_layout()
-fig.savefig('allsky_recovery.png')
-
-
-total_number_steps = 6 * 50.
-fig, ax = plt.subplots()
-ax.hist(df.runTime/total_number_steps, bins=50)
-ax.set_xlabel('run-time per step [s]')
-fig.tight_layout()
-fig.savefig('runTimeHist.png')
-

Paper/AllSkyMC/submitfile

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

Paper/AllSkyMCNoiseOnly/AllSkyMCNoiseOnly_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
-
-for ((n=0;n<90;n++))
-do
-/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/NoiseOnlyMCResults_"$1".txt $(pwd)/CollectedOutput

Paper/AllSkyMCNoiseOnly/generate_data.py

-import pyfstat
-import numpy as np
-import os
-import sys
-import time
-
-ID = sys.argv[1]
-outdir = sys.argv[2]
-
-label = 'run_{}'.format(ID)
-data_label = '{}_data'.format(label)
-results_file_name = '{}/NoiseOnlyMCResults_{}.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_center = 30
-F1_center = -1e-10
-F2 = 0
-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)
-
-DeltaAlpha = 0.02
-DeltaDelta = 0.02
-
-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 = 0
-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(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)
-
-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,
-    detector='H1,L1')
-data.make_data()
-
-startTime = time.time()
-theta_prior = {'F0': {'type': 'unif',
-                      'lower': F0_center-DeltaF0,
-                      'upper': F0_center+DeltaF0},
-               'F1': {'type': 'unif',
-                      'lower': F1_center-DeltaF1,
-                      'upper': F1_center+DeltaF1},
-               'F2': F2,
-               'Alpha': {'type': 'unif',
-                         'lower': Alpha_center-DeltaAlpha,
-                         'upper': Alpha_center+DeltaAlpha},
-               'Delta': {'type': 'unif',
-                         'lower': Delta_center-DeltaDelta,
-                         'upper': Delta_center+DeltaDelta},
-               }
-
-ntemps = 2
-log10temperature_min = -1
-nwalkers = 100
-
-mcmc = pyfstat.MCMCFollowUpSearch(
-    label=label, outdir=outdir,
-    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.run(run_setup=run_setup, create_plots=False, log_table=False,
-         gen_tex_table=False)
-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} {}\n'
-            .format(dF0, dF1, maxtwoF, runTime))
-os.system('rm {}/*{}*'.format(outdir, label))

Paper/AllSkyMCNoiseOnly/generate_table.py

-import pyfstat
-import numpy as np
-
-outdir = 'data'
-
-label = 'allsky_setup'
-data_label = '{}_data'.format(label)
-
-# Properties of the GW data
-sqrtSX = 2e-23
-tstart = 1000000000
-Tspan = 100*86400
-tend = tstart + Tspan
-
-# Fixed properties of the signal
-F0_center = 30
-F1_center = 1e-10
-F2 = 0
-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)
-
-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)
-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
-
-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,
-    detector='H1,L1')
-data.make_data()
-predicted_twoF = data.predict_fstat()
-
-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/2.,
-                         'upper': Alpha+DeltaAlpha/2.},
-               'Delta': {'type': 'unif',
-                         'lower': Delta-DeltaDelta/2.,
-                         'upper': Delta+DeltaDelta/2.},
-               }
-
-ntemps = 1
-log10temperature_min = -1
-nwalkers = 100
-
-mcmc = pyfstat.MCMCFollowUpSearch(
-    label=label, outdir=outdir,
-    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.run(run_setup)

Paper/AllSkyMCNoiseOnly/plot_data.py

-import pyfstat
-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
-from scipy.stats import rv_continuous, chi2
-
-filenames = glob.glob("CollectedOutput/*.txt")
-
-plt.style.use('paper')
-
-Tspan = 100 * 86400
-
-
-class maxtwoFinNoise_gen(rv_continuous):
-    def _pdf(self, twoF, Ntrials):
-        F = twoF/2.0
-        alpha = (1 + F)*np.exp(-F)
-        a = Ntrials/2.0*F*np.exp(-F)
-        b = (1 - alpha)**(Ntrials-1)
-        return a*b
-
-
-df_list = []
-for fn in filenames:
-    df = pd.read_csv(
-        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,
-        label='Monte-Carlo histogram')
-
-maxtwoFinNoise = maxtwoFinNoise_gen(a=0)
-Ntrials_effective, loc, scale = maxtwoFinNoise.fit(df.twoF.values, floc=0, fscale=1)
-print 'Ntrials effective = {:1.2e}'.format(Ntrials_effective)
-twoFsmooth = np.linspace(0, df.twoF.max(), 1000)
-best_fit_pdf = maxtwoFinNoise.pdf(twoFsmooth, Ntrials_effective)
-ax.plot(twoFsmooth, best_fit_pdf, '-r',
-        label=r'$p(2\mathcal{{F}}_{{\rm max}})$ for {} $N_{{\rm trials}}$'
-        .format(pyfstat.texify_float(Ntrials_effective, d=2)))
-
-pval = 1e-6
-twoFsmooth_HD = np.linspace(
-    twoFsmooth[np.argmax(best_fit_pdf)], df.twoF.max(), 100000)
-best_fit_pdf_HD = maxtwoFinNoise.pdf(twoFsmooth_HD, Ntrials_effective)
-spacing = twoFsmooth_HD[1]-twoFsmooth_HD[0]
-print twoFsmooth_HD[np.argmin(np.abs(best_fit_pdf_HD - pval))], spacing
-
-ax.set_xlabel('$\widetilde{2\mathcal{F}}$')
-ax.set_xlim(0, 60)
-ax.legend(frameon=False, fontsize=6, loc=2)
-fig.tight_layout()
-fig.savefig('allsky_noise_twoF_histogram.png')
-
-from latex_macro_generator import write_to_macro
-write_to_macro('AllSkyMCNoiseOnlyMaximum', '{:1.1f}'.format(np.max(df.twoF)),
-               '../macros.tex')
-write_to_macro('AllSkyMCNoiseN', len(df), '../macros.tex')

Paper/AllSkyMCNoiseOnly/submitfile

-Executable= AllSkyMCNoiseOnly_repeat.sh
-Arguments= $(Cluster)_$(Process)
-Universe=vanilla
-Input=/dev/null
-accounting_group = ligo.dev.o2.cw.explore.test
-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/DirectedMC/DirectedMC_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
-
-for ((n=0;n<10;n++))
-do
-/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_"$1".txt $(pwd)/CollectedOutput

Paper/DirectedMC/generate_data.py

-import pyfstat
-import numpy as np
-import os
-import sys
-import time
-
-ID = sys.argv[1]
-outdir = sys.argv[2]
-
-label = 'run_{}'.format(ID)
-data_label = '{}_data'.format(label)
-results_file_name = '{}/MCResults_{}.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_center = 30
-F1_center = -1e-10
-F2 = 0
-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)
-
-depths = np.linspace(100, 400, 9)
-depths = [118.75, 156.25]
-
-nsteps = 25
-run_setup = [((nsteps, 0), 20, False),
-             ((nsteps, 0), 7, False),
-             ((nsteps, 0), 2, False),
-             ((nsteps, nsteps), 1, False)]
-
-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
-
-    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,
-        detector='H1,L1')
-    data.make_data()
-
-    startTime = time.time()
-    theta_prior = {'F0': {'type': 'unif',
-                          'lower': F0_center-DeltaF0,
-                          'upper': F0_center+DeltaF0},
-                   'F1': {'type': 'unif',
-                          'lower': F1_center-DeltaF1,
-                          'upper': F1_center+DeltaF1},
-                   'F2': F2,
-                   'Alpha': Alpha,
-                   'Delta': Delta
-                   }
-
-    ntemps = 2
-    log10temperature_min = -1
-    nwalkers = 100
-
-    mcmc = pyfstat.MCMCFollowUpSearch(
-        label=label, outdir=outdir,
-        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.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(depth, h0, dF0, dF1, maxtwoF, runTime))
-    os.system('rm {}/*{}*'.format(outdir, label))

Paper/DirectedMC/generate_table.py

-import pyfstat
-import numpy as np
-
-outdir = 'data'
-
-label = 'directed_setup'
-data_label = '{}_data'.format(label)
-
-# Properties of the GW data
-sqrtSX = 1e-23
-tstart = 1000000000
-Tspan = 100*86400
-tend = tstart + Tspan
-
-# Fixed properties of the signal
-F0_center = 30
-F1_center = 1e-10
-F2 = 0
-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)
-
-depth = 100
-
-nsteps = 25
-run_setup = [((nsteps, 0), 20, False),
-             ((nsteps, 0), 7, False),
-             ((nsteps, 0), 2, False),
-             ((nsteps, nsteps), 1, False)]
-
-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)
-
-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,
-    detector='H1,L1')
-data.make_data()
-predicted_twoF = data.predict_fstat()
-
-theta_prior = {'F0': {'type': 'unif',
-                      'lower': F0-DeltaF0,
-                      'upper': F0+DeltaF0},
-               'F1': {'type': 'unif',
-                      'lower': F1-DeltaF1,
-                      'upper': F1+DeltaF1},
-               'F2': F2,
-               'Alpha': Alpha,
-               'Delta': Delta
-               }
-
-ntemps = 1
-log10temperature_min = -1
-nwalkers = 100
-
-mcmc = pyfstat.MCMCFollowUpSearch(
-    label=label, outdir=outdir,
-    sftfilepath='{}/*{}*sft'.format(outdir, data_label),
-    theta_prior=theta_prior,
-    tref=tref, minStartTime=tstart, maxStartTime=tend,
-    nwalkers=nwalkers, ntemps=ntemps, nsteps=[nsteps, nsteps],
-    log10temperature_min=log10temperature_min)
-#mcmc.run(Nsegs0=20, R=10)
-mcmc.run(run_setup)

Paper/DirectedMC/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
-
-filenames = glob.glob("CollectedOutput/*.txt")
-
-plt.style.use('paper')
-
-Tspan = 100 * 86400
-
-
-def Recovery(Tspan, Depth, twoFstar=60, detectors='H1,L1'):
-    numDetectors = len(detectors.split(','))
-    cmd = ("DetectionProbabilityStackSlide('Nseg', 1, 'Tdata', {},"
-           "'misHist', createDeltaHist(0), 'avg2Fth', {}, 'detectors', '{}',"
-           "'Depth', {})"
-           ).format(numDetectors*Tspan, twoFstar, detectors, Depth)
-    return octave.eval(cmd, verbose=False)
-
-
-def binomialConfidenceInterval(N, K, confidence=0.95):
-    cmd = '[fLow, fUpper] = binomialConfidenceInterval({}, {}, {})'.format(
-        N, K, confidence)
-    [l, u] = octave.eval(cmd, verbose=False, return_both=True)[0].split('\n')
-    return float(l.split('=')[1]), float(u.split('=')[1])
-
-df_list = []
-for fn in filenames:
-    df = pd.read_csv(
-        fn, sep=' ', names=['depth', 'h0', 'dF0', 'dF1', '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)
-recovery_fraction = []
-recovery_fraction_CI = []
-for d in depths:
-    twoFs = df[df.depth == d].twoF.values
-    N = len(twoFs)
-    K = np.sum(twoFs > twoFstar)
-    print d, N, K
-    recovery_fraction.append(K/float(N))
-    [fLower, fUpper] = binomialConfidenceInterval(N, K)
-    recovery_fraction_CI.append([fLower, fUpper])
-
-yerr = np.abs(recovery_fraction - np.array(recovery_fraction_CI).T)
-fig, ax = plt.subplots()
-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', zorder=10)
-
-fname = 'analytic_data_{}.txt'.format(twoFstar)
-if os.path.isfile(fname):
-    depths_smooth, recovery_analytic = np.loadtxt(fname)
-else:
-    depths_smooth = np.linspace(10, 550, 100)
-    recovery_analytic = []
-    for d in tqdm(depths_smooth):
-        recovery_analytic.append(Recovery(Tspan, d, twoFstar))
-    np.savetxt(fname, np.array([depths_smooth, recovery_analytic]))
-depths_smooth = np.concatenate(([0], depths_smooth))
-recovery_analytic = np.concatenate(([1], recovery_analytic))
-ax.plot(depths_smooth, recovery_analytic, '-k', label='Theoretical maximum')
-
-
-ax.set_ylim(0, 1.05)
-ax.set_xlabel(r'Sensitivity depth', size=10)
-ax.set_ylabel(r'Recovered fraction', size=10)
-ax.legend(loc=1, frameon=False)
-
-fig.tight_layout()
-fig.savefig('directed_recovery.png')
-
-
-total_number_steps = 5*25.
-fig, ax = plt.subplots()
-ax.hist(df.runTime/total_number_steps, bins=50)
-ax.set_xlabel('run-time per step [s]')
-fig.tight_layout()
-fig.savefig('runTimeHist.png')
-

Paper/DirectedMC/plot_recovery.py

-import matplotlib.pyplot as plt
-import numpy as np
-import scipy.stats
-
-
-def Recovery(Tspan, Depth, twoFstar=60):
-    rho2 = 4*Tspan/25./Depth**2
-    twoF_Hs = scipy.stats.distributions.ncx2(df=4, nc=rho2)
-    return 1 - twoF_Hs.cdf(twoFstar)
-
-N = 500
-Tspan = np.linspace(0.1, 365*86400, N)
-Depth = np.linspace(10, 300, N)
-
-X, Y = np.meshgrid(Tspan, Depth)
-X = X / 86400
-
-Z = [[Recovery(t, d) for t in Tspan] for d in Depth]
-
-fig, ax = plt.subplots()
-pax = ax.pcolormesh(X, Y, Z, cmap=plt.cm.viridis)
-CS = ax.contour(X, Y, Z, [0.95])
-plt.clabel(CS, inline=1, fontsize=12, fmt='%s', manual=[(200, 180)])
-plt.colorbar(pax, label='Recovery fraction')
-ax.set_xlabel(r'$T_{\rm span}$ [days]', size=16)
-ax.set_ylabel(r'Depth=$\frac{\sqrt{S_{\rm n}}}{h_0}$', size=14)
-ax.set_xlim(min(Tspan)/86400., max(Tspan)/86400.)
-ax.set_ylim(min(Depth), max(Depth))
-
-fig.savefig('recovery.png')

Paper/DirectedMC/submitfile

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

Paper/DirectedMCNoiseOnly/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
-
-for ((n=0;n<100;n++))
-do
-/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/NoiseOnlyMCResults_"$1".txt $(pwd)/CollectedOutput