home/diploma/lws_root/thesis.py

import logging
from collections import defaultdict

from matplotlib import pyplot as plt
from matplotlib.ticker import FuncFormatter
import scipy

import numpy as np

import lws.localize as localizer
import utils.accelerated as speed
from utils import path

log = logging.getLogger('lws')

def currentApSurfacePlot(gui):
    num_contours = gui.surface_number_of_contours
    gui.surface_view.visible = False
    
    gui.reloadAP('108')
    
    gui.ipw_z.ipw.slice_index = 23
    
    for so in gui.scene_object_groups:
        so.visible = so.name in ('OG 1', 'Stuff', 'Stairs')
        
    gui.showScalarBar(title='RSSI in dbm')
    
    gui.scene.background = (1.0,1.0,1.0)
    gui.ipw_z.visible = False
    
    gui.surface_number_of_contours = 150
    gui.scene.mlab.view(51.2, 13.49, 40.13, [ 38.82, 6.77 , -0.58])
    gui.scene.set_size((1200, 680))
    gui.surface_view.visible = True
    
    # gui.surface_view.contour.auto_contours = False
    #~ _debug()
    
    gui.scene.mlab.savefig(r'D:\loco-dev\dirk\thesis\figures\current_ap_surface.png')
    gui.scene.mlab.savefig(r'D:\loco-dev\dirk\thesis\figures\current_ap_surface.jpg')
    #~ gui.surface_number_of_contours = num_contours


def pruningPrepare(gui, cube_width=3, prune_to=3000):
    gui.device_name = 'iconia'
    gui.doLoadPaths()
    
    for s in gui.localization_paths:
        if s.id == 'og1_eg_right':
            gui.selected_path_sequence = s
    
    gui.measurement_runid = '4'
    gui.loadMeasurementRun()
    
    gui.cube_width = cube_width
    
    gui.doLocalize(threaded=False)
    
def pruningPlot(gui):
    gui.surface_number_of_contours = 5
    
    pruningPrepare(gui)
    
    for so in gui.scene_object_groups:
        so.visible = so.name in ('OG 1', 'Stuff', 'Stairs')
    
    gui.ipw_z.visible = False
    
    gui.pruning_idx = 6
    gui.pruning_visible_count = 300
    gui.viewPruningSingle()
    
    gui.surface_number_of_contours = 100
    
    gui.scene.mlab.view(-57.603641999035005, 25.219052810285703, 27.428001533423213, [ 32.94838163,   9.41659425,   1.45038033])
    gui.scene.set_size((1200, 680))
    
    gui.surface_view.visible = True
    
    gui.scene.background = (0.8,0.8,0.8)
    
    gui.scene.mlab.savefig(r'D:\loco-dev\dirk\thesis\figures\thesis_pruning3d_plot.png')
    gui.scene.mlab.savefig(r'D:\loco-dev\dirk\thesis\figures\thesis_pruning3d_plot.jpg')
    
    gui.scene.background = (0.5,0.5,0.5)

def localizationPlot(gui):
    
    
    pruningPrepare(gui, cube_width=2)
    
    for so in gui.scene_object_groups:
        so.visible = so.name in ('EG', 'OG 1', 'Stairs')
    
    gui.showScalarBar(title='Error in m')
    
    gui.ipw_z.visible = False
    
    #~ gui.pruning_idx = 6
    #~ gui.pruning_visible_count = 300
    #~ gui.viewPruningSingle()
    
    gui.surface_number_of_contours = 50
    
    gui.scene.mlab.view(-61.482491732159787, 22.479891467276911, 35.710743801652669, [40.16922802,   8.7447422 ,  -2.92382367])
    gui.scene.set_size((1100, 570))
    
    gui.surface_view.visible = True
    
    gui.scene.background = (0.8,0.8,0.8)
    
    gui.scene.mlab.savefig(r'D:\loco-dev\dirk\thesis\figures\thesis_localization_result3d.png')
    gui.scene.mlab.savefig(r'D:\loco-dev\dirk\thesis\figures\thesis_localization_result3d.jpg')
    
    gui.scene.background = (0.5,0.5,0.5)
    
