1 | #region License Information
|
---|
2 | /*
|
---|
3 | * This file is part of SimSharp which is licensed under the MIT license.
|
---|
4 | * See the LICENSE file in the project root for more information.
|
---|
5 | */
|
---|
6 | #endregion
|
---|
7 |
|
---|
8 | using System;
|
---|
9 | using System.Collections.Generic;
|
---|
10 | using System.Linq;
|
---|
11 | using System.Text;
|
---|
12 |
|
---|
13 | namespace SimSharp {
|
---|
14 | /// <summary>
|
---|
15 | /// This class calculates weighted statistics of a time series variable.
|
---|
16 | /// The weight is given as the duration of the variable's value.
|
---|
17 | ///
|
---|
18 | /// It is typically used to calculate utilization of resources or inventory levels.
|
---|
19 | /// </summary>
|
---|
20 | /// <remarks>
|
---|
21 | /// The monitor updates the statistics online, except for <see cref="GetMedian"/>
|
---|
22 | /// which is calculated given the collected data whenever it is called.
|
---|
23 | /// When the monitor was initialized with collect = false in the constructor,
|
---|
24 | /// median and other percentiles cannot be computed (double.NaN is returned).
|
---|
25 | /// Also to print a histogram requires that the monitor was initialized to collect
|
---|
26 | /// all the changes to the variable's value.
|
---|
27 | ///
|
---|
28 | /// Collecting the data naturally incurs some memory overhead.
|
---|
29 | /// </remarks>
|
---|
30 | public sealed class TimeSeriesMonitor : ITimeSeriesMonitor {
|
---|
31 | private readonly Simulation env;
|
---|
32 | /// <summary>
|
---|
33 | /// Can only be set in the constructor.
|
---|
34 | /// When it is true, median and percentiles can be computed and a
|
---|
35 | /// histogram can be printed. In addition <see cref="Series"/>
|
---|
36 | /// may return all the remembered values for further processing.
|
---|
37 | /// </summary>
|
---|
38 | public bool Collect { get; }
|
---|
39 |
|
---|
40 | private bool active;
|
---|
41 | /// <summary>
|
---|
42 | /// The monitor can be set to suppress updates. When it is set
|
---|
43 | /// to false, the statistics, except for <see cref="Current"/>
|
---|
44 | /// will not be updated and new samples are ignored.
|
---|
45 | /// When set to true it will resume updates from the current
|
---|
46 | /// simulation time onward.
|
---|
47 | /// </summary>
|
---|
48 | public bool Active {
|
---|
49 | get { return active; }
|
---|
50 | set {
|
---|
51 | if (active == value) return;
|
---|
52 | if (firstSampleObserved && lastUpdateTime < env.NowD) {
|
---|
53 | if (value) {
|
---|
54 | // enable it
|
---|
55 | if (Current < Min) Min = Current;
|
---|
56 | if (Current > Max) Max = Current;
|
---|
57 |
|
---|
58 | lastUpdateTime = env.NowD;
|
---|
59 | series.Add(new Entry { Date = env.NowD, Level = Current });
|
---|
60 | } else if (!value) {
|
---|
61 | // disable it
|
---|
62 | OnlineUpdate();
|
---|
63 | }
|
---|
64 | active = value;
|
---|
65 | OnUpdated();
|
---|
66 | } else active = value;
|
---|
67 | }
|
---|
68 | }
|
---|
69 |
|
---|
70 | /// <summary>
|
---|
71 | /// The name of the variable that is being monitored.
|
---|
72 | /// Used for output in <see cref="Summarize(bool, int, double?, double?)"/>.
|
---|
73 | /// </summary>
|
---|
74 | public string Name { get; set; }
|
---|
75 |
|
---|
76 | public double TotalTimeD { get; private set; }
|
---|
77 | public TimeSpan TotalTime { get { return env.ToTimeSpan(TotalTimeD); } }
|
---|
78 |
|
---|
79 | public double Min { get; private set; }
|
---|
80 | public double Max { get; private set; }
|
---|
81 | public double Area {
|
---|
82 | get {
|
---|
83 | if (!UpToDate) OnlineUpdate();
|
---|
84 | return area;
|
---|
85 | }
|
---|
86 | private set => area = value;
|
---|
87 | }
|
---|
88 | double INumericMonitor.Sum { get { return Area; } }
|
---|
89 | public double Mean {
|
---|
90 | get {
|
---|
91 | if (!UpToDate) OnlineUpdate();
|
---|
92 | return mean;
|
---|
93 | }
|
---|
94 | private set => mean = value;
|
---|
95 | }
|
---|
96 | public double StdDev { get { return Math.Sqrt(Variance); } }
|
---|
97 | public double Variance {
|
---|
98 | get {
|
---|
99 | if (!UpToDate) OnlineUpdate();
|
---|
100 | return (TotalTimeD > 0) ? variance / TotalTimeD : 0.0;
|
---|
101 | }
|
---|
102 | }
|
---|
103 | public double Current { get; private set; }
|
---|
104 | double INumericMonitor.Last { get { return Current; } }
|
---|
105 |
|
---|
106 | private bool UpToDate { get { return !Active || env.NowD == lastUpdateTime; } }
|
---|
107 |
|
---|
108 | private double lastUpdateTime;
|
---|
109 | private double variance;
|
---|
110 |
|
---|
111 | private bool firstSampleObserved;
|
---|
112 | private double area;
|
---|
113 | private double mean;
|
---|
114 |
|
---|
115 | private List<Entry> series;
|
---|
116 | /// <summary>
|
---|
117 | /// Returns the list of collected values, or an empty enumerable
|
---|
118 | /// when <see cref="Collect"/> was initialized to false.
|
---|
119 | /// </summary>
|
---|
120 | public IEnumerable<Entry> Series {
|
---|
121 | get {
|
---|
122 | if (!UpToDate) OnlineUpdate();
|
---|
123 | return series != null ? series.AsEnumerable() : Enumerable.Empty<Entry>();
|
---|
124 | }
|
---|
125 | }
|
---|
126 |
|
---|
127 | /// <summary>
|
---|
128 | /// Calls <see cref="GetPercentile(double)"/>.
|
---|
129 | /// </summary>
|
---|
130 | /// <remarks>
|
---|
131 | /// Median can only be computed when the monitor was initialized to collect the data.
|
---|
132 | ///
|
---|
133 | /// The data is preprocessed on every call, the runtime complexity of this method is therefore O(n * log(n)).
|
---|
134 | /// </remarks>
|
---|
135 | /// <returns>The median (50th percentile) of the time series.</returns>
|
---|
136 | public double GetMedian() {
|
---|
137 | return GetPercentile(0.5);
|
---|
138 | }
|
---|
139 |
|
---|
140 | /// <summary>
|
---|
141 | /// Calculates the weighted p-percentile of the sampled levels (duration is the weight).
|
---|
142 | /// </summary>
|
---|
143 | /// <remarks>
|
---|
144 | /// Percentiles can only be computed when the monitor was initialized to collect the data.
|
---|
145 | ///
|
---|
146 | /// The data is preprocessed on every call, the runtime complexity of this method is therefore O(n * log(n)).
|
---|
147 | /// </remarks>
|
---|
148 | /// <exception cref="ArgumentException">Thrown when <paramref name="p"/> is outside the valid range.</exception>
|
---|
149 | /// <param name="p">The percentile has to be in the range [0;1].</param>
|
---|
150 | /// <returns>The respective percentile of the time series.</returns>
|
---|
151 | public double GetPercentile(double p) {
|
---|
152 | if (p < 0 || p > 1) throw new ArgumentException("Percentile must be between 0 and 1", "p");
|
---|
153 | if (!Collect) return double.NaN;
|
---|
154 | return GetPercentile(series, p);
|
---|
155 | }
|
---|
156 | private static double GetPercentile(IList<Entry> s, double p) {
|
---|
157 | var seq = Cumulate(s).ToList();
|
---|
158 | if (seq.Count == 0) return double.NaN;
|
---|
159 | var total = seq.Last().CumulatedDuration;
|
---|
160 | var n = total * p;
|
---|
161 | var ilower = 0;
|
---|
162 | while (ilower < seq.Count) {
|
---|
163 | if (seq[ilower].CumulatedDuration >= n) break;
|
---|
164 | ilower++;
|
---|
165 | }
|
---|
166 | if (seq[ilower].CumulatedDuration == n) {
|
---|
167 | if (ilower < seq.Count - 1)
|
---|
168 | return (seq[ilower].Level + seq[ilower + 1].Level) / 2.0;
|
---|
169 | }
|
---|
170 | return seq[ilower].Level;
|
---|
171 | }
|
---|
172 |
|
---|
173 | private static IEnumerable<LevelDuration> Cumulate(IList<Entry> s) {
|
---|
174 | var totalDuration = 0.0;
|
---|
175 | // The last entry will always be skipped as it has Duration == 0
|
---|
176 | foreach (var g in s.Where(x => x.Duration > 0).GroupBy(x => x.Level, x => x.Duration).OrderBy(x => x.Key)) {
|
---|
177 | var duration = g.Sum();
|
---|
178 | totalDuration += duration;
|
---|
179 | yield return new LevelDuration { Level = g.Key, Duration = duration, CumulatedDuration = totalDuration };
|
---|
180 | }
|
---|
181 | }
|
---|
182 |
|
---|
183 | public TimeSeriesMonitor(Simulation env, string name = null, bool collect = false) {
|
---|
184 | this.env = env;
|
---|
185 | active = true;
|
---|
186 | Name = name;
|
---|
187 | Collect = collect;
|
---|
188 | lastUpdateTime = env.NowD;
|
---|
189 | if (Collect) series = new List<Entry>(64);
|
---|
190 | }
|
---|
191 | public TimeSeriesMonitor(Simulation env, double initial, string name = null, bool collect = false) {
|
---|
192 | this.env = env;
|
---|
193 | active = true;
|
---|
194 | Name = name;
|
---|
195 | Collect = collect;
|
---|
196 | lastUpdateTime = env.NowD;
|
---|
197 | firstSampleObserved = true;
|
---|
198 | Current = Min = Max = mean = initial;
|
---|
199 | if (Collect) series = new List<Entry>(64) { new Entry { Date = env.NowD, Level = initial } };
|
---|
200 | }
|
---|
201 |
|
---|
202 | public void Reset() {
|
---|
203 | TotalTimeD = 0;
|
---|
204 | Current = Min = Max = area = mean = 0;
|
---|
205 | if (Collect) series.Clear();
|
---|
206 | variance = 0;
|
---|
207 | firstSampleObserved = false;
|
---|
208 | lastUpdateTime = env.NowD;
|
---|
209 | }
|
---|
210 |
|
---|
211 | public void Reset(double initial) {
|
---|
212 | TotalTimeD = 0;
|
---|
213 | Current = Min = Max = mean = initial;
|
---|
214 | if (Collect) {
|
---|
215 | series.Clear();
|
---|
216 | series.Add(new Entry { Date = env.NowD, Level = initial });
|
---|
217 | }
|
---|
218 | area = 0;
|
---|
219 | variance = 0;
|
---|
220 | firstSampleObserved = true;
|
---|
221 | lastUpdateTime = env.NowD;
|
---|
222 | }
|
---|
223 |
|
---|
224 | public void Increase(double value = 1) {
|
---|
225 | UpdateTo(Current + value);
|
---|
226 | }
|
---|
227 |
|
---|
228 | public void Decrease(double value = 1) {
|
---|
229 | UpdateTo(Current - value);
|
---|
230 | }
|
---|
231 |
|
---|
232 | public void UpdateTo(double value) {
|
---|
233 | if (!Active) {
|
---|
234 | Current = value;
|
---|
235 | return;
|
---|
236 | }
|
---|
237 |
|
---|
238 | if (!firstSampleObserved) {
|
---|
239 | Min = Max = mean = value;
|
---|
240 | firstSampleObserved = true;
|
---|
241 | lastUpdateTime = env.NowD;
|
---|
242 | if (Collect) series.Add(new Entry { Date = env.NowD, Level = value });
|
---|
243 | Current = value;
|
---|
244 | } else {
|
---|
245 | if (value < Min) Min = value;
|
---|
246 | if (value > Max) Max = value;
|
---|
247 |
|
---|
248 | OnlineUpdate();
|
---|
249 |
|
---|
250 | if (Current != value) {
|
---|
251 | if (Collect) series.Add(new Entry { Date = env.NowD, Level = value });
|
---|
252 | Current = value;
|
---|
253 | }
|
---|
254 | }
|
---|
255 |
|
---|
256 | OnUpdated();
|
---|
257 | }
|
---|
258 |
|
---|
259 | private void OnlineUpdate() {
|
---|
260 | var duration = env.NowD - lastUpdateTime;
|
---|
261 | if (duration > 0) {
|
---|
262 | if (Collect) {
|
---|
263 | var prevIdx = series.Count - 1;
|
---|
264 | var prev = series[prevIdx];
|
---|
265 | prev.Duration = env.NowD - prev.Date;
|
---|
266 | series[prevIdx] = prev;
|
---|
267 | }
|
---|
268 | area += (Current * duration);
|
---|
269 | var oldMean = mean;
|
---|
270 | mean = oldMean + (Current - oldMean) * duration / (duration + TotalTimeD);
|
---|
271 | variance = variance + (Current - oldMean) * (Current - mean) * duration;
|
---|
272 | TotalTimeD += duration;
|
---|
273 | }
|
---|
274 | lastUpdateTime = env.NowD;
|
---|
275 | }
|
---|
276 |
|
---|
277 | public event EventHandler Updated;
|
---|
278 | private void OnUpdated() {
|
---|
279 | Updated?.Invoke(this, EventArgs.Empty);
|
---|
280 | }
|
---|
281 |
|
---|
282 | string IMonitor.Summarize() {
|
---|
283 | return Summarize();
|
---|
284 | }
|
---|
285 |
|
---|
286 | /// <summary>
|
---|
287 | /// Provides a summary of the statistics in a certain format.
|
---|
288 | /// If the monitor is configured to collect data, it may also print a histogram.
|
---|
289 | /// </summary>
|
---|
290 | /// <param name="withHistogram">Whether to suppress the histogram.
|
---|
291 | /// This is only effective if <see cref="Collect"/> was set to true, otherwise
|
---|
292 | /// the data to produce the histogram is not available in the first place.</param>
|
---|
293 | /// <param name="maxBins">The maximum number of bins that should be used.
|
---|
294 | /// Note that the bin width and thus the number of bins is also governed by
|
---|
295 | /// <paramref name="binWidth"/> if it is defined.
|
---|
296 | /// This is only effective if <see cref="Collect"/> and <paramref name="withHistogram"/>
|
---|
297 | /// was set to true, otherwise the data to produce the histogram is not available
|
---|
298 | /// in the first place.</param>
|
---|
299 | /// <param name="histMin">The minimum for the histogram to start at or the sample
|
---|
300 | /// minimum in case the default (null) is given.
|
---|
301 | /// This is only effective if <see cref="Collect"/> and <paramref name="withHistogram"/>
|
---|
302 | /// was set to true, otherwise the data to produce the histogram is not available
|
---|
303 | /// in the first place.</param>
|
---|
304 | /// <param name="binWidth">The interval for the bins of the histogram or the
|
---|
305 | /// range (<see cref="Max"/> - <see cref="Min"/>) divided by the number of bins
|
---|
306 | /// (<paramref name="maxBins"/>) in case the default value (null) is given.
|
---|
307 | /// This is only effective if <see cref="Collect"/> and <paramref name="withHistogram"/>
|
---|
308 | /// was set to true, otherwise the data to produce the histogram is not available
|
---|
309 | /// in the first place.</param>
|
---|
310 | /// <returns>A formatted string that provides a summary of the statistics.</returns>
|
---|
311 | public string Summarize(bool withHistogram = true, int maxBins = 20, double? histMin = null, double? binWidth = null) {
|
---|
312 | var nozero = Collect ? series.Where(x => x.Level != 0 && x.Duration > 0).ToList() : new List<Entry>();
|
---|
313 | var nozeromin = nozero.Count > 0 ? nozero.Min(x => x.Level) : double.NaN;
|
---|
314 | var nozeromax = nozero.Count > 0 ? nozero.Max(x => x.Level) : double.NaN;
|
---|
315 | var nozeroduration = Collect ? nozero.Sum(x => x.Duration) : double.NaN;
|
---|
316 | var nozeromean = nozero.Count > 1 ? nozero.Sum(x => x.Level * x.Duration / nozeroduration) : double.NaN;
|
---|
317 | var nozerostdev = nozero.Count > 2 ? Math.Sqrt(nozero.Sum(x => x.Duration * (x.Level - nozeromean) * (x.Level - nozeromean)) / nozeroduration) : double.NaN;
|
---|
318 | var sb = new StringBuilder();
|
---|
319 | sb.Append("Time series statistics");
|
---|
320 | if (!string.IsNullOrEmpty(Name))
|
---|
321 | sb.Append(" of " + Name);
|
---|
322 | sb.AppendLine();
|
---|
323 | sb.AppendLine(" all excl.zero zero ");
|
---|
324 | sb.AppendLine("--------------- --------------- --------------- ---------------");
|
---|
325 | sb.AppendLine(string.Format("{0,15} {1,15} {2,15} {3,15}", "Duration", Formatter.Format15(TotalTimeD), Formatter.Format15(nozeroduration), Formatter.Format15(TotalTimeD - nozeroduration)));
|
---|
326 | sb.AppendLine(string.Format("{0,15} {1,15} {2,15}", "Mean", Formatter.Format15(Mean), Formatter.Format15(nozeromean)));
|
---|
327 | sb.AppendLine(string.Format("{0,15} {1,15} {2,15}", "Std.dev", Formatter.Format15(StdDev), Formatter.Format15(nozerostdev)));
|
---|
328 | sb.AppendLine();
|
---|
329 | sb.AppendLine(string.Format("{0,15} {1,15} {2,15}", "Minimum", Formatter.Format15(Min), Formatter.Format15(nozeromin)));
|
---|
330 | if (Collect) {
|
---|
331 | sb.AppendLine(string.Format("{0,15} {1,15} {2,15}", "Percentile-5%", Formatter.Format15(GetPercentile(0.05)), Formatter.Format15(GetPercentile(nozero, 0.05))));
|
---|
332 | sb.AppendLine(string.Format("{0,15} {1,15} {2,15}", "Median", Formatter.Format15(GetMedian()), Formatter.Format15(GetPercentile(nozero, 0.5))));
|
---|
333 | sb.AppendLine(string.Format("{0,15} {1,15} {2,15}", "Percentile-95%", Formatter.Format15(GetPercentile(0.95)), Formatter.Format15(GetPercentile(nozero, 0.95))));
|
---|
334 | }
|
---|
335 | sb.AppendLine(string.Format("{0,15} {1,15} {2,15}", "Maximum", Formatter.Format15(Max), Formatter.Format15(nozeromax)));
|
---|
336 |
|
---|
337 | if (Collect && withHistogram) {
|
---|
338 | var histData = Cumulate(series);
|
---|
339 | sb.AppendLine();
|
---|
340 | sb.AppendLine("Histogram");
|
---|
341 | sb.AppendLine("<= duration % cum% ");
|
---|
342 | sb.AppendLine("--------------- --------------- ----- ------");
|
---|
343 | var totStars = 0;
|
---|
344 | var iter = histData.GetEnumerator();
|
---|
345 | var cumul = 0.0;
|
---|
346 | if (!iter.MoveNext()) {
|
---|
347 | sb.AppendLine("no data");
|
---|
348 | } else {
|
---|
349 | var kvp = iter.Current;
|
---|
350 | var moreData = true;
|
---|
351 | for (var bin = 0; bin <= maxBins; bin++) {
|
---|
352 | var next = (histMin ?? Min) + bin * (binWidth ?? (Max - Min) / 20.0);
|
---|
353 | var dur = 0.0;
|
---|
354 | var prob = 0.0;
|
---|
355 | while (moreData && (kvp.Level <= next || bin == maxBins)) {
|
---|
356 | dur += kvp.Duration;
|
---|
357 | prob += kvp.Duration / TotalTimeD;
|
---|
358 | cumul = kvp.CumulatedDuration / TotalTimeD;
|
---|
359 | moreData = iter.MoveNext();
|
---|
360 | kvp = iter.Current;
|
---|
361 | }
|
---|
362 | var probstars = (int)Math.Round(100 * prob / 2);
|
---|
363 | var cumulstars = (int)Math.Round(100 * cumul / 2);
|
---|
364 | var numstars = probstars;
|
---|
365 | if (numstars + totStars < cumulstars) numstars++;
|
---|
366 | var stars = string.Join("", Enumerable.Repeat("*", numstars));
|
---|
367 | totStars += numstars;
|
---|
368 | var cumulbar = "|".PadLeft(totStars + 1 - numstars);
|
---|
369 | sb.AppendLine(string.Format("{0,15} {1,15} {2,5:F1} {3,5:F1} {4}{5}",
|
---|
370 | (!moreData && next < Max) ? "inf" : Formatter.Format15(next),
|
---|
371 | Formatter.Format15(dur), prob * 100, cumul * 100, stars, cumulbar));
|
---|
372 | if (!moreData) break;
|
---|
373 | }
|
---|
374 | }
|
---|
375 | }
|
---|
376 | return sb.ToString();
|
---|
377 | }
|
---|
378 |
|
---|
379 | public override string ToString() {
|
---|
380 | return Summarize();
|
---|
381 | }
|
---|
382 |
|
---|
383 | public struct Entry {
|
---|
384 | public double Date;
|
---|
385 | public double Duration;
|
---|
386 | public double Level;
|
---|
387 | }
|
---|
388 |
|
---|
389 | private struct LevelDuration {
|
---|
390 | public double Level;
|
---|
391 | public double Duration;
|
---|
392 | public double CumulatedDuration;
|
---|
393 | }
|
---|
394 | }
|
---|
395 |
|
---|
396 | [Obsolete("Use the class TimeSeriesMonitor instead.")]
|
---|
397 | public sealed class ContinuousStatistics {
|
---|
398 | private readonly Simulation env;
|
---|
399 |
|
---|
400 | public int Count { get; private set; }
|
---|
401 | public double TotalTimeD { get; private set; }
|
---|
402 | public TimeSpan TotalTime { get { return env.ToTimeSpan(TotalTimeD); } }
|
---|
403 |
|
---|
404 | public double Min { get; private set; }
|
---|
405 | public double Max { get; private set; }
|
---|
406 | public double Area { get; private set; }
|
---|
407 | public double Mean { get; private set; }
|
---|
408 | public double StdDev { get { return Math.Sqrt(Variance); } }
|
---|
409 | public double Variance { get { return (TotalTimeD > 0) ? variance / TotalTimeD : 0.0; } }
|
---|
410 |
|
---|
411 | private double lastUpdateTime;
|
---|
412 | private double lastValue;
|
---|
413 | private double variance;
|
---|
414 |
|
---|
415 | private bool firstSample;
|
---|
416 |
|
---|
417 |
|
---|
418 | public ContinuousStatistics(Simulation env) {
|
---|
419 | this.env = env;
|
---|
420 | lastUpdateTime = env.NowD;
|
---|
421 | }
|
---|
422 |
|
---|
423 | public void Update(double value) {
|
---|
424 | Count++;
|
---|
425 |
|
---|
426 | if (!firstSample) {
|
---|
427 | Min = Max = Mean = value;
|
---|
428 | firstSample = true;
|
---|
429 | } else {
|
---|
430 | if (value < Min) Min = value;
|
---|
431 | if (value > Max) Max = value;
|
---|
432 |
|
---|
433 | var duration = env.NowD - lastUpdateTime;
|
---|
434 | if (duration > 0) {
|
---|
435 | Area += (lastValue * duration);
|
---|
436 | var oldMean = Mean;
|
---|
437 | Mean = oldMean + (lastValue - oldMean) * duration / (duration + TotalTimeD);
|
---|
438 | variance = variance + (lastValue - oldMean) * (lastValue - Mean) * duration;
|
---|
439 | TotalTimeD += duration;
|
---|
440 | }
|
---|
441 | }
|
---|
442 |
|
---|
443 | lastUpdateTime = env.NowD;
|
---|
444 | lastValue = value;
|
---|
445 | }
|
---|
446 | }
|
---|
447 | }
|
---|