Index: /trunk/bin/analyze
===================================================================
--- /trunk/bin/analyze	(revision 4)
+++ /trunk/bin/analyze	(revision 4)
@@ -0,0 +1,48 @@
+#!/usr/bin/env python3
+#-*- mode: Python;-*-
+
+import sys
+import os
+import time
+import random
+import tempfile
+import argparse
+import socket
+import json
+try:
+    import requests
+except:
+    sys.stderr.write('ERROR: Could not import requests module.  Ensure it is installed.\n')
+    sys.stderr.write('       Under Debian, the package name is "python3-requests"\n.')
+    sys.exit(1)
+
+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 *
+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)
+
+start = time.time()
+num_probes = analyzeProbes(db)
+end = time.time()
+print("analyzed %d probes' packets in: %f" % (num_probes, end-start))
Index: /trunk/bin/sampler
===================================================================
--- /trunk/bin/sampler	(revision 4)
+++ /trunk/bin/sampler	(revision 4)
@@ -0,0 +1,188 @@
+#!/usr/bin/env python3
+#-*- mode: Python;-*-
+
+import sys
+import os
+import time
+import random
+import tempfile
+import argparse
+import socket
+import json
+try:
+    import requests
+except:
+    sys.stderr.write('ERROR: Could not import requests module.  Ensure it is installed.\n')
+    sys.stderr.write('       Under Debian, the package name is "python3-requests"\n.')
+    sys.exit(1)
+
+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 *
+import nanownlib.storage
+
+parser = argparse.ArgumentParser(
+    description="")
+parser.add_argument('-c', dest='cases', type=str, default='{"short":10000,"long":13000}',
+                    help='JSON representation of echo timing cases.')
+parser.add_argument('--tcpts', action='store_true', help='Enable TCP timestamp profiling')
+parser.add_argument('session_name', default=None,
+                    help='Name for the sampler session (used in output filenames)')
+parser.add_argument('sample_count', type=int, default=None,
+                    help='Number of samples to collect')
+parser.add_argument('host', default=None,
+                    help='IP address or host name of server')
+parser.add_argument('port', nargs='?', type=int, default=80,
+                    help='TCP port number of HTTP service (default: 80)')
+options = parser.parse_args()
+
+
+session_name = options.session_name
+num_samples=options.sample_count
+hostname = options.host
+port = options.port
+protocol = 'http'
+
+cases = json.loads(options.cases)
+
+def extractReportedRuntime(headers, body):
+    if body.startswith("waited: "):
+        return int(body.split("waited: ")[1], 10)
+
+
+def sendRequest(data=None):
+    method = 'GET'
+    path = '/data/clamav-audit-results.txt?t=' + data
+    url = "%s://%s:%d%s" % (protocol,hostname,port,path)
+    headers = {}
+    body = (b'')
+    req = requests.Request('GET', url).prepare()
+
+    retry = True
+    while retry:
+        try:
+            session = requests.Session()
+            response = session.send(req)
+            #print(repr(response.raw._original_response.local_address))
+            reported = extractReportedRuntime(response.headers, response.text)
+            retry = False
+        except Exception as e:
+            sys.stderr.write("ERROR: HTTP request problem: %s\n" % repr(e))
+            time.sleep(1.0)
+            sys.stderr.write("ERROR: retrying...\n")
+    
+    return {'userspace_rtt':response.elapsed.microseconds*1000,
+            'reported':reported,
+            'local_port':response.raw._original_response.local_address[1]}
+
+
+def fetch(probedata, data):
+    #   http://docs.python-requests.org/en/latest/api/#requests.Response
+    result = sendRequest("%d" % data)
+    result.update(probedata)
+    
+    return result
+
+
+setCPUAffinity()
+setTCPTimestamps()
+host_ip = socket.gethostbyname(hostname) #XXX: what about multiple A records?
+
+meta = {'tcpts_mean':None,'tcpts_stddev':None,'tcpts_slopes':None}
+if options.tcpts:
+    sys.stderr.write("INFO: Probing target for TCP timestamp precision...\n")
+    sniffer_fp = tempfile.NamedTemporaryFile('w+t')
+    sniffer = startSniffer(host_ip, port, sniffer_fp.name)
+    time.sleep(1.0)
+    ports = runTimestampProbes(host_ip, port, hostname, 12)
+    time.sleep(1.0)
+    stopSniffer(sniffer)
+    sniffer_fp.seek(0)
+    mean,stddev,slopes = computeTimestampPrecision(sniffer_fp, ports)
+    meta = {'tcpts_mean':mean,'tcpts_stddev':stddev,'tcpts_slopes':json.dumps(slopes)}
+    
+if meta['tcpts_mean'] == None:
+    sys.stderr.write("INFO: TCP timestamps not supported.\n")
+elif meta['tcpts_stddev'] == None:
+    sys.stderr.write("INFO: Not enough TCP timestamp samples.\n")
+elif meta['tcpts_stddev']/meta['tcpts_mean'] > 0.05:
+    sys.stderr.write("INFO: TCP timestamp precision is inconsistent.\n")
+    sys.stderr.write("INFO:     mean=%(tcpts_mean)f,stddev=%(tcpts_stddev)f,samples=%(tcpts_slopes)s\n"
+                     % meta)
+else:
+    sys.stderr.write("INFO: Estimated TCP timestamp precision: %f (stddev: %f, %f%%)\n"
+                     % (meta['tcpts_mean'], meta['tcpts_stddev'], 100*meta['tcpts_stddev']/meta['tcpts_mean']))
+
+
+sniffer_fp = tempfile.NamedTemporaryFile('w+t')
+db_file = "%s.db" % session_name
+
+sniffer = startSniffer(host_ip, port, sniffer_fp.name)
+db = nanownlib.storage.db(db_file)
+db.addMeta(meta)
+time.sleep(0.5) # ensure sniffer is fully ready and our process is migrated
+
+sid = 0
+sample_types = [('train',int(num_samples*2/7)),
+                ('train_null',int(num_samples*3/7)+1),
+                ('test',int(num_samples*2/7)),
+                ('test_null',0)]
+report_interval = 20
+start = time.time()
+next_report = start+report_interval
+for st,count in sample_types:
+    if sniffer.poll() != None:
+        sys.stderr.write('ERROR: Sniffer quit, sender exiting...\n')
+        break
+        
+    for k in range(0,count):
+        sample_order = list(cases.items())
+        random.shuffle(sample_order)
+        if st.endswith('null'):
+            for i in range(1,len(sample_order)):
+                sample_order[i] = (sample_order[i][0],sample_order[0][1])
+            random.shuffle(sample_order)
+            #print('after', sample_order)
+            
+        results = []
+        now = int(time.time()*1000000000)
+        for i in range(len(sample_order)):
+            results.append(fetch({'sample':sid, 'test_case':sample_order[i][0],
+                                  'type':st, 'tc_order':i, 'time_of_day':now},
+                                 sample_order[i][1]))
+
+        #print(results)
+        db.addProbes(results)
+        db.conn.commit()
+        sid += 1
+
+        if (time.time() > next_report):
+            #s = time.time()
+            reportProgress(db, sample_types, start)
+            #print("reportProgress time:", time.time()-s)
+            next_report += report_interval
+
+print("probes complete in %f" % (time.time()-start))
+time.sleep(2.0) # Give sniffer a chance to collect remaining packets
+stopSniffer(sniffer)
+
+start = time.time()
+associatePackets(sniffer_fp, db)
+sniffer_fp.close()
+end = time.time()
+print("associate time:", end-start)
+
+start = time.time()
+num_probes = analyzeProbes(db)
+end = time.time()
+print("analyzed %d probes' packets in: %f" % (num_probes, end-start))
Index: /trunk/bin/train
===================================================================
--- /trunk/bin/train	(revision 4)
+++ /trunk/bin/train	(revision 4)
@@ -0,0 +1,442 @@
+#!/usr/bin/env python3
+#-*- mode: Python;-*-
+
+import sys
+import os
+import time
+import random
+import statistics
+import functools
+import argparse
+import threading
+import queue
+import pprint
+import json
+
+
+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 *
+import nanownlib.storage
+from nanownlib.stats import boxTest,multiBoxTest,subsample,bootstrap,bootstrap2,trimean,midhinge,midhingeTest,samples2Distributions,samples2MeanDiffs
+
+parser = argparse.ArgumentParser(
+    description="")
+#parser.add_argument('-c', dest='cases', type=str, default='{"short":10000,"long":1010000}',
+#                    help='JSON representation of echo timing cases. Default: {"short":10000,"long":1010000}')
+parser.add_argument('session_data', default=None,
+                    help='Database file storing session information')
+options = parser.parse_args()
+
+
+
+class WorkerThreads(object):
+    workq = None
+    resultq = None
+    target = None
+    
+    def __init__(self, num_workers, target):
+        self.workq = queue.Queue()
+        self.resultq = queue.Queue()
+        self.target = target
+        
+        self.workers = []
+        for i in range(num_workers):
+            t = threading.Thread(target=self._worker)
+            t.daemon = True
+            t.start()
+            self.workers.append(t)
+
+    def _worker(self):
+        while True:
+            item = self.workq.get()
+            if item == None:
+                self.workq.task_done()
+                break
+
+            job_id,args = item
+            self.resultq.put((job_id, self.target(*args)))
+            self.workq.task_done()
+
+    def addJob(self, job_id, args):
+        self.workq.put((job_id, args))
+            
+    def wait(self):
+        self.workq.join()
+
+    def stop(self):
+        for i in range(0,len(workers)):
+            self.workq.put(None)
+        for w in self.workers:
+            w.join()
+
+            
+def trainBoxTest(db, test_cases, longest, subsample_size):
+
+    def trainAux(low,high,num_trials):
+        estimator = functools.partial(boxTest, {'low':low, 'high':high})
+        estimates = bootstrap(estimator, db, 'train', test_cases, subsample_size, num_trials)
+        null_estimates = bootstrap(estimator, db, 'train_null', test_cases, subsample_size, num_trials)
+
+        #XXX: need to have a configurable policy on what we're looking for.
+        #     which is longest or which is shortest?
+        bad_estimates = len([e for e in estimates if e != longest])
+        bad_null_estimates = len([e for e in null_estimates if e != None])
+        
+        false_negatives = 100.0*bad_estimates/num_trials
+        false_positives = 100.0*bad_null_estimates/num_trials
+        return false_positives,false_negatives
+
+    start = time.time()
+    wt = WorkerThreads(2, trainAux)
+    
+    width = 2.0
+    performance = []
+    percentiles = list(range(0,50,2))
+    for low in percentiles:
+        wt.addJob(low, (low,low+width,200))
+    wt.wait()
+    while not wt.resultq.empty():
+        job_id,errors = wt.resultq.get()
+        fp,fn = errors
+        performance.append(((fp+fn)/2.0, job_id, fn, fp))
+    performance.sort()
+    pprint.pprint(performance)
+    print(time.time()-start)
+    
+    lows = [p[1] for p in performance[0:5]]
+    widths = [w/10.0 for w in range(0,105,5)]
+    performance = []
+    for width in widths:
+        false_positives = []
+        false_negatives = []
+        for low in lows:
+            wt.addJob(low,(low,low+width,150))
+        wt.wait()
+        while not wt.resultq.empty():
+            job_id,errors = wt.resultq.get()
+            fp,fn = errors
+            false_negatives.append(fn)
+            false_positives.append(fp)
+
+        #print(width, false_negatives)
+        #print(width, false_positives)
+        performance.append(((statistics.mean(false_positives)+statistics.mean(false_negatives))/2.0,
+                            width, statistics.mean(false_negatives), statistics.mean(false_positives)))
+    performance.sort()
+    pprint.pprint(performance)
+    good_width = performance[0][1]
+    print("good_width:",good_width)
+
+
+    lc = {}
+    for low in lows:
+        if low-1 > 0:
+            lc[low-1] = None
+        lc[low] = None
+        lc[low+1] = None
+    lows = lc.keys()
+    
+    performance = []
+    for low in lows:
+        wt.addJob(low, (low,low+good_width,300))
+    wt.wait()
+    while not wt.resultq.empty():
+        job_id,errors = wt.resultq.get()
+        fp,fn = errors
+        performance.append(((fp+fn)/2.0, job_id, fn, fp))
+    performance.sort()
+    pprint.pprint(performance)
+    best_low = performance[0][1]
+    print("best_low:", best_low)
+
+    
+    widths = [good_width-0.4,good_width-0.3,good_width-0.2,good_width-0.1,
+              good_width,good_width+0.1,good_width+0.2,good_width+0.3,good_width+0.4]
+    performance = []
+    for width in widths:
+        wt.addJob(width, (best_low,best_low+width,200))
+    wt.wait()
+    while not wt.resultq.empty():
+        job_id,errors = wt.resultq.get()
+        fp,fn = errors
+        performance.append(((fp+fn)/2.0, job_id, fn, fp))
+    performance.sort()
+    pprint.pprint(performance)
+    best_width=performance[0][1]
+    print("best_width:",best_width)
+    print("final_performance:", performance[0][0])
+
+    return {"low":best_low,"high":best_low+good_width}
+
+
+def trainBoxTest2(db, unusual_case, greater, subsample_size):
+
+    def trainAux(low,high,num_trials):
+        estimator = functools.partial(multiBoxTest, {'low':low, 'high':high}, unusual_case, greater)
+        estimates = bootstrap2(estimator, db, 'train', subsample_size, num_trials)
+        null_estimates = bootstrap2(estimator, db, 'train_null', subsample_size, num_trials)
+
+        bad_estimates = len([e for e in estimates if e != 1])
+        bad_null_estimates = len([e for e in null_estimates if e != 0])
+        
+        false_negatives = 100.0*bad_estimates/num_trials
+        false_positives = 100.0*bad_null_estimates/num_trials
+        return false_positives,false_negatives
+
+    start = time.time()
+    wt = WorkerThreads(2, trainAux)
+    
+    num_trials = 200
+    width = 2.0
+    performance = []
+    percentiles = list(range(0,50,2))
+    for low in percentiles:
+        wt.addJob(low, (low,low+width,num_trials))
+    wt.wait()
+    while not wt.resultq.empty():
+        job_id,errors = wt.resultq.get()
+        fp,fn = errors
+        performance.append(((fp+fn)/2.0, job_id, fn, fp))
+    performance.sort()
+    pprint.pprint(performance)
+    print(time.time()-start)
+    
+    num_trials = 150
+    lows = [p[1] for p in performance[0:5]]
+    widths = [w/10.0 for w in range(0,105,5)]
+    performance = []
+    for width in widths:
+        false_positives = []
+        false_negatives = []
+        for low in lows:
+            wt.addJob(low,(low,low+width,num_trials))
+        wt.wait()
+        while not wt.resultq.empty():
+            job_id,errors = wt.resultq.get()
+            fp,fn = errors
+            false_negatives.append(fn)
+            false_positives.append(fp)
+
+        #print(width, false_negatives)
+        #print(width, false_positives)
+        performance.append(((statistics.mean(false_positives)+statistics.mean(false_negatives))/2.0,
+                            width, statistics.mean(false_negatives), statistics.mean(false_positives)))
+    performance.sort()
+    pprint.pprint(performance)
+    good_width = performance[0][1]
+    print("good_width:",good_width)
+
+
+    lc = {}
+    for low in lows:
+        if low-1 >= 0:
+            lc[low-1] = None
+        lc[low] = None
+        lc[low+1] = None
+    lows = lc.keys()
+    print("candidate lows:")
+    pprint.pprint(lows)
+    
+    num_trials = 300
+    performance = []
+    for low in lows:
+        wt.addJob(low, (low,low+good_width,num_trials))
+    wt.wait()
+    while not wt.resultq.empty():
+        job_id,errors = wt.resultq.get()
+        fp,fn = errors
+        performance.append(((fp+fn)/2.0, job_id, fn, fp))
+    performance.sort()
+    pprint.pprint(performance)
+    best_low = performance[0][1]
+    print("best_low:", best_low)
+
+    num_trials = 200
+    widths = [good_width-0.4,good_width-0.3,good_width-0.2,good_width-0.1,
+              good_width,good_width+0.1,good_width+0.2,good_width+0.3,good_width+0.4]
+    performance = []
+    for width in widths:
+        wt.addJob(width, (best_low,best_low+width,num_trials))
+    wt.wait()
+    while not wt.resultq.empty():
+        job_id,errors = wt.resultq.get()
+        fp,fn = errors
+        performance.append(((fp+fn)/2.0, job_id, fn, fp))
+    performance.sort()
+    pprint.pprint(performance)
+    best_width=performance[0][1]
+    print("best_width:",best_width)
+    print("final_performance:", performance[0][0])
+    
+    params = json.dumps({"low":best_low,"high":best_low+good_width})
+    return {'algorithm':"boxtest",
+            'params':params,
+            'sample_size':subsample_size,
+            'num_trials':num_trials,
+            'trial_type':"train",
+            'false_positives':performance[0][3],
+            'false_negatives':performance[0][2]}
+
+
+def trainMidhinge(db, unusual_case, greater, subsample_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)
+
+        bad_estimates = len([e for e in estimates if e != 1])
+        bad_null_estimates = len([e for e in null_estimates if e != 0])
+        
+        false_negatives = 100.0*bad_estimates/num_trials
+        false_positives = 100.0*bad_null_estimates/num_trials
+        return false_positives,false_negatives
+
+    #determine expected delta based on differences
+    start = time.time()
+    mean_diffs = list(samples2MeanDiffs(subsample(db, 'train'), 'packet_rtt', unusual_case))
+    threshold = trimean(mean_diffs)/2.0
+    print("initial threshold:", threshold)
+    print("median threshold:", statistics.median(mean_diffs)/2.0)
+    print("midhinge threshold:", midhinge(mean_diffs)/2.0)
+    print("trimean threshold:", trimean(mean_diffs)/2.0)
+    
+    mean_diffs = list(samples2MeanDiffs(subsample(db, 'train_null'), 'packet_rtt', unusual_case))
+    print(len(mean_diffs))
+    print("null mean:", statistics.mean(mean_diffs))
+    print("null median:", statistics.median(mean_diffs))
+    print("null midhinge:", midhinge(mean_diffs))
+    print("null trimean:", trimean(mean_diffs))
+    print(time.time()-start)
+
+    
+    start = time.time()
+    wt = WorkerThreads(1, trainAux)
+    
+    num_trials = 200
+    performance = []
+    #for distance in range(1,46,4):
+    for distance in range(25,46,4):
+        wt.addJob(distance, (distance,threshold,num_trials))
+    wt.wait()
+    while not wt.resultq.empty():
+        job_id,errors = wt.resultq.get()
+        fp,fn = errors
+        performance.append(((fp+fn)/2.0, job_id, fn, fp))
+    performance.sort()
+    pprint.pprint(performance)
+    print(time.time()-start)
+    good_distance = performance[0][1]
+    print("good_distance:",good_distance)
+
+    
+    num_trials = 200
+    start = time.time()
+    performance = []
+    for t in range(80,125,5):
+        wt.addJob(threshold*(t/100.0), (good_distance,threshold*(t/100.0),num_trials))
+    wt.wait()
+    while not wt.resultq.empty():
+        job_id,errors = wt.resultq.get()
+        fp,fn = errors
+        performance.append(((fp+fn)/2.0, job_id, fn, fp))
+    performance.sort()
+    pprint.pprint(performance)
+    print(time.time()-start)
+    good_threshold = performance[0][1]
+    print("good_threshold:", good_threshold)
+
+    
+    num_trials = 200
+    start = time.time()
+    performance = []
+    for d in range(-4,5):
+        wt.addJob(good_distance+d, (good_distance+d,good_threshold,num_trials))
+    wt.wait()
+    while not wt.resultq.empty():
+        job_id,errors = wt.resultq.get()
+        fp,fn = errors
+        performance.append(((fp+fn)/2.0, job_id, fn, fp))
+    performance.sort()
+    pprint.pprint(performance)
+    print(time.time()-start)
+    best_distance = performance[0][1]
+    print("best_distance:",best_distance)
+
+    num_trials = 200
+    start = time.time()
+    performance = []
+    for t in range(95,106):
+        wt.addJob(good_threshold*(t/100.0), (best_distance,good_threshold*(t/100.0),num_trials))
+    wt.wait()
+    while not wt.resultq.empty():
+        job_id,errors = wt.resultq.get()
+        fp,fn = errors
+        performance.append(((fp+fn)/2.0, job_id, fn, fp))
+    performance.sort()
+    pprint.pprint(performance)
+    print(time.time()-start)
+    best_threshold = performance[0][1]
+    print("best_threshold:", best_threshold)
+
+    params = json.dumps({'distance':best_distance,'threshold':best_threshold})
+    return {'algorithm':"midhinge",
+            'params':params,
+            'sample_size':subsample_size,
+            'num_trials':num_trials,
+            'trial_type':"train",
+            'false_positives':performance[0][3],
+            'false_negatives':performance[0][2]}
+
+
+#classifiers = {'boxtest':{'train':trainBoxTest2, 'test':multiBoxTest},
+#               'midhinge':{'train':trainMidhinge, 'test':midhinge}}
+
+
+db = nanownlib.storage.db(options.session_data)
+#cursor = db.cursor()
+#cursor.execute("SELECT min(sample) min, max(sample) max FROM probes")
+#train_start,test_end = cursor.fetchone()
+#train_end = int(test_end-train_start)
+#test_start = train_end+1
+#subsample_size = min(10000,(train_end-train_start+1)/4)
+
+start = time.time()
+unusual_case,unusual_diff = findUnusualTestCase(db)
+greater = (unusual_diff > 0)
+print("unusual_case:", unusual_case)
+print("unusual_diff:", unusual_diff)
+end = time.time()
+print(":", end-start)
+
+start = time.time()
+results = trainMidhinge(db, unusual_case, greater, 6000)
+db.addClassifierResults(results)
+print("midhinge result:", results)
+end = time.time()
+print(":", end-start)
+
+start = time.time()
+results = trainBoxTest2(db, unusual_case, greater, 6000)
+db.addClassifierResults(results)
+print("multi box test result:", results)
+end = time.time()
+print(":", end-start)
+
+#start = time.time()
+#print("box test params:", trainBoxTest(db, test_cases, 'long', 100))
+#end = time.time()
+#print(":", end-start)
+
+
Index: /trunk/lib/nanownlib/__init__.py
===================================================================
--- /trunk/lib/nanownlib/__init__.py	(revision 4)
+++ /trunk/lib/nanownlib/__init__.py	(revision 4)
@@ -0,0 +1,456 @@
+#!/usr/bin/env python3
+#-*- mode: Python;-*-
+
+import sys
+import time
+import random
+import argparse
+import socket
+import datetime
+import http.client
+import threading
+import queue
+import subprocess
+import multiprocessing
+import csv
+import json
+import gzip
+import statistics
+import numpy
+import netifaces
+try:
+    import requests
+except:
+    sys.stderr.write('ERROR: Could not import requests module.  Ensure it is installed.\n')
+    sys.stderr.write('       Under Debian, the package name is "python3-requests"\n.')
+    sys.exit(1)
+
+from .stats import *
+
+
+def getLocalIP(remote_host, remote_port):
+    connection = socket.create_connection((remote_host, remote_port))
+    ret_val = connection.getsockname()[0]
+    connection.close()
+
+    return ret_val
+
+
+def getIfaceForIP(ip):
+    for iface in netifaces.interfaces():
+        addrs = netifaces.ifaddresses(iface).get(netifaces.AF_INET, None)
+        if addrs:
+            for a in addrs:
+                if a.get('addr', None) == ip:
+                    return iface
+
+
+def setTCPTimestamps(enabled=True):
+    fh = open('/proc/sys/net/ipv4/tcp_timestamps', 'r+b')
+    ret_val = False
+    if fh.read(1) == b'1':
+        ret_val = True
+
+    fh.seek(0)
+    if enabled:
+        fh.write(b'1')
+    else:
+        fh.write(b'0')
+    fh.close()
+    
+    return ret_val
+
+
+def trickleHTTPRequest(ip,port,hostname):
+    my_port = None
+    try:
+        sock = socket.create_connection((ip, port))
+        my_port = sock.getsockname()[1]
+        
+        #print('.')
+        sock.sendall(b'GET / HTTP/1.1\r\n')
+        time.sleep(0.5)
+        rest = b'''Host: '''+hostname.encode('utf-8')+b'''\r\nUser-Agent: Secret Agent Man\r\nX-Extra: extra read all about it!\r\nConnection: close\r\n'''
+        for r in rest:
+            sock.sendall(bytearray([r]))
+            time.sleep(0.05)
+
+        time.sleep(0.5)
+        sock.sendall('\r\n')
+
+        r = None
+        while r != b'':
+            r = sock.recv(16)
+
+        sock.close()
+    except Exception as e:
+        pass
+
+    return my_port
+
+
+def runTimestampProbes(host_ip, port, hostname, num_trials, concurrency=4): 
+    myq = queue.Queue()
+    def threadWrapper(*args):
+        try:
+            myq.put(trickleHTTPRequest(*args))
+        except Exception as e:
+            sys.stderr.write("ERROR from trickleHTTPRequest: %s\n" % repr(e))
+            myq.put(None)
+
+    threads = []
+    ports = []
+    for i in range(num_trials):
+        if len(threads) >= concurrency:
+            ports.append(myq.get())
+        t = threading.Thread(target=threadWrapper, args=(host_ip, port, hostname))
+        t.start()
+        threads.append(t)
+
+    for t in threads:
+        t.join()
+
+    while myq.qsize() > 0:
+        ports.append(myq.get())
+
+    return ports
+
+
+def computeTimestampPrecision(sniffer_fp, ports):
+    rcvd = []
+    for line in sniffer_fp:
+        p = json.loads(line)
+        if p['sent']==0:
+            rcvd.append((p['observed'],p['tsval'],int(p['local_port'])))
+
+    slopes = []
+    for port in ports:
+        trcvd = [tr for tr in rcvd if tr[2]==port and tr[1]!=0]
+
+        if len(trcvd) < 2:
+            sys.stderr.write("WARN: Inadequate data points.\n")
+            continue
+        
+        if trcvd[0][1] > trcvd[-1][1]:
+            sys.stderr.write("WARN: TSval wrap.\n")
+            continue
+
+        x = [tr[1] for tr in trcvd]
+        y = [tr[0] for tr in trcvd]
+
+        slope,intercept = OLSRegression(x, y)
+        slopes.append(slope)
+
+    if len(slopes) == 0:
+        return None,None,None
+
+    m = statistics.mean(slopes)
+    if len(slopes) == 1:
+        return (m, None, slopes)
+    else:
+        return (m, statistics.stdev(slopes), slopes)
+
+    
+def OLSRegression(x,y):
+    #print(x,y)
+    x = numpy.array(x)
+    y = numpy.array(y)
+    #A = numpy.vstack([x, numpy.ones(len(x))]).T
+    #m, c = numpy.linalg.lstsq(A, y)[0] # broken
+    #c,m = numpy.polynomial.polynomial.polyfit(x, y, 1) # less accurate
+    c,m = numpy.polynomial.Polynomial.fit(x,y,1).convert().coef
+
+    #print(m,c)
+
+    #import matplotlib.pyplot as plt
+    #plt.clf()
+    #plt.scatter(x, y)
+    #plt.plot(x, m*x + c, 'r', label='Fitted line')
+    #plt.show()
+    
+    return (m,c)
+
+
+def startSniffer(target_ip, target_port, output_file):
+    my_ip = getLocalIP(target_ip, target_port)
+    my_iface = getIfaceForIP(my_ip)
+    return subprocess.Popen(['chrt', '-r', '99', './bin/csamp', my_iface, my_ip,
+                             target_ip, "%d" % target_port, output_file, '0'])
+
+def stopSniffer(sniffer):
+    sniffer.terminate()
+    sniffer.wait(2)
+    if sniffer.poll() == None:
+        sniffer.kill()
+        sniffer.wait(1)
+
+        
+def setCPUAffinity():
+    import ctypes
+    from ctypes import cdll,c_int,byref
+    cpus = multiprocessing.cpu_count()
+    
+    libc = cdll.LoadLibrary("libc.so.6")
+    #libc.sched_setaffinity(os.getpid(), 1, ctypes.byref(ctypes.c_int(0x01)))
+    return libc.sched_setaffinity(0, 4, byref(c_int(0x00000001<<(cpus-1))))
+
+
+# Monkey patching that instruments the HTTPResponse to collect connection source port info
+class MonitoredHTTPResponse(http.client.HTTPResponse):
+    local_address = None
+
+    def __init__(self, sock, *args, **kwargs):
+        self.local_address = sock.getsockname()
+        super(MonitoredHTTPResponse, self).__init__(sock,*args,**kwargs)
+            
+requests.packages.urllib3.connection.HTTPConnection.response_class = MonitoredHTTPResponse
+
+
+def removeDuplicatePackets(packets):
+    #return packets
+    suspect = None
+    seen = {}
+    # XXX: Need to review this deduplication algorithm and make sure it is correct
+    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']):
+            seen[key] = p
+    
+    if len(seen) < len(packets):
+        suspect = 'd'
+        #sys.stderr.write("INFO: removed %d duplicate packets.\n" % (len(packets) - len(seen)))
+
+    return suspect,seen.values()
+
+
+def analyzePackets(packets, timestamp_precision, trim_sent=0, trim_rcvd=0):
+    suspect,packets = removeDuplicatePackets(packets)
+
+    #sort_key = lambda d: (d['tcpseq'],d['tcpack'])
+    sort_key = lambda d: (d['observed'],d['tcpseq'])
+    sent = sorted((p for p in packets if p['sent']==1 and p['payload_len']>0), key=sort_key)
+    rcvd = sorted((p for p in packets if p['sent']==0 and p['payload_len']>0), key=sort_key)
+
+    if len(sent) <= trim_sent:
+        last_sent = sent[-1]
+    else:
+        last_sent = sent[trim_sent]
+
+    if len(rcvd) <= trim_rcvd:
+        last_rcvd = rcvd[0]
+    else:
+        last_rcvd = rcvd[len(rcvd)-1-trim_rcvd]
+    
+    packet_rtt = last_rcvd['observed'] - last_sent['observed']
+    if packet_rtt < 0:
+        sys.stderr.write("WARN: Negative packet_rtt. last_rcvd=%s,last_sent=%s\n" % (last_rcvd, last_sent))
+
+    last_sent_ack = None
+    try:
+        last_sent_ack = min(((p['observed'],p) for p in packets
+                             if p['sent']==0 and p['payload_len']+last_sent['tcpseq']==p['tcpack']))[1]
+    except Exception as e:
+        sys.stderr.write("WARN: Could not find last_sent_ack.\n")
+
+    tsval_rtt = None
+    if None not in (timestamp_precision, last_sent_ack):
+        tsval_rtt = int(round((last_rcvd['tsval'] - last_sent_ack['tsval'])*timestamp_precision))
+
+    return {'packet_rtt':packet_rtt,
+            'tsval_rtt':tsval_rtt,
+            'suspect':suspect,
+            'sent_trimmed':trim_sent,
+            'rcvd_trimmed':trim_rcvd},len(sent),len(rcvd)
+
+
+# trimean and mad for each dist of differences
+def evaluateTrim(db, unusual_case, strim, rtrim):
+    cursor = db.conn.cursor()
+    query="""
+      SELECT packet_rtt-(SELECT avg(packet_rtt) FROM probes,trim_analysis 
+                         WHERE sent_trimmed=:strim AND rcvd_trimmed=:rtrim AND trim_analysis.probe_id=probes.id AND probes.test_case!=:unusual_case AND sample=u.sample AND probes.type in ('train','test'))
+      FROM (SELECT probes.sample,packet_rtt FROM probes,trim_analysis WHERE sent_trimmed=:strim AND rcvd_trimmed=:rtrim AND trim_analysis.probe_id=probes.id AND probes.test_case=:unusual_case AND probes.type in ('train','test')) u
+    """
+
+    params = {"strim":strim,"rtrim":rtrim,"unusual_case":unusual_case}
+    cursor.execute(query, params)
+    differences = [row[0] for row in cursor]
+    
+    return trimean(differences),mad(differences)
+
+
+
+def analyzeProbes(db):
+    db.conn.execute("CREATE INDEX IF NOT EXISTS packets_probe ON packets (probe_id)")
+    pcursor = db.conn.cursor()
+    kcursor = db.conn.cursor()
+
+    pcursor.execute("SELECT tcpts_mean FROM meta")
+    try:
+        timestamp_precision = pcursor.fetchone()[0]
+    except:
+        timestamp_precision = None
+    
+    pcursor.execute("DELETE FROM trim_analysis")
+    db.conn.commit()
+    
+    count = 0
+    sent_tally = []
+    rcvd_tally = []
+    for pid, in pcursor.execute("SELECT id FROM probes"):
+        kcursor.execute("SELECT * FROM packets WHERE probe_id=?", (pid,))
+        try:
+            analysis,s,r = analyzePackets(kcursor.fetchall(), timestamp_precision)
+            analysis['probe_id'] = pid
+            sent_tally.append(s)
+            rcvd_tally.append(r)
+        except Exception as e:
+            print(e)
+            sys.stderr.write("WARN: couldn't find enough packets for probe_id=%s\n" % pid)
+        
+        #print(pid,analysis)
+        db.addTrimAnalyses([analysis])
+        count += 1
+    db.conn.commit()
+    num_sent = statistics.mode(sent_tally)
+    num_rcvd = statistics.mode(rcvd_tally)
+    sent_tally = None
+    rcvd_tally = None
+    print("num_sent: %d, num_rcvd: %d" % (num_sent,num_rcvd))
+    
+    for strim in range(0,num_sent):
+        for rtrim in range(0,num_rcvd):
+            if strim == 0 and rtrim == 0:
+                continue # no point in doing 0,0 again
+            for pid, in pcursor.execute("SELECT id FROM probes"):
+                kcursor.execute("SELECT * FROM packets WHERE probe_id=?", (pid,))
+                try:
+                    analysis,s,r = analyzePackets(kcursor.fetchall(), timestamp_precision, strim, rtrim)
+                    analysis['probe_id'] = pid
+                except Exception as e:
+                    print(e)
+                    sys.stderr.write("WARN: couldn't find enough packets for probe_id=%s\n" % pid)
+                    
+                db.addTrimAnalyses([analysis])
+            db.conn.commit()
+
+    # Populate analysis table so findUnusualTestCase can give us a starting point
+    pcursor.execute("DELETE FROM analysis")
+    db.conn.commit()
+    pcursor.execute("INSERT INTO analysis SELECT id,probe_id,suspect,packet_rtt,tsval_rtt FROM trim_analysis WHERE sent_trimmed=0 AND rcvd_trimmed=0")
+    
+    unusual_case,delta = findUnusualTestCase(db)
+    evaluations = {}
+    for strim in range(0,num_sent):
+        for rtrim in range(0,num_rcvd):
+            evaluations[(strim,rtrim)] = evaluateTrim(db, unusual_case, strim, rtrim)
+
+    import pprint
+    pprint.pprint(evaluations)
+
+    delta_margin = 0.1
+    best_strim = 0
+    best_rtrim = 0
+    good_delta,good_mad = evaluations[(0,0)]
+    
+    for strim in range(1,num_sent):
+        delta,mad = evaluations[(strim,0)]
+        if abs(good_delta - delta) < abs(delta_margin*good_delta) and mad < good_mad:
+            best_strim = strim
+        else:
+            break
+
+    good_delta,good_mad = evaluations[(best_strim,0)]
+    for rtrim in range(1,num_rcvd):
+        delta,mad = evaluations[(best_strim,rtrim)]
+        if (abs(delta) > abs(good_delta) or abs(good_delta - delta) < abs(delta_margin*good_delta)) and mad < good_mad:
+            best_rtrim = rtrim
+        else:
+            break
+
+    print("selected trim parameters:",(best_strim,best_rtrim))
+    
+    if best_strim != 0 or best_rtrim !=0:
+        pcursor.execute("DELETE FROM analysis")
+        db.conn.commit()
+        pcursor.execute("INSERT INTO analysis SELECT id,probe_id,suspect,packet_rtt,tsval_rtt FROM trim_analysis WHERE sent_trimmed=? AND rcvd_trimmed=?",
+                        (best_strim,best_rtrim))
+    
+    return count
+
+
+        
+def parseJSONLines(fp):
+    for line in fp:
+        yield json.loads(line)
+
+
+def associatePackets(sniffer_fp, db):
+    sniffer_fp.seek(0)
+
+    # now combine sampler data with packet data
+    buffered = []
+
+    cursor = db.conn.cursor()
+    cursor.execute("SELECT count(*) count,min(time_of_day) start,max(time_of_day+userspace_rtt) end from probes")
+    ptimes = cursor.fetchone()
+    window_size = 100*int((ptimes['end']-ptimes['start'])/ptimes['count'])
+    print("associate window_size:", window_size)
+
+    db.addPackets(parseJSONLines(sniffer_fp), window_size)
+
+    cursor.execute("SELECT count(*) count FROM packets WHERE probe_id is NULL")
+    unmatched = cursor.fetchone()['count']
+    if unmatched > 0:
+        sys.stderr.write("WARNING: %d observed packets didn't find a home...\n" % unmatched)
+ 
+    return None
+
+
+def enumStoredTestCases(db):
+    cursor = db.conn.cursor()
+    cursor.execute("SELECT test_case FROM probes GROUP BY test_case")
+    return [tc[0] for tc in cursor]
+
+
+def findUnusualTestCase(db):
+    test_cases = enumStoredTestCases(db)
+
+    cursor = db.conn.cursor()
+    cursor.execute("SELECT packet_rtt FROM probes,analysis WHERE probes.id=analysis.probe_id AND probes.type in ('train','test')")
+    global_tm = trimean([row['packet_rtt'] for row in cursor])
+
+    tm_abs = []
+    tm_map = {}
+    # XXX: if more speed needed, percentile extension to sqlite might be handy...
+    for tc in test_cases:
+        cursor.execute("SELECT packet_rtt FROM probes,analysis WHERE probes.id=analysis.probe_id AND probes.type in ('train','test') AND probes.test_case=?", (tc,))
+        tm_map[tc] = trimean([row['packet_rtt'] for row in cursor])
+        tm_abs.append((abs(tm_map[tc]-global_tm), tc))
+
+    magnitude,tc = max(tm_abs)
+    cursor.execute("SELECT packet_rtt FROM probes,analysis WHERE probes.id=analysis.probe_id AND probes.type in ('train','test') AND probes.test_case<>?", (tc,))
+    remaining_tm = trimean([row['packet_rtt'] for row in cursor])
+
+    return (tc, tm_map[tc]-remaining_tm)
+
+
+def reportProgress(db, sample_types, start_time):
+    cursor = db.conn.cursor()
+    output = ''
+    total_completed = 0
+    total_requested = 0
+    for st in sample_types:
+        cursor.execute("SELECT count(id) c FROM (SELECT id FROM probes WHERE type=? AND time_of_day>? GROUP BY sample)", (st[0],int(start_time*1000000000)))
+        count = cursor.fetchone()[0]
+        output += " | %s remaining: %d" % (st[0], st[1]-count)
+        total_completed += count
+        total_requested += st[1]
+
+    rate = total_completed / (time.time() - start_time)
+    total_time = total_requested / rate        
+    eta = datetime.datetime.fromtimestamp(start_time+total_time)
+    print("STATUS:",output[3:],"| est. total_time: %s | est. ETA: %s" % (str(datetime.timedelta(seconds=total_time)), str(eta)))
Index: /trunk/lib/nanownlib/stats.py
===================================================================
--- /trunk/lib/nanownlib/stats.py	(revision 4)
+++ /trunk/lib/nanownlib/stats.py	(revision 4)
@@ -0,0 +1,352 @@
+
+import sys
+import os
+import math
+import statistics
+import gzip
+import random
+import scipy
+import scipy.stats
+import numpy
+
+# Don't trust numpy's seeding
+numpy.random.seed(random.SystemRandom().randint(0,2**32-1))
+
+
+def mad(arr):
+    """ Median Absolute Deviation: a "Robust" version of standard deviation.
+        Indices variabililty of the sample.
+        https://en.wikipedia.org/wiki/Median_absolute_deviation 
+    """
+    arr = numpy.ma.array(arr).compressed() # should be faster to not use masked arrays.
+    med = numpy.median(arr)
+    return numpy.median(numpy.abs(arr - med))
+
+
+def cov(x,y):
+    mx = statistics.mean(x)
+    my = statistics.mean(y)
+    products = []
+    for i in range(0,len(x)):
+        products.append((x[i] - mx)*(y[i] - my))
+
+    return statistics.mean(products)
+
+
+def difference(ls):
+    return ls[0]-ls[1]
+
+def product(ls):
+    return ls[0]*ls[1]
+
+def hypotenuse(ls):
+    return math.hypot(ls[0],ls[1])
+    
+def trustValues(derived, trustFunc):
+    ret_val = []
+    for k,v in derived.items():
+        ret_val.append((trustFunc((v['long'],v['short'])), k))
+
+    ret_val.sort()
+    return ret_val
+
+
+def prunedWeights(derived, trust, alpha):
+    weights = {}
+
+    threshold = len(trust)*(1.0-alpha)
+    for i in range(0,len(trust)):
+        if i < threshold:
+            weights[trust[i][1]] = 1.0
+        else:
+            weights[trust[i][1]] = 0.0
+
+    return weights
+
+
+def linearWeights(derived, trust, alpha):
+    x1 = trust[0][0]
+    y1 = 1.0 + (alpha*10)
+    x2 = trust[(len(trust)-1)//3][0]
+    y2 = 1.0
+    m = (y1-y2)/(x1-x2)
+    b = y1 - m*x1
+
+    weights = {}
+    for t,k in trust:
+        weights[k] = m*t+b
+        if weights[k] < 0.0:
+            weights[k] = 0.0
+
+    return weights
+
+
+def invertedWeights(derived,trust,alpha):
+    # (x+1-first_sample)^(-alpha)
+    #scale = trust[0][0]
+
+    #weights = {}
+    #for t,k in trust:
+    #    weights[k] = (t+1-scale)**(-1.0*alpha)
+    #    if weights[k] < 0.0:
+    #        weights[k] = 0.0
+
+    weights = {}
+    for i in range(len(trust)):
+        w = 10.0/(i+2.0)-0.2
+        if w < 0.0:
+            w = 0.0
+        weights[trust[i][1]] = w
+        
+    
+    return weights
+
+
+
+def arctanWeights(derived,trust,alpha):
+    shift = trust[int((len(trust)-1)*(1.0-alpha))][0]
+    minimum = trust[0][0]
+    
+    weights = {}
+    for i in range(len(trust)):
+        w = math.pi/2.0 - math.atan(2*(trust[i][0] - shift)/(shift-minimum))
+        if w < 0.0:
+            w = 0.0
+        weights[trust[i][1]] = w
+
+    return weights
+
+
+def arctanWeights2(derived,trust,alpha):
+    shift = trust[int((len(trust)-1)*(1.0-alpha))][0]
+    minimum = trust[0][0]
+    stretch = trust[int((len(trust)-1)*0.5)][0] - minimum # near median
+    
+    weights = {}
+    for i in range(len(trust)):
+        w = math.pi/2.0 - math.atan(3*(trust[i][0] - shift)/(shift-minimum))
+        if w < 0.0:
+            w = 0.0
+        weights[trust[i][1]] = w
+
+    return weights
+
+
+def midhinge(values, distance=25):
+    return (numpy.percentile(values, 50-distance) + numpy.percentile(values, 50+distance))/2.0
+
+def trimean(values, distance=25):
+    return (midhinge(values, distance) + statistics.median(values))/2
+
+def weightedMean(derived, weights):
+    normalizer = sum(weights.values())/len(weights)
+    return statistics.mean([w*(derived[k]['long']-derived[k]['short'])/normalizer for k,w in weights.items()])
+
+def weightedMeanTsval(derived, weights):
+    normalizer = sum(weights.values())/len(weights)
+    return statistics.mean([w*(derived[k]['long_tsval']-derived[k]['short_tsval'])/normalizer for k,w in weights.items()])
+
+
+def estimateMean(trustFunc, weightFunc, alpha, derived):
+    trust = trustValues(derived, trustFunc)
+    weights = weightFunc(derived, trust, alpha)
+    return weightedMean(derived, weights)
+
+
+def estimateMeanTsval(trustFunc, weightFunc, alpha, derived):
+    trust = trustValues(derived, trustFunc)
+    weights = weightFunc(derived, trust, alpha)
+    return weightedMeanTsval(derived, weights)
+
+
+#def estimateMedian(trustFunc, weightFunc, alpha, derived):
+#    trust = trustValues(derived, trustFunc)
+#    weights = weightFunc(derived, trust, alpha)
+
+#    return statistics.median([(derived[k]['long']-derived[k]['short']) for k,w in weights.items() if w > 0.0])
+
+def estimateMedian(derived):
+    return statistics.median([(d['long']-d['short']) for d in derived.values()])
+
+
+def estimateMidhinge(derived):
+    return midhinge([(d['long']-d['short']) for d in derived.values()])
+
+
+def estimateTrimean(derived):
+    return trimean([(d['long']-d['short']) for d in derived.values()])
+
+
+def tTest(expected_mean, derived):
+    diffs = [(d['long']-d['short']) for d in derived.values()]
+    null_tval, null_pval = scipy.stats.ttest_1samp(diffs, 0.0)
+    tval, pval = scipy.stats.ttest_1samp(diffs, expected_mean)
+
+    if pval < null_pval:
+        return 1
+    else:
+        return 0
+
+    
+def diffMedian(derived):
+    l = [tc['long']-tc['short'] for s,tc in derived.items()]
+    return statistics.median(l)
+
+
+def subsample_ids(db, probe_type, subsample_size=None):
+    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]
+    #print("population_size:", population_size)
+    if subsample_size == None or subsample_size > population_size:
+        subsample_size = population_size
+    
+    start = numpy.random.random_integers(0,population_size-1)
+    cursor.execute("SELECT sample FROM probes WHERE type=? GROUP BY sample ORDER BY sample LIMIT ? OFFSET ?", (probe_type,subsample_size,start))
+    for row in cursor:
+        subsample_size -= 1
+        yield row['sample']
+
+    if subsample_size > 0:
+        cursor.execute("SELECT sample FROM probes WHERE type=? GROUP BY sample ORDER BY sample LIMIT ?", (probe_type,subsample_size))
+        for row in cursor:
+            yield row['sample']
+    
+
+def subsample(db, probe_type, subsample_size=None):
+    cursor = db.conn.cursor()
+    cursor.execute("SELECT count(test_case) FROM (SELECT test_case FROM probes GROUP BY test_case)")
+    num_test_cases = cursor.fetchone()[0]
+    
+    for sid in subsample_ids(db, probe_type, subsample_size):
+        cursor.execute("SELECT test_case,tc_order,time_of_day,reported,userspace_rtt,suspect,packet_rtt,tsval_rtt FROM probes p,analysis a WHERE p.sample=? and a.probe_id=p.id", (sid,))
+        probes = cursor.fetchall()
+        if len(probes) != num_test_cases:
+            sys.stderr.write("WARN: sample %d had %d probes, but %d expected!  Discarding...\n" % (sid, len(probes), num_test_cases))
+            continue
+        yield (sid,[dict(r) for r in probes])
+
+# if test_cases=None, include all of them.  Otherwise, include only the specified test cases.
+def samples2Distributions(samples, field, test_cases=None):
+    ret_val = {}
+    
+    for sid,probes in samples:
+        for p in probes:
+            if p['test_case'] in ret_val:
+                ret_val[p['test_case']].append(p[field])
+            elif test_cases == None or p['test_case'] in test_cases:
+                ret_val[p['test_case']] = [p[field]]
+                
+                
+    return ret_val
+
+
+def samples2MeanDiffs(samples, field, unusual_case):
+    ret_val = {}
+    
+    for sid,probes in samples:
+        unusual_value = None
+        for p in probes:
+            if p['test_case'] == unusual_case:
+                unusual_value = p[field]
+                break
+        yield statistics.mean([unusual_value-p[field] for p in probes if p['test_case'] != unusual_case])
+
+
+def bootstrap(estimator, db, probe_type, test_cases, subsample_size, num_trials):
+    ret_val = []
+    for t in range(num_trials):
+        ret_val.append(estimator(test_cases, subsample(db, probe_type, subsample_size)))
+
+    return ret_val
+
+
+def bootstrap2(estimator, db, probe_type, subsample_size, num_trials):
+    ret_val = []
+    for t in range(num_trials):
+        ret_val.append(estimator(subsample(db, probe_type, subsample_size)))
+
+    return ret_val
+
+
+# Returns the test case name that clearly has higher RTT; otherwise, returns None
+def boxTest(params, test_cases, samples):
+    if len(test_cases) != 2:
+        # XXX: somehow generalize the box test to handle more than 2 cases
+        raise Exception()
+    dists = samples2Distributions(samples,'packet_rtt', test_cases) #XXX: field from params
+
+    tmp1,tmp2 = dists.items()
+    test_case1,dist1 = tmp1
+    test_case2,dist2 = tmp2
+    
+    d1_high = numpy.percentile(dist1, params['high'])
+    d2_low = numpy.percentile(dist2, params['low'])
+    if d1_high < d2_low:
+        return test_case2
+
+    d1_low = numpy.percentile(dist1, params['low'])
+    d2_high = numpy.percentile(dist2, params['high'])
+
+    if d2_high < d1_low:
+        return test_case1
+    
+    return None
+
+
+# Returns 1 if unusual_case is unusual in the expected direction
+#         0 if it isn't unusual
+#        -1 if it is unusual in the wrong direction
+def multiBoxTest(params, unusual_case, greater, samples):
+    #XXX: packet_rtt field from params
+    dists = samples2Distributions(samples, 'packet_rtt')
+    
+    uc = dists[unusual_case]
+    rest = []
+    for tc,d in dists.items():
+        if tc != unusual_case:
+            rest.extend(d)
+
+    uc_high = numpy.percentile(uc, params['high'])
+    rest_low = numpy.percentile(rest, params['low'])
+    if uc_high < rest_low:
+        if greater:
+            return -1
+        else:
+            return 1
+
+    uc_low = numpy.percentile(uc, params['low'])
+    rest_high = numpy.percentile(rest, params['high'])
+    if rest_high < uc_low:
+        if greater:
+            return 1
+        else:
+            return -1
+        
+    return 0
+
+
+# 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))
+
+    mh = midhinge(diffs, params['distance'])
+    if greater:
+        if mh > params['threshold']:
+            return 1
+        else:
+            return 0
+    else:
+        if mh < params['threshold']:
+            return 1
+        else:
+            return 0
+
+
+def rmse(expected, measurements):
+    s = sum([(expected-m)**2 for m in measurements])/len(measurements)
+    return math.sqrt(s)
+
+def nrmse(expected, measurements):
+    return rmse(expected, measurements)/(max(measurements)-min(measurements))
Index: /trunk/lib/nanownlib/storage.py
===================================================================
--- /trunk/lib/nanownlib/storage.py	(revision 4)
+++ /trunk/lib/nanownlib/storage.py	(revision 4)
@@ -0,0 +1,130 @@
+#-*- mode: Python;-*-
+
+import sys
+import os
+import uuid
+import threading
+import sqlite3
+
+def _newid():
+    return uuid.uuid4().hex
+
+
+class db(threading.local):
+    conn = None
+    cursor = None
+    def __init__(self, path):
+        exists = os.path.exists(path)
+        self.conn = sqlite3.connect(path)
+        self.conn.execute("PRAGMA foreign_keys = ON;")
+        self.conn.row_factory = sqlite3.Row
+        
+        if not exists:
+            self.conn.execute(
+                """CREATE TABLE meta (id BLOB PRIMARY KEY,
+                                      tcpts_mean REAL,
+                                      tcpts_stddev REAL,
+                                      tcpts_slopes TEXT)
+                """)
+
+            self.conn.execute(
+                """CREATE TABLE probes (id BLOB PRIMARY KEY,
+                                        sample INTEGER,
+                                        test_case TEXT,
+                                        type TEXT,
+                                        tc_order INTEGER,
+                                        time_of_day INTEGER,
+                                        local_port INTEGER,
+                                        reported INTEGER,
+                                        userspace_rtt INTEGER,
+                                        UNIQUE (sample, test_case))
+                """)
+
+            self.conn.execute(
+                """CREATE TABLE packets (id BLOB PRIMARY KEY,
+                                         probe_id REFERENCES probes(id) ON DELETE CASCADE,
+                                         sent INTEGER,
+                                         observed INTEGER,
+                                         tsval INTEGER,
+                                         payload_len INTEGER,
+                                         tcpseq INTEGER,
+                                         tcpack INTEGER)
+                """)
+
+            self.conn.execute(
+                """CREATE TABLE analysis (id BLOB PRIMARY KEY,
+                                          probe_id UNIQUE REFERENCES probes(id) ON DELETE CASCADE,
+                                          suspect TEXT,
+                                          packet_rtt INTEGER,
+                                          tsval_rtt INTEGER)
+                """)
+
+            self.conn.execute(
+                """CREATE TABLE trim_analysis (id BLOB PRIMARY KEY,
+                                               probe_id REFERENCES probes(id) ON DELETE CASCADE,
+                                               suspect TEXT,
+                                               packet_rtt INTEGER,
+                                               tsval_rtt INTEGER,
+                                               sent_trimmed INTEGER,
+                                               rcvd_trimmed INTEGER)
+                """)
+
+            self.conn.execute(
+                """CREATE TABLE classifier_results (id BLOB PRIMARY KEY,
+                                                    algorithm TEXT,
+                                                    params TEXT,
+                                                    sample_size INTEGER,
+                                                    num_trials INTEGER,
+                                                    trial_type TEXT,
+                                                    false_positives REAL,
+                                                    false_negatives REAL)
+                """)
+
+    def __del__(self):
+        if self.conn:
+            self.conn.commit()
+            self.conn.close()
+
+    def _insert(self, table, row):
+        rid = _newid()
+        keys = row.keys()
+        columns = ','.join(keys)
+        placeholders = ':'+', :'.join(keys)
+        query = "INSERT INTO %s (id,%s) VALUES ('%s',%s)" % (table, columns, rid, placeholders)
+        #print(row)
+        self.conn.execute(query, row)
+        return rid
+
+    def addMeta(self, meta):
+        ret_val = self._insert('meta', meta)
+        self.conn.commit()
+        return ret_val
+    
+    def addProbes(self, p):
+        return [self._insert('probes', row) for row in p]
+
+    def addPackets(self, pkts, window_size):
+        query = ("INSERT INTO packets (id,probe_id,sent,observed,tsval,payload_len,tcpseq,tcpack)"
+                 " VALUES(randomblob(16),"
+                 "(SELECT id FROM probes WHERE local_port=:local_port AND :observed>time_of_day"
+                 " AND :observed<time_of_day+userspace_rtt+%d" 
+                 " ORDER BY time_of_day ASC LIMIT 1),"
+                 ":sent,:observed,:tsval,:payload_len,:tcpseq,:tcpack)") % window_size
+        self.conn.execute("PRAGMA foreign_keys = OFF;")
+        self.conn.execute("CREATE INDEX IF NOT EXISTS probes_port ON probes (local_port)")
+        cursor = self.conn.cursor()
+        #print(query, list(pkts)[0:3])
+        cursor.executemany(query, pkts)
+        self.conn.commit()
+        self.conn.execute("PRAGMA foreign_keys = ON;")
+
+    def addAnalyses(self, analyses):
+        return [self._insert('analysis', row) for row in analyses]
+
+    def addTrimAnalyses(self, analyses):
+        return [self._insert('trim_analysis', row) for row in analyses]
+
+    def addClassifierResults(self, results):
+        ret_val = self._insert('classifier_results', results)
+        self.conn.commit()
+        return ret_val
Index: /trunk/src/compile.sh
===================================================================
--- /trunk/src/compile.sh	(revision 4)
+++ /trunk/src/compile.sh	(revision 4)
@@ -0,0 +1,3 @@
+#!/bin/sh
+
+gcc -ggdb -Wl,-z,relro,-z,now -fstack-protector-strong -Wformat -Werror=format-security -D_FORTIFY_SOURCE=2 csamp.c -lpcap -o ../bin/csamp
Index: /trunk/src/csamp.c
===================================================================
--- /trunk/src/csamp.c	(revision 4)
+++ /trunk/src/csamp.c	(revision 4)
@@ -0,0 +1,392 @@
+/*
+Portions are Copyright 2002 Tim Carstens. All rights reserved. 
+Redistribution and use, with or without modification, are permitted provided
+that the following conditions are met:
+
+- Redistribution must retain the above copyright notice and this list of 
+  conditions.
+
+- The name of Tim Carstens may not be used to endorse or promote products
+  derived from this document without specific prior written permission.
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <inttypes.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <errno.h>
+#include <netinet/in.h>
+#include <pcap.h>
+
+#define DEBUG 0
+
+/* ethernet headers are always exactly 14 bytes */
+/* XXX: ^^^ really? */
+#define SIZE_ETHERNET 14
+
+/* Ethernet addresses are 6 bytes */
+#define ETHER_ADDR_LEN	6
+
+/* Ethernet header */
+/*
+struct sniff_ethernet
+{
+  u_char ether_dhost[ETHER_ADDR_LEN]; / * Destination host address * /
+  u_char ether_shost[ETHER_ADDR_LEN]; / * Source host address * /
+  u_short ether_type; / * IP? ARP? RARP? etc * /
+};
+*/
+
+/* IP header */
+struct sniff_ip
+{
+  u_char ip_vhl;		/* version << 4 | header length >> 2 */
+  u_char ip_tos;		/* type of service */
+  u_short ip_len;		/* total length */
+  u_short ip_id;		/* identification */
+  u_short ip_off;		/* fragment offset field */
+#define IP_RF 0x8000		/* reserved fragment flag */
+#define IP_DF 0x4000		/* dont fragment flag */
+#define IP_MF 0x2000		/* more fragments flag */
+#define IP_OFFMASK 0x1fff	/* mask for fragmenting bits */
+  u_char ip_ttl;		/* time to live */
+  u_char ip_p;		/* protocol */
+  u_short ip_sum;		/* checksum */
+  struct in_addr ip_src,ip_dst; /* source and dest address */
+};
+#define IP_HL(ip)		(((ip)->ip_vhl) & 0x0f)
+#define IP_V(ip)		(((ip)->ip_vhl) >> 4)
+
+/* TCP header */
+typedef u_int tcp_seq;
+
+struct sniff_tcp
+{
+  u_short th_sport;	/* source port */
+  u_short th_dport;	/* destination port */
+  tcp_seq th_seq;		/* sequence number */
+  tcp_seq th_ack;		/* acknowledgement number */
+  u_char th_offx2;	/* data offset, rsvd */
+#define TH_OFF(th)	(((th)->th_offx2 & 0xf0) >> 4)
+  u_char th_flags;
+#define TH_FIN 0x01
+#define TH_SYN 0x02
+#define TH_RST 0x04
+#define TH_PUSH 0x08
+#define TH_ACK 0x10
+#define TH_URG 0x20
+#define TH_ECE 0x40
+#define TH_CWR 0x80
+#define TH_FLAGS (TH_FIN|TH_SYN|TH_RST|TH_ACK|TH_URG|TH_ECE|TH_CWR)
+  u_short th_win;		/* window */
+  u_short th_sum;		/* checksum */
+  u_short th_urp;		/* urgent pointer */
+};
+
+typedef struct {
+  uint8_t kind;
+  uint8_t size;
+} tcp_option_t;
+
+
+
+struct in_addr my_ip_raw;
+struct in_addr target_ip_raw;
+uint16_t target_port;
+
+FILE* output;
+
+struct useful_info
+{
+  uint32_t tcpseq;
+  uint32_t tcpack;
+  uint32_t tsval;
+  uint16_t local_port;
+  uint16_t payload_len;
+  uint8_t sent;
+};
+
+
+char* output_template = "{\"local_port\":%u,\"sent\":%u,\"payload_len\":%u,\"tcpseq\":%u,\"tcpack\":%u,\"observed\":%lu%06lu000,\"tsval\":%u}\n";
+uint8_t payloads_only = 1;
+
+
+pcap_t* create_listener(const char* dev, int snaplen, int promisc, int to_ms, char* errbuf)
+{
+  pcap_t* ret_val;
+  int status;
+  int tstcount, i, besttst;
+  int* tstypes;
+
+  ret_val = pcap_create(dev, errbuf);
+  if (ret_val == NULL)
+    return (NULL);
+  
+  status = pcap_set_snaplen(ret_val, snaplen);
+  if (status < 0)
+    goto fail;
+  
+  status = pcap_set_promisc(ret_val, promisc);
+  if (status < 0)
+    goto fail;
+  
+  status = pcap_set_timeout(ret_val, to_ms);
+  if (status < 0)
+    goto fail;
+
+
+  /* Try to select the best timestamp source for packet capture */
+  tstcount = pcap_list_tstamp_types(ret_val, &tstypes);
+  if(tstcount > 0)
+  {
+#if DEBUG
+    fprintf(stderr, "INFO: Available Packet Timers: ");
+#endif
+    besttst = -1;
+    for(i=0; i < tstcount; i++)
+    {
+#if DEBUG
+        fprintf(stderr, " %s", pcap_tstamp_type_val_to_name(tstypes[i]));
+#endif
+      switch (tstypes[i])
+      {
+      case PCAP_TSTAMP_HOST:
+        break;
+      case PCAP_TSTAMP_HOST_LOWPREC:
+        break;
+      case PCAP_TSTAMP_HOST_HIPREC:
+        {
+          if(besttst != PCAP_TSTAMP_ADAPTER_UNSYNCED
+             && besttst != PCAP_TSTAMP_ADAPTER)
+            besttst = PCAP_TSTAMP_HOST_HIPREC;
+          break;
+        }
+      case PCAP_TSTAMP_ADAPTER:
+        {
+          if(besttst != PCAP_TSTAMP_ADAPTER_UNSYNCED)
+            besttst = PCAP_TSTAMP_ADAPTER;
+          break;
+        }
+      case PCAP_TSTAMP_ADAPTER_UNSYNCED:
+          /*besttst = PCAP_TSTAMP_ADAPTER_UNSYNCED;*/
+          break;
+      default:
+#if DEBUG
+          fprintf(stderr, " unknown_type_%d", tstypes[i]);
+#endif
+        break;
+      }
+    }
+    pcap_free_tstamp_types(tstypes);
+#if DEBUG
+    fprintf(stderr, "\n");
+#endif
+    /*besttst = PCAP_TSTAMP_HOST;*/
+    if(besttst != -1)
+    {
+      fprintf(stderr, "INFO: Attempting to set the timestamp source to: %s\n",
+              pcap_tstamp_type_val_to_name(besttst));
+      if(pcap_set_tstamp_type(ret_val, besttst) != 0)
+        fprintf(stderr, "WARN: Failed to set preferred timestamp source.\n");
+    }
+  }
+  
+  /* Attempt to set nanosecond timestamp precision */    
+  if(pcap_set_tstamp_precision(ret_val, PCAP_TSTAMP_PRECISION_NANO) != 0)
+    fprintf(stderr, "INFO: Failed to set packet capture nanosecond precision.\n");
+  else
+    output_template = "{\"local_port\":%u,\"sent\":%u,\"payload_len\":%u,\"tcpseq\":%u,\"tcpack\":%u,\"observed\":%lu%09lu,\"tsval\":%u}\n";
+
+
+  status = pcap_activate(ret_val);
+  if (status < 0)
+    goto fail;
+  
+  return ret_val;
+
+ fail:
+  pcap_close(ret_val);
+  return NULL;
+}
+
+
+static int extract_packet_fields(const struct pcap_pkthdr *header,
+                                 const u_char* packet,
+                                 struct useful_info* fields)
+{
+  const struct sniff_ip* ip;               /* The IP header */
+  const struct sniff_tcp* tcp;             /* The TCP header */
+  u_int iphdr_size;
+  u_int tcphdr_size;
+  uint16_t ip_len;
+  uint16_t tsval;
+  uint8_t* opt;
+  tcp_option_t* _opt;
+  
+  /*printf("snaplen: %u\n", header->caplen);*/
+  
+  ip = (struct sniff_ip*)(packet + SIZE_ETHERNET);
+  iphdr_size = IP_HL(ip)*4;
+  if (iphdr_size < 20)
+    return 0;
+  ip_len = ntohs(ip->ip_len);
+  
+  tcp = (struct sniff_tcp*)(packet + SIZE_ETHERNET + iphdr_size);
+  tcphdr_size = TH_OFF(tcp)*4;
+  if (tcphdr_size < 20)
+    return 0;
+
+  /*fprintf(stderr, "%d\n", (ip_len - iphdr_size - tcphdr_size));*/
+  if ((ip_len - iphdr_size - tcphdr_size < 0)
+      || (payloads_only && (ip_len - iphdr_size - tcphdr_size == 0)))
+  { return 0; }
+
+  fields->tcpseq = ntohl(tcp->th_seq);
+  fields->tcpack = ntohl(tcp->th_ack);
+  fields->payload_len = ip_len - iphdr_size - tcphdr_size;
+  if(ip->ip_src.s_addr == target_ip_raw.s_addr && ntohs(tcp->th_sport) == target_port)
+  {
+    fields->sent = 0;
+    fields->local_port = ntohs(tcp->th_dport);
+  }
+  else
+  {
+    fields->sent = 1;
+    fields->local_port = ntohs(tcp->th_sport);
+  }
+  /*printf("src: %lX / dst: %lX\n", ip->ip_src, ip->ip_dst);*/
+
+  fields->tsval = 0;
+  if(tcphdr_size > 20)
+  {
+    opt = (uint8_t*)(packet + SIZE_ETHERNET + iphdr_size + 20);
+    while((*opt != 0) && (opt - packet) < header->caplen)
+    {
+      _opt = (tcp_option_t*)opt;
+      if(_opt->kind == 1) /* NOP */
+      {
+        ++opt;  // NOP is one byte;
+        continue;
+      }
+      
+      if(_opt->kind == 8) /* timestamp */
+      {
+        fields->tsval = ntohl(*(uint32_t*)(opt+2));
+        break;
+      }
+      
+      opt += _opt->size;
+    }
+  }
+
+  return 1;
+}
+
+
+void process_packet(u_char *args, const struct pcap_pkthdr *header, const u_char *packet)
+{
+  struct useful_info fields;
+
+  if(extract_packet_fields(header, packet, &fields))
+  {
+    fprintf(output, output_template, fields.local_port, fields.sent,
+            fields.payload_len, fields.tcpseq, fields.tcpack,
+            header->ts.tv_sec, header->ts.tv_usec, fields.tsval);
+    fflush(output);
+  }
+}
+
+
+
+int main(int argc, char** argv)
+{
+  char* dev;                      /* The device to sniff on */
+  char* my_ip;
+  char* target_ip;
+  pcap_t* handle;                 /* Session handle */
+  char bpf[255];                  /* The filter expression (length of 168 should be enough)*/
+  char errbuf[PCAP_ERRBUF_SIZE];  /* Error string */
+  struct bpf_program fp;          /* The compiled filter */
+  bpf_u_int32 mask;               /* Our netmask */
+  bpf_u_int32 net;                /* Our IP */
+  struct pcap_pkthdr header;      /* The header that pcap gives us */
+  const u_char *packet;           /* The actual packet */
+    
+  if(argc < 6)
+  {
+    fprintf(stderr, "USAGE:\n  %s {interface} {my_ip} {target_ip} {target_port} {output_file} [{payloads_only?}]\n", argv[0]);
+    return 1;
+  }
+
+  dev = argv[1];
+  my_ip = argv[2];
+  target_ip = argv[3];
+  target_port = atoi(argv[4]);
+  if(argc == 7 && argv[6][0] == '0')
+    payloads_only = 0;
+  
+  if(!(output = fopen(argv[5], "w+")))
+  {
+    fprintf(stderr, "ERROR: could not open output file due to: %s\n", strerror(errno));
+    return 2;
+  }
+  
+  snprintf(bpf, 255, "(src host %s and dst host %s and tcp and src port %u) or (dst host %s and src host %s and tcp and dst port %u)",
+           target_ip, my_ip, target_port, target_ip, my_ip, target_port);
+
+
+  if(!inet_aton(my_ip, &my_ip_raw))
+  {
+    fprintf(stderr, "Couldn't parse my_ip.\n");
+    return 1;
+  }
+
+  if(!inet_aton(target_ip, &target_ip_raw))
+  {
+    fprintf(stderr, "Couldn't parse target_ip.\n");
+    return 1;
+  }
+
+  memset(errbuf, 0, PCAP_ERRBUF_SIZE);
+
+  /* Find the properties for the device */
+  if (pcap_lookupnet(dev, &net, &mask, errbuf) == -1)
+  {
+    fprintf(stderr, "Couldn't get netmask for device %s: %s\n", dev, errbuf);
+    net = 0;
+    mask = 0;
+  }
+  
+  /* Open the session in promiscuous mode */
+  /* XXX: does to_ms timeout (param 4) matter? */
+  handle = create_listener(dev, BUFSIZ, 0, 1000, errbuf);
+  if (handle == NULL)
+  {
+    fprintf(stderr, "Couldn't open device %s: %s\n", dev, errbuf);
+    return 2;
+  }
+
+  /* Compile and apply the filter */
+  if (pcap_compile(handle, &fp, bpf, 0, net) == -1)
+  {
+    fprintf(stderr, "Couldn't parse filter %s: %s\n", bpf, pcap_geterr(handle));
+    return 2;
+  }
+  
+  if (pcap_setfilter(handle, &fp) == -1)
+  {
+    fprintf(stderr, "Couldn't install filter %s: %s\n", bpf, pcap_geterr(handle));
+    return 2;
+  }
+
+  /* XXX: report errors */
+  pcap_loop(handle, 0, process_packet, NULL);
+  
+  /* And close the session */
+  pcap_close(handle);
+
+  return 0;
+}