def pruningScatter(gui):
    destfile = r'D:\loco-dev\dirk\thesis\figures\pruning_scatter.png'
    
    pruningPrepare(gui, cube_width=3, prune_to=500)
    
    
    POINT_SIZE = 1.0
    size_x, size_y, size_z = gui.localizer.cubed_data_shape
    ravel = lambda x, y, z: size_z * size_y * x + size_z * y + z
    
    aspectRatio = 0.25
    dpi = 90
    width = 1000
    figsize = width / float(dpi), width / float(dpi) * aspectRatio

    fig = plt.figure(figsize=figsize)
    ax = plt.axes([0.08, 0.16, 0.91, 0.81], xlabel='time $t$', ylabel=r'states $s$')
    ax.set_ylim(0, size_x * size_y * size_z)
    ax.set_xlim(0, len(gui.localizer.decoder.pruning_history))
    
    for i, p_data in enumerate(gui.localizer.decoder.pruning_history):
        xs = []
        ys = []
        
        log.info('processing idx %s of %s' % (i, len(gui.localizer.decoder.pruning_history)))
        curr_ys = []
        for tx, ty, tz, cost in p_data:
            r = ravel(tx, ty, tz)
            ys.append(r)
            curr_ys.append(r)
            xs.append(i)
            
        INTERPOLATE = 20
        for j in range(1, INTERPOLATE):
            for y in curr_ys:
                xs.append(i + j / float(INTERPOLATE))
                ys.append(y)
        log.info('scattering...')
        ax.scatter(xs, ys, s=POINT_SIZE, linewidth=0, color='black', marker='.',  antialiased=False, alpha=0.002)
        
        
    
    #~ leg = plt.legend(apids)
    #~ ltext  = leg.get_texts()
    #~ llines = leg.get_lines()
    #~ plt.setp(ltext, fontsize='x-small')
    #~ plt.setp(llines, linewidth=0.7)
    #~ leg.draw_frame(False)
    
    #~ if fname is None:
        #~ plt.show()
    #~ else:
    
    log.info('saving to %s...' % destfile)
    fig.savefig(destfile, dpi=dpi)

def latex3dmodel(gui, mesh):
    print mesh
    matnames =  ['MatConcrete',	'MatLightWalls', 'MatDoors',
    'MatGlassDoor', 'MatIronDoor', 'MatFascade', 'MatGlassWindow',
    'MatCupboard', 'MatTable', 'MatHardware', 'MatRailing']
    
    mat2avgthickness = defaultdict(lambda : '')
    mat2volumeratio = defaultdict(lambda : 0.3)
    
    mat2avgthickness['MatDoors'] = 0.03
    mat2volumeratio['MatDoors'] = 0.08
    
    mat2avgthickness['MatGlassDoor'] = 0.02
    mat2volumeratio['MatGlassDoor'] = 0.05
    
    mat2avgthickness['MatIronDoor'] = 0.05
    mat2volumeratio['MatIronDoor'] = 0.13
    
    mat2avgthickness['MatGlassWindow'] = 0.05
    mat2volumeratio['MatGlassWindow'] = 0.13
    
    
    mat2avgthickness['MatConcrete'] = 0.3
    mat2volumeratio['MatConcrete'] = 0.75
    
    mat2avgthickness['MatLightWalls'] = 0.15
    mat2volumeratio['MatLightWalls'] = 0.40
    
    mat2avgthickness['MatFascade'] = 0.4
    mat2volumeratio['MatFascade'] = 0.8
    
    mat2triangles = defaultdict(set)
    mat2vertices = defaultdict(set)
    
    for objname, triangles in mesh.objects.items():
        matname = mesh.materials[objname]
        mat2triangles[matname].update(triangles)
        for v1, v2, v3 in triangles:
            mat2vertices[matname].update([v1, v2, v3])
    
   
    env = localizer.Environment(gui.objfile, aps=localizer.getAPsFromGui(gui))
    
    shape = (gui.resolution.x, gui.resolution.y, gui.resolution.z)
    size = gui.resolution.x * gui.resolution.y * gui.resolution.z
    
    cubesize_m3 = env.di.sx * env.di.sy * env.di.sz
    print cubesize_m3
    
    mat2volume = defaultdict(list)
    for objname, triangles in mesh.objects.items():
        
        
        np_ts = np.zeros((len(triangles), 3, 3))
        vv = mesh.vertices
        for i, (v1, v2, v3) in enumerate(triangles):
            
            p1 = (vv[0][v1], vv[1][v1], vv[2][v1])
            p2 = (vv[0][v2], vv[1][v2], vv[2][v2])
            p3 = (vv[0][v3], vv[1][v3], vv[2][v3])
            
            np_ts[i, 0, 0] = p1[0]
            np_ts[i, 0, 1] = p1[1]
            np_ts[i, 0, 2] = p1[2]
            np_ts[i, 1, 0] = p2[0]
            np_ts[i, 1, 1] = p2[1]
            np_ts[i, 1, 2] = p2[2]
            np_ts[i, 2, 0] = p3[0]
            np_ts[i, 2, 1] = p3[1]
            np_ts[i, 2, 2] = p3[2]
        
        print 'building %s' % objname
        
        CACHE = True
        cachefile = path('r:/%s.npy' % objname)
        if CACHE and cachefile.exists():
            blockedCubes = np.load(cachefile)
        else:
            blockedCubes = speed.buildTriangleCubeIntersection(np_ts, np.zeros((0, 3, 3)), np.zeros((0, 3, 3)), cubed_data_shape=shape, cube_width=1, vi_dimensions=env.di)    
            np.save(cachefile, blockedCubes)
        
        numblocked = len(np.where(blockedCubes > 0)[0])
        
        matname = mesh.materials[objname]
        mat2volume[matname].append(numblocked * cubesize_m3)
    
    print '----'
    for matname in matnames:
        
        avgt = mat2avgthickness[matname]
        ratio = mat2volumeratio[matname]
        vol = sum(mat2volume[matname]) * ratio
        
        print '    %s & %s & %s & %d & %s \\\\ \\hline' % (
                matname[3:], len(mat2vertices[matname]), len(mat2triangles[matname]), vol, avgt)
    print '----'

