analysis.py

Go to the documentation of this file.

#!/usr/bin/env python
#  tvaLib Copyright (c) 2012-2016 Paul G. St-Aubin
#  Ecole Polytechnique de Montreal, McGill University
#  Python 2.7; (dt) Spyder Windows 10 64-bit; ipython Ubuntu 15.04 64-bit

import os, sys, math
from copy    import deepcopy
from inspect import getfile      as inspect_getfile
from inspect import currentframe as inspect_currentframe

try:              import numpy as np   
except Exception: raise Exception, [0101, 'Numpy is not installed.']

if __name__ == '__main__':
    print('Analysis library loaded directly.')
    sys.path.insert(0,os.path.dirname(os.path.dirname(os.path.abspath(inspect_getfile(inspect_currentframe()))))) 
import lib.tools       as tvaLib
     

def rateString(obs, total):
    ''' Take observations and total sample and format a string like "10/100 (10%)"  '''
    if(total == 0): return 'None'
    if(type(obs) is float): obs_ = round(obs,1)
    else:                   obs_ = obs
    return str(obs_)+'/'+str(total)+' ('+str(round(float(obs)/float(total)*100,1))+'%)'


class Measures():
    ''' The most basic analysis container contains individual values, e.g. one 
        value for each observation, or each site analysis. This could contain, 
        for example the mean speed at a site, or, the total flow, etc. '''
    def __init__(self, prototype=None):
        self.values    = []
        self.prototype = prototype
        return
        
    def __len__(self):               return len(self.values)
    def __getitem__(self, i):        return self.values[i]
        
    def add(self, *args):
        if(self.prototype == None): return False
        self.values.append(deepcopy(self.prototype))
        self.values[-1].addInstructions(*args)
        return True
        
    def increment(self, *args, **kwargs):
        value = self.values[-1].increment(*args, **kwargs)
        return value

    def getNonZeroValueIndeces(self):
        return_ = []
        for i in range(len(self.values)):
            if(self.values[i].value): return_.append(i)
        return return_
    def getMeasures(self, round_=-1):
        ''' Return measure of all analyses '''
        indeces = self.getNonZeroValueIndeces()
        measures = []
        for i in indeces:
            measures.append(self.values[i].value)
        if(round_>=0 and isinstance(round_, int)):
            measures = [round(x, round_) for x in measures]
        return measures
    def getDepths(self, round_=-1):
        ''' Return depths of all analyses '''
        indeces = self.getNonZeroValueIndeces()
        depths = []
        for i in indeces:
            depths.append(self.values[i].depth)
        if(round_>=0 and isinstance(round_, int)):
            depths = [round(x, round_) for x in depths]
        return depths
    def getMeasuresAggregatedBySite(self, round_=-1):
        ''' Return measure of all sites '''
        indeces = self.getNonZeroValueIndeces()
        siteIds = self.getSiteIds()
        measures = []
        for s in siteIds:
            localmeasure = 0
            localDepth = 0
            for i in indeces:
                if(self.values[i].metadata.idx == s): localmeasure = tvaLib.Math.combineMean([localmeasure,localDepth], [self.values[i].value, self.values[i].depth])
            if(localmeasure): measures.append(localmeasure)
        if(round_>=0 and isinstance(round_, int)):
            measures = [round(x, round_) for x in measures]
        return measures
    def getSiteIds(self, NonZeroOnly=False):
        ''' Return list of site IDs '''
        ids = []
        for i in self.values:
            if(NonZeroOnly and i.value == 0): continue
            ids.append(i.metadata.idx)
        ids = tvaLib.unique(ids)
        return ids
    
    def getSiteRankBySaIx(self, saIx):
        try:    return self.getSiteIdRanking().index(saIx)
        except: return '-'

    def getSiteRankingByIndex(self):
        ''' Return ranked id aggregated by sites '''
        measures = self.getMeasuresAggregatedBySite()
        return sorted(range(len(measures)), key=measures.__getitem__)
    def getSiteIdRanking(self):
        ''' Return ranked id aggregated by sites '''
        indeces = self.getSiteRankingByIndex()
        siteIds = self.getSiteIds(NonZeroOnly=True)
        return [siteIds[x] for x in indeces]       
    def getSiteRankingByMeasure(self, round_=-1):
        ''' Return sorted measure aggregated by sites '''
        return sorted(self.getMeasuresAggregatedBySite(round_))
    
    def getMeasureRankingByIndex(self, style='verbose'):
        ''' Return ranked index of all analyses '''
        measures = self.getMeasures()
        siteIndeces = self.getNonZeroValueIndeces()
        if(measures and siteIndeces):
            measures, siteIndeces = zip(*sorted(zip(measures, siteIndeces)))
        if(style=='verbose'): return [self.values[x].metadata.name for x in siteIndeces]
        else:                 return [x for x in siteIndeces]  
    def getMeasureRankingByMeasure(self, round_=-1):
        ''' Return sorted measure of all analyses '''
        return sorted(self.getMeasures(round_))
        
    def printResult(self,label='default', round_=-1):
        print '============================================'
        print 'Ranking by '+label+' aggregated by site'
        print self.getSiteIdRanking()
        print self.getSiteRankingByMeasure(round_)
        print 'Ranking by '+label+' for all analyses'
        print self.getMeasureRankingByIndex()
        print self.getMeasureRankingByMeasure(round_)        

