Simple Simulation in Racket

I updated the simple simulation example from the simulation collection to Racket. The simple simulation example shows how a basic discrete-event simulation engine can be written in pure Racket (specifically, the racket/base language). Fully commented it is only about two and a half pages of code and does not contain any references to external packages. If you're curious about the usage of continuations, this might be a good example to look at. [If you're an expert in continuations, any critique would be welcome - particularly on the event loop in the start-simulation procedure.]

The code uses parameters for its global variables, which may be unfamiliar to some people. For example, the simulation clock is maintained by the current-time parameter. This parameter is created by (make-parameter current-time 0.0). A call to (current-time) will return its current value. The call (current-time (event-time (current-event))) in start-simulation sets its value to the time of the next event. Finally, the parameterize call in run-simulation creates new instances of the parameters and (I think) is cleaner than just resetting their values.

And, here is the actual code.

#lang racket/base
;;; Simplified Simulation System

;;; Global simulation control variables

;;; future-event-list : (parameter/c (list? event?))
;;; current-time : (parameter/c real?)
;;; current-event : (parameter/c (or/c false/c event?))
;;; event-loop-exit : (parameter/c (or/c false/c continuation?))
;;; event-loop-next : (parameter/c (or/c false/c continuation?))
(define future-event-list (make-parameter '()))
(define current-time (make-parameter 0.0))
(define current-event (make-parameter #f))
(define event-loop-exit (make-parameter #f))
(define event-loop-next (make-parameter #f))

;;; Event definition and scheduling

;;; (struct event (time function arguments))
;;;   time : (>=/c 0.0)
;;;   function : (or/c false/c procedure?)
;;;   arguments : list?
;;; Each event has a time the event is to be executed, the function to
;;; be executed, and the (evaluated) arguments to the function.
(struct event (time function arguments))

;;; (schedule event) -> any
;;;   event : event?
;;; Add an event to the event list.
(define (schedule event)
 (future-event-list (event-schedule event (future-event-list))))

;;; (event-schedule event event-list) -> (list-of event?)
;;;   event : event?
;;;   event-list : (list-of event?)
;;; Return a new list of events corresponding to the given event added
;;; to the given list of events.
(define (event-schedule event event-list)
 (cond ((null? event-list)
        (list event))
       ((< (event-time event)
           (event-time (car event-list)))
        (cons event event-list))
       (else
        (cons (car event-list)
              (event-schedule event (cdr event-list))))))
;;; Simulation control routines

;;; (wait/work delay) -> any
;;;   delay : (>=/c 0.0)
;;; Simulate the delay while work is being done.  Add an event to
;;; execute the current continuation to the event list.
(define (wait/work delay)
 (let/cc continue
   ;; Add an event to execute the current continuation
   (schedule (event (+ (current-time) delay) continue '()))
   ;; Return to the main loop
   ((event-loop-next))))

;;; (start-simulation) -> any
;;; This is the main simulation loop.  As long as there are events to
;;; be executed (or until the simulation is explicitly stopped), remove
;;; the next event from the event list, advance the clock to the time
;;; of the event, and apply the event's functions to its arguments.
(define (start-simulation)
 (let/ec exit
   ;; Save the event loop exit continuation
   (event-loop-exit exit)
   ;; Event loop
   (let loop ()
     ;; Exit if no more events
     (when (null? (future-event-list))
       ((event-loop-exit)))
     (let/cc next
       ;; Save the event loop next continuation
       (event-loop-next next)
       ;; Execute the next event
       (current-event (car (future-event-list)))
       (future-event-list (cdr (future-event-list)))
       (current-time (event-time (current-event)))
       (apply (event-function (current-event))
              (event-arguments (current-event))))
     (loop))))

;;; (stop-simulation) -> any
;;; Stop the execution of the current simulation (by jumping to its
;;; exit continuation).
(define (stop-simulation)
 ((event-loop-exit)))

;;; Random Distributions (to remove external dependencies)

;;; (random-flat a b) -> inexact-real?
;;;   a : real?
;;;   b : real?
;;; Returns a random real number from a uniform distribution between a
;;; and b.
(define (random-flat a b)
 (+ a (* (random) (- b a))))

;;; (random-exponential mu) -> inexact-real?
;;;   mu : real?
;;; Returns a random real number from an exponential distribution with
;;; mean mu.
(define (random-exponential mu)
 (* (- mu) (log (random))))

;;; Example Simulation Model

;;; (generator n) -> any
;;;   n : exact-positive-integer?
;;; Process to generate n customers arriving into the system.
(define (generator n)
 (do ((i 0 (+ i 1)))
     ((= i n) (void))
   (wait/work (random-exponential 4.0))
   (schedule (event (current-time) customer (list i)))))

;;; (customer i) -> any
;;;   i : exact-nonnegative-integer?
;;; The ith customer into the system.  The customer is in the system
;;; 2 to 10 minutes and then leaves.
(define (customer i)
 (printf "~a: customer ~a enters~n" (current-time) i)
 (wait/work (random-flat 2.0 10.0))
 (printf "~a: customer ~a leaves~n" (current-time) i))

;;; (run-simulation n) -> any
;;;   n : exact-positive-integer?
;;; Run the simulation for n customers (or until explicitly stopped at
;;; some specified time).
(define (run-simulation n)
 ;; Create new global values
 (parameterize ((future-event-list '())
                (current-time 0.0)
                (current-event #f)
                (event-loop-exit #f)
                (event-loop-next #f))
   ;; Schedule the customer generator
   (schedule (event 0.0 generator (list n)))
   ;; Stop the simulation at the specified time (optional)
   ;(schedule (event 50.0 stop-simulation '()))
   ;; Start the simulation main loop
   (start-simulation)))

;;; Run the simulation for 10 customers.
(run-simulation 10)

An example of the output looks like:

2.3985738870908615: customer 0 enters
5.395876334125138: customer 1 enters
7.716207787604494: customer 2 enters
8.062394508850671: customer 1 leaves
9.204092461998275: customer 0 leaves
14.431739507072376: customer 3 enters
15.15611565215575: customer 2 leaves
19.107487138771972: customer 4 enters
19.67563344212178: customer 3 leaves
26.51109574171283: customer 5 enters
27.060305975366514: customer 4 leaves
31.04777263526701: customer 6 enters
32.83444054122948: customer 5 leaves
36.31758748274069: customer 7 enters
39.27687483880874: customer 6 leaves
40.643350655011744: customer 8 enters
41.76876758081738: customer 7 leaves
42.82235216823591: customer 8 leaves
47.63479457664656: customer 9 enters
57.31661407095231: customer 9 leaves

Your output will vary slightly because of the random numbers. That is, this is a stochastic model.

Status Update

I haven't posted for a while, so I just wanted to provide an update on the status and (hopeful) direction of my various PLaneT packages. This post will cover the science, simulation, and inference collections, which work together to provide a knowledge-based simulation capability.

We at SET Corporation are using the science, simulation, and inference collections as part of our Multi-Agent, Dynamic NEtwork Simulation System (MADNESS). The inference collection is also used in our COntent GENeration from Templates (COGENT) system. These are used on a daily basis in our agent-based simulation work.

Science Collection

The current release of the science collection is version 3.10. I only expect to do bug fixes to the version 3 code. The next release should be version 4.0. For version 4.0 I plan to convert the code to Typed Racket, which has gotten to the point that typed code generally runs as fast or faster than my hand-optimized numeric code.

The only major new functionality planned for version 4.0 are fast Fourier transforms. This code is already written and ready for release, but the documentation has not been updated. Vincent St-Amour at Northeastern University converted the FFT code to Typed Racket to test the optimized code generation for numeric processing. The results were extremely impressive and I need to make a separate post on this.

Simulation Collection

The current version of the simulation collection is version 3.4. I only expect to do bug fixes to the version 3 code. At some point I will release a version 4.0, which may (or may not) be in Typed Racket. The decision to use Typed Racket will be based on my experience in converting the science collection. I will also rewrite the simulation language using syntax-parse to provide better syntax error handling.

Inference Collection

The current version of the inference collection is version 2.7. There was a flurry of releases earlier this year as we really began to use the inference collection in MADNESS and COGENT. Unfortunately, the documentation has not yet caught up with the code changes.

I made earlier posts on the upcoming version 3 rewrite. As with the simulation collection, this may (or may not) be in Typed Scheme. I will definitely use syntax-parse to implement a rule compiler for the rule language.