def plotAPCoverage(gui):
    destfile = r'D:\loco-dev\dirk\thesis\figures\ap_coverage_hist.png'
    
    to_be_summed = []
    
    for ap in gui.aps:
        print 'processing: %s' % ap
        if not ap.used:
            continue
        
        gui.lazyLoadAP(ap)
        a = np.zeros(ap.data_in_db.shape, dtype=np.int)
        a[np.where(ap.data_in_db > -100)] = 1
        to_be_summed.append(a)
        
    print 'combining...'
    coverage_map = np.add.reduce(to_be_summed)
    
    print 'plotting'
    
    aspectRatio = 0.3
    dpi = 90
    width = 1000
    figsize = width / float(dpi), width / float(dpi) * aspectRatio

    fig = plt.figure(figsize=figsize)
    ax = plt.axes([0.095, 0.16, 0.90, 0.80], xlabel='number of visible APs', ylabel=r'covered Volume')
    
    b = coverage_map.ravel()
    avg_numaps = sum(b) / float(len(b))
    
    plt.axvline(avg_numaps, color='black', linewidth=2.0)
    print 'avg: %.1f' % avg_numaps
    
    ax.hist(b, bins=20, color='#9999FF')
    #~ ax.set_ylim(0, size_x * size_y * size_z)
    ax.set_xlim(0, 21)
    
    ax.yaxis.set_major_formatter(FuncFormatter(lambda x, pos=None: '%d $m^3$' % (x * 0.2**3)))
    ax.grid()
    
    log.info('saving to %s...' % destfile)
    
    fig.savefig(destfile, dpi=dpi)


def plotReadingCounts(counts):
    destfile = r'D:\loco-dev\dirk\thesis\figures\reading_counts_hist.png'
    
    print 'plotting'
    
    aspectRatio = 0.3
    dpi = 90
    width = 1000
    figsize = width / float(dpi), width / float(dpi) * aspectRatio

    fig = plt.figure(figsize=figsize)
    ax = plt.axes([0.095, 0.16, 0.90, 0.80], xlabel='XX', ylabel=r'YY')
    
    plt.axvline(scipy.average(counts), color='black', linewidth=2.0)
    print 'avg: %.1f' % scipy.average(counts)
    
    ax.hist(counts, color='#9999FF')
    #~ ax.set_ylim(0, size_x * size_y * size_z)
    #~ ax.set_xlim(0, 21)
    
    #~ ax.yaxis.set_major_formatter(FuncFormatter(lambda x, pos=None: '%d $m^3$' % (x * 0.2**3)))
    ax.grid()
    
    log.info('saving to %s...' % destfile)
    
    fig.savefig(destfile, dpi=dpi)
    
