Index: trunk/bin/graph
===================================================================
--- trunk/bin/graph	(revision 6)
+++ trunk/bin/graph	(revision 6)
@@ -0,0 +1,341 @@
+#!/usr/bin/env python3
+
+import sys
+import os
+import time
+import random
+import tempfile
+import argparse
+import socket
+import json
+
+import matplotlib.mlab as mlab
+import matplotlib.pyplot as plt
+
+
+VERSION = "{DEVELOPMENT}"
+if VERSION == "{DEVELOPMENT}":
+    script_dir = '.'
+    try:
+        script_dir = os.path.dirname(os.path.realpath(__file__))
+    except:
+        try:
+            script_dir = os.path.dirname(os.path.abspath(sys.argv[0]))
+        except:
+            pass
+    sys.path.append("%s/../lib" % script_dir)
+
+from nanownlib import *
+from nanownlib.stats import *
+import nanownlib.storage
+
+
+parser = argparse.ArgumentParser(
+    description="")
+parser.add_argument('db_file', default=None,
+                    help='')
+options = parser.parse_args()
+db = nanownlib.storage.db(options.db_file)
+
+
+def differences(db, unusual_case, column='packet_rtt'):
+    cursor = db.conn.cursor()
+    query="""
+      SELECT %(column)s-(SELECT avg(%(column)s) FROM probes,analysis
+                         WHERE analysis.probe_id=probes.id AND probes.test_case!=:unusual_case AND probes.type in ('train','test') AND sample=u.sample)
+      FROM (SELECT probes.sample,%(column)s FROM probes,analysis 
+                         WHERE analysis.probe_id=probes.id AND probes.test_case =:unusual_case AND probes.type in ('train','test')) u
+      """ % {"column":column}
+    params = {"unusual_case":unusual_case}
+    cursor.execute(query, params)
+    for row in cursor:
+        yield row[0]
+
+
+def timeSeries(db, probe_type, unusual_case):
+    cursor = db.conn.cursor()
+    query="""
+      SELECT time_of_day,packet_rtt AS uc,(SELECT avg(packet_rtt) FROM probes,analysis
+                                           WHERE analysis.probe_id=probes.id AND probes.test_case!=:unusual_case AND probes.type=:probe_type AND sample=u.sample) AS oc
+      FROM (SELECT time_of_day,probes.sample,packet_rtt FROM probes,analysis 
+                                           WHERE analysis.probe_id=probes.id AND probes.test_case =:unusual_case AND probes.type=:probe_type) u
+    """
+    
+    params = {"probe_type":probe_type,"unusual_case":unusual_case}
+    cursor.execute(query, params)
+    for row in cursor:
+        yield {'time_of_day':row['time_of_day'],unusual_case:row['uc'],'other_cases':row['oc']}
+#samples,derived,null_derived = parse_data(input1)
+
+#trust = trustValues(derived, sum)
+#weights = linearWeights(derived, trust, 0.25)
+#print('(test): %f' % weightedMean(derived,weights))
+
+diffs = list(differences(db, 'long'))
+reported_diffs = list(differences(db, 'long', 'reported'))
+#shorts = [s['packet_rtt'] for s in samples.values() if s['test_case']=='short']
+#longs = [s['packet_rtt'] for s in samples.values() if s['test_case']=='long']
+
+short_overtime = [(sample['time_of_day'],sample['short']) for sample in timeSeries(db,'train','short')]
+long_overtime = [(sample['time_of_day'],sample['long']) for sample in timeSeries(db,'train','long')]
+diff_overtime = [(sample['time_of_day'],sample['long']-sample['other_cases']) for sample in timeSeries(db,'train','long')]
+short_overtime.sort()
+long_overtime.sort()
+diff_overtime.sort()
+
+print('packet_rtt diff median: %f' % statistics.median(diffs))
+print('packet_rtt diff midhinge: %f' % midhinge(diffs))
+print('packet_rtt diff trimean: %f' % trimean(diffs))
+print('packet_rtt diff MAD: %f' % mad(diffs))
+print('reported diff trimean: %f' % trimean(reported_diffs))
+print('reported diff MAD: %f' % mad(reported_diffs))
+
+
+#all_data = longs+shorts
+#all_data.sort()
+#cut_off_low = all_data[0]
+#cut_off_high = all_data[int(len(all_data)*0.997)]
+
+
+plt.clf()
+plt.title("Packet RTT over time")
+plt.xlabel('Time of Day')
+plt.ylabel('RTT')
+s = plt.scatter([t for t,rtt in short_overtime], [rtt for t,rtt in short_overtime], s=1, color='red', alpha=0.6)
+l = plt.scatter([t for t,rtt in long_overtime], [rtt for t,rtt in long_overtime], s=1, color='blue', alpha=0.6)
+d = plt.scatter([t for t,rtt in diff_overtime], [rtt for t,rtt in diff_overtime], s=1, color='purple', alpha=0.6)
+plt.legend((s,l,d), ('short','long','difference'), scatterpoints=1)
+#plt.savefig('paper/figures/comcast-powerboost1.png')
+plt.show()
+
+short_overtime,long_overtime,diff_overtime = None,None,None
+
+
+num_bins = 300
+reported_diffs.sort()
+cut_off_low = reported_diffs[int(len(diffs)*0.003)]
+cut_off_high = reported_diffs[int(len(diffs)*0.997)]
+
+plt.clf()
+# the histogram of the data
+n, bins, patches = plt.hist(reported_diffs, num_bins, normed=1, color='black', histtype='step', alpha=0.8,
+                            range=(cut_off_low,cut_off_high))
+plt.xlabel('RTT Difference')
+plt.ylabel('Probability')
+plt.title(r'Histogram - distribution of differences')
+
+# Tweak spacing to prevent clipping of ylabel
+plt.subplots_adjust(left=0.15)
+#plt.legend()
+plt.show()
+#plt.savefig('paper/graphs/dists-vs-dist-of-diffs2.svg')
+
+
+
+
+num_bins = 300
+diffs.sort()
+cut_off_low = diffs[int(len(diffs)*0.003)]
+cut_off_high = diffs[int(len(diffs)*0.997)]
+
+plt.clf()
+# the histogram of the data
+n, bins, patches = plt.hist(diffs, num_bins, normed=1, color='purple', histtype='step', alpha=0.8,
+                            range=(cut_off_low,cut_off_high))
+plt.xlabel('RTT Difference')
+plt.ylabel('Probability')
+plt.title(r'Histogram - distribution of differences')
+
+# Tweak spacing to prevent clipping of ylabel
+plt.subplots_adjust(left=0.15)
+#plt.legend()
+plt.show()
+#plt.savefig('paper/graphs/dists-vs-dist-of-diffs2.svg')
+
+sys.exit(0)
+
+
+
+num_bins = 150
+# the histogram of the data
+n, bins, patches = plt.hist((shorts,longs), num_bins, normed=1, label=['short', 'long'], color=['red','blue'], histtype='step', alpha=0.8,
+                            range=(cut_off_low,cut_off_high))
+#n, bins, patches = plt.hist(shorts2+longs2, num_bins, normed=1, facecolor='blue', histtype='step', alpha=0.3)
+# add a 'best fit' line
+#y = mlab.normpdf(bins, mu, sigma)
+#plt.plot(bins, y, 'r--')
+plt.xlabel('packet_rtt')
+plt.ylabel('Probability')
+plt.title(r'Histogram - RTT short and long')
+
+# Tweak spacing to prevent clipping of ylabel
+plt.subplots_adjust(left=0.15)
+plt.legend()
+#plt.show()
+plt.savefig('paper/figures/comcast-powerboost2.svg')
+
+
+
+
+num_trials = 200
+
+
+subsample_sizes = (50,150,300,500,700,1000,2000,3000,5000,7000,10000,15000,20000)
+estimator = functools.partial(boxTest, 0.07, 0.08)
+performance = []
+for subsample_size in subsample_sizes:
+    estimates = bootstrap(derived, subsample_size, num_trials, estimator)
+    performance.append(100.0*len([e for e in estimates if e == 1])/num_trials)
+
+null_performance = []
+for subsample_size in subsample_sizes:
+    null_estimates = bootstrap(null_derived, subsample_size, num_trials, estimator)
+    null_performance.append(100.0*len([e for e in null_estimates if e == 0])/num_trials)
+
+plt.clf()
+plt.title("boxTest bootstrap")
+plt.xlabel('sample size')
+plt.ylabel('performance')
+plt.scatter(subsample_sizes, performance, s=2, color='red', alpha=0.6)
+plt.scatter(subsample_sizes, null_performance, s=2, color='blue', alpha=0.6)
+plt.show()
+
+
+
+subsample_sizes = (50,150,300,400,500,700,1000,2000,3000,4000,5000,7000,10000)
+estimator = diffMedian
+performance = []
+for subsample_size in subsample_sizes:
+    estimates = bootstrap(derived, subsample_size, num_trials, estimator)
+    performance.append(100.0*len([e for e in estimates if e > expected_mean*0.9 and e < expected_mean*1.1])/num_trials)
+
+plt.clf()
+plt.title("diff median bootstrap")
+plt.xlabel('sample size')
+plt.ylabel('performance')
+plt.scatter(subsample_sizes, performance, s=1, color='red', alpha=0.6)
+plt.show()
+
+
+
+
+subsample_sizes = (50,150,300,400,500,700,1000,2000,3000,4000,5000,7000,10000)
+weight_funcs = (linearWeights, prunedWeights)
+for wf in weight_funcs:
+    estimator = functools.partial(estimateMean, hypotenuse, wf, 0.40)
+    performance = []
+    for subsample_size in subsample_sizes:
+        estimates = bootstrap(derived, subsample_size, num_trials, estimator)
+        performance.append(100.0*len([e for e in estimates if e > expected_mean*0.9 and e < expected_mean*1.1])/num_trials)
+
+    plt.clf()
+    plt.title(repr(wf))
+    plt.xlabel('sample size')
+    plt.ylabel('performance')
+    plt.scatter(subsample_sizes, performance, s=1, color='red', alpha=0.6)
+    plt.show()
+
+
+
+num_bins = 300
+# the histogram of the data
+n, bins, patches = plt.hist((tsshorts,tslongs), num_bins, normed=1, label=['short', 'long'], color=['red','blue'], histtype='step', alpha=0.8)
+#n, bins, patches = plt.hist(shorts2+longs2, num_bins, normed=1, facecolor='blue', histtype='step', alpha=0.3)
+# add a 'best fit' line
+#y = mlab.normpdf(bins, mu, sigma)
+#plt.plot(bins, y, 'r--')
+plt.xlabel('packet_rtt')
+plt.ylabel('Probability')
+plt.title(r'Histogram - tsval_rtt short vs long')
+
+# Tweak spacing to prevent clipping of ylabel
+plt.subplots_adjust(left=0.15)
+plt.legend()
+plt.show()
+
+
+
+    
+####
+#trust_methods = [min,max,sum,difference,product]
+trust_methods = [sum,product,hypotenuse]
+colors = ['red','blue','green','purple','orange','black']
+weight_methods = [prunedWeights, linearWeights]
+alphas = [i/100.0 for i in range(0,100,2)]
+
+
+
+
+plt.clf()
+plt.title(r'Trust Method Comparison - Linear')
+plt.xlabel('Alpha')
+plt.ylabel('Mean error')
+paths = []
+for tm in trust_methods:
+    trust = trustValues(derived, tm)
+    series = []
+    for alpha in alphas:
+        weights = linearWeights(derived, trust, alpha)
+        series.append(weightedMean(derived, weights) - expected_mean)
+
+    paths.append(plt.scatter(alphas, series, s=1, color=colors[len(paths)],alpha=0.6))
+
+plt.legend(paths, [repr(tm) for tm in trust_methods], scatterpoints=1)
+plt.show()
+
+
+
+plt.clf()
+plt.title(r'Trust Method Comparison - Pruned')
+plt.xlabel('Alpha')
+plt.ylabel('Mean error')
+paths = []
+for tm in trust_methods:
+    trust = trustValues(derived, tm)
+    series = []
+    for alpha in alphas:
+        weights = prunedWeights(derived, trust, alpha)
+        series.append(weightedMean(derived, weights) - expected_mean)
+
+    paths.append(plt.scatter(alphas, series, s=1, color=colors[len(paths)],alpha=0.6))
+
+plt.legend(paths, [repr(tm) for tm in trust_methods], scatterpoints=1)
+plt.show()
+
+
+sys.exit(0)
+
+plt.clf()
+plt.title(r'Trust Method Comparison - Inverted')
+plt.xlabel('Alpha')
+plt.ylabel('Mean error')
+paths = []
+for tm in trust_methods:
+    trust = trustValues(derived, tm)
+    series = []
+    for alpha in alphas:
+        weights = invertedWeights(derived, trust, alpha)
+        series.append(weightedMean(derived, weights) - expected_mean)
+
+    paths.append(plt.scatter(alphas, series, s=1, color=colors[len(paths)],alpha=0.6))
+
+plt.legend(paths, [repr(tm) for tm in trust_methods], scatterpoints=1)
+plt.show()
+
+
+plt.clf()
+plt.title(r'Trust Method Comparison - Arctangent')
+plt.xlabel('Alpha')
+plt.ylabel('Mean error')
+paths = []
+for tm in trust_methods:
+    trust = trustValues(derived, tm)
+    series = []
+    for alpha in alphas:
+        weights = arctanWeights(derived, trust, alpha)
+        series.append(weightedMean(derived, weights) - expected_mean)
+
+    paths.append(plt.scatter(alphas, series, s=1, color=colors[len(paths)],alpha=0.6))
+
+plt.legend(paths, [repr(tm) for tm in trust_methods], scatterpoints=1)
+plt.show()
Index: trunk/bin/train
===================================================================
--- trunk/bin/train	(revision 5)
+++ trunk/bin/train	(revision 6)
@@ -290,10 +290,10 @@
 
 
-def trainMidhinge(db, unusual_case, greater, subsample_size):
+def trainMidhinge(db, unusual_case, greater, subseries_size):
 
     def trainAux(distance, threshold, num_trials):
