Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/HeuristicLab.ExtLibs/HeuristicLab.SimSharp/3.3.1/SimSharp-3.3.1/Core/Environment.cs @ 17487

Last change on this file since 17487 was 17487, checked in by abeham, 4 years ago

#3065: update Sim# to 3.3.1

File size: 48.6 KB
Line 
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
8using System;
9using System.Collections.Generic;
10using System.Diagnostics;
11using System.IO;
12using System.Threading;
13using System.Threading.Tasks;
14
15namespace SimSharp {
16  /// <summary>
17  /// Simulation hold the event queues, schedule and process events.
18  /// </summary>
19  /// <remarks>
20  /// This class is not thread-safe against manipulation of the event queue. If you supply a termination
21  /// event that is set outside the simulation thread, please use the <see cref="ThreadSafeSimulation"/> environment.
22  ///
23  /// For most purposes <see cref="Simulation"/> is however the better and faster choice.
24  /// </remarks>
25  public class Simulation {
26    private const int InitialMaxEvents = 1024;
27
28    /// <summary>
29    /// Describes the number of seconds that a logical step of 1 in the *D-API takes.
30    /// </summary>
31    protected double DefaultTimeStepSeconds { get; private set; }
32
33    /// <summary>
34    /// Calculates the logical date of the simulation by the amount of default steps
35    /// that have passed.
36    /// </summary>
37    public double NowD {
38      get { return (Now - StartDate).TotalSeconds / DefaultTimeStepSeconds; }
39    }
40
41    private DateTime now;
42    /// <summary>
43    /// The current simulation time as a calendar date.
44    /// </summary>
45    public virtual DateTime Now { get => now; protected set => now = value; }
46
47    /// <summary>
48    /// The calendar date when the simulation started. This defaults to 1970-1-1 if
49    /// no other date has been specified in the overloaded constructor.
50    /// </summary>
51    public DateTime StartDate { get; protected set; }
52
53    /// <summary>
54    /// The random number generator that is to be used in all events in
55    /// order to produce reproducible results.
56    /// </summary>
57    protected IRandom Random { get; set; }
58
59    protected EventQueue ScheduleQ;
60    public Process ActiveProcess { get; set; }
61
62    public TextWriter Logger { get; set; }
63    public int ProcessedEvents { get; protected set; }
64
65    public Simulation() : this(new DateTime(1970, 1, 1)) { }
66    public Simulation(TimeSpan? defaultStep) : this(new DateTime(1970, 1, 1), defaultStep) { }
67    public Simulation(int randomSeed, TimeSpan? defaultStep = null) : this(new DateTime(1970, 1, 1), randomSeed, defaultStep) { }
68    public Simulation(DateTime initialDateTime, TimeSpan? defaultStep = null) : this(new PcgRandom(), initialDateTime, defaultStep) { }
69    public Simulation(DateTime initialDateTime, int randomSeed, TimeSpan? defaultStep = null) : this(new PcgRandom(randomSeed), initialDateTime, defaultStep) { }
70    public Simulation(IRandom random, DateTime initialDateTime, TimeSpan? defaultStep = null) {
71      DefaultTimeStepSeconds = (defaultStep ?? TimeSpan.FromSeconds(1)).Duration().TotalSeconds;
72      StartDate = initialDateTime;
73      Now = initialDateTime;
74      Random = random;
75      ScheduleQ = new EventQueue(InitialMaxEvents);
76      Logger = Console.Out;
77    }
78
79    public double ToDouble(TimeSpan span) {
80      return span.TotalSeconds / DefaultTimeStepSeconds;
81    }
82
83    public TimeSpan ToTimeSpan(double span) {
84      return TimeSpan.FromSeconds(DefaultTimeStepSeconds * span);
85    }
86
87    /// <summary>
88    /// Creates a new process from an event generator. The process is automatically
89    /// scheduled to be started at the current simulation time.
90    /// </summary>
91    /// <param name="generator">The generator function that represents the process.</param>
92    /// <param name="priority">The priority to rank events at the same time (smaller value = higher priority).</param>
93    /// <returns>The scheduled process that was created.</returns>
94    public Process Process(IEnumerable<Event> generator, int priority = 0) {
95      return new Process(this, generator, priority);
96    }
97
98    /// <summary>
99    /// Creates and returns a new timeout.
100    /// </summary>
101    /// <param name="delay">The time after which the timeout is fired.</param>
102    /// <param name="priority">The priority to rank events at the same time (smaller value = higher priority).</param>
103    /// <returns>The scheduled timeout event that was created.</returns>
104    public Timeout TimeoutD(double delay, int priority = 0) {
105      return Timeout(TimeSpan.FromSeconds(DefaultTimeStepSeconds * delay), priority);
106    }
107
108    /// <summary>
109    /// Creates and returns a new timeout.
110    /// </summary>
111    /// <param name="delay">The time after which the timeout is fired.</param>
112    /// <param name="priority">The priority to rank events at the same time (smaller value = higher priority).</param>
113    /// <returns>The scheduled timeout event that was created.</returns>
114    public Timeout Timeout(TimeSpan delay, int priority = 0) {
115      return new Timeout(this, delay, priority: priority);
116    }
117
118    public virtual void Reset(int randomSeed) {
119      ProcessedEvents = 0;
120      Now = StartDate;
121      Random = new PcgRandom(randomSeed);
122      ScheduleQ = new EventQueue(InitialMaxEvents);
123      useSpareNormal = false;
124    }
125
126    public virtual void ScheduleD(double delay, Event @event) {
127      Schedule(TimeSpan.FromSeconds(DefaultTimeStepSeconds * delay), @event);
128    }
129
130    /// <summary>
131    /// Schedules an event to occur at the same simulation time as the call was made.
132    /// </summary>
133    /// <param name="event">The event that should be scheduled.</param>
134    /// <param name="priority">The priority to rank events at the same time (smaller value = higher priority).</param>
135    public virtual void Schedule(Event @event, int priority = 0) {
136      DoSchedule(Now, @event, priority);
137    }
138
139    /// <summary>
140    /// Schedules an event to occur after a certain (positive) delay.
141    /// </summary>
142    /// <exception cref="ArgumentException">
143    /// Thrown when <paramref name="delay"/> is negative.
144    /// </exception>
145    /// <param name="delay">The (positive) delay after which the event should be fired.</param>
146    /// <param name="event">The event that should be scheduled.</param>
147    /// <param name="priority">The priority to rank events at the same time (smaller value = higher priority).</param>
148    public virtual void Schedule(TimeSpan delay, Event @event, int priority = 0) {
149      if (delay < TimeSpan.Zero)
150        throw new ArgumentException("Negative delays are not allowed in Schedule(TimeSpan, Event).");
151      var eventTime = Now + delay;
152      DoSchedule(eventTime, @event, priority);
153    }
154
155    protected virtual EventQueueNode DoSchedule(DateTime date, Event @event, int priority = 0) {
156      if (ScheduleQ.MaxSize == ScheduleQ.Count) {
157        // the capacity has to be adjusted, there are more events in the queue than anticipated
158        var oldSchedule = ScheduleQ;
159        ScheduleQ = new EventQueue(ScheduleQ.MaxSize * 2);
160        foreach (var e in oldSchedule) ScheduleQ.Enqueue(e.PrimaryPriority, e.Event, e.SecondaryPriority);
161      }
162      return ScheduleQ.Enqueue(date, @event, priority);
163    }
164
165    public virtual object RunD(double? until = null) {
166      if (!until.HasValue) return Run();
167      return Run(Now + TimeSpan.FromSeconds(DefaultTimeStepSeconds * until.Value));
168    }
169
170    public virtual object Run(TimeSpan span) {
171      return Run(Now + span);
172    }
173
174    public virtual object Run(DateTime until) {
175      if (until <= Now) throw new InvalidOperationException("Simulation end date must lie in the future.");
176      var stopEvent = new Event(this);
177      var node = DoSchedule(until, stopEvent);
178      // stop event is always the first to execute at the given time
179      node.InsertionIndex = -1;
180      ScheduleQ.OnNodeUpdated(node);
181      return Run(stopEvent);
182    }
183
184    protected CancellationTokenSource _stop = null;
185    /// <summary>
186    /// Run until a certain event is processed.
187    /// </summary>
188    /// <remarks>
189    /// This simulation environment is not thread-safe, thus triggering this event outside the environment
190    /// leads to potential race conditions. Please use the <see cref="ThreadSafeSimulation"/> environment in case you
191    /// require this functionality. Note that the performance of <see cref="ThreadSafeSimulation"/> is lower due to locking.
192    ///
193    /// For real-time based termination, you can also call <see cref="StopAsync"/> which sets a flag indicating the simulation
194    /// to stop before processing the next event.
195    /// </remarks>
196    /// <param name="stopEvent">The event that stops the simulation.</param>
197    /// <returns></returns>
198    public virtual object Run(Event stopEvent = null) {
199      _stop = new CancellationTokenSource();
200      if (stopEvent != null) {
201        if (stopEvent.IsProcessed) {
202          return stopEvent.Value;
203        }
204        stopEvent.AddCallback(StopSimulation);
205      }
206      OnRunStarted();
207      try {
208        var stop = ScheduleQ.Count == 0 || _stop.IsCancellationRequested;
209        while (!stop) {
210          Step();
211          stop = ScheduleQ.Count == 0 || _stop.IsCancellationRequested;
212        }
213      } catch (StopSimulationException e) { OnRunFinished(); return e.Value; }
214      OnRunFinished();
215      if (stopEvent == null) return null;
216      if (!_stop.IsCancellationRequested && !stopEvent.IsTriggered) throw new InvalidOperationException("No scheduled events left but \"until\" event was not triggered.");
217      return stopEvent.Value;
218    }
219
220    public virtual void StopAsync() {
221      _stop?.Cancel();
222    }
223
224    public event EventHandler RunStarted;
225    protected void OnRunStarted() {
226      RunStarted?.Invoke(this, EventArgs.Empty);
227    }
228
229    public event EventHandler RunFinished;
230    protected void OnRunFinished() {
231      RunFinished?.Invoke(this, EventArgs.Empty);
232    }
233
234    /// <summary>
235    /// Performs a single step of the simulation, i.e. process a single event
236    /// </summary>
237    /// <remarks>
238    /// This method is not thread-safe
239    /// </remarks>
240    public virtual void Step() {
241      Event evt;
242      var next = ScheduleQ.Dequeue();
243      Now = next.PrimaryPriority;
244      evt = next.Event;
245      evt.Process();
246      ProcessedEvents++;
247    }
248
249    /// <summary>
250    /// Peeks at the time of the next event in terms of the defined step
251    /// </summary>
252    /// <remarks>
253    /// This method is not thread-safe
254    /// </remarks>
255    public virtual double PeekD() {
256      if (ScheduleQ.Count == 0) return double.MaxValue;
257      return (Peek() - StartDate).TotalSeconds / DefaultTimeStepSeconds;
258    }
259
260    /// <summary>
261    /// Peeks at the time of the next event
262    /// </summary>
263    /// <remarks>
264    /// This method is not thread-safe
265    /// </remarks>
266    public virtual DateTime Peek() {
267      return ScheduleQ.Count > 0 ? ScheduleQ.First.PrimaryPriority : DateTime.MaxValue;
268    }
269
270    protected virtual void StopSimulation(Event @event) {
271      throw new StopSimulationException(@event.Value);
272    }
273
274    public virtual void Log(string message, params object[] args) {
275      if (Logger != null)
276        Logger.WriteLine(message, args);
277    }
278
279    #region Random number distributions
280    public double RandUniform(IRandom random, double a, double b) {
281      return a + (b - a) * random.NextDouble();
282    }
283    public double RandUniform(double a, double b) {
284      return RandUniform(Random, a, b);
285    }
286
287    public TimeSpan RandUniform(IRandom random, TimeSpan a, TimeSpan b) {
288      return TimeSpan.FromSeconds(RandUniform(random, a.TotalSeconds, b.TotalSeconds));
289    }
290    public TimeSpan RandUniform(TimeSpan a, TimeSpan b) {
291      return RandUniform(Random, a, b);
292    }
293    public double RandTriangular(IRandom random, double low, double high) {
294      var u = random.NextDouble();
295      if (u > 0.5)
296        return high + (low - high) * Math.Sqrt(((1.0 - u) / 2));
297      return low + (high - low) * Math.Sqrt(u / 2);
298    }
299    public double RandTriangular(double low, double high) {
300      return RandTriangular(Random, low, high);
301    }
302
303    public TimeSpan RandTriangular(IRandom random, TimeSpan low, TimeSpan high) {
304      return TimeSpan.FromSeconds(RandTriangular(random, low.TotalSeconds, high.TotalSeconds));
305    }
306    public TimeSpan RandTriangular(TimeSpan low, TimeSpan high) {
307      return RandTriangular(Random, low, high);
308    }
309
310    public double RandTriangular(IRandom random, double low, double high, double mode) {
311      var u = random.NextDouble();
312      var c = (mode - low) / (high - low);
313      if (u > c)
314        return high + (low - high) * Math.Sqrt(((1.0 - u) * (1.0 - c)));
315      return low + (high - low) * Math.Sqrt(u * c);
316    }
317    public double RandTriangular(double low, double high, double mode) {
318      return RandTriangular(Random, low, high, mode);
319    }
320
321    public TimeSpan RandTriangular(IRandom random, TimeSpan low, TimeSpan high, TimeSpan mode) {
322      return TimeSpan.FromSeconds(RandTriangular(random, low.TotalSeconds, high.TotalSeconds, mode.TotalSeconds));
323    }
324    public TimeSpan RandTriangular(TimeSpan low, TimeSpan high, TimeSpan mode) {
325      return RandTriangular(Random, low, high, mode);
326    }
327
328    /// <summary>
329    /// Returns a number that is exponentially distributed given a certain mean.
330    /// </summary>
331    /// <remarks>
332    /// Unlike in other APIs here the mean should be given and not the lambda parameter.
333    /// </remarks>
334    /// <param name="random">The random number generator to use.</param>
335    /// <param name="mean">The mean(!) of the distribution is 1 / lambda.</param>
336    /// <returns>A number that is exponentially distributed</returns>
337    public double RandExponential(IRandom random, double mean) {
338      return -Math.Log(1 - random.NextDouble()) * mean;
339    }
340    /// <summary>
341    /// Returns a number that is exponentially distributed given a certain mean.
342    /// </summary>
343    /// <remarks>
344    /// Unlike in other APIs here the mean should be given and not the lambda parameter.
345    /// </remarks>
346    /// <param name="mean">The mean(!) of the distribution is 1 / lambda.</param>
347    /// <returns>A number that is exponentially distributed</returns>
348    public double RandExponential(double mean) {
349      return RandExponential(Random, mean);
350    }
351
352    /// <summary>
353    /// Returns a timespan that is exponentially distributed given a certain mean.
354    /// </summary>
355    /// <remarks>
356    /// Unlike in other APIs here the mean should be given and not the lambda parameter.
357    /// </remarks>
358    /// <param name="random">The random number generator to use.</param>
359    /// <param name="mean">The mean(!) of the distribution is 1 / lambda.</param>
360    /// <returns>A number that is exponentially distributed</returns>
361    public TimeSpan RandExponential(IRandom random, TimeSpan mean) {
362      return TimeSpan.FromSeconds(RandExponential(random, mean.TotalSeconds));
363    }
364    /// <summary>
365    /// Returns a timespan that is exponentially distributed given a certain mean.
366    /// </summary>
367    /// <remarks>
368    /// Unlike in other APIs here the mean should be given and not the lambda parameter.
369    /// </remarks>
370    /// <param name="mean">The mean(!) of the distribution is 1 / lambda.</param>
371    /// <returns>A number that is exponentially distributed</returns>
372    public TimeSpan RandExponential(TimeSpan mean) {
373      return RandExponential(Random, mean);
374    }
375
376    private bool useSpareNormal = false;
377    private double spareNormal = double.NaN;
378
379    /// <summary>
380    /// Uses the Marsaglia polar method to generate a random variable
381    /// from two uniform random distributed values.
382    /// </summary>
383    /// <remarks>
384    /// Unlike <see cref="RandNormal(double, double)"/> this method does not
385    /// make use of a spare random variable. It discards the spare and thus
386    /// requires twice the number of calls to the underlying IRandom instance.
387    /// </remarks>
388    /// <param name="random">The random number generator to use.</param>
389    /// <param name="mu">The mean of the normal distribution.</param>
390    /// <param name="sigma">The standard deviation of the normal distribution.</param>
391    /// <returns>A number that is normal distributed.</returns>
392    public virtual double RandNormal(IRandom random, double mu, double sigma) {
393      return MarsagliaPolar(random, mu, sigma, out _); // do not reuse the spare normal in this case, because it could be from a different RNG
394    }
395    /// <summary>
396    /// Uses the Marsaglia polar method to generate a random variable
397    /// from two uniform random distributed values.
398    /// </summary>
399    /// <remarks>
400    /// A spare random variable is generated from the second uniformly
401    /// distributed value. Thus, the two calls to the uniform random number
402    /// generator will be made only every second call.
403    /// </remarks>
404    /// <param name="mu">The mean of the normal distribution.</param>
405    /// <param name="sigma">The standard deviation of the normal distribution.</param>
406    /// <returns>A number that is normal distributed.</returns>
407    public virtual double RandNormal(double mu, double sigma) {
408      if (useSpareNormal) {
409        useSpareNormal = false;
410        return spareNormal * sigma + mu;
411      } else {
412        useSpareNormal = true;
413        return MarsagliaPolar(Random, mu, sigma, out spareNormal);
414      }
415    }
416    private double MarsagliaPolar(IRandom random, double mu, double sigma, out double spare) {
417      double u, v, s;
418      do {
419        u = random.NextDouble() * 2 - 1;
420        v = random.NextDouble() * 2 - 1;
421        s = u * u + v * v;
422      } while (s > 1 || s == 0);
423      var mul = Math.Sqrt(-2.0 * Math.Log(s) / s);
424      spare = v * mul;
425      return mu + sigma * u * mul;
426    }
427
428    /// <summary>
429    /// Uses the Marsaglia polar method to generate a random variable
430    /// from two uniform random distributed values.
431    /// </summary>
432    /// <remarks>
433    /// A spare random variable is generated from the second uniformly
434    /// distributed value. Thus, the two calls to the uniform random number
435    /// generator will be made only every second call.
436    /// </remarks>
437    /// <param name="random">The random number generator to use.</param>
438    /// <param name="mu">The mean of the normal distribution.</param>
439    /// <param name="sigma">The standard deviation of the normal distribution.</param>
440    /// <returns>A number that is normal distributed.</returns>
441    public TimeSpan RandNormal(IRandom random, TimeSpan mu, TimeSpan sigma) {
442      return TimeSpan.FromSeconds(RandNormal(random, mu.TotalSeconds, sigma.TotalSeconds));
443    }
444    /// <summary>
445    /// Uses the Marsaglia polar method to generate a random variable
446    /// from two uniform random distributed values.
447    /// </summary>
448    /// <remarks>
449    /// A spare random variable is generated from the second uniformly
450    /// distributed value. Thus, the two calls to the uniform random number
451    /// generator will be made only every second call.
452    /// </remarks>
453    /// <param name="mu">The mean of the normal distribution.</param>
454    /// <param name="sigma">The standard deviation of the normal distribution.</param>
455    /// <returns>A number that is normal distributed.</returns>
456    public TimeSpan RandNormal(TimeSpan mu, TimeSpan sigma) {
457      return RandNormal(Random, mu, sigma);
458    }
459
460    public double RandNormalPositive(IRandom random, double mu, double sigma) {
461      double val;
462      do {
463        val = RandNormal(random, mu, sigma);
464      } while (val <= 0);
465      return val;
466    }
467    public double RandNormalPositive(double mu, double sigma) {
468      return RandNormalPositive(Random, mu, sigma);
469    }
470
471    public TimeSpan RandNormalPositive(IRandom random, TimeSpan mu, TimeSpan sigma) {
472      return TimeSpan.FromSeconds(RandNormalPositive(random, mu.TotalSeconds, sigma.TotalSeconds));
473    }
474    public TimeSpan RandNormalPositive(TimeSpan mu, TimeSpan sigma) {
475      return RandNormalPositive(Random, mu, sigma);
476    }
477
478    public double RandNormalNegative(IRandom random, double mu, double sigma) {
479      double val;
480      do {
481        val = RandNormal(random, mu, sigma);
482      } while (val >= 0);
483      return val;
484    }
485    public double RandNormalNegative(double mu, double sigma) {
486      return RandNormalNegative(Random, mu, sigma);
487    }
488
489    public TimeSpan RandNormalNegative(IRandom random, TimeSpan mu, TimeSpan sigma) {
490      return TimeSpan.FromSeconds(RandNormalNegative(random, mu.TotalSeconds, sigma.TotalSeconds));
491    }
492    public TimeSpan RandNormalNegative(TimeSpan mu, TimeSpan sigma) {
493      return RandNormalNegative(Random, mu, sigma);
494    }
495
496    /// <summary>
497    /// Returns values from a log-normal distribution with the mean
498    /// exp(mu + sigma^2 / 2)
499    /// and the standard deviation
500    /// sqrt([exp(sigma^2)-1] * exp(2 * mu + sigma^2))
501    /// </summary>
502    /// <param name="random">The random number generator to use.</param>
503    /// <param name="mu">The mu parameter of the log-normal distribution (not the mean).</param>
504    /// <param name="sigma">The sigma parameter of the log-normal distribution (not the standard deviation).</param>
505    /// <returns>A log-normal distributed random value.</returns>
506    public double RandLogNormal(IRandom random, double mu, double sigma) {
507      return Math.Exp(RandNormal(random, mu, sigma));
508    }
509    /// <summary>
510    /// Returns values from a log-normal distribution with the mean
511    /// exp(mu + sigma^2 / 2)
512    /// and the standard deviation
513    /// sqrt([exp(sigma^2)-1] * exp(2 * mu + sigma^2))
514    /// </summary>
515    /// <param name="mu">The mu parameter of the log-normal distribution (not the mean).</param>
516    /// <param name="sigma">The sigma parameter of the log-normal distribution (not the standard deviation).</param>
517    /// <returns>A log-normal distributed random value.</returns>
518    public double RandLogNormal(double mu, double sigma) {
519      return RandLogNormal(Random, mu, sigma);
520    }
521
522    /// <summary>
523    /// Returns values from a log-normal distribution with
524    /// the mean <paramref name="mean"/> and standard deviation <paramref name="stdev"/>.
525    /// </summary>
526    /// <param name="random">The random number generator to use.</param>
527    /// <param name="mean">The distribution mean.</param>
528    /// <param name="stdev">The distribution standard deviation.</param>
529    /// <returns>A log-normal distributed random value.</returns>
530    public double RandLogNormal2(IRandom random, double mean, double stdev) {
531      if (stdev == 0) return mean;
532      var sigma = Math.Sqrt(Math.Log(stdev * stdev / (mean * mean) + 1));
533      var mu = Math.Log(mean) - 0.5 * sigma * sigma;
534      return Math.Exp(RandNormal(random, mu, sigma));
535    }
536    /// <summary>
537    /// Returns values from a log-normal distribution with
538    /// the mean <paramref name="mean"/> and standard deviation <paramref name="stdev"/>.
539    /// </summary>
540    /// <param name="mean">The distribution mean.</param>
541    /// <param name="stdev">The distribution standard deviation.</param>
542    /// <returns>A log-normal distributed random value.</returns>
543    public double RandLogNormal2(double mean, double stdev) {
544      return RandLogNormal2(Random, mean, stdev);
545    }
546
547    /// <summary>
548    /// Returns a timespan value from a log-normal distribution with the mean
549    /// exp(mu + sigma^2 / 2)
550    /// and the standard deviation
551    /// sqrt([exp(sigma^2)-1] * exp(2 * mu + sigma^2))
552    /// </summary>
553    /// <param name="random">The random number generator to use.</param>
554    /// <param name="mu">The mu parameter of the log-normal distribution (not the mean).</param>
555    /// <param name="sigma">The sigma parameter of the log-normal distribution (not the standard deviation).</param>
556    /// <returns>A log-normal distributed random timespan.</returns>
557    public TimeSpan RandLogNormal(IRandom random, TimeSpan mu, TimeSpan sigma) {
558      return TimeSpan.FromSeconds(RandLogNormal(random, mu.TotalSeconds, sigma.TotalSeconds));
559    }
560    /// <summary>
561    /// Returns a timespan value from a log-normal distribution with the mean
562    /// exp(mu + sigma^2 / 2)
563    /// and the standard deviation
564    /// sqrt([exp(sigma^2)-1] * exp(2 * mu + sigma^2))
565    /// </summary>
566    /// <param name="mu">The mu parameter of the log-normal distribution (not the mean).</param>
567    /// <param name="sigma">The sigma parameter of the log-normal distribution (not the standard deviation).</param>
568    /// <returns>A log-normal distributed random timespan.</returns>
569    public TimeSpan RandLogNormal(TimeSpan mu, TimeSpan sigma) {
570      return RandLogNormal(Random, mu, sigma);
571    }
572
573    /// <summary>
574    /// Returns a timespan value from a log-normal distribution with
575    /// the mean <paramref name="mean"/> and standard deviation <paramref name="stdev"/>.
576    /// </summary>
577    /// <param name="random">The random number generator to use.</param>
578    /// <param name="mean">The distribution mean.</param>
579    /// <param name="stdev">The distribution standard deviation.</param>
580    /// <returns>A log-normal distributed random timespan.</returns>
581    public TimeSpan RandLogNormal2(IRandom random, TimeSpan mean, TimeSpan stdev) {
582      return TimeSpan.FromSeconds(RandLogNormal2(random, mean.TotalSeconds, stdev.TotalSeconds));
583    }
584    /// <summary>
585    /// Returns a timespan value from a log-normal distribution with
586    /// the mean <paramref name="mean"/> and standard deviation <paramref name="stdev"/>.
587    /// </summary>
588    /// <param name="mean">The distribution mean.</param>
589    /// <param name="stdev">The distribution standard deviation.</param>
590    /// <returns>A log-normal distributed random timespan.</returns>
591    public TimeSpan RandLogNormal2(TimeSpan mean, TimeSpan stdev) {
592      return RandLogNormal2(Random, mean, stdev);
593    }
594
595    public double RandCauchy(IRandom random, double x0, double gamma) {
596      return x0 + gamma * Math.Tan(Math.PI * (random.NextDouble() - 0.5));
597    }
598    public double RandCauchy(double x0, double gamma) {
599      return RandCauchy(Random, x0, gamma);
600    }
601
602    public TimeSpan RandCauchy(IRandom random, TimeSpan x0, TimeSpan gamma) {
603      return TimeSpan.FromSeconds(RandCauchy(random, x0.TotalSeconds, gamma.TotalSeconds));
604    }
605    public TimeSpan RandCauchy(TimeSpan x0, TimeSpan gamma) {
606      return RandCauchy(Random, x0, gamma);
607    }
608
609    public double RandWeibull(IRandom random, double alpha, double beta) {
610      return alpha * Math.Pow(-Math.Log(1 - random.NextDouble()), 1 / beta);
611    }
612    public double RandWeibull(double alpha, double beta) {
613      return RandWeibull(Random, alpha, beta);
614    }
615
616    public TimeSpan RandWeibull(IRandom random, TimeSpan alpha, TimeSpan beta) {
617      return TimeSpan.FromSeconds(RandWeibull(random, alpha.TotalSeconds, beta.TotalSeconds));
618    }
619    public TimeSpan RandWeibull(TimeSpan alpha, TimeSpan beta) {
620      return RandWeibull(Random, alpha, beta);
621    }
622
623
624    /// <summary>
625    /// Generates a random sample from a given source
626    /// </summary>
627    /// <typeparam name="T">The type of the element in parameter source</typeparam>
628    /// <exception cref="ArgumentException">
629    /// Thrown when <paramref name="source"/> and <paramref name="weights"/> have different size.
630    /// or when <paramref name="weights"/> contains an invalid or negative value.
631    /// or when <paramref name="weights"/> sum equals zero or an invalid value.
632    /// </exception>
633    /// <param name="random">The random number generator to use.</param>
634    /// <param name="source">a random sample is generated from its elements.</param>
635    /// <param name="weights">The weight associated with each entry in source.</param>
636    /// <returns>The generated random samples</returns>
637    public T RandChoice<T>(IRandom random, IList<T> source, IList<double> weights) {
638      if (source.Count != weights.Count) {
639        throw new ArgumentException("source and weights must have same size");
640      }
641
642      double totalW = 0;
643      foreach (var w in weights) {
644        if (w < 0) {
645          throw new ArgumentException("weight values must be non-negative", nameof(weights));
646        }
647        totalW += w;
648      }
649
650      if (double.IsNaN(totalW) || double.IsInfinity(totalW))
651        throw new ArgumentException("Not a valid weight", nameof(weights));
652      if (totalW == 0)
653        throw new ArgumentException("total weight must be greater than 0", nameof(weights));
654
655      var rnd = random.NextDouble();
656      double aggWeight = 0;
657      int idx = 0;
658      foreach (var w in weights) {
659        if (w > 0) {
660          aggWeight += (w / totalW);
661          if (rnd <= aggWeight) {
662            break;
663          }
664        }
665        idx++;
666      }
667      return source[idx];
668    }
669    /// <summary>
670    /// Generates a random sample from a given source
671    /// </summary>
672    /// <typeparam name="T">The type of the element in parameter source</typeparam>
673    /// <exception cref="ArgumentException">
674    /// Thrown when <paramref name="source"/> and <paramref name="weights"/> have different size.
675    /// or when <paramref name="weights"/> contains an invalid or negative value.
676    /// or when <paramref name="weights"/> sum equals zero
677    /// </exception>
678    /// <param name="source">a random sample is generated from its elements.</param>
679    /// <param name="weights">The weight associated with each entry in source.</param>
680    /// <returns>The generated random samples</returns>
681    public T RandChoice<T>(IList<T> source, IList<double> weights) {
682      return RandChoice(Random, source, weights);
683    }
684
685    #endregion
686
687    #region Random timeouts
688    public Timeout TimeoutUniformD(IRandom random, double a, double b) {
689      return new Timeout(this, ToTimeSpan(RandUniform(random, a, b)));
690    }
691    public Timeout TimeoutUniformD(double a, double b) {
692      return TimeoutUniformD(Random, a, b);
693    }
694
695    public Timeout TimeoutUniform(IRandom random, TimeSpan a, TimeSpan b) {
696      return new Timeout(this, RandUniform(random, a, b));
697    }
698    public Timeout TimeoutUniform(TimeSpan a, TimeSpan b) {
699      return TimeoutUniform(Random, a, b);
700    }
701
702    public Timeout TimeoutTriangularD(IRandom random, double low, double high) {
703      return new Timeout(this, ToTimeSpan(RandTriangular(random, low, high)));
704    }
705    public Timeout TimeoutTriangularD(double low, double high) {
706      return TimeoutTriangularD(Random, low, high);
707    }
708
709    public Timeout TimeoutTriangular(IRandom random, TimeSpan low, TimeSpan high) {
710      return new Timeout(this, RandTriangular(random, low, high));
711    }
712    public Timeout TimeoutTriangular(TimeSpan low, TimeSpan high) {
713      return TimeoutTriangular(Random, low, high);
714    }
715
716    public Timeout TimeoutTriangularD(IRandom random, double low, double high, double mode) {
717      return new Timeout(this, ToTimeSpan(RandTriangular(random, low, high, mode)));
718    }
719    public Timeout TimeoutTriangularD(double low, double high, double mode) {
720      return TimeoutTriangularD(Random, low, high, mode);
721    }
722
723    public Timeout TimeoutTriangular(IRandom random, TimeSpan low, TimeSpan high, TimeSpan mode) {
724      return new Timeout(this, RandTriangular(random, low, high, mode));
725    }
726    public Timeout TimeoutTriangular(TimeSpan low, TimeSpan high, TimeSpan mode) {
727      return TimeoutTriangular(Random, low, high, mode);
728    }
729
730    public Timeout TimeoutExponentialD(IRandom random, double mean) {
731      return new Timeout(this, ToTimeSpan(RandExponential(random, mean)));
732    }
733    public Timeout TimeoutExponentialD(double mean) {
734      return TimeoutExponentialD(Random, mean);
735    }
736
737    public Timeout TimeoutExponential(IRandom random, TimeSpan mean) {
738      return new Timeout(this, RandExponential(random, mean));
739    }
740    public Timeout TimeoutExponential(TimeSpan mean) {
741      return TimeoutExponential(Random, mean);
742    }
743
744    public Timeout TimeoutNormalPositiveD(IRandom random, double mu, double sigma) {
745      return new Timeout(this, ToTimeSpan(RandNormalPositive(random, mu, sigma)));
746    }
747    public Timeout TimeoutNormalPositiveD(double mu, double sigma) {
748      return TimeoutNormalPositiveD(Random, mu, sigma);
749    }
750
751    public Timeout TimeoutNormalPositive(IRandom random, TimeSpan mu, TimeSpan sigma) {
752      return new Timeout(this, RandNormalPositive(random, mu, sigma));
753    }
754    public Timeout TimeoutNormalPositive(TimeSpan mu, TimeSpan sigma) {
755      return TimeoutNormalPositive(Random, mu, sigma);
756    }
757
758    public Timeout TimeoutLogNormalD(IRandom random, double mu, double sigma) {
759      return new Timeout(this, ToTimeSpan(RandLogNormal(random, mu, sigma)));
760    }
761    public Timeout TimeoutLogNormalD(double mu, double sigma) {
762      return TimeoutLogNormalD(Random, mu, sigma);
763    }
764
765    public Timeout TimeoutLogNormal2D(IRandom random, double mean, double stdev) {
766      return new Timeout(this, ToTimeSpan(RandLogNormal2(random, mean, stdev)));
767    }
768    public Timeout TimeoutLogNormal2D(double mean, double stdev) {
769      return TimeoutLogNormal2D(Random, mean, stdev);
770    }
771
772    public Timeout TimeoutLogNormal(IRandom random, TimeSpan mu, TimeSpan sigma) {
773      return new Timeout(this, RandLogNormal(random, mu, sigma));
774    }
775    public Timeout TimeoutLogNormal(TimeSpan mu, TimeSpan sigma) {
776      return TimeoutLogNormal(Random, mu, sigma);
777    }
778
779    public Timeout TimeoutLogNormal2(IRandom random, TimeSpan mean, TimeSpan stdev) {
780      return new Timeout(this, RandLogNormal2(random, mean, stdev));
781    }
782    public Timeout TimeoutLogNormal2(TimeSpan mean, TimeSpan stdev) {
783      return TimeoutLogNormal2(Random, mean, stdev);
784    }
785    #endregion
786  }
787
788  /// <summary>
789  /// Provides a simulation environment that is thread-safe against manipulations of the event queue.
790  /// Its performance is somewhat lower than the non-thread-safe environment (cf. <see cref="Simulation"/>)
791  /// due to the locking involved.
792  /// </summary>
793  /// <remarks>
794  /// Please carefully consider if you must really schedule the stop event in a separate thread. You can also
795  /// call <see cref="Simulation.StopAsync"/> to request termination after the current event has been processed.
796  ///
797  /// The simulation will still run in only one thread and execute all events sequentially.
798  /// </remarks>
799  public class ThreadSafeSimulation : Simulation {
800    protected object _locker;
801
802    public ThreadSafeSimulation() : this(new DateTime(1970, 1, 1)) { }
803    public ThreadSafeSimulation(TimeSpan? defaultStep) : this(new DateTime(1970, 1, 1), defaultStep) { }
804    public ThreadSafeSimulation(DateTime initialDateTime, TimeSpan? defaultStep = null) : this(new PcgRandom(), initialDateTime, defaultStep) { }
805    public ThreadSafeSimulation(int randomSeed, TimeSpan? defaultStep = null) : this(new DateTime(1970, 1, 1), randomSeed, defaultStep) { }
806    public ThreadSafeSimulation(DateTime initialDateTime, int randomSeed, TimeSpan? defaultStep = null) : this(new PcgRandom(randomSeed), initialDateTime, defaultStep) { }
807    public ThreadSafeSimulation(IRandom random, DateTime initialDateTime, TimeSpan? defaultStep = null) : base(random, initialDateTime, defaultStep) {
808      _locker = new object();
809    }
810
811
812    /// <summary>
813    /// Schedules an event to occur at the same simulation time as the call was made.
814    /// </summary>
815    /// <remarks>
816    /// This method is thread-safe against manipulations of the event queue
817    /// </remarks>
818    /// <param name="event">The event that should be scheduled.</param>
819    /// <param name="priority">The priority to rank events at the same time (smaller value = higher priority).</param>
820    public override void Schedule(Event @event, int priority = 0) {
821      lock (_locker) {
822        DoSchedule(Now, @event, priority);
823      }
824    }
825
826    /// <summary>
827    /// Schedules an event to occur after a certain (positive) delay.
828    /// </summary>
829    /// <remarks>
830    /// This method is thread-safe against manipulations of the event queue
831    /// </remarks>
832    /// <exception cref="ArgumentException">
833    /// Thrown when <paramref name="delay"/> is negative.
834    /// </exception>
835    /// <param name="delay">The (positive) delay after which the event should be fired.</param>
836    /// <param name="event">The event that should be scheduled.</param>
837    /// <param name="priority">The priority to rank events at the same time (smaller value = higher priority).</param>
838    public override void Schedule(TimeSpan delay, Event @event, int priority = 0) {
839      if (delay < TimeSpan.Zero)
840        throw new ArgumentException("Negative delays are not allowed in Schedule(TimeSpan, Event).");
841      lock (_locker) {
842        var eventTime = Now + delay;
843        DoSchedule(eventTime, @event, priority);
844      }
845    }
846
847    /// <summary>
848    /// Run until a certain event is processed.
849    /// </summary>
850    /// <remarks>
851    /// This method is thread-safe against manipulations of the event queue
852    /// </remarks>
853    /// <param name="stopEvent">The event that stops the simulation.</param>
854    /// <returns></returns>
855    public override object Run(Event stopEvent = null) {
856      _stop = new CancellationTokenSource();
857      if (stopEvent != null) {
858        if (stopEvent.IsProcessed) {
859          return stopEvent.Value;
860        }
861        stopEvent.AddCallback(StopSimulation);
862      }
863      OnRunStarted();
864      try {
865        var stop = false;
866        lock (_locker) {
867          stop = ScheduleQ.Count == 0 || _stop.IsCancellationRequested;
868        }
869        while (!stop) {
870          Step();
871          lock (_locker) {
872            stop = ScheduleQ.Count == 0 || _stop.IsCancellationRequested;
873          }
874        }
875      } catch (StopSimulationException e) { OnRunFinished(); return e.Value; }
876      OnRunFinished();
877      if (stopEvent == null) return null;
878      if (!_stop.IsCancellationRequested && !stopEvent.IsTriggered) throw new InvalidOperationException("No scheduled events left but \"until\" event was not triggered.");
879      return stopEvent.Value;
880    }
881
882    public Task<object> RunAsync(TimeSpan duration) {
883      return Task.Run(() => Run(duration));
884    }
885
886    public Task<object> RunAsync(DateTime until) {
887      return Task.Run(() => Run(until));
888    }
889
890    /// <summary>
891    /// Run until a certain event is processed, but does not block.
892    /// </summary>
893    /// <param name="stopEvent">The event that stops the simulation.</param>
894    /// <returns></returns>
895    public Task<object> RunAsync(Event stopEvent = null) {
896      return Task.Run(() => Run(stopEvent));
897    }
898
899    /// <summary>
900    /// Performs a single step of the simulation, i.e. process a single event
901    /// </summary>
902    /// <remarks>
903    /// This method is thread-safe against manipulations of the event queue
904    /// </remarks>
905    public override void Step() {
906      Event evt;
907      lock (_locker) {
908        var next = ScheduleQ.Dequeue();
909        Now = next.PrimaryPriority;
910        evt = next.Event;
911      }
912      evt.Process();
913      ProcessedEvents++;
914    }
915
916    /// <summary>
917    /// Peeks at the time of the next event in terms of the defined step
918    /// </summary>
919    /// <remarks>
920    /// This method is thread-safe against manipulations of the event queue
921    /// </remarks>
922    public override double PeekD() {
923      lock (_locker) {
924        if (ScheduleQ.Count == 0) return double.MaxValue;
925        return (Peek() - StartDate).TotalSeconds / DefaultTimeStepSeconds;
926      }
927    }
928
929    /// <summary>
930    /// Peeks at the time of the next event
931    /// </summary>
932    /// <remarks>
933    /// This method is thread-safe against manipulations of the event queue
934    /// </remarks>
935    public override DateTime Peek() {
936      lock (_locker) {
937        return ScheduleQ.Count > 0 ? ScheduleQ.First.PrimaryPriority : DateTime.MaxValue;
938      }
939    }
940  }
941
942  /// <summary>
943  /// Provides a simulation environment where delays in simulation time may result in a similar
944  /// delay in wall-clock time. The environment is not an actual realtime simulation environment
945  /// in that there is no guarantee that 3 seconds in model time are also exactly 3 seconds in
946  /// observed wall-clock time. This simulation environment is a bit slower, as the overhead of
947  /// the simulation kernel (event creation, queuing, processing, etc.) is not accounted for.
948  ///
949  /// However, it features a switch between virtual and realtime, thus allowing it to be used
950  /// in contexts where realtime is only necessary sometimes (e.g. during interaction with
951  /// long-running co-processes). Such use cases may arise in simulation control problems.
952  /// </summary>
953  public class PseudoRealtimeSimulation : ThreadSafeSimulation {
954    public const double DefaultRealtimeScale = 1;
955
956    /// <summary>
957    /// The scale at which the simulation runs in comparison to realtime. A value smaller
958    /// than 1 results in longer-than-realtime delays, while a value larger than 1 results
959    /// in shorter-than-realtime delays. A value of exactly 1 is realtime.
960    /// </summary>
961    public double? RealtimeScale { get; protected set; } = DefaultRealtimeScale;
962    /// <summary>
963    /// Whether a non-null <see cref="RealtimeScale"/> has been set.
964    /// </summary>
965    public bool IsRunningInRealtime => RealtimeScale.HasValue;
966
967    private object _timeLocker = new object();
968    /// <summary>
969    /// The current model time. Note that, while in realtime, this may continuously change.
970    /// </summary>
971    public override DateTime Now {
972      get { lock (_timeLocker) { return base.Now + _rtDelayTime.Elapsed; } }
973      protected set => base.Now = value;
974    }
975
976    protected CancellationTokenSource _rtDelayCtrl = null;
977    protected Stopwatch _rtDelayTime = new Stopwatch();
978
979
980    public PseudoRealtimeSimulation() : this(new DateTime(1970, 1, 1)) { }
981    public PseudoRealtimeSimulation(TimeSpan? defaultStep) : this(new DateTime(1970, 1, 1), defaultStep) { }
982    public PseudoRealtimeSimulation(DateTime initialDateTime, TimeSpan? defaultStep = null) : this(new PcgRandom(), initialDateTime, defaultStep) { }
983    public PseudoRealtimeSimulation(int randomSeed, TimeSpan? defaultStep = null) : this(new DateTime(1970, 1, 1), randomSeed, defaultStep) { }
984    public PseudoRealtimeSimulation(DateTime initialDateTime, int randomSeed, TimeSpan? defaultStep = null) : this(new PcgRandom(randomSeed), initialDateTime, defaultStep) { }
985    public PseudoRealtimeSimulation(IRandom random, DateTime initialDateTime, TimeSpan? defaultStep = null) : base(random, initialDateTime, defaultStep) { }
986   
987    protected override EventQueueNode DoSchedule(DateTime date, Event @event, int priority = 0) {
988      if (ScheduleQ.Count > 0 && date < ScheduleQ.First.PrimaryPriority) _rtDelayCtrl?.Cancel();
989      return base.DoSchedule(date, @event, priority);
990    }
991
992    public override void Step() {
993      var delay = TimeSpan.Zero;
994      double? rtScale = null;
995      lock (_locker) {
996        if (IsRunningInRealtime) {
997          rtScale = RealtimeScale;
998          var next = ScheduleQ.First.PrimaryPriority;
999          delay = next - base.Now;
1000          if (rtScale.Value != 1.0) delay = TimeSpan.FromMilliseconds(delay.TotalMilliseconds / rtScale.Value);
1001          _rtDelayCtrl = CancellationTokenSource.CreateLinkedTokenSource(_stop.Token);
1002        }
1003      }
1004
1005      if (delay > TimeSpan.Zero) {
1006        _rtDelayTime.Start();
1007        Task.Delay(delay, _rtDelayCtrl.Token).ContinueWith(_ => { }).Wait();
1008        _rtDelayTime.Stop();
1009        var observed = _rtDelayTime.Elapsed;
1010
1011        lock (_locker) {
1012          if (rtScale.Value != 1.0) observed = TimeSpan.FromMilliseconds(observed.TotalMilliseconds / rtScale.Value);
1013          if (_rtDelayCtrl.IsCancellationRequested && observed < delay) {
1014            lock (_timeLocker) {
1015              Now = base.Now + observed;
1016              _rtDelayTime.Reset();
1017            }
1018            return; // next event is not processed, step is not actually completed
1019          }
1020        }
1021      }
1022
1023      Event evt;
1024      lock (_locker) {
1025        var next = ScheduleQ.Dequeue();
1026        lock (_timeLocker) {
1027          _rtDelayTime.Reset();
1028          Now = next.PrimaryPriority;
1029        }
1030        evt = next.Event;
1031      }
1032      evt.Process();
1033      ProcessedEvents++;
1034    }
1035
1036    /// <summary>
1037    /// Switches the simulation to virtual time mode, i.e., running as fast as possible.
1038    /// In this mode, events are processed without delay just like in a <see cref="ThreadSafeSimulation"/>.
1039    /// </summary>
1040    /// <remarks>
1041    /// An ongoing real-time delay is being canceled when this method is called. Usually, this
1042    /// is only the case when this method is called from a thread other than the main simulation thread.
1043    ///
1044    /// If the simulation is already in virtual time mode, this method has no effect.
1045    /// </remarks>
1046    public virtual void SetVirtualtime() {
1047      lock (_locker) {
1048        if (!IsRunningInRealtime) return;
1049        RealtimeScale = null;
1050        _rtDelayCtrl?.Cancel();
1051      }
1052    }
1053
1054    /// <summary>
1055    /// Switches the simulation to real time mode. The real time factor of
1056    /// this default mode is configurable.
1057    /// </summary>
1058    /// <remarks>
1059    /// If this method is called while running in real-time mode, but given a different
1060    /// <paramref name="realtimeScale"/>, the current delay is canceled and the remaining
1061    /// time is delayed using the new time factor.
1062    ///
1063    /// The default factor is 1, i.e., real time - a timeout of 5 seconds would cause
1064    /// a wall-clock delay of 5 seconds. With a factor of 2, the delay as measured by
1065    /// a wall clock would be 2.5 seconds, whereas a factor of 0.5, a wall-clock delay of
1066    /// 10 seconds would be observed.
1067    /// </remarks>
1068    /// <param name="realtimeScale">A value strictly greater than 0 used to scale real time events.</param>
1069    public virtual void SetRealtime(double realtimeScale = DefaultRealtimeScale) {
1070      lock (_locker) {
1071        if (realtimeScale <= 0.0) throw new ArgumentException("The simulation speed scaling factor must be strictly positive.", nameof(realtimeScale));
1072        if (IsRunningInRealtime && realtimeScale != RealtimeScale) _rtDelayCtrl?.Cancel();
1073        RealtimeScale = realtimeScale;
1074      }
1075    }
1076
1077    /// <summary>
1078    /// This is only a convenience for mixed real- and virtual time simulations.
1079    /// It creates a new pseudo realtime process which will set the simulation
1080    /// to realtime every time it continues (e.g., if it has been set to virtual time).
1081    /// The process is automatically scheduled to be started at the current simulation time.
1082    /// </summary>
1083    /// <param name="generator">The generator function that represents the process.</param>
1084    /// <param name="priority">The priority to rank events at the same time (smaller value = higher priority).</param>
1085    /// <param name="realtimeScale">A value strictly greater than 0 used to scale real time events (1 = realtime).</param>
1086    /// <returns>The scheduled process that was created.</returns>
1087    public Process PseudoRealtimeProcess(IEnumerable<Event> generator, int priority = 0, double realtimeScale = DefaultRealtimeScale) {
1088      return new PseudoRealtimeProcess(this, generator, priority, realtimeScale);
1089    }
1090  }
1091
1092  /// <summary>
1093  /// Environments hold the event queues, schedule and process events.
1094  /// </summary>
1095  [Obsolete("Use class Simulation or ThreadSafeSimulation instead. Due to name clashes with System.Environment the class SimSharp.Environment is being outphased.")]
1096  public class Environment : ThreadSafeSimulation {
1097    public Environment()
1098      : base() {
1099      Random = new SystemRandom();
1100    }
1101    public Environment(TimeSpan? defaultStep)
1102      : base(defaultStep) {
1103      Random = new SystemRandom();
1104    }
1105    public Environment(int randomSeed, TimeSpan? defaultStep = null)
1106      : base(randomSeed, defaultStep) {
1107      Random = new SystemRandom(randomSeed);
1108    }
1109    public Environment(DateTime initialDateTime, TimeSpan? defaultStep = null)
1110      : base(initialDateTime, defaultStep) {
1111      Random = new SystemRandom();
1112    }
1113    public Environment(DateTime initialDateTime, int randomSeed, TimeSpan? defaultStep = null)
1114      : base(initialDateTime, randomSeed, defaultStep) {
1115      Random = new SystemRandom(randomSeed);
1116    }
1117
1118    protected static readonly double NormalMagicConst = 4 * Math.Exp(-0.5) / Math.Sqrt(2.0);
1119    public override double RandNormal(double mu, double sigma) {
1120      return RandNormal(Random, mu, sigma);
1121    }
1122    public override double RandNormal(IRandom random, double mu, double sigma) {
1123      double z, zz, u1, u2;
1124      do {
1125        u1 = random.NextDouble();
1126        u2 = 1 - random.NextDouble();
1127        z = NormalMagicConst * (u1 - 0.5) / u2;
1128        zz = z * z / 4.0;
1129      } while (zz > -Math.Log(u2));
1130      return mu + z * sigma;
1131    }
1132  }
1133}
Note: See TracBrowser for help on using the repository browser.