1#!/usr/local/bin/python3.8 2# Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.org/sumo 3# Copyright (C) 2007-2019 German Aerospace Center (DLR) and others. 4# This program and the accompanying materials 5# are made available under the terms of the Eclipse Public License v2.0 6# which accompanies this distribution, and is available at 7# http://www.eclipse.org/legal/epl-v20.html 8# SPDX-License-Identifier: EPL-2.0 9 10# @file flowFromRoutes.py 11# @author Daniel Krajzewicz 12# @author Jakob Erdmann 13# @author Michael Behrisch 14# @date 2007-06-28 15# @version $Id$ 16 17from __future__ import absolute_import 18from __future__ import print_function 19import math 20import sys 21import os 22 23from xml.sax import make_parser, handler 24from optparse import OptionParser 25 26SUMO_HOME = os.environ.get('SUMO_HOME', 27 os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', '..')) 28sys.path.append(os.path.join(SUMO_HOME, 'tools')) 29 30from sumolib.miscutils import uMin, uMax # noqa 31import detector # noqa 32from detector import relError # noqa 33 34 35def make_geh(interval_minutes): 36 def geh(m, c): 37 # GEH expects hourly flow 38 m = m * 60 / interval_minutes 39 c = c * 60 / interval_minutes 40 """Error function for hourly traffic flow measures after Geoffrey E. Havers""" 41 if m + c == 0: 42 return 0 43 else: 44 return math.sqrt(2 * (m - c) * (m - c) / float(m + c)) 45 return geh 46 47 48SEP = ";" 49 50 51def print_record(*args, **kwargs): 52 comment = '#' + SEP if kwargs.get('comment', False) else '' 53 print(comment + SEP.join([repr(a) if type(a) is float else str(a) for a in args])) 54 55 56class LaneMap: 57 58 def get(self, key, default): 59 return key[0:-2] 60 61 62class DetectorRouteEmitterReader(handler.ContentHandler): 63 64 def __init__(self, detFile): 65 self._routes = {} 66 self._detReader = detector.DetectorReader(detFile, LaneMap()) 67 self._edgeFlow = {} 68 self._parser = make_parser() 69 self._parser.setContentHandler(self) 70 self._begin = None 71 self._end = None 72 self.minTime = uMax 73 self.maxTime = uMin 74 75 def reset(self, start, end): 76 self._routes = {} 77 self._edgeFlow = {} 78 self._begin = 60 * start 79 self._end = 60 * end 80 81 def addRouteFlow(self, route, flow): 82 for edge in self._routes[route]: 83 if edge not in self._edgeFlow: 84 self._edgeFlow[edge] = 0 85 self._edgeFlow[edge] += flow 86 87 def startElement(self, name, attrs): 88 if name == 'route': 89 if 'id' in attrs: 90 self._routes[attrs['id']] = attrs['edges'].split() 91 if name == 'vehicle': 92 if self._begin is None or float(attrs['depart']) >= self._begin: 93 self.addRouteFlow(attrs['route'], 1) 94 self.minTime = min(self.minTime, float(attrs['depart'])) 95 self.maxTime = max(self.maxTime, float(attrs['depart'])) 96 if name == 'flow': 97 if 'route' in attrs: 98 if self._begin is None or float(attrs['begin']) >= self._begin and float(attrs['end']) <= self._end: 99 self.addRouteFlow(attrs['route'], float(attrs['number'])) 100 self.minTime = min(self.minTime, float(attrs['begin'])) 101 self.maxTime = max(self.maxTime, float(attrs['end'])) 102 if name == 'routeDistribution': 103 if 'routes' in attrs: 104 routes = attrs['routes'].split() 105 nums = attrs['probabilities'].split() 106 for r, n in zip(routes, nums): 107 self.addRouteFlow(r, float(n)) 108 109 def readDetFlows(self, flowFile, flowCol): 110 if self._begin is None: 111 return self._detReader.readFlows(flowFile, flow=flowCol) 112 else: 113 return self._detReader.readFlows(flowFile, flow=options.flowcol, time="Time", 114 timeVal=self._begin / 60, timeMax=self._end / 60) 115 116 def getDataIntervalMinutes(self, fallback=60): 117 if self.maxTime == uMin: 118 return fallback 119 interval = (self.maxTime - self.minTime) / 60 120 if interval == 0: 121 interval = fallback 122 return interval 123 124 def clearFlows(self): 125 self._detReader.clearFlows() 126 127 def calcStatistics(self, interval, geh_threshold): 128 rSum = 0. 129 dSum = 0. 130 sumAbsDev = 0. 131 sumSquaredDev = 0. 132 sumSquaredPercent = 0. 133 sumGEH = 0. 134 nGEHthresh = 0. 135 n = 0 136 geh = make_geh(interval) 137 for edge, detData in self._detReader._edge2DetData.items(): 138 rFlow = self._edgeFlow.get(edge, 0) 139 for group in detData: 140 if group.isValid: 141 dFlow = group.totalFlow 142 if dFlow > 0 or options.respectzero: 143 rSum += rFlow 144 dSum += dFlow 145 dev = float(abs(rFlow - dFlow)) 146 sumAbsDev += dev 147 sumSquaredDev += dev * dev 148 if dFlow > 0: 149 sumSquaredPercent += dev * dev / dFlow / dFlow 150 sumGEH += geh(rFlow, dFlow) 151 if geh(rFlow, dFlow) < geh_threshold: 152 nGEHthresh += 1 153 n += 1 154 if self._begin is not None: 155 print_record('interval', self._begin, comment=True) 156 print_record('avgRouteFlow', 'avgDetFlow', 'avgDev', 'RMSE', 'RMSPE', 'GEH', 'GEH%', comment=True) 157 if n == 0: 158 # avoid division by zero 159 n = -1 160 print_record( 161 rSum / n, 162 dSum / n, 163 sumAbsDev / n, 164 math.sqrt(sumSquaredDev / n), 165 math.sqrt(sumSquaredPercent / n), 166 sumGEH / n, 167 100. * nGEHthresh / n, 168 comment=True) 169 170 def printFlows(self, includeDets, useGEH, interval): 171 edgeIDCol = ["edge"] if options.edgenames else [] 172 timeCol = ["Time"] if options.writetime else [] 173 measureCol = "ratio" 174 measure = relError 175 if useGEH: 176 measureCol = "GEH" 177 measure = make_geh(interval) 178 179 if includeDets: 180 cols = ['Detector'] + edgeIDCol + timeCol + ['RouteFlow', 'DetFlow', measureCol] 181 else: 182 cols = ['Detector'] + edgeIDCol + timeCol + ['qPkw'] 183 print_record(*cols) 184 output = [] 185 time = self._begin if self._begin is not None else 0 186 for edge, detData in self._detReader._edge2DetData.items(): 187 detString = [] 188 dFlow = [] 189 for group in detData: 190 if group.isValid: 191 detString.append(group.getName( 192 options.longnames, options.firstname)) 193 dFlow.append(group.totalFlow) 194 rFlow = len(detString) * [self._edgeFlow.get(edge, 0)] 195 edges = len(detString) * [edge] 196 if includeDets: 197 output.extend(zip(detString, edges, rFlow, dFlow)) 198 else: 199 output.extend(zip(detString, edges, rFlow)) 200 if includeDets: 201 for group, edge, rflow, dflow in sorted(output): 202 if dflow > 0 or options.respectzero: 203 edgeCol = [edge] if options.edgenames else [] 204 timeCol = [time] if options.writetime else [] 205 cols = [group] + edgeCol + timeCol + [rflow, dflow, measure(rflow, dflow)] 206 print_record(*cols) 207 else: 208 for group, edge, flow in sorted(output): 209 edgeCol = [edge] if options.edgenames else [] 210 timeCol = [time] if options.writetime else [] 211 cols = [group] + edgeCol + timeCol + [flow] 212 print_record(*cols) 213 214 215optParser = OptionParser() 216optParser.add_option("-d", "--detector-file", dest="detfile", 217 help="read detectors from FILE (mandatory)", metavar="FILE") 218optParser.add_option("-r", "--routes", dest="routefile", 219 help="read routes from FILE (mandatory)", metavar="FILE") 220optParser.add_option("-e", "--emitters", dest="emitfile", 221 help="read emitters from FILE (mandatory)", metavar="FILE") 222optParser.add_option("-f", "--detector-flow-file", dest="flowfile", 223 help="read detector flows to compare to from FILE", metavar="FILE") 224optParser.add_option("-c", "--flow-column", dest="flowcol", default="qPKW", 225 help="which column contains flows", metavar="STRING") 226optParser.add_option("-z", "--respect-zero", action="store_true", dest="respectzero", 227 default=False, help="respect detectors without data (or with permanent zero) with zero flow") 228optParser.add_option("-D", "--dfrouter-style", action="store_true", dest="dfrstyle", 229 default=False, help="emitter files in dfrouter style (explicit routes)") 230optParser.add_option("-i", "--interval", type="int", help="aggregation interval in minutes") 231optParser.add_option("--long-names", action="store_true", dest="longnames", 232 default=False, help="do not use abbreviated names for detector groups") 233optParser.add_option("--first-name", action="store_true", dest="firstname", 234 default=False, help="use first id in group as representative") 235optParser.add_option("--edge-names", action="store_true", dest="edgenames", 236 default=False, help="include detector group edge name in output") 237optParser.add_option("-b", "--begin", type="float", default=0, help="begin time in minutes") 238optParser.add_option("--end", type="float", default=None, help="end time in minutes") 239optParser.add_option("--geh", action="store_true", dest="geh", 240 default=False, help="compare flows using GEH measure") 241optParser.add_option("--geh-treshold", type="float", default=5, dest="geh_threshold", 242 help="report percentage of detectors below threshold") 243optParser.add_option("--write-time", action="store_true", dest="writetime", 244 default=False, help="write time in output") 245optParser.add_option("-v", "--verbose", action="store_true", dest="verbose", 246 default=False, help="tell me what you are doing") 247(options, args) = optParser.parse_args() 248if not options.detfile or not options.routefile or not options.emitfile: 249 optParser.print_help() 250 sys.exit() 251parser = make_parser() 252if options.verbose: 253 print("Reading detectors") 254reader = DetectorRouteEmitterReader(options.detfile) 255parser.setContentHandler(reader) 256if options.interval: 257 haveFlows = True 258 start = options.begin # minutes 259 while ((options.end is None and haveFlows) or 260 (options.end is not None and start < options.end)): 261 end = start + options.interval 262 if options.end is not None: 263 end = min(end, options.end) 264 reader.reset(start, end) 265 if options.verbose: 266 print("Reading routes") 267 parser.parse(options.routefile) 268 if options.verbose: 269 print("Reading emitters") 270 parser.parse(options.emitfile) 271 if options.flowfile: 272 if options.verbose: 273 print("Reading flows") 274 haveFlows = reader.readDetFlows(options.flowfile, options.flowcol) 275 if haveFlows: 276 reader.printFlows(bool(options.flowfile), options.geh, options.interval) 277 if options.flowfile: 278 reader.calcStatistics(options.interval, options.geh_threshold) 279 reader.clearFlows() 280 start += options.interval 281else: 282 if options.verbose: 283 print("Reading routes") 284 parser.parse(options.routefile) 285 if options.verbose: 286 print("Reading emitters") 287 parser.parse(options.emitfile) 288 if options.flowfile: 289 if options.verbose: 290 print("Reading flows") 291 reader.readDetFlows(options.flowfile, options.flowcol) 292 fallbackInterval = 60 293 if options.begin and options.end: 294 fallbackInterval = options.end - options.begin 295 dataInterval = reader.getDataIntervalMinutes(fallbackInterval) 296 reader.printFlows(bool(options.flowfile), options.geh, dataInterval) 297 if options.flowfile: 298 reader.calcStatistics(dataInterval, options.geh_threshold) 299