home/diploma/lws_root/measure.py

from collections import namedtuple, defaultdict
import sqlite3
from datetime import datetime
import time

import rpyc

import numpy as np
from scipy.misc import pilutil
import scipy.ndimage.filters as image_filters
import scipy.ndimage as ndimage

from utils.path import path
from utils import getPyPlot

fmtts = lambda ts: str(datetime.fromtimestamp(ts)).split('.')[0]
UTC_OFFSET = 3600 * 2 # 2 hours summertime

WISPY_OFFSET = 0.7

CHANNELS = range(1, 15)
BURST_PERIOD = 6.0
ZERO_PERIOD = 6.0
CHANNELS_B = {
    #~ 0: (2401, 2495),
    1: (2401, 2423), 2: (2406, 2428), 3: (2411, 2433), 
    4: (2416, 2438), 5: (2421, 2443), 6: (2426, 2448), 
    7: (2431, 2453), 8: (2436, 2458), 9: (2441, 2463), 
    10: (2446, 2468), 11: (2451, 2473), 12: (2456, 2478), 
    13: (2461, 2483), 14: (2473, 2495)
}
#~ CHANNELS_B = dict([(e[0], (e[1][0] + 7, e[1][1] - 7)) for e in CHANNELS_B.items()])

BURST_STATIONS = {
    '107': ('127.0.0.1', 18107),
    '108': ('127.0.0.1', 18108),
    '110': ('127.0.0.1', 18110),
}

def runBurst(station):
    ''' connect to AP and send bursts over different channels'''
    host, port = BURST_STATIONS[station]
    
    channel2ts = {}
    conn = rpyc.classic.connect(host, port=port)
    conn.modules.sys.path.append('/sbin')
    chaburst = conn.modules.chaburst

    try:
        chaburst.stopDefaultWifi()
        chaburst.startBurstDevice()

        for chan in CHANNELS:
            print 'measure channel %s' % chan
            chaburst.setChannel(chan)
            ts1 = time.time()
            chaburst.sendBurst(BURST_PERIOD)
            channel2ts[chan] = (ts1, time.time())
        
    finally:
        chaburst.startDefaultWifi()
        conn.close()
    
    print 'measure zero channel'
    # store timestamps for a period without traffic at fake channel 0
    #~ ts1 = time.time()
    #~ time.sleep(ZERO_PERIOD)
    #~ channel2ts[0] = (ts1, time.time())
            
    return channel2ts
    
def gatherTimestampsForPositions(stations):
    ''' collect measurements at different locations interactivly
    '''
    print 'Enter first location index:', 

    loc2ts = {} # key (station, location), value {channel: (ts_start, ts_end)}
    while True:
        try:
            loc = raw_input()
            loc = int(loc)
        except ValueError:
            if loc == 'q':
                print 'done...'
                break
            else:
                print 'need a int value for location!'
                continue
            
        ts1 = time.time()
        for station in stations:
            print 'start measure %s for location %s at station %s' % (fmtts(ts1), loc, station)
            #~ print 'end measurement with [return]',
            #~ raw_input()
            try:
                channel2ts = runBurst(station)
            except Exception, KeyboardInterrupt:
                print 'error during running burst'
                continue
        
            #~ ts2 = time.time()
            
            loc2ts[(station, loc)] = channel2ts
            
        print 'continue with next location: ',

    print 'got the following locations'
    for (loc, station), ts in loc2ts.items():
        print loc, station, ts
    
    return loc2ts
    
def getSignalStrength(data):
    '''estimate the paramters of a 1d GMM with 2 modes'''
    from scikits.learn.mixture import GMM
    # assume 2 modes - one gaussian for the noise, and one for the signal
    clf = GMM(n_states=2, cvtype='full')
    clf.fit(np.array(data).reshape(len(data), 1))
    # assume the larger mean is our signal
    
    signalStrengthMax = max(clf.means)
    idx = list(clf.means).index(signalStrengthMax)
    varianceMax = list(clf.covars)[idx]
    
    signalStrengthMin = min(clf.means)
    return signalStrengthMax, varianceMax, signalStrengthMin