def latexDevice2Measurements(gui):
    device2measuredloc = defaultdict(set)
    device2measuredap = defaultdict(set)
    device2readingcount = defaultdict(int)
    device2readingcount_all = defaultdict(int)
    device2api = {'iconia': 'android', 'nexus':'android', 'galaxyp':'android', 'nexus1':'android', 'messbook':'libpcap', 'wispy':'custom'}
    device2class = {'iconia': 'tablet', 'nexus':'phone', 'galaxyp':'phone', 'nexus1':'phone', 'messbook':'usb-adapter', 'wispy':'spectrum analyzer'}
    
    device2ap2count = defaultdict(lambda: defaultdict(int))
    
    device2apgroup2stddev = defaultdict(lambda: defaultdict(list))
    device2apgroup2count = defaultdict(lambda: defaultdict(int))
    device2stddevs = defaultdict(list)
    locid2readingcount = defaultdict(int)
    
    dnames = ['iconia', 'nexus', 'galaxyp', 'nexus1', 'messbook', 'wispy']
    for dname in dnames:
        gui.retriveMeasurements(device_name=dname, autoreload=False, overlay=False)
        gui.loadLocations(dname)
        
        for ap in gui.aps:
            gui.loadMeasurements(ap)
            for m in gui.measurements:
                if m.available:
                    device2measuredloc[dname].add(m.locid)
                    device2readingcount[dname] += 1
                    device2readingcount_all[dname] += len(m.rssis)
                    
                    device2ap2count[dname][ap.id] +=1
                    #~ print scipy.std(m.rssis)
                    #~ print m.rssis
                    device2stddevs[dname].append(scipy.std(m.rssis))
                    device2apgroup2stddev[dname][ap.powid].append(scipy.std(m.rssis))
                    device2apgroup2count[dname][ap.powid] += 1
                    locid2readingcount[m.locid] += len(m.rssis)
                    device2measuredap[dname].add(ap.id)
    
    # a little bit ugly - but the AP has been moved is therefore not included anymore
    device2measuredap['wispy'].add('110')
    device2readingcount['wispy'] += 6
    device2readingcount_all['wispy'] += 6
    
    print '---'
    for dname in dnames:
        clocs = len(device2measuredloc[dname])
        rc = device2readingcount[dname]
        arc = device2readingcount_all[dname]
        numaps = len(device2measuredap[dname])
        print '    %s & %s & %s & %s & %s & %s & %s\\\\ \\hline' % (dname, device2class[dname], device2api[dname], clocs, numaps, rc, arc)
    
    print '---'
    
    _fmt = lambda powid: powid.replace('wrt54g', 'wrt').replace('eduroam', 'edu')
    powids = device2apgroup2stddev['iconia'].keys()
    
    print r'    \begin{tabularx}{0.98\textwidth}{|X|c|c|' + ''.join(['c|'] * len(powids) *2) + '}\hline'
    print r'    \rowcolor[gray]{.92}'
    print r'    \textbf{Device} & \textbf{$N_{all}$} & \textbf{$\sigma_{all}$} &', 
    print ' & '.join(r'\textbf{$N_{%s}$} & \textbf{$\sigma_{%s}$}' % (_fmt(e), _fmt(e)) for e in powids),
    print r'\\ \hline\hline'
    
    _avg = lambda l: sum(l) / float(len(l))
    
    totals_powids = []
    for dname in dnames[:5]:
        print '    %s &' % dname,
        print '%s &' % device2readingcount[dname],
        
        print ' %.1f &' % _avg(device2stddevs[dname]),
        
        print ' & '.join('%s & %.1f' % (device2apgroup2count[dname][e], _avg(device2apgroup2stddev[dname][e])) for e in powids),
            
        #~ print ' & '.join(str(e) for e in device2ap2count[dname]values()), 
        print '\\\\ \\hline'
    
    print 'total/avg & ',
    print sum(device2readingcount[dname] for dname in dnames[:5]),
    print ' & ',
    print '%.1f' % _avg(reduce(lambda a, b: a + b, (device2stddevs[dname] for dname in dnames[:4]))),
    print ' & ',
    for i, e in enumerate(powids):
        print sum(device2apgroup2count[dname][e] for dname in dnames[:5]),
        print ' & ' ,
        print '%.1f' % _avg(reduce(lambda a, b: a + b, (device2apgroup2stddev[dname][e] for dname in dnames[:4]))),
        
        if i < len(powids)-1:
            print '& ',
            
    print '\\\\ \\hline'
    #~ print sum(device2apgroup2count[dname][e] for e in powid for dname in dnames[:4]),
        
    print r'    \end{tabularx}'
    
    print '----'
    print 'avg reading count: %d' % scipy.average(locid2readingcount.values())
    plotReadingCounts(locid2readingcount.values())
    