class Frequencies(Measures):
    ''' This analysis contains a large number of observations.  '''
    def getNonZeroValueIndeces(self):
        return_ = []
        for i in range(len(self.values)):
            #print self.values[i].binData
            if(self.values[i].binData): return_.append(i)
        return return_
    def getBins(self):
        ''' Return bins of all analyses '''
        indeces = self.getNonZeroValueIndeces()
        bins = []
        for i in indeces:
            bins.append(self.values[i].binStops)
        return bins
    def getFrequencies(self):
        ''' Return frequencies of all analyses '''
        indeces = self.getNonZeroValueIndeces()
        frequencies = []
        for i in indeces:
            frequencies.append(self.values[i].binData)
        return frequencies
    def getSumFrequency(self):
        freq = self.getFrequencies()
        for i in range(len(freq)):
            if(i==0): sumFreq = freq[i]
            else:     sumFreq = [x+y for x,y in zip(sumFreq,freq[i])]
        return sumFreq
    def getSTDEVofBins(self):
        ''' This function returns the standard deviation on each bin from the
            sample (freqSamples).  '''
        return_ = []
        for data in self.getFrequencies():
            if(not data or type(data) is not list): return False
            mean = sum(data)/len(data)
            squareOfResult=[]
            for obs in data:
                squareOfResult.append(math.pow(obs-mean,2))
            return_.append(math.sqrt(sum(squareOfResult)/len(squareOfResult)))
        return return_
    def getMAXIMAofBins(self, mode='max'):
        ''' This function returns the min/max outliers on each bin from the
            sample (freqSamples).  '''
        return_ = []
        for data in self.getFrequencies():
            if(not data or type(data) is not list): return False
            if(mode=='min'): return_.append(min(data))
            else:            return_.append(max(data))
        return return_
    
    def clusterByKStest(self, nklusters):
        ''' Cluster distributions by k-test distances into nklusters groups
        
            Input:
            ======
            self:      should be a correctly instantiated, non-empty list of
                       frequency objects
            nklusters: number of cluster centroids (equivalent to k in k-means)
            
            Output:
            =======
            returns a one-dimensional array of corresponding cluster associated
            with each frequency object contained by this object
            Also overwrites any self.clusters
            
            '''
        if(not self.values): return False  
        if(sum([1 for x in self.values if x]) < nklusters): nklusters = sum([1 for x in self.values if x]) 
        
        
        meanKsDistances = []
        for distroNeedleIx in range(len(self.values)):
            if(not self.values[distroNeedleIx].getCDF()): 
                meanKsDistances.append(None)
                continue
            ksDistances = []
            for distroStackIx in range(len(self.values)):
                if(distroNeedleIx == distroStackIx or not self.values[distroStackIx].getCDF()): continue
                ksDistances.append(tvaLib.Math.ksTest(self.values[distroNeedleIx].getCDF(),self.values[distroStackIx].getCDF()))
            if(len(ksDistances) == 0): meanKsDistances.append(0)
            else:                      meanKsDistances.append(sum(ksDistances)/float(len(ksDistances)))
        self.centroidIxs = sorted(range(len(meanKsDistances)), key=lambda i: meanKsDistances[i])[meanKsDistances.count(None):nklusters+meanKsDistances.count(None)]
        
        
        self.clusters = []
        
        for distroNeedleIx in range(len(self.values)):
            if(not self.values[distroNeedleIx].getCDF()): 
                self.clusters.append(None)
                continue
            ksDistances = []
            isCentroid  = False
            for klusterStackIx in range(nklusters):
                if(distroNeedleIx == self.centroidIxs[klusterStackIx]): 
                    isCentroid = True
                    break
                ksDistances.append(tvaLib.Math.ksTest(self.values[distroNeedleIx].getCDF(),self.values[self.centroidIxs[klusterStackIx]].getCDF()))
            if(isCentroid): self.clusters.append(klusterStackIx)
            else:           self.clusters.append(ksDistances.index(min(ksDistances)))
        
        return self.clusters

