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