plot_data.py 2.78 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
import os
from tqdm import tqdm
from oct2py import octave
Gregory Ashton's avatar
Gregory Ashton committed
7
8
9
import glob

filenames = glob.glob("CollectedOutput/*.txt")
Gregory Ashton's avatar
Gregory Ashton committed
10
11

plt.style.use('paper')
12
13
14
15

Tspan = 100 * 86400


Gregory Ashton's avatar
Gregory Ashton committed
16
17
18
19
20
21
22
23
24
25
26
27
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)
Gregory Ashton's avatar
Gregory Ashton committed
28
    [l, u] = octave.eval(cmd, verbose=False, return_both=True)[0].split('\n')
Gregory Ashton's avatar
Gregory Ashton committed
29
    return float(l.split('=')[1]), float(u.split('=')[1])
30

Gregory Ashton's avatar
Gregory Ashton committed
31
32
33
34
35
36
37
38
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)
39
40
41
42

twoFstar = 60
depths = np.unique(df.depth.values)
recovery_fraction = []
Gregory Ashton's avatar
Gregory Ashton committed
43
recovery_fraction_CI = []
44
45
for d in depths:
    twoFs = df[df.depth == d].twoF.values
Gregory Ashton's avatar
Gregory Ashton committed
46
47
48
49
50
51
    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])
52

Gregory Ashton's avatar
Gregory Ashton committed
53
yerr = np.abs(recovery_fraction - np.array(recovery_fraction_CI).T)
54
fig, ax = plt.subplots()
Gregory Ashton's avatar
Gregory Ashton committed
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
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')

72

Gregory Ashton's avatar
Gregory Ashton committed
73
74
75
76
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)
77

Gregory Ashton's avatar
Gregory Ashton committed
78
fig.tight_layout()
79
fig.savefig('directed_recovery.png')
Gregory Ashton's avatar
Gregory Ashton committed
80
81


Gregory Ashton's avatar
Gregory Ashton committed
82
83
total_number_steps = 5*20.
df_clean = df[df.CLUSTER_ID == '969049']  # Hack due to a change in the code
Gregory Ashton's avatar
Gregory Ashton committed
84
fig, ax = plt.subplots()
Gregory Ashton's avatar
Gregory Ashton committed
85
86
87
ax.hist(df_clean.runTime/total_number_steps, bins=50)
ax.set_xlabel('run-time per step [s]')
fig.tight_layout()
Gregory Ashton's avatar
Gregory Ashton committed
88
89
fig.savefig('runTimeHist.png')