def generateSyntheticPaths(gui):
    aps = localizer.getAPsFromGui(gui)
    #~ for ap in aps:
        #~ gui.lazyLoadAP(ap)
    env = localizer.Environment(gui.objfile, aps=aps)
    SPEED = 1.0
    GRANULARITY = 0.75
    
    destpath = path(r'D:\loco-dev\dirk\tmp\tracked_path_synthetic')
    #~ gui.cube_width = 3
    
    assert len(gui.localization_paths) > 0
    
    for NOISE in np.arange(0, 18, 0.5):
        for p in gui.localization_paths:
            if p.id.startswith('still_') or p.id == 'custom':
                continue
            print p.id, 'noise: %s' % NOISE
                
            d = destpath.joinpath('sigma_%.1f/%s' % (NOISE, p.id))
            if d.exists():
                d.rmtree()
            d.makedirs()
            
            gui._loadPathSequence(p)
            gui._localization_path_changed(refreshCubes=False)
            
                    
            #~ half_cube = (gui.res_dx * gui.cube_width) / 2.0
            #~ positions = [np.array(gui.translateToMeter(*(xyz * gui.cube_width))) + half_cube for xyz in cubes]
            for x in range(20):
                corners = [(c.x, c.y, c.z) for c in gui.localization_path]
                measurements = env.generateSyntheticMeasurementsFromCorners(corners=corners, 
                                speed=SPEED, 
                                noise=NOISE,
                                granularity=GRANULARITY)
        
                #~ measurements = env.generateSyntheticMeasurementsAtPositions(positions, speed=SPEED, noise=NOISE, sample_environment=half_cube/2.0)
                f = d.joinpath('measurements_%02d.txt' % x)
                f.write_text(str(measurements))
                print 'writing %s' % f.abspath()
                
            #~ print measurements
            
    
def pathInfo(gui):
    assert len(gui.localization_paths) > 0
    print '---'
    devices = 'iconia', 'nexus'
    
    stats  = defaultdict(dict)
    
    device2pathid2collected = defaultdict(lambda: defaultdict(int))
    for device in devices:
        gui.device_name = device
        gui.doLoadPaths()

        
        for p in gui.localization_paths:
           
            if p.id.startswith('still_') or p.id == 'custom':
                continue
            gui._loadPathSequence(p)
            gui._localization_path_changed(refreshCubes=False)
            
            pid = p.id
            if pid.endswith('_r'):
                pid = p.id[:-2]
            
            device2pathid2collected[device][p.id] = p.collect_runs            
            
            corners = gui.localization_path
            dist = 0
            for c1, c2 in zip(corners, corners[1:]):
                dist += ((c2.x - c1.x)**2 + (c2.y - c1.y)**2 + (c2.z - c1.z)**2)**0.5
            
            stats[pid]['dist'] = dist
            stats[pid]['turns'] = len(gui.localization_path)
            
            num_aps = []
            num_signals = []
            durations = []
            
            for r in gui.measurement_runs:
                run = gui.runid2run.get(r)
                
                measurements = localizer.Measurements()
                measurements.loadFromString(run.measurements)
                #~ measurements.interpolateSignals()
                
                for m in measurements:
                    num_aps.append(len(m.apdata))
                
                durations.append((measurements[-1].timestamp - measurements[0].timestamp).total_seconds())
                
                num_signals.append(len(measurements))
            
            if not 'num_aps' in stats[pid]:
                stats[pid]['num_aps'] = []
            stats[pid]['num_aps'].extend(num_aps)
            
            if not 'num_signals' in stats[pid]:
                stats[pid]['num_signals'] = []
            stats[pid]['num_signals'].extend(num_signals)
            
            if not 'durations' in stats[pid]:
                stats[pid]['durations'] = []
            stats[pid]['durations'].extend(durations)
        
            #~ gui.measurement_runid = r
            #~ gui.loadMeasurementRun()
    
    for pid, data in sorted(stats.items()):
        avg_signals = sum(data['num_signals']) / float(len(data['num_signals']))
        avg_aps = sum(data['num_aps']) / float(len(data['num_aps']))
        p = pid.replace('_', '-')
        print '    %s & $%.1fm$ & %d & %.1f & %.1f & %s\\\\ \\hline' % (
                p, data['dist'], data['turns'],
                avg_signals, avg_aps, r'\ref{fig:path_%s}' % pid)
    
    print '---'
    pids = stats.keys()
    for pid in pids:
        p = pid.replace('_', '-')
        a = device2pathid2collected['iconia'][pid]
        b = device2pathid2collected['iconia'][pid + '_r']
        c = device2pathid2collected['nexus'][pid]
        d = device2pathid2collected['nexus'][pid + '_r']
        fig = r'\ref{fig:path_%s}' % pid
        speed = duration = 0
        
        dd = stats[pid]['durations']
        avg_duration = sum(dd) / float(len(dd))
        speed = stats[pid]['dist'] / avg_duration
        
        print '    %s & $%ds$ & $%.1fm/s$ & %d/%d & %d/%d\\\\ \\hline' % (p, avg_duration, speed, a, b, c, d)