class Profiles(Measures):
    def getBins(self):
        ''' Return profile bins. '''
        for value in self.values:
            try:    return value.bins
            except: continue
    
    
    def getMean(self, min_up_coverage=0, min_down_coverage=0):
        mean = []
        for bIx in range(len(self.getBins())):
            try: mean.append(tvaLib.Math.combineMean([value.binData[bIx] for value in self.values if value.binData and value.up_coverage > min_up_coverage and value.down_coverage > min_down_coverage], [value.binWeights[bIx] for value in self.values if value.binData and value.up_coverage > min_up_coverage and value.down_coverage > min_down_coverage]))
            except ZeroDivisionError: mean.append(0)
        return mean
        
    def getStdDev(self, min_up_coverage=0, min_down_coverage=0):
        stdev = []
        for bIx in range(len(self.getBins())):
            try: stdev.append(tvaLib.Math.combineStdDev([value.binStdDev[bIx] for value in self.values if value.binData and value.up_coverage > min_up_coverage and value.down_coverage > min_down_coverage], [value.binData[bIx] for value in self.values if value.binData], [value.binWeights[bIx] for value in self.values if value.binData and value.up_coverage > min_up_coverage and value.down_coverage > min_down_coverage]))
            except ZeroDivisionError: stdev.append(0)
        return stdev
    



class AnalysisInstrument():
    metadata = None

    def addInstructions(self, metadata=None):
        self.metadata = metadata
        return True
        
class Measure(AnalysisInstrument):
    def __init__(self):
        self.value = 0.0
        self.depth = 0
        return
    
    def increment(self, data, weights=1):
        ''' data is a list of observations '''
        if(type(data) is not list):    data    = [data]
        if(type(weights) is not list): weights = [weights]
        if(len(data) <= 0): return False
        value = self.incrementByMean(data, weights)
        return value   
    
    def incrementByMean(self, data, weights=1):
        ''' data is a simple list '''
        if(type(data) is not list):                data    = [data]
        if(weights and type(weights) is not list): weights = [weights]
        if(len(data) <= 0): return False
        if(weights and len(weights) == len(data)):
            sum_val = 0
            len_val = 0
            for i in range(len(data)):
                sum_val += data[i]*weights[i]
                len_val += weights[i]
            mean_val = sum_val/float(len_val)
        else:
            len_val = len(data)
            mean_val = sum(data)/float(len_val)
        self.value = tvaLib.Math.combineMean([self.value, mean_val], [self.depth, len_val])
        self.depth += len(data)
        return mean_val


class MeasureByHour(AnalysisInstrument):
    def __init__(self):
        ''' To get hourly values, use getCurve()[hour], etc. '''
        self.values = [[] for x in range(1440)]
        self.depths = [0 for x in range(1440)]
        return
    def increment(self, data, startTimes, startTime, duration, framerate):
        ''' data:       a list of measures. 
            startTimes: a list of starttimes corresponding to each measure in
                        data (lists should be the same size)
            startTime:  start time of measure recording  
            duration:   duration of measure recording  
            '''
        if(len(data) <= 0): return False
        
        for measure,start in zip(data,startTimes):
            timestampBin = int((startTime.hour*60 + startTime.minute + start/60.0/framerate))
            self.values[timestampBin].append(measure)
        for bin in range(startTime.hour*60 + startTime.minute, int(startTime.hour*60 + startTime.minute + duration/60.0)):
            self.depths[bin] += 1
        return True
    def getCurve(self, bins=24, type='rate', emptyVal=None, divideBy=1, divideByValueVariety=False, strictDepthCompleteness=False):
        ''' 2 bins is am/pm, 24 bins is hourly, 96 bins is quarter hourly, etc.
            
            divideByValueVariety might be handy if the measures are
            categorical.
            
            If strictDepthCompleteness is set to True, only return values for
            bins with a complete depth coverage (in other words the entire
            period has had data collected over that period). '''
        if(divideByValueVariety): divideBy = len(list(set(tvaLib.flatten_list(self.values))))
        returnList  = []
        returnDepth = []
        for bin in range(bins):
            if(strictDepthCompleteness and [] in self.values[(bin)*int(1440/bins):(bin+1)*int(1440/bins)]): 
                returnList.append(emptyVal)
                returnDepth.append(0)
            elif(self.values[(bin)*int(1440/bins):(bin+1)*int(1440/bins)].count([]) == len(self.values[(bin)*int(1440/bins):(bin+1)*int(1440/bins)])): 
                returnList.append(emptyVal)
                returnDepth.append(0)
            else:
                returnList.append(tvaLib.flatten_list(self.values[(bin)*int(1440/bins):(bin+1)*int(1440/bins)]))
                returnDepth.append(sum(self.depths[(bin)*int(1440/bins):(bin+1)*int(1440/bins)]))
        for i,depth in zip(range(len(returnList)),returnDepth):
            if(returnList[i] and returnList[i] != emptyVal):
                if(type=='mean'):  returnList[i] = sum(returnList[i])/float(len(returnList[i]))
                elif(type=='std'): returnList[i] = np.std(returnList[i])
                else:
                    if(depth): returnList[i] = len(returnList[i])/float(depth)*int(1440/bins)/float(divideBy)
                    else:      returnList[i] = emptyVal
        return returnList