-        estimator = functools.partial(midhingeTest, {'distance':distance,'threshold':threshold}, unusual_case, greater)
-        estimates = bootstrap2(estimator, db, 'train', subsample_size, num_trials)
-        null_estimates = bootstrap2(estimator, db, 'train_null', subsample_size, num_trials)
+        estimator = functools.partial(midhingeTest, {'distance':distance,'threshold':threshold}, greater)
+        estimates = bootstrap3(estimator, db, 'train', unusual_case, subseries_size, num_trials)
+        null_estimates = bootstrap3(estimator, db, 'train_null', unusual_case, subseries_size, num_trials)
 
         bad_estimates = len([e for e in estimates if e != 1])
@@ -306,5 +306,5 @@
     #determine expected delta based on differences
     start = time.time()
-    mean_diffs = list(samples2MeanDiffs(subsample(db, 'train'), 'packet_rtt', unusual_case))
+    mean_diffs = [s['unusual_case']-s['other_cases'] for s in subseries(db, 'train', unusual_case)]
     threshold = trimean(mean_diffs)/2.0
     print("initial threshold:", threshold)
@@ -313,5 +313,5 @@
     print("trimean threshold:", trimean(mean_diffs)/2.0)
     
-    mean_diffs = list(samples2MeanDiffs(subsample(db, 'train_null'), 'packet_rtt', unusual_case))
+    mean_diffs = [s['unusual_case']-s['other_cases'] for s in subseries(db, 'train_null', unusual_case)]
     print(len(mean_diffs))
     print("null mean:", statistics.mean(mean_diffs))