def updateSignalData(locations, ssid, wispy_db, loc2ts):
    ''' update Location() data in locations with measured rssi values for given ssid
    the values are extracted from the wispy_db file by using timestamps from measurements_file
    '''
    # currently only first ts is used
    for locid, (ts1, ts2) in loc2ts.items():
        loc2ts[locid] = ts1
        
    # now use these timestamps to extract signals from wispy-db
    firstts, lastts = min(loc2ts.values()), max(loc2ts.values())
    
    res = {}
    with sqlite3.connect(wispy_db) as conn:
        c = conn.cursor()
        sql = '''select 
                    rssi, 
                    milliseconds_since_epoch 
                from 
                    network_data 
                where 
                    milliseconds_since_epoch > ?
                    and milliseconds_since_epoch < ?
                    and ssid = ?'''
        
        c.execute(sql, ((firstts - 10 + UTC_OFFSET) * 1000, (lastts + 10 + UTC_OFFSET) * 1000, ssid))
        ts2rssi = {}
        for row in c:
            ts = row[1] / 1000 - UTC_OFFSET
            ts2rssi[ts] = row[0]
            #~ print datetime.fromtimestamp(ts)
        
        for loc, ts in loc2ts.items():
            # iterate over all stored rssi values
            # and pick that one with the minimum timestamp delta (thats not fast but should be ok)
            delttats, key = min([(abs(ts - e), e) for e in ts2rssi.keys()])
            #~ print '%s at location %s with rssi: %s [delta %.2f sec]' % (fmtts(ts), loc, ts2rssi[key], delttats)
            locations[loc] = locations[loc]._replace(rssi=ts2rssi[key], ts=ts)

def ts2dbts(ts):
    ''' convert ts to milliseconds_since_epoch '''
    return (ts + UTC_OFFSET) * 1000

def dbts2ts(dbts):
    return dbts / 1000.0 - UTC_OFFSET
    
def rasterize(xs, ys, zs, readings_per_sweep, destfile):
    first_ts = np.min(xs)
    last_ts = np.max(xs)
    delta = last_ts - first_ts
    #~ print delta
    resolution_x = len(set(xs))
    resolution_y = readings_per_sweep
    data = np.zeros((resolution_x, resolution_y))
    
    for x, y, z in zip(xs, ys, zs):
        _x = int(round((x - first_ts) / float(delta) * (resolution_x - 1)))
        data[_x][y] = z
    
    data = pilutil.imresize(data, (1000, resolution_y))
    
    data = data.transpose()
    pilutil.imsave(destfile, data)
    
        
Location = namedtuple('Location', 'id x y z rssi ts')
def loadLocations(f):
    print 'loading locations from %s' % f.abspath()
    lines = [l for l in f.lines() if l.strip() and not l.strip().startswith('#')]
    res = {}
    for l in lines:
        i, x, y, z = [float(e.strip()) for e in l.split()]
        i = int(i)
        res[i] = Location(i, x, y, z, None, None)
    return res

    
def loadMeasurements():
    ''' return a list of locations with measured rssi signal strengths '''
    p = path(r'D:\loco-dev\maps\whs19\experiment\measure_locs.txt')
    locs = loadLocations(p)
    wispy_db = path('d:/loco-dev/maps/whs19/experiment/measure.wsx')
    run_file = path('d:/loco-dev/maps/whs19/experiment/run1.txt')
    
    updateSignalData(locs, 'CS31', wispy_db, run_file)
    return locs


