Changeset 10

Timestamp:

07/13/15 19:16:30 (9 years ago)

Author:

tim

Message:

.

Location:

trunk

Files:

• bin/analyze_packets (modified) (1 diff)
• bin/graph (modified) (2 diffs)
• bin/train (modified) (12 diffs)
• lib/nanownlib/__init__.py (modified) (13 diffs)
• lib/nanownlib/stats.py (modified) (8 diffs)
• lib/nanownlib/storage.py (modified) (7 diffs)

trunk/bin/analyze_packets

@@ -44 +44 @@
 
 start = time.time()
+import cProfile
+#cProfile.run('num_probes = analyzeProbes(db)')
 num_probes = analyzeProbes(db)
 end = time.time()

trunk/bin/graph

@@ -10 +10 @@
 import json
 
+import numpy
 import matplotlib.mlab as mlab
 import matplotlib.pyplot as plt
@@ -85 +86 @@
 
 print('packet_rtt diff median: %f' % statistics.median(diffs))
-print('packet_rtt diff midhinge: %f' % midhinge(diffs))
+print('packet_rtt diff midhinge: %f' % midsummary(diffs))
 print('packet_rtt diff trimean: %f' % trimean(diffs))
+print('packet_rtt diff quadsummary: %f' % quadsummary(diffs))
+print('packet_rtt diff ubersummary: %f' % ubersummary(diffs))
 print('packet_rtt diff MAD: %f' % mad(diffs))
 print('reported diff trimean: %f' % trimean(reported_diffs))
+print('reported diff quadsummary: %f' % quadsummary(reported_diffs))
+print('reported diff ubersummary: %f' % ubersummary(reported_diffs))
 print('reported diff MAD: %f' % mad(reported_diffs))
+
+import cProfile
+kresults = kfilter({},diffs)
+#print('packet_rtt diff kfilter: ', numpy.mean(kresults['est']), kresults['var'])
+print('packet_rtt diff kfilter: ', kresults['est'][-1], kresults['var'][-1])
+kresults = kfilter({},reported_diffs)
+#print('reported diff kfilter: ', numpy.mean(kresults['est']), kresults['var'][-1])
+print('reported diff kfilter: ', kresults['est'][-1], kresults['var'][-1])
 
 

trunk/bin/train

@@ -26 +26 @@
 
 from nanownlib import *
+from nanownlib.stats import *
+from nanownlib.parallel import WorkerThreads
 import nanownlib.storage
-from nanownlib.stats import boxTest,multiBoxTest,subsample,bootstrap,bootstrap2,trimean,midhinge,midhingeTest,samples2Distributions,samples2MeanDiffs
-from nanownlib.parallel import WorkerThreads
+
 
 
@@ -35 +36 @@
 #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('--retrain', action='append', default=[], help='Force a classifier to be retrained.  May be specified multiple times.')
+parser.add_argument('--retest', action='append', default=[], help='Force a classifier to be retested.  May be specified multiple times.')
 parser.add_argument('session_data', default=None,
                     help='Database file storing session information')
 options = parser.parse_args()
 
-            
-
-def trainBoxTest(db, unusual_case, greater, subseries_size):
-
+
+def trainBoxTest(db, unusual_case, greater, num_observations):
+    db.resetOffsets()
+    
     def trainAux(low,high,num_trials):
         estimator = functools.partial(multiBoxTest, {'low':low, 'high':high}, 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)
+        estimates = bootstrap3(estimator, db, 'train', unusual_case, num_observations, num_trials)
+        null_estimates = bootstrap3(estimator, db, 'train_null', unusual_case, num_observations, num_trials)
 
         bad_estimates = len([e for e in estimates if e != 1])
@@ -116 +119 @@
     
     num_trials = 500
-    widths = [good_width+(x/100.0) for x in range(-60,75,5) if good_width+(x/100.0) > 0.0]
+    widths = [good_width+(x/100.0) for x in range(-70,75,5) if good_width+(x/100.0) > 0.0]
     performance = []
     for width in widths:
@@ -124 +127 @@
         job_id,errors = wt.resultq.get()
         fp,fn = errors
-        performance.append(((fp+fn)/2.0, job_id, fn, fp))
+        #performance.append(((fp+fn)/2.0, job_id, fn, fp))
+        performance.append((abs(fp-fn), job_id, fn, fp))
     performance.sort()
     #pprint.pprint(performance)
@@ -133 +137 @@
     wt.stop()
     params = json.dumps({"low":best_low,"high":best_low+best_width})
-    return {'algorithm':"boxtest",
+    return {'trial_type':"train",
+            'num_observations':num_observations,
+            'num_trials':num_trials,
             'params':params,
-            'sample_size':subseries_size,
-            'num_trials':num_trials,
-            'trial_type':"train",
             'false_positives':performance[0][3],
             'false_negatives':performance[0][2]}
 
 
-def trainMidhinge(db, unusual_case, greater, subseries_size):
-
+def trainSummary(summaryFunc, db, unusual_case, greater, num_observations):
+    db.resetOffsets()
+    stest = functools.partial(summaryTest, summaryFunc)
+    
     def trainAux(distance, threshold, 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)
+        estimator = functools.partial(stest, {'distance':distance,'threshold':threshold}, greater)
+        estimates = bootstrap3(estimator, db, 'train', unusual_case, num_observations, num_trials)
+        null_estimates = bootstrap3(estimator, db, 'train_null', unusual_case, num_observations, num_trials)
 
         bad_estimates = len([e for e in estimates if e != 1])
@@ -158 +163 @@
     #determine expected delta based on differences
     mean_diffs = [s['unusual_case']-s['other_cases'] for s in db.subseries('train', unusual_case)]
-    threshold = trimean(mean_diffs)/2.0
+    threshold = summaryFunc(mean_diffs)/2.0
     #print("init_threshold:", threshold)
     
@@ -181 +186 @@
     num_trials = 500
     performance = []
-    for t in range(50,154,4):
+    for t in range(80,122,2):
         wt.addJob(threshold*(t/100.0), (good_distance,threshold*(t/100.0),num_trials))
     wt.wait()
@@ -187 +192 @@
         job_id,errors = wt.resultq.get()
         fp,fn = errors
-        performance.append(((fp+fn)/2.0, job_id, fn, fp))
+        #performance.append(((fp+fn)/2.0, job_id, fn, fp))
+        performance.append((abs(fp-fn), job_id, fn, fp))
     performance.sort()
     #pprint.pprint(performance)
@@ -217 +223 @@
         job_id,errors = wt.resultq.get()
         fp,fn = errors
-        performance.append(((fp+fn)/2.0, job_id, fn, fp))
+        #performance.append(((fp+fn)/2.0, job_id, fn, fp))
+        performance.append((abs(fp-fn), job_id, fn, fp))
     performance.sort()
     #pprint.pprint(performance)
@@ -225 +232 @@
     wt.stop()
     params = json.dumps({'distance':best_distance,'threshold':best_threshold})
-    return {'algorithm':"midhinge",
+    return {'trial_type':"train",
+            'num_observations':num_observations,
+            'num_trials':num_trials,
             'params':params,
-            'sample_size':subseries_size,
-            'num_trials':num_trials,
-            'trial_type':"train",
             'false_positives':performance[0][3],
             'false_negatives':performance[0][2]}
 
 
-classifiers = {'boxtest':{'train':trainBoxTest, 'test':multiBoxTest},
-               'midhinge':{'train':trainMidhinge, 'test':midhinge}}
+def trainKalman(db, unusual_case, greater, num_observations):
+    db.resetOffsets()
+
+    def trainAux(params, num_trials):
+        estimator = functools.partial(kalmanTest, params, greater)
+        estimates = bootstrap3(estimator, db, 'train', unusual_case, num_observations, num_trials)
+        null_estimates = bootstrap3(estimator, db, 'train_null', unusual_case, num_observations, 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
+    
+    mean_diffs = [s['unusual_case']-s['other_cases'] for s in db.subseries('train', unusual_case)]
+    good_threshold = kfilter({},mean_diffs)['est'][-1]/2.0
+
+    wt = WorkerThreads(2, trainAux)
+    num_trials = 200
+    performance = []
+    for t in range(90,111):
+        params = {'threshold':good_threshold*(t/100.0)}
+        wt.addJob(good_threshold*(t/100.0), (params,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.append((abs(fp-fn), job_id, fn, fp))
+    performance.sort()
+    #pprint.pprint(performance)
+    best_threshold = performance[0][1]
+    #print("best_threshold:", best_threshold)
+    params = {'threshold':best_threshold}
+
+    wt.stop()
+    
+    return {'trial_type':"train",
+            'num_observations':num_observations,
+            'num_trials':num_trials,
+            'params':json.dumps(params),
+            'false_positives':performance[0][3],
+            'false_negatives':performance[0][2]}
+
+    
+    #determine expected delta based on differences
+classifiers = {'boxtest':{'train':trainBoxTest, 'test':multiBoxTest, 'train_results':[]},
+               'midsummary':{'train':functools.partial(trainSummary, midsummary), 'test':midsummaryTest, 'train_results':[]},
+               #'ubersummary':{'train':functools.partial(trainSummary, ubersummary), 'test':ubersummaryTest, 'train_results':[]},
+               'quadsummary':{'train':functools.partial(trainSummary, quadsummary), 'test':quadsummaryTest, 'train_results':[]},
+               'kalman':{'train':trainKalman, 'test':kalmanTest, 'train_results':[]},
+               #'_trimean':{'train':None, 'test':trimeanTest, 'train_results':[]},
+              }
 
 
@@ -242 +300 @@
 import cProfile
 
-def trainClassifier(db, unusual_case, greater, trainer):
+def trainClassifier(db, unusual_case, greater, classifier, retrain=False):
+    if retrain:
+        print("Dropping stored training results...")
+        db.deleteClassifierResults(classifier, 'train')
+    
+    trainer = classifiers[classifier]['train']
     threshold = 5.0 # in percent
-    size = 4000
+    num_obs = 1000
+    max_obs = int(db.populationSize('train')/5)
     result = None
-    while size < db.populationSize('train')/5:
-        size = min(size*2, int(db.populationSize('train')/5))
-        result = trainer(db,unusual_case,greater,size)
+    while num_obs < max_obs:
+        num_obs = min(int(num_obs*1.5), max_obs)
+        result = db.fetchClassifierResult(classifier, 'train', num_obs)
+        if result != None:
+            train_time = "(stored)"
+        else:
+            start = time.time()
+            result = trainer(db,unusual_case,greater,num_obs)
+            result['classifier'] = classifier
+            train_time = "%f" % (time.time()-start)
+            
         error = statistics.mean([result['false_positives'],result['false_negatives']])
-        print("subseries size: %d | error: %f | false_positives: %f | false_negatives: %f"
-              % (size,error,result['false_positives'],result['false_negatives']))
+        print("number of observations: %d | error: %f | false_positives: %f | false_negatives: %f | train time: %s | params: %s"
+              % (num_obs, error, result['false_positives'],result['false_negatives'], train_time, result['params']))
+        db.addClassifierResults(result)
+        classifiers[classifier]['train_results'].append(result)
+
         if error < threshold:
             break
-    if result != None:
-        db.addClassifierResults(result)
 
     return result
+
+
+
+def testClassifier(db, unusual_case, greater, classifier, retest=False):
+    target_error = 5.0 # in percent
+    num_trials = 1000
+    max_obs = int(db.populationSize('test')/5)
+
+    tester = classifiers[classifier]['test']
+    
+    def testAux(params, num_trials, num_observations):
+        estimator = functools.partial(tester, params, greater)
+        estimates = bootstrap3(estimator, db, 'test', unusual_case, num_observations, num_trials)
+        null_estimates = bootstrap3(estimator, db, 'train_null', unusual_case, num_observations, 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
+        print("testAux:", num_observations, false_positives, false_negatives, params)
+        return false_positives,false_negatives
+
+
+    if retest:
+        print("Dropping stored test results...")
+        db.deleteClassifierResults(classifier, 'test')
+
+
+    test_results = []
+    lte = math.log(target_error/100.0)
+    for tr in classifiers[classifier]['train_results']:
+        db.resetOffsets()
+        params = json.loads(tr['params'])
+        num_obs = tr['num_observations']
+    
+        print("initial test")
+        fp,fn = testAux(params, num_trials, num_obs)
+        error = (fp+fn)/2.0
+        print("walking up")
+        while (error > target_error) and (num_obs < max_obs):
+            increase_factor = 1.5 * lte/math.log(error/100.0) # don't ask how I came up with this
+            #print("increase_factor:", increase_factor)
+            num_obs = min(int(increase_factor*num_obs), max_obs)
+            fp,fn = testAux(params, num_trials, num_obs)
+            error = (fp+fn)/2.0
+
+        print("walking down")
+        while (num_obs > 0):
+            current_best = (num_obs,error,params,fp,fn)
+            num_obs = int(0.95*num_obs)
+            fp,fn = testAux(params, num_trials, num_obs)
+            error = (fp+fn)/2.0
+            if error > target_error:
+                break
+        
+        test_results.append(current_best)
+
+    test_results.sort()
+    best_obs,error,best_params,fp,fn = test_results[0]
+    
+    return {'classifier':classifier,
+            'trial_type':"test",
+            'num_observations':best_obs,
+            'num_trials':num_trials,
+            'params':best_params,
+            'false_positives':fp,
+            'false_negatives':fn}
 
 
@@ -268 +409 @@
 print(":", end-start)
 
-
-for c,funcs in classifiers.items():
+for c in sorted(classifiers.keys()):
+    if classifiers[c]['train'] == None:
+        continue
     start = time.time()
     print("Training %s..." % c)
-    result = trainClassifier(db, unusual_case, greater, funcs['train'])
+    result = trainClassifier(db, unusual_case, greater, c, c in options.retrain)
     print("%s result:" % c)
     pprint.pprint(result)
     print("completed in:", time.time()-start)
 
-sys.exit(0)
-
-start = time.time()
-results = trainBoxTest(db, unusual_case, greater, 6000)
-#db.addClassifierResults(results)
-print("multi box test result:")
-pprint.pprint(results)
-print(":", time.time()-start)
+db.clearCache()
+
+for c in sorted(classifiers.keys()):
+    start = time.time()
+    print("Testing %s..." % c)
+    result = testClassifier(db, unusual_case, greater, c, c in options.retest)
+    print("%s result:" % c)
+    pprint.pprint(result)
+    classifiers[c]['test_error'] = (result['false_positives']+result['false_negatives'])/2.0
+    print("completed in:", time.time()-start)

trunk/lib/nanownlib/__init__.py

@@ -4 +4 @@
 import sys
 import time
+import traceback
 import random
 import argparse
@@ -209 +210 @@
 def removeDuplicatePackets(packets):
     #return packets
-    suspect = None
+    suspect = ''
     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']):
+        if (key not in seen):
             seen[key] = p
-    
-    if len(seen) < len(packets):
-        suspect = 'd'
-        #sys.stderr.write("INFO: removed %d duplicate packets.\n" % (len(packets) - len(seen)))
+            continue
+        if p['sent']==1 and (seen[key]['observed'] > p['observed']): #earliest sent
+            seen[key] = p
+            suspect += 's'
+            continue 
+        if p['sent']==0 and (seen[key]['observed'] > p['observed']): #earliest rcvd
+            seen[key] = p
+            suspect += 'r'
+            continue
+    
+    #if len(seen) < len(packets):
+    #   sys.stderr.write("INFO: removed %d duplicate packets.\n" % (len(packets) - len(seen)))
 
     return suspect,seen.values()
@@ -230 +236 @@
     suspect,packets = removeDuplicatePackets(packets)
 
-    #sort_key = lambda d: (d['tcpseq'],d['tcpack'])
-    sort_key = lambda d: (d['observed'],d['tcpseq'])
+    sort_key = lambda d: (d['tcpseq'],d['observed'])
     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]
+    alt_key = lambda d: (d['observed'],d['tcpseq'])
+    rcvd_alt = sorted((p for p in packets if p['sent']==0 and p['payload_len']>0), key=alt_key)
+
+    s_off = trim_sent
+    if s_off >= len(sent):
+        s_off = -1
+    last_sent = sent[s_off]
+
+    r_off = len(rcvd) - trim_rcvd - 1
+    if r_off <= 0:
+        r_off = 0
+    last_rcvd = rcvd[r_off]
+    if last_rcvd != rcvd_alt[r_off]:
+        suspect += 'R'
     
     packet_rtt = last_rcvd['observed'] - last_sent['observed']
@@ -272 +282 @@
     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
+                         WHERE sent_trimmed=:strim AND rcvd_trimmed=:rtrim AND trim_analysis.probe_id=probes.id AND probes.test_case!=:unusual_case AND sample=u.s AND probes.type in ('train','test'))
+      FROM (SELECT probes.sample s,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') AND 1 NOT IN (select 1 from probes p,trim_analysis t WHERE p.sample=s AND t.probe_id=p.id AND t.suspect LIKE '%R%')) u
+    """
+    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.s AND probes.type in ('train','test'))
+      FROM (SELECT probes.sample s,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
     """
 
@@ -280 +295 @@
     differences = [row[0] for row in cursor]
     
-    return trimean(differences),mad(differences)
+    return ubersummary(differences),mad(differences)
 
 
@@ -287 +302 @@
     db.conn.execute("CREATE INDEX IF NOT EXISTS packets_probe ON packets (probe_id)")
     pcursor = db.conn.cursor()
-    kcursor = db.conn.cursor()
+    db.conn.commit()
 
     pcursor.execute("SELECT tcpts_mean FROM meta")
@@ -297 +312 @@
     pcursor.execute("DELETE FROM trim_analysis")
     db.conn.commit()
+
+    def loadPackets(db):
+        cursor = db.conn.cursor()
+        cursor.execute("SELECT * FROM packets ORDER BY probe_id")
+
+        probe_id = None
+        entry = []
+        ret_val = []
+        for p in cursor:
+            if probe_id == None:
+                probe_id = p['probe_id']
+            if p['probe_id'] != probe_id:
+                ret_val.append((probe_id,entry))
+                probe_id = p['probe_id']
+                entry = []
+            entry.append(dict(p))
+        ret_val.append((probe_id,entry))
+        return ret_val
+    
+    start = time.time()
+    packet_cache = loadPackets(db)
+    print("packets loaded in: %f" % (time.time()-start))
     
     count = 0
     sent_tally = []
     rcvd_tally = []
-    for pid, in pcursor.execute("SELECT id FROM probes"):
-        kcursor.execute("SELECT * FROM packets WHERE probe_id=?", (pid,))
+    for probe_id,packets in packet_cache:
         try:
-            analysis,s,r = analyzePackets(kcursor.fetchall(), timestamp_precision)
-            analysis['probe_id'] = pid
+            analysis,s,r = analyzePackets(packets, timestamp_precision)
+            analysis['probe_id'] = probe_id
             sent_tally.append(s)
             rcvd_tally.append(r)
+            db.addTrimAnalyses([analysis])
         except Exception as e:
-            print(e)
-            sys.stderr.write("WARN: couldn't find enough packets for probe_id=%s\n" % pid)
+            traceback.print_exc()
+            sys.stderr.write("WARN: couldn't find enough packets for probe_id=%s\n" % probe_id)
         
         #print(pid,analysis)
-        db.addTrimAnalyses([analysis])
         count += 1
     db.conn.commit()
@@ -326 +362 @@
             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,))
+            for probe_id,packets in packet_cache:
                 try:
-                    analysis,s,r = analyzePackets(kcursor.fetchall(), timestamp_precision, strim, rtrim)
-                    analysis['probe_id'] = pid
+                    analysis,s,r = analyzePackets(packets, timestamp_precision, strim, rtrim)
+                    analysis['probe_id'] = probe_id
                 except Exception as e:
                     print(e)
@@ -336 +371 @@
                     
                 db.addTrimAnalyses([analysis])
-            db.conn.commit()
+    db.conn.commit()
 
     # Populate analysis table so findUnusualTestCase can give us a starting point
@@ -359 +394 @@
     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:
+        if delta*good_delta > 0.0 and (abs(good_delta) - abs(delta)) < abs(delta_margin*good_delta) and mad < good_mad:
             best_strim = strim
         else:
@@ -367 +402 @@
     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:
+        if delta*good_delta > 0.0 and (abs(good_delta) - abs(delta)) < abs(delta_margin*good_delta) and mad < good_mad:
             best_rtrim = rtrim
         else:
@@ -425 +460 @@
     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])
+    global_tm = quadsummary([row['packet_rtt'] for row in cursor])
 
     tm_abs = []
@@ -432 +467 @@
     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_map[tc] = quadsummary([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])
+    remaining_tm = quadsummary([row['packet_rtt'] for row in cursor])
 
     ret_val = (tc, tm_map[tc]-remaining_tm)

trunk/lib/nanownlib/stats.py

@@ -2 +2 @@
 import sys
 import os
+import functools
 import math
 import statistics
@@ -133 +134 @@
 
 
-def midhinge(values, distance=25):
-    return (numpy.percentile(values, 50-distance) + numpy.percentile(values, 50+distance))/2.0
+def midsummary(values, distance=25):
+    #return (numpy.percentile(values, 50-distance) + numpy.percentile(values, 50+distance))/2.0
+    l,h = numpy.percentile(values, (50-distance,50+distance))
+    return (l+h)/2.0
 
 def trimean(values, distance=25):
-    return (midhinge(values, distance) + statistics.median(values))/2
-
+    return (midsummary(values, distance) + statistics.median(values))/2
+
+def ubersummary(values, distance=25):
+    left2 = 50-distance
+    left1 = left2/2.0
+    left3 = (left2+50)/2.0
+    right2 = 50+distance
+    right3 = (right2+50)/2.0
+    right1 = (right2+100)/2.0
+    l1,l2,l3,r3,r2,r1 = numpy.percentile(values, (left1,left2,left3,right3,right2,right1))
+    #print(left1,left2,left3,50,right3,right2,right1)
+    #print(l1,l2,l3,m,r3,r2,r1)
+    return (l1+l2*4+l3+r3+r2*4+r1)/12.0
+    #return statistics.mean((l1,l2,l3,m,r3,r2,r1))
+
+def quadsummary(values, distance=25):
+    left2 = 50-distance
+    left1 = left2/2.0
+    right2 = 50+distance
+    right1 = (right2+100)/2.0
+    l1,l2,r2,r1 = numpy.percentile(values, (left1,left2,right2,right1))
+    #print(left1,left2,left3,50,right3,right2,right1)
+    #print(l1,l2,l3,m,r3,r2,r1)
+    return (l1+l2+r2+r1)/4.0
+    #return statistics.mean((l1,l2,l3,m,r3,r2,r1))
+
+def quadsummary(values, distance=25):
+    left1 = 50-distance
+    left2 = (left1+50)/2.0
+    right1 = 50+distance
+    right2 = (right1+50)/2.0
+    l1,l2,r2,r1 = numpy.percentile(values, (left1,left2,right2,right1))
+    #print(left1,left2,left3,50,right3,right2,right1)
+    #print(l1,l2,l3,m,r3,r2,r1)
+    return (l1+l2+r2+r1)/4.0
+    #return statistics.mean((l1,l2,l3,m,r3,r2,r1))
+
+    
 def weightedMean(derived, weights):
     normalizer = sum(weights.values())/len(weights)
@@ -170 +209 @@
 
 
-def estimateMidhinge(derived):
-    return midhinge([(d['long']-d['short']) for d in derived.values()])
+def estimateMidsummary(derived):
+    return midsummary([(d['long']-d['short']) for d in derived.values()])
 
 
@@ -348 +387 @@
     rest = [s['other_cases'] for s in samples]
     
-    uc_high = numpy.percentile(uc, params['high'])
-    rest_low = numpy.percentile(rest, params['low'])
+    uc_high,uc_low = numpy.percentile(uc, (params['high'],params['low']))
+    rest_high,rest_low = numpy.percentile(rest, (params['high'],params['low']))
     if uc_high < rest_low:
         if greater:
@@ -356 +395 @@
             return 1
 
-    uc_low = numpy.percentile(uc, params['low'])
-    rest_high = numpy.percentile(rest, params['high'])
     if rest_high < uc_low:
         if greater:
@@ -369 +406 @@
 # Returns 1 if unusual_case is unusual in the expected direction
 #         0 otherwise
-def midhingeTest(params, greater, samples):
+def summaryTest(f, params, greater, samples):
     diffs = [s['unusual_case']-s['other_cases'] for s in samples]
 
-    mh = midhinge(diffs, params['distance'])
-    #mh = trimean(diffs, params['distance'])
+    mh = f(diffs, params['distance'])
     if greater:
         if mh > params['threshold']:
@@ -385 +421 @@
             return 0
 
+midsummaryTest = functools.partial(summaryTest, midsummary)
+trimeanTest = functools.partial(summaryTest, trimean)
+ubersummaryTest = functools.partial(summaryTest, ubersummary)
+quadsummaryTest = functools.partial(summaryTest, quadsummary)
 
 def rmse(expected, measurements):
@@ -392 +432 @@
 def nrmse(expected, measurements):
     return rmse(expected, measurements)/(max(measurements)-min(measurements))
+
+
+class KalmanFilter1D:
+    def __init__(self, x0, P, R, Q):
+        self.x = x0
+        self.P = P
+        self.R = R
+        self.Q = Q
+
+    def update(self, z):
+        self.x = (self.P * z + self.x * self.R) / (self.P + self.R)
+        self.P = 1. / (1./self.P + 1./self.R)
+
+    def predict(self, u=0.0):
+        self.x += u
+        self.P += self.Q
+
+
+def kfilter(params, observations):
+    x = numpy.array(observations)
+    movement = 0
+    est = []    
+    var = []
+    kf = KalmanFilter1D(x0 = quadsummary(x), # initial state
+                        #P  = 10000,          # initial variance
+                        P  = 10,          # initial variance
+                        R  = numpy.std(x),   # msensor noise
+                        Q  = 0)              # movement noise
+    for round in range(1):
+        for d in x:
+            kf.predict(movement)
+            kf.update(d)
+            est.append(kf.x)
+            var.append(kf.P)
+
+    return({'est':est, 'var':var})
+
+
+def kalmanTest(params, greater, samples):
+    diffs = [s['unusual_case']-s['other_cases'] for s in samples]
+
+    m = kfilter(params, diffs)['est'][-1]
+    if greater:
+        if m > params['threshold']:
+            return 1
+        else:
+            return 0
+    else:
+        if m < params['threshold']:
+            return 1
+        else:
+            return 0
+
+
+def kalmanTest2(params, greater, samples):
+    diffs = [s['unusual_case']-s['other_cases'] for s in samples]
+
+    estimates = []
+    size = 500
+    for i in range(100):
+        off = random.randrange(0,len(diffs))
+        sub = diffs[off:size]
+        if len(sub) < size:
+            sub += diffs[0:size-len(sub)]
+        estimates.append(kfilter(params, sub)['est'][-1])
+            
+    m = quadsummary(estimates)
+    if greater:
+        if m > params['threshold']:
+            return 1
+        else:
+            return 0
+    else:
+        if m < params['threshold']:
+            return 1
+        else:
+            return 0
+

trunk/lib/nanownlib/storage.py

@@ -4 +4 @@
 import os
 import uuid
+import random
 import threading
 import sqlite3
 
 import numpy
+# Don't trust numpy's seeding
+numpy.random.seed(random.SystemRandom().randint(0,2**32-1))
 
 def _newid():
@@ -18 +21 @@
     _population_sizes = None
     _population_cache = None
+    _offset_cache = None
+    _cur_offsets = None
     
     def __init__(self, path):
@@ -26 +31 @@
         self._population_sizes = {}
         self._population_cache = {}
+        self._offset_cache = {}
+        self._cur_offsets = {}
         
         if not exists:
@@ -79 +86 @@
             self.conn.execute(
                 """CREATE TABLE classifier_results (id BLOB PRIMARY KEY,
-                                                    algorithm TEXT,
+                                                    classifier TEXT,
+                                                    trial_type TEXT,
+                                                    num_observations INTEGER,
+                                                    num_trials INTEGER,
                                                     params TEXT,
-                                                    sample_size INTEGER,
-                                                    num_trials INTEGER,
-                                                    trial_type TEXT,
                                                     false_positives REAL,
                                                     false_negatives REAL)
@@ -109 +116 @@
 
     def subseries(self, probe_type, unusual_case, size=None, offset=None, field='packet_rtt'):
-        if (probe_type,unusual_case,field) not in self._population_cache:
+        cache_key = (probe_type,unusual_case,field)
+        
+        if cache_key not in self._population_cache:
             query="""
             SELECT %(field)s AS unusual_case,
@@ -121 +130 @@
             cursor = self.conn.cursor()
             cursor.execute(query, params)
-            self._population_cache[(probe_type,unusual_case,field)] = [dict(row) for row in cursor.fetchall()]
-
-        population = self._population_cache[(probe_type,unusual_case,field)]
+            p = [dict(row) for row in cursor.fetchall()]
+            self._population_cache[cache_key] = p
+            self._offset_cache[cache_key] = tuple(numpy.random.random_integers(0,len(p)-1, len(p)/5))
+            self._cur_offsets[cache_key] = 0
+
+        population = self._population_cache[cache_key]
 
         if size == None or size > len(population):
             size = len(population)
         if offset == None or offset >= len(population) or offset < 0:
-            offset = numpy.random.random_integers(0,len(population)-1)
-
+            offset = self._offset_cache[cache_key][self._cur_offsets[cache_key]]
+            self._cur_offsets[cache_key] = (offset + 1) % len(self._offset_cache[cache_key])
+        
         try:
-            ret_val = population[offset:offset+size]
+            offset = int(offset)
+            size = int(size)
         except Exception as e:
             print(e, offset, size)
-            
+            return None
+        
+        ret_val = population[offset:offset+size]
         if len(ret_val) < size:
             ret_val += population[0:size-len(ret_val)]
         
         return ret_val
-
-    
+    
+    
+    def resetOffsets(self):
+        for k in self._cur_offsets.keys():
+            self._cur_offsets[k] = 0
+
+            
     def clearCache(self):
         self._population_cache = {}
+        self._offset_cache = {}
+        self._cur_offsets = {}
 
         
@@ -188 +211 @@
         self.conn.commit()
         return ret_val
+
+    def fetchClassifierResult(self, classifier, trial_type, num_observations):
+        query = """
+          SELECT * FROM classifier_results
+          WHERE classifier=? AND trial_type=? AND num_observations=? 
+          ORDER BY false_positives+false_negatives
+          LIMIT 1;
+        """
+        cursor = self.conn.cursor()
+        cursor.execute(query, (classifier, trial_type, num_observations))
+        ret_val = cursor.fetchone()
+
+        if ret_val != None:
+            ret_val = dict(ret_val)
+        return ret_val
+
+    def deleteClassifierResults(self, classifier, trial_type, num_observations=None):
+        params = {"classifier":classifier,"trial_type":trial_type,"num_observations":num_observations}
+        query = """
+          DELETE FROM classifier_results
+          WHERE classifier=:classifier AND trial_type=:trial_type
+        """
+        if num_observations != None:
+            query += " AND num_observations=:num_observations"
+        
+        self.conn.execute(query, params)
+        self.conn.commit()
+