sumolib.statistics

  1# Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
  2# Copyright (C) 2012-2026 German Aerospace Center (DLR) and others.
  3# This program and the accompanying materials are made available under the
  4# terms of the Eclipse Public License 2.0 which is available at
  5# https://www.eclipse.org/legal/epl-2.0/
  6# This Source Code may also be made available under the following Secondary
  7# Licenses when the conditions for such availability set forth in the Eclipse
  8# Public License 2.0 are satisfied: GNU General Public License, version 2
  9# or later which is available at
 10# https://www.gnu.org/licenses/old-licenses/gpl-2.0-standalone.html
 11# SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
 12
 13# @file    statistics.py
 14# @author  Jakob Erdmann
 15# @author  Michael Behrisch
 16# @date    2012-05-08
 17
 18from __future__ import absolute_import
 19from __future__ import print_function
 20from __future__ import division
 21import math
 22import warnings
 23from collections import defaultdict
 24try:
 25    from numpy import sqrt, set_printoptions
 26except ImportError:
 27    from math import sqrt
 28
 29
 30def round(value):  # to round in Python 3 like in Python 2
 31    if value < 0:
 32        return math.ceil(value - 0.5)
 33    else:
 34        return math.floor(value + 0.5)
 35
 36
 37def identity(value):
 38    return value
 39
 40
 41class _ExtremeType(object):
 42    """
 43    see http://www.python.org/dev/peps/pep-0326/
 44    """
 45
 46    def __init__(self, isMax, rep):
 47        object.__init__(self)
 48        self._isMax = isMax
 49        self._rep = rep
 50
 51    def __eq__(self, other):
 52        return isinstance(other, self.__class__) and other._isMax == self._isMax
 53
 54    def __ne__(self, other):
 55        return not self == other
 56
 57    def __gt__(self, other):
 58        return self._isMax and not self == other
 59
 60    def __ge__(self, other):
 61        return self._isMax
 62
 63    def __lt__(self, other):
 64        return not self._isMax and not self == other
 65
 66    def __le__(self, other):
 67        return not self._isMax
 68
 69    def __repr__(self):
 70        return self._rep
 71
 72
 73uMax = _ExtremeType(True, "uMax")
 74uMin = _ExtremeType(False, "uMin")
 75
 76
 77def setPrecision(formatstr, precision, isArray=False, preformatted=False):
 78    if isArray:
 79        set_printoptions(precision=2)
 80        return formatstr.replace('%.2f', '%s')
 81    elif preformatted:
 82        return formatstr.replace('%.2f', '%s')
 83    else:
 84        return formatstr.replace('%.2f', '%.' + str(int(precision)) + 'f')
 85
 86
 87class Statistics:
 88
 89    def __init__(self, label=None, abs=False, histogram=False, printMin=True, scale=1, printDev=False, printSum=False):
 90        self.label = label
 91        self.min = uMax
 92        self.min_label = None
 93        self.max = uMin
 94        self.max_label = None
 95        self.values = []
 96        self.abs = abs
 97        self.printMin = printMin
 98        self.scale = scale
 99        self.printDev = printDev