Sweep = namedtuple('Sweep', 'dt data setting channel')
Setting = namedtuple('Setting', 'id starting_frequency frequency_resolution readings_per_sweep amplitude_offset amplitude_resolution rssi_max')
class Measure(object):
    def __init__(self, location_file, wispy_db, plotpath, loc2tsfile):
        self.location_file = location_file
        self.wispy_db = wispy_db
        self.loc2tsfile = loc2tsfile
        if not self.loc2tsfile.exists():
            self.loc2tsfile.touch()
        
        # contents wispy_db setting table
        self.settings = {} # key: setting id, value: Setting()
        # load data into self.settings
        self.loadWiSpySettings()
        
        self.loc2ts = defaultdict(list) # key: location id, values: list of dicts containing timestamps for each channel
        # load data into self.loc2ts
        self.loadMeasurements()
        
        self.plotpath = plotpath
        if not plotpath.exists():
            plotpath.makedirs()
        
    def loadWiSpySettings(self):
        self.settings.clear()
        
        with sqlite3.connect(self.wispy_db) as conn:
            c = conn.cursor()
            
            # load all settings
            sql = '''select 
                        id, 
                        starting_frequency,
                        frequency_resolution,
                        readings_per_sweep,
                        amplitude_offset,
                        amplitude_resolution,
                        rssi_max
                    from 
                        setting
                    '''
            c.execute(sql)
            for row in c:
                self.settings[row[0]] = Setting(*row)
    
    
    def loadSweepDataForLocation(self, locid, station):
        
        sweeps = defaultdict(list) # keys: measure_idx, values: [Sweep(), ]
        with sqlite3.connect(self.wispy_db) as conn:
            c = conn.cursor()
            
            channeldata = self.loc2ts[(locid, station)]
            
            for measure_idx, channel2ts in enumerate(channeldata):
                
                allts = []
                for ts1, ts2 in channel2ts.values():
                    allts.extend([ts1, ts2])
                
                fromts, untilts = min(allts), max(allts)
                
                # load sweep data from firsts to lastts
                sql = '''select 
                            sweep_data, 
                            milliseconds_since_epoch,
                            setting_id
                        from 
                            sweep 
                        where 
                            milliseconds_since_epoch > ?
                            and milliseconds_since_epoch < ?
                        order by
                            milliseconds_since_epoch
                        '''
                MORE_SECONDS = 10 # fetch more seconds at the left/right side of intervall
                c.execute(sql, ((fromts - MORE_SECONDS + UTC_OFFSET) * 1000, (untilts + MORE_SECONDS + UTC_OFFSET) * 1000))
                
                for row in c:
                    ts = dbts2ts(row[1])
                    for channel, (ts1, ts2) in channel2ts.items():
                        if ts1 <= ts <= ts2:
                            sweep = Sweep(datetime.fromtimestamp(ts), row[0], row[2], channel)
                            sweeps[measure_idx].append(sweep)
                
                # check if we found any data in wispy db
                if not measure_idx in sweeps:
                    print datetime.fromtimestamp(fromts - MORE_SECONDS)
                    print 'got no sweep data from wispy db for measurement %s at location %s' % (measure_idx, locid)
        #~ print sweeps
        return sweeps
                
    
    def combineWispyDBs(self, source, dest, clear=False):
        if clear:
            with sqlite3.connect(dest) as conn:
                c = conn.cursor()
                c.execute('delete from network_data')
                c.execute('delete from sweep')
                c.execute('delete from setting')
                
        with sqlite3.connect(source) as conn:
                
            c = conn.cursor()
            # select last setting
            sql = '''select 
                        id,
                        device,
                        start_time_msepoch,
                        starting_frequency,
                        frequency_resolution,
                        readings_per_sweep,
                        amplitude_offset,
                        amplitude_resolution,
                        rssi_max
                    from 
                        setting 
                    order by 
                        id desc'''
            setting = list(c.execute(sql))[-1]
            sql = '''select
                        milliseconds_since_epoch,
                        sweep_data
                    from
                        sweep
                    where
                        setting_id = ?'''
            
            sweeps = list(c.execute(sql, str(setting[0])))
            print 'got %s sweeps' % len(sweeps)
            sql = '''select
                        ssid,
                        channel,
                        encryption,
                        mac,
                        rssi,
                        supported_rates,
                        signal_quality,
                        mode,
                        milliseconds_since_epoch
                    from
                        network_data'''
            network_data = list(c.execute(sql))
        
        print 'inserting data into %s' % dest
        with sqlite3.connect(dest) as conn:
            c = conn.cursor()
            sql = '''select max(id) from setting'''
            try:
                setting_id = list(c.execute(sql))[0][0] + 1
            except Exception:
                setting_id = 0
            print 'using new setting id %s' % setting_id
            
            sql = '''insert into setting (id, device, start_time_msepoch, starting_frequency, 
                            frequency_resolution, readings_per_sweep, amplitude_offset, amplitude_resolution, rssi_max)
                    values(?, ?, ?, ?, ?, ?, ?, ?, ?)'''
            
            c.execute(sql, [setting_id] + list(setting[1:]))
            
            sql = '''select max(id) from sweep'''
            try:
                next_sweep_id = list(c.execute(sql))[0][0] + 1
                
            except Exception:
                next_sweep_id = 0
            sql = '''insert into sweep (id, setting_id, milliseconds_since_epoch, sweep_data)
                    values(?, ?, ?, ?)'''
            
            for s_id, sweep in enumerate(sweeps, next_sweep_id):
                c.execute(sql, [s_id, setting_id] + list(sweep))
            
            sql = '''select max(id) from network_data'''
            try:
                next_nd_id = list(c.execute(sql))[0][0] + 1
            except Exception:
                next_nd_id = 0
            sql = '''insert into network_data (id, ssid, channel, encryption, mac, rssi, supported_rates, 
                        signal_quality, mode, milliseconds_since_epoch)
                    values(?, ?, ?, ?, ?, ? , ?, ?, ?, ?)'''
            
            for nd_id, nd in enumerate(network_data, next_nd_id):
                c.execute(sql, [nd_id] + list(nd))
            
            
            
    def runMeasurement(self, stations):
        ## check if stations are reachable
        for s in stations:
            try:
                conn = rpyc.classic.connect(BURST_STATIONS[s][0], port=BURST_STATIONS[s][1])
            except:
                raise RuntimeError('station %s is not reachable' % s)
            else:
                print 'station %s is reachable' % s
                conn.close()
            
        loc2ts = gatherTimestampsForPositions(stations)
        s = '\n'.join('%s/%s %r' % (station, locid, channel2ts) for (station, locid), channel2ts in loc2ts.items())
        print 'storing measurements result to %s' % self.loc2tsfile.abspath()
        self.loc2tsfile.write_text('run from %s at locations %s\n%s\n\n' % (datetime.now(), ','.join(str(e) for e in loc2ts), s), append=True)
        
        # refresh self.loc2ts
        self.loadMeasurements()
        self.loadWiSpySettings()
        
        print 'please save wispy db now!',
        raw_input()

    def loadMeasurements(self):
        ''' load measurements from loc2tsfile (next to location_file) and return
        dict with keys: location id, values: list of dicts that contain (start/end) timestamps for each channel
        '''
        self.loc2ts.clear()
        
        print 'loading measurements from %s' % self.loc2tsfile.abspath()
        for line in [l for l in self.loc2tsfile.lines() if l.strip()]:
            if line.startswith('run from'):
                continue
            #line format: 107/1 {1: (1315234649.809, 1315234651.031), 2: (1315234651 ...
            l = line.split()
            
            station, locid = l[0].split('/')
            locid = int(locid)
            
            channel2ts = eval(' '.join(l[1:]))
            
            self.loc2ts[(locid, station)].append(channel2ts)
            
        # adjust wispy timestamp offset (empirically determined!)
        for channeldata in self.loc2ts.values():
            for channel2ts in channeldata:
                for chanid, (ts1, ts2) in channel2ts.items():
                    channel2ts[chanid] = (ts1 + WISPY_OFFSET, ts2 + WISPY_OFFSET)
    
    def analyzeMeasurements(self, locid, measure_idx, station):
        #~ locs = loadLocations(self.location_file)
        sweep = self.loadSweepDataForLocation(locid, station)
        
        self.drawMeasure(sweep[measure_idx], locid, station, measure_idx)
        
        avg_sig = self.analyzeBurst(sweep[measure_idx], locid, measure_idx, draw=True)
        print 'got avg sig %.2f' % avg_sig
        
    
    def buildOptimizeFile(self, optimizefile, station):
        measure_idx = 0
        locs = loadLocations(self.location_file)
        loc2sig = []
        for locid, loc in locs.items():
            sweep = self.loadSweepDataForLocation(locid, station)
            avg_sig = self.analyzeBurst(sweep[measure_idx], locid, measure_idx)
            if avg_sig is not None:
                loc2sig.append((loc, avg_sig))
        
        s = '\n'.join('%s %s %s %s %.2f' % (loc.id, loc.x, loc.y, loc.z, avg_sig) for loc, avg_sig in loc2sig)
        optimizefile = path(optimizefile % station)
        print 'storing optimize file to %s' % optimizefile.abspath()
        optimizefile.write_text(s)
            
        
    def drawMeasure(self, sweeps, locid, station, measure_idx):
        channel2ts = self.loc2ts[(locid, station)][measure_idx] 
        
        plt, mlab, dates, font_manager, ticker = getPyPlot()
        
        width = 1200
        aspectRatio = 0.5
        dpi = 70
        
        figsize = width / float(dpi), width / float(dpi) * aspectRatio
        
        fig = plt.figure(figsize=figsize)
        ax1 = plt.axes([0.08, 0.21, 0.85, 0.75])
        
        # assert all Sweep() instances are from the same wispy setting
        assert len(set([s.setting for s in sweeps])) == 1
        setting = self.settings[sweeps[0].setting]
        
        starting_frequency = setting.starting_frequency
        frequency_resolution = setting.frequency_resolution / 1000
        amplitude_offset = setting.amplitude_offset
        amplitude_resolution = setting.amplitude_resolution
        readings_per_sweep = setting.readings_per_sweep
        
        
        xs, ys, zs = [], [], []
        for sweep in sweeps:
            assert len(sweep.data) == readings_per_sweep
            for i, c in enumerate(sweep.data):
                # TODO: is this correct?
                amplitude = amplitude_offset + ord(c) * amplitude_resolution
                
                #~ print amplitude
                
                if amplitude > -100:
                    xs.append(sweep.dt)
                    
                    y = (starting_frequency + (i * frequency_resolution)) / 1000.0
                    #~ print starting_frequency, i, frequency_resolution
                    ys.append(i)
                    zs.append(ord(c))
                    
        xs = mlab.date2num(xs)
        _ys = np.array(ys)
        _ys = ((_ys * frequency_resolution) + starting_frequency) / 1000
        ax1.scatter(xs, _ys , s=2, c='#AAAACC', linewidths=0)
        
        for channel, (ts1, ts2) in channel2ts.items():
            _xs = mlab.date2num([datetime.fromtimestamp(ts) for ts in (ts1, ts2)])
            ax1.plot(_xs, [CHANNELS_B[channel]] * 2, c='#7766FF')
        
        ax1.xaxis_date()
        ax1.xaxis.set_major_formatter(dates.DateFormatter('%H:%M:%S'))
        
        
        fig.savefig(self.plotpath.joinpath('full_%s_%s.png' % (locid, measure_idx)), format='png', dpi=dpi)
        plt.close(fig)
        
        rasterize(xs, ys, zs, readings_per_sweep, self.plotpath.joinpath('rast_%s_%s.png' % (locid, measure_idx)))
    
    def default_analyzer(self, channeldata):
        sigstrenghts = {}
        for c, cdata in channeldata.items():
            if c == 0:
                continue
                
            sigstrenghts[c] = getSignalStrength(cdata)
        
        MIN_STRENGTH = 3 # ignore signals that are too near at sig0
        all_sigs = [sig for sig, variance, sig0 in sigstrenghts.values() if sig > sig0 + MIN_STRENGTH and variance < 50]
        if len(all_sigs) > 0:
            avg_sig = sum(all_sigs) / float(len(all_sigs))
        else:
            avg_sig = -100
        
        return sigstrenghts, avg_sig
        
    def analyzeBurst(self, sweeps, locid, measure_idx, draw=False):
        c = len(set([s.setting for s in sweeps]))
        if c != 1:
            print 'anticipated exactly one setting for locid %s, measure_idx %s, but got %s' % (locid, measure_idx, c)
            return
            
        setting = self.settings[sweeps[0].setting]
        
        #~ y = (starting_frequency + (i * frequency_resolution)) / 1000.0
        channeldata = defaultdict(list)
        for sweep in sweeps:
            
            minfreq , maxfreq = CHANNELS_B[sweep.channel]
            
            _toidx = lambda f: int((f - (setting.starting_frequency / 1000.0)) / (setting.frequency_resolution / 10.0**6))
            minidx = _toidx(minfreq)
            maxidx = _toidx(maxfreq)
            
            for c in sweep.data[minidx:maxidx]:
                ampl = setting.amplitude_offset + ord(c) * setting.amplitude_resolution
                channeldata[sweep.channel].append(ampl)
                
        
        sigstrenghts, avg_sig = default_analyzer(channeldata)
        
        
        if draw:
            plt, mlab, dates, font_manager, ticker = getPyPlot()
            for c, cdata in channeldata.items():
                if c == 0:
                    continue
                    
                fig = plt.figure(figsize=(12, 6))
                ax = plt.axes()
                
                r, BIN_WIDTH = (-120, 0), 2
                bins = [e for e in range(r[0], r[1], BIN_WIDTH)]
                
                ax.hist(cdata, bins=bins)
                ax.set_xlim(r)
                
                sig, variance, sig0 = sigstrenghts[c]
                
                ax.axvline(sig)
                ax.axvline(sig0, color='red')
                ax.axvline(avg_sig, color='green', linewidth=2.0)
                
                
                ax1 = ax.twinx()
                ax1.plot(range(-120, 50), [mlab.normpdf(x, sig, variance)[0][0] for x in range(-120, 50)])
                ax1.set_xlim(r)
                
                plt.title('loc %s, run %s, chan %s [mu: %.1f, sigma: %.1f, mu0: %.1f, mu_total: %.2f]' % (locid, measure_idx, c, sig, variance, sig0, avg_sig))
                
                
                fig.savefig(self.plotpath.joinpath('channel_%02d_%02d_%02d.png' % (locid, measure_idx, c)), dpi=80)
                plt.close(fig)
    
        return avg_sig
        