class MeanSpeedByHour(MeasureByHour):
    def increment(self, data, startTime, duration, framerate, mps_kmh=3.6):
        ''' data is a list of objects (trajectories). '''
        if(len(data) <= 0): return False
        
        for datum in data:
            timestampBin = int((startTime.hour*60 + startTime.minute + datum.getFirstInstant()/60.0/framerate))
            self.values[timestampBin].append(sum(datum.velocities.positions[2])/float(len(datum.velocities.positions[2]))*mps_kmh*framerate)
        return True

class Pool(AnalysisInstrument):
    ''' Pool together raw data. Usefull if the data is structured using a custom scheme. '''
    def __init__(self):
        self.value = []
        return
    
    def increment(self, data):
        ''' Data is a list of observations '''
        self.value.append(data)
        return self.value   

class Frequency(AnalysisInstrument):
    def __init__(self, binSize=0.25, xstart=0, xend=10):
        self.binStops = [xstart]
        self.binData  = []
        self.binSize  = binSize
        self.depth    = 0
        self.source   = []
        
        while(True):
            xstart += self.binSize
            self.binStops.append(xstart)
            self.binData.append(0)
            if(xstart > xend): break
        return
    
    def increment(self, data, weights=None, method='byData', keepSource=False):
        ''' Increment data.
            
            Input:
            ======
            method=='byData':
            Data is a simple list of individual observations
            method=='byMeanHisto':
            Use this method if the histogram has been pre calculated '''
        if(len(data) <= 0): return False
        
        if(method=='byMeanHisto'):
            (histo,xs) = data
            if(self.depth == 0):
                self.binStops    = xs
                self.binSize     = xs[1]-xs[0]
                self.binData     = list(histo)
                self.depth       = sum(histo)
                return True
            if(len(histo) != len(self.binData)): return False
            for histIx in range(len(histo)):
                self.binData[histIx]     = (self.binData[histIx]*self.depth+histo[histIx]*sum(histo))/(self.depth+sum(histo))
            self.depth   += sum(histo)
        else:
            if(keepSource): self.source += data
            if(weights): result = list(np.histogram(data, bins=self.binStops, weights=weights)[0])
            else:        result = list(np.histogram(data, bins=self.binStops)[0])
            for i in range(len(result)):
                self.binData[i] += result[i]
            self.depth += len(data)
        return True

        
    def verifyBinStops(self, data):
        ''' This method is called to build/expand'''
        data = sorted(data)
        if(data[0] < self.binStops[-1] or data[-1] > self.binStops[-1]): return False
        else:                                                            return True
    
    def getFreq(self):
        if(not sum(self.binData)): return False
        return self.binData
    
    def getPDF(self):
        if(not sum(self.binData)): return False
        return [x/float(sum(self.binData)) for x in self.binData]
    
    def getCDF(self):
        frequency = self.getPDF()
        if(not frequency): return False
        return_ = []
        for i in range(len(frequency)):
            return_.append(sum(frequency[0:i+1]))
        return return_
        
    def ksTestWith(self, refData):
        thisData = self.getCDF()
        if(isinstance(refData, Frequency)): refData = refData.getCDF()
        if(not thisData or type(refData) != list or len(refData) != len(thisData)): return False
        supremum = 0
        for i in range(len(thisData)):
            if(math.fabs(thisData[i] - refData[i]) > supremum): supremum = math.fabs(thisData[i] - refData[i])
        return supremum
        