100        self.printSum = printSum
101        self.isArray = False
102        if histogram:
103            self.counts = defaultdict(int)
104        else:
105            self.counts = None
106
107    def add(self, v, label=None):
108        self.values.append(v)
109        isMin = v < self.min
110        if type(isMin) is not bool:
111            self.isArray = True
112            v = tuple(v)
113        if v < self.min:
114            self.min = v
115            self.min_label = label
116        if v > self.max:
117            self.max = v
118            self.max_label = label
119        if self.counts is not None:
120            self.counts[int(round(v / self.scale))] += 1
121
122    def update(self, other):
123        for v in other.values:
124            self.add(v)
125
126    def clear(self):
127        self.min = uMax
128        self.min_label = None
129        self.max = uMin
130        self.max_label = None
131        self.values = []
132        if self.counts:
133            self.counts.clear()
134
135    def count(self):
136        return len(self.values)
137
138    def avg(self):
139        """return the mean value"""
140        # XXX rename this method
141        if len(self.values) > 0:
142            return sum(self.values) / float(len(self.values))
143        else:
144            return None
145
146    def avg_abs(self):
147        """return the mean of absolute values"""
148        # XXX rename this method
149        if len(self.values) > 0:
150            return sum(map(abs, self.values)) / float(len(self.values))
151        else:
152            return None
153
154    def meanAndStdDev(self, limit=None):
155        """return the mean and the standard deviation optionally limited to the last limit values"""
156        if limit is None or len(self.values) < limit:
157            limit = len(self.values)
158        if limit > 0:
159            mean = sum(self.values[-limit:]) / float(limit)
160            sumSq = 0.
161            for v in self.values[-limit:]:
162                sumSq += (v - mean) * (v - mean)
163            return mean, sqrt(sumSq / limit)
164        else:
165            return None
166
167    def relStdDev(self, limit=None):
168        """return the relative standard deviation optionally limited to the last limit values"""
169        moments = self.meanAndStdDev(limit)
170        if moments is None:
171            return None
172        return moments[1] / moments[0]
173
174    def mean(self):
175        warnings.warn("mean() is deprecated, because the name is misleading, use median() instead")
176        return self.median()
177
178    def mean_abs(self):
179        warnings.warn("mean_abs() is deprecated, because the name is misleading, use median_abs() instead")
180        return self.median_abs()
181
182    def average_absolute_deviation_from_mean(self):
183        if len(self.values) > 0:
184            m = self.avg()
185            return sum([abs(v - m) for v in self.values]) / len(self.values)
186        else:
187            return None
188
189    def median(self):
190        """return the median value"""
191        if len(self.values) > 0:
192            return sorted(self.values)[len(self.values) // 2]
193        else:
194            return None
195
196    def median_abs(self):
197        """return the median of absolute values"""
198        if len(self.values) > 0:
199            return sorted(map(abs, self.values))[len(self.values) // 2]
200        else:
201            return None
202
203    def quartiles(self):
204        s = sorted(self.values, key=lambda v: tuple(v) if self.isArray else v)
205        return s[len(self.values) // 4], s[len(self.values) // 2], s[3 * len(self.values) // 4]
206
207    def rank(self, fraction):
208        if len(self.values) > 0:
209            return sorted(self.values)[int(round(len(self.values) * fraction + 0.5))]
210        else:
211            return None
212
213    def histogram(self):
214        if self.counts is not None:
215            return [(k * self.scale, self.counts[k]) for k in sorted(self.counts.keys())]
216        else:
217            return "Histogramm is deactivated"
218
219    def toString(self, precision=2, histStyle=1, fmt=identity):
220        """histStyle
221            0 : not shown
222            1 : one line
223            2 : fancy
224            """
225        if len(self.values) > 0:
226            pre = fmt != identity
227            min = ''
228            if self.printMin:
229                min = setPrecision('min %.2f%s, ', precision, self.isArray, pre) % (
230                    fmt(self.min), ('' if self.min_label is None else ' (%s)' % (self.min_label,)))
231            result = setPrecision('%s: count %s, %smax %.2f%s, mean %.2f', precision, self.isArray, pre) % (
232                self.label, len(self.values), min,
233                fmt(self.max),
234                ('' if self.max_label is None else ' (%s)' %
235                 (self.max_label,)),
236                fmt(self.avg()))
237            result += setPrecision(', Q1 %.2f, median %.2f, Q3 %.2f', precision, self.isArray, pre) % (
238                tuple(map(fmt, self.quartiles())))
239            if self.abs:
240                result += setPrecision(', mean_abs %.2f, median_abs %.2f', precision, self.isArray, pre) % (
241                    fmt(self.avg_abs()), fmt(self.median_abs()))
242            if self.printDev:
243                result += (setPrecision(', stdDev %.2f', precision, self.isArray) % (self.meanAndStdDev()[1]))
244            if self.printSum:
245                result += (setPrecision(', sum %.2f', precision, self.isArray) % sum(self.values))
246            if self.counts is not None:
247                if histStyle == 1:
248                    result += '\n histogram: %s' % self.histogram()
249                elif histStyle == 2:
250                    keylen = len("%.0f" % (self.scale * max(self.counts.keys())))
251                    formatStr = "%%%i.0f: %%s" % keylen
252                    result = 'histogram of %s:\n%s\n%s' % (self.label,
253                                                           '\n'.join([formatStr % x for x in self.histogram()]),
254                                                           result)
255            return result
256        else:
257            return '%s: no values' % self.label
258
259    def toXML(self, precision=2, tag="statistic", indent=4, label=None, fmt=identity, extraAttributes={}):
260        pre = fmt != identity
261        if label is None:
262            label = self.label
263        description = ' description="%s"' % label if label != '' else ''
264
265        result = ' ' * indent + '<%s%s' % (tag, description)
266        for k, v in extraAttributes.items():
267            result += ' %s="%s"' % (k, v)
268        if self.count() > 0:
269            result += ' count="%i"' % self.count()
270            result += (setPrecision(' min="%.2f" minLabel="%s" max="%.2f" maxLabel="%s" mean="%.2f"',
271                                    precision, self.isArray, pre) %
272                       (fmt(self.min), self.min_label, fmt(self.max), self.max_label, fmt(self.avg())))
273            result += setPrecision(' Q1="%.2f" median="%.2f" Q3="%.2f"', precision, self.isArray, pre) % (
274                tuple(map(fmt, self.quartiles())))
275            result += (setPrecision(' meanAbs="%.2f" medianAbs="%.2f"', precision, self.isArray, pre) %
276                       (fmt(self.avg_abs()), fmt(self.median_abs())))
277            if self.printDev:
278                result += (setPrecision(' stdDev="%.2f"', precision, self.isArray) %
279                           (self.meanAndStdDev()[1]))
280        if self.counts is not None:
281            result += '>\n'
282            for kv in self.histogram():
283                result += setPrecision(8 * ' ' + '<hist key="%.2f" value="%i"/>\n', precision, self.isArray) % kv
284            result += ' ' * indent + '</%s>\n' % tag
285        else:
286            result += '/>\n'
287        return result
288
289    def __str__(self):
290        return self.toString()
291
292    def normalise_to_range(self, n_min=0, n_max=1):
293        """Normalises the stored list of values between n_min and n_max, Default: [0,1]"""
294        ret = []
295        range_length = n_max - n_min
296        values_diff = max(self.values) - min(self.values)
297        for val in self.values:
298            temp = (((val - min(self.values))*range_length)/values_diff) + n_min
299            ret.append(temp)
300        return ret
301
302
303def geh(m, c):
304    """Error function for hourly traffic flow measures after Geoffrey E. Havers"""
305    if m + c == 0:
306        return 0
307    else:
308        return math.sqrt(2 * (m - c) * (m - c) / (m + c))
309
310
311def sqv(self, m, c, scaling_factor=1000):
312    """Scaling Quality Value Calculation, Ref: https://journals.sagepub.com/doi/10.1177/0361198119838849
313        scaling_factor:
314        Number of person trips per day (total, per mode, per purpose) : 1
315        Mean trip distance in kilometers : 10
316        Duration of all trips per person per day in minutes : 100
317        Traffic volume per hour : 1000
318        Traffic volume per day : 10000
319    """
320    return 1/(1 + math.sqrt(((m-c)*(m-c))/(scaling_factor*c)))