@@ -323,7 +323,7 @@
     
     start = time.time()
-    wt = WorkerThreads(1, trainAux)
-    
-    num_trials = 200
+    wt = WorkerThreads(2, trainAux)
+    
+    num_trials = 20
     performance = []
     #for distance in range(1,46,4):
@@ -335,4 +335,9 @@
         fp,fn = errors
         performance.append(((fp+fn)/2.0, job_id, fn, fp))
+    #for distance in range(25,46,4):
+    #    job_id = distance
+    #    fp,fn = trainAux(distance, threshold, num_trials)
+    #    performance.append(((fp+fn)/2.0, job_id, fn, fp))
+    
     performance.sort()
     pprint.pprint(performance)
@@ -342,5 +347,5 @@
 
     
-    num_trials = 200
+    num_trials = 20
     start = time.time()
     performance = []
@@ -359,5 +364,5 @@
 
     
-    num_trials = 200
+    num_trials = 20
     start = time.time()
     performance = []
@@ -374,6 +379,6 @@
     best_distance = performance[0][1]
     print("best_distance:",best_distance)
-
-    num_trials = 200
+    
+    num_trials = 20
     start = time.time()
     performance = []
@@ -394,5 +399,5 @@
     return {'algorithm':"midhinge",
             'params':params,
-            'sample_size':subsample_size,
+            'sample_size':subseries_size,
             'num_trials':num_trials,
             'trial_type':"train",
@@ -421,10 +426,15 @@
 print(":", end-start)
 
+import cProfile
+
 start = time.time()
-results = trainMidhinge(db, unusual_case, greater, 6000)
-db.addClassifierResults(results)
+#cProfile.run('results = trainMidhinge(db, unusual_case, greater, 100)')
+results = trainMidhinge(db, unusual_case, greater, 100)
+#db.addClassifierResults(results)
 print("midhinge result:", results)
 end = time.time()
 print(":", end-start)
+
+sys.exit(0)
 
 start = time.time()
Index: trunk/lib/nanownlib/__init__.py
===================================================================
--- trunk/lib/nanownlib/__init__.py	(revision 5)
+++ trunk/lib/nanownlib/__init__.py	(revision 6)
@@ -214,8 +214,8 @@
     for p in packets:
         key = (p['sent'],p['tcpseq'],p['tcpack'],p['payload_len'])
-        #if (key not in seen)\
-        #   or p['sent']==1 and (seen[key]['observed'] < p['observed'])\
-        #   or p['sent']==0 and (seen[key]['observed'] > p['observed']):
-        if (key not in seen) or (seen[key]['observed'] > p['observed']):
+        if (key not in seen)\
+           or p['sent']==1 and (seen[key]['observed'] < p['observed'])\
+           or p['sent']==0 and (seen[key]['observed'] > p['observed']):
+            #if (key not in seen) or (seen[key]['observed'] > p['observed']):
             seen[key] = p
     
@@ -352,5 +352,5 @@
     pprint.pprint(evaluations)
 
-    delta_margin = 0.1
+    delta_margin = 0.15
     best_strim = 0
     best_rtrim = 0
@@ -439,5 +439,7 @@
     remaining_tm = trimean([row['packet_rtt'] for row in cursor])
 
-    return (tc, tm_map[tc]-remaining_tm)
+    ret_val = (tc, tm_map[tc]-remaining_tm)
+    print("unusual_case: %s, delta: %f" % ret_val)
+    return ret_val
 
 
Index: trunk/lib/nanownlib/stats.py
===================================================================
--- trunk/lib/nanownlib/stats.py	(revision 5)
+++ trunk/lib/nanownlib/stats.py	(revision 6)
@@ -227,4 +227,38 @@
         yield (sid,[dict(r) for r in probes])
 
+
+def subseries(db, probe_type, unusual_case, size=None, offset=None, field='packet_rtt'):
+    cursor = db.conn.cursor()
+    cursor.execute("SELECT max(c) FROM (SELECT count(sample) c FROM probes WHERE type=? GROUP BY test_case)", (probe_type,))
+    population_size = cursor.fetchone()[0]
+
+    if size == None or size > population_size:
+        size = population_size
+    if offset == None or offset >= population_size or offset < 0:
+        offset = numpy.random.random_integers(0,population_size-1)
+
+    query="""
+      SELECT %(field)s AS unusual_case,
+             (SELECT avg(%(field)s) FROM probes,analysis
+              WHERE analysis.probe_id=probes.id AND probes.test_case!=:unusual_case AND probes.type=:probe_type AND sample=u.sample) AS other_cases
+      FROM   (SELECT probes.sample,%(field)s FROM probes,analysis 
+              WHERE analysis.probe_id=probes.id AND probes.test_case =:unusual_case AND probes.type=:probe_type) u
+      LIMIT :size OFFSET :offset
+    """ % {"field":field}
+    
+    params = {"probe_type":probe_type, "unusual_case":unusual_case, "offset":offset, "size":size}
+    cursor.execute(query, params)
+    for row in cursor:
+        size -= 1
+        yield dict(row)
+
+    if size > 0:
+        params['offset'] = 0
+        params['size'] = size
+        cursor.execute(query, params)
+        for row in cursor:
+            yield dict(row)
+    
+
 # if test_cases=None, include all of them.  Otherwise, include only the specified test cases.
 def samples2Distributions(samples, field, test_cases=None):
@@ -270,4 +304,12 @@
 
 
+def bootstrap3(estimator, db, probe_type, unusual_case, subseries_size, num_trials):
+    ret_val = []
+    for t in range(num_trials):
+        ret_val.append(estimator(subseries(db, probe_type, unusual_case, subseries_size)))
+
+    return ret_val
+
+
 # Returns the test case name that clearly has higher RTT; otherwise, returns None
 def boxTest(params, test_cases, samples):
@@ -329,6 +371,6 @@
 # Returns 1 if unusual_case is unusual in the expected direction
 #         0 otherwise
-def midhingeTest(params, unusual_case, greater, samples):
-    diffs = list(samples2MeanDiffs(samples, 'packet_rtt', unusual_case))
+def midhingeTest(params, greater, samples):
+    diffs = [s['unusual_case']-s['other_cases'] for s in samples]
 
     mh = midhinge(diffs, params['distance'])