class Profile(AnalysisInstrument):
    def __init__(self):
        self.binData       = []
        self.bins          = []
        self.binWeights    = []
        self.binStdDev     = []
        self.binMins       = []
        self.binMaxs       = []
        self.up_coverage   = 0
        self.down_coverage = 0
        return
    
    def setBins(self, values): self.bins = values; return True
    
    def increment(self, data, weights, stdDevs=None, mins=None, maxs=None):
        ''' Increment data. '''
        if(len(data) <= 0): return False
        
        if(not self.binData):
            self.binData    = data
            self.binWeights = weights
            if(stdDevs): self.binStdDev = stdDevs
            if(mins):    self.binMins = mins
            if(maxs):    self.binMaxs = maxs
        else:
            if(len(data) != len(self.binData) or len(data) != len(weights)): return False
            if(stdDevs): self.binStdDev = [tvaLib.Math.combineStdDev([stdDevs[ix], self.binStdDev[ix]], [data[ix], self.binData[ix]], [weights[ix], self.binWeights[ix]]) for ix in range(len(data))]
            self.binData                = [tvaLib.Math.combineMean([data[ix], self.binData[ix]], [weights[ix], self.binWeights[ix]]) for ix in range(len(data))]
            if(mins):    self.binMins   = [m1 if m1 < m2 else m2 for m1,m2 in zip(mins, self.binMins)]
            if(maxs):    self.binMaxs   = [m1 if m1 > m2 else m2 for m1,m2 in zip(maxs, self.binMaxs)]
            for ix in range(len(data)): self.binWeights[ix] += weights[ix]
        return True
        



class Exposure(Measure):
    def increment(self, data, duration=3600):
        ''' data is a generic rate (int) '''
        if(not isinstance(data, list)): data = [data]
        self.incrementByMean([x*float(3600/duration) for x in data])
        return True
        
class Total(Measure):
    def increment(self, data):
        ''' data is a generic total (int). '''
        self.value += data
        return True
        
class MeanSpeed(Measure):
    def increment(self, data, speed_conv):
        ''' data is a list of objects (trajectories). '''
        if(len(data) <= 0): return False
        
        data_ = []
        for datum in data:
            data_.append(sum(datum.velocities.positions[2])/float(len(datum.velocities.positions[2]))*speed_conv)
        mean_val = self.incrementByMean(data_)
        return mean_val
        
class FlowRatePerLanePerHour(Measure):
    def increment(self, data, duration=3600, location=0):
        ''' data is a list of objects (trajectories). '''
        if(len(data) <= 0): return False
        
        data_ = [0 for i in range(len(self.metadata.site.alignments))]
        for datum in data:
            data_[datum.curvilinearPositions.getLanes()[location]] += 1
        for i in range(len(data_)):
            data_[i] = data_[i] * 3600/duration
        if(len(filter(None, data_))): mean_val = self.incrementByMean(filter(None, data_))
        else:                         mean_val = self.incrementByMean([0])
        return mean_val
        
class VehKmTraveled(Measure):
    def increment(self, data):
        ''' data is a list of objects (trajectories). '''
        if(len(data) <= 0): return False
        value = 0
        for datum in data:
            value += sum(datum.velocities.positions[2])/1000.0
            self.depth += 1
        self.value += value
        return value

class IndicatorThresholds(Measure):
    def __init__(self, thresh_ranges):
        self.thresh_ranges = thresh_ranges
        self.values  = [0 for i in thresh_ranges]
        self.weights = [[] for i in thresh_ranges]
        self.total   = 0
        return
    def increment(self, data, weights=None):
        ''' data is a a list of indicators 
            weights is a list of weights (0-1) of equal size to data
            '''
        self.total += len(data)
        for thresh in range(len(self.thresh_ranges)):
            for datum in data:
                if(datum <= self.thresh_ranges[thresh]):
                    self.values[thresh] += 1
            if(weights):
                for weight,datum in zip(weights,data):
                    if(datum <= self.thresh_ranges[thresh]):
                        self.weights[thresh].append(weight)
        return True
        
    def getValue(self, thresh):         return self.values[thresh]
    def getWeightedValue(self, thresh): return self.values[thresh]*sum(self.weights[thresh])/float(len(self.weights[thresh]))
        
#K-S test