plot_data.py 2.83 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
df_list = []
for fn in filenames:
    df = pd.read_csv(
Gregory Ashton's avatar
Gregory Ashton committed
34
        fn, sep=' ', names=['depth', 'h0', 'dF0', 'dF1', 'twoF', 'runTime'])
Gregory Ashton's avatar
Gregory Ashton committed
35
    df['CLUSTER_ID'] = fn.split('_')[1]
Gregory Ashton's avatar
Gregory Ashton committed
36
37
38
39
    if len(df) != 9:
        print len(df), fn
    else:
        df_list.append(df)
Gregory Ashton's avatar
Gregory Ashton committed
40
df = pd.concat(df_list)
41
42
43
44

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

Gregory Ashton's avatar
Gregory Ashton committed
55
yerr = np.abs(recovery_fraction - np.array(recovery_fraction_CI).T)
56
fig, ax = plt.subplots()
Gregory Ashton's avatar
Gregory Ashton committed
57
ax.errorbar(depths, recovery_fraction, yerr=yerr, fmt='sr', marker='s', ms=2,
Gregory Ashton's avatar
Gregory Ashton committed
58
            capsize=1, capthick=0.5, elinewidth=0.5,
Gregory Ashton's avatar
Gregory Ashton committed
59
            label='Monte-Carlo result', zorder=10)
Gregory Ashton's avatar
Gregory Ashton committed
60

61
fname = 'analytic_data_{}.txt'.format(twoFstar)
Gregory Ashton's avatar
Gregory Ashton committed
62
63
64
65
66
67
68
69
70
71
72
73
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')

74

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

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


Gregory Ashton's avatar
Gregory Ashton committed
84
85
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
86
fig, ax = plt.subplots()
Gregory Ashton's avatar
Gregory Ashton committed
87
88
89
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
90
91
fig.savefig('runTimeHist.png')