plot_data.py 2.52 KB
Newer Older
1
2
3
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
Gregory Ashton's avatar
Gregory Ashton committed
4
5
6
7
8
import os
from tqdm import tqdm
from oct2py import octave

plt.style.use('paper')
9
10
11
12

Tspan = 100 * 86400


Gregory Ashton's avatar
Gregory Ashton committed
13
14
15
16
17
18
19
20
21
22
23
24
25
26
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])
27
28
29
30
31

results_file_name = 'MCResults.txt'

df = pd.read_csv(
    results_file_name, sep=' ', names=['depth', 'h0', 'dF0', 'dF1',
32
                                       'twoF_predicted', 'twoF', 'runTime'])
33
34
35
36

twoFstar = 60
depths = np.unique(df.depth.values)
recovery_fraction = []
Gregory Ashton's avatar
Gregory Ashton committed
37
recovery_fraction_CI = []
38
39
for d in depths:
    twoFs = df[df.depth == d].twoF.values
Gregory Ashton's avatar
Gregory Ashton committed
40
41
42
43
44
45
    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])
46

Gregory Ashton's avatar
Gregory Ashton committed
47
yerr = np.abs(recovery_fraction - np.array(recovery_fraction_CI).T)
48
fig, ax = plt.subplots()
Gregory Ashton's avatar
Gregory Ashton committed
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
ax.errorbar(depths, recovery_fraction, yerr=yerr, fmt='sk', marker='s', ms=2,
            capsize=1, capthick=0.5, elinewidth=0.5,
            label='Monte-Carlo result')

fname = 'analytic_data.txt'
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')

66

Gregory Ashton's avatar
Gregory Ashton committed
67
68
69
70
ax.set_ylim(0, 1.05)
ax.set_xlabel(r'Signal depth', size=10)
ax.set_ylabel(r'Recovered fraction', size=10)
ax.legend(loc=1, frameon=False)
71

Gregory Ashton's avatar
Gregory Ashton committed
72
fig.tight_layout()
73
fig.savefig('directed_recovery.png')
Gregory Ashton's avatar
Gregory Ashton committed
74
75
76
77
78
79
80


fig, ax = plt.subplots()
ax.hist(df.runTime, bins=50)
ax.set_xlabel('runTime per follow-up [s]')
fig.savefig('runTimeHist.png')