if __name__ == '__main__':
    #python /usr/lib/python2.6/site-packages/rpyc/scripts/rpyc_classic.py -m forking -p 20001
    
    #~ wispy_db = path('d:/loco-dev/maps/whs19/experiment/measure.wsx')
    #~ location_file = path(r'D:\loco-dev\maps\whs19\experiment\measure_locs.txt')
    
    location_file = path(r'D:\loco-dev\maps\umic\experiment\locations.txt')
    wispy_db = path('d:/loco-dev/maps/umic/experiment/measure.wsx')
    loc2tsfile = path(r'D:\loco-dev\maps\umic\experiment\loc2ts.txt')
    plotpath = path('r:/')
    optimizefile = path('d:/loco-dev/maps/umic/experiment/measure_%s.txt')
    
    m = Measure(location_file, wispy_db, plotpath, loc2tsfile)
    
    stations = ['107', '108', '110']
    #~ stations = ['110']
    #~ stations = ['110']
    
    #~ m.runMeasurement(stations)
    m.combineWispyDBs('d:/loco-dev/maps/umic/experiment/measure1.wsx', 'd:/loco-dev/maps/umic/experiment/measure.wsx', clear=True)
    m.combineWispyDBs('d:/loco-dev/maps/umic/experiment/measure2.wsx', 'd:/loco-dev/maps/umic/experiment/measure.wsx')
    m.combineWispyDBs('d:/loco-dev/maps/umic/experiment/measure3.wsx', 'd:/loco-dev/maps/umic/experiment/measure.wsx')
    m.combineWispyDBs('d:/loco-dev/maps/umic/experiment/measure5.wsx', 'd:/loco-dev/maps/umic/experiment/measure.wsx',)
    #~ m.combineWispyDBs('d:/loco-dev/maps/umic/experiment/measure4.wsx', 'd:/loco-dev/maps/umic/experiment/measure.wsx')
    #~ m.combineWispyDBs('d:/loco-dev/maps/umic/experiment/measure6.wsx', 'd:/loco-dev/maps/umic/experiment/measure.wsx')
    
    #~ for i in range(1, 12):
        #~ m.analyzeMeasurements(i, 0)
    #~ m.analyzeMeasurements(3, 0, '110')
    #~ m.analyzeMeasurements(4, 0, '110')
    #~ m.analyzeMeasurements(15, 0, '110')
    for s in stations:
        m.buildOptimizeFile(optimizefile, s)