15. Simulation and Frameworks

Where Are We

Up to this point, we have developed code that can read a model into memory and write it back out again, detecting a variety of errors in the model. The only point of writing the model out is to verify that the model has indeed been read correctly.

This highway network model could be used for many purposes:

• We could use the model (augmented with map coordinates) to display a map, the way Google Maps does.
• We could use the model to provide driving directions by running a shortest path algorithm over the model.
• We could use the same model to run a simulation of traffic loading to allow testing of the consequences of various proposed road improvements.

The logic circuit model we have discussed could also serve many purposes:

• We could use the model (augmented with 3-d coordinates) to display the layout of the electronic components of a computer system. This might be an intermediate point in chip design.
• We could use the model to run a simulation of a system of logic gates in order to test the logical design of a new chip.

A neuron network model could also serve many purposes:

• We could use the model (augmented with 3-d coordinates) to display a map of the physical anatomy of the neurons in an organism. A neuroanatomist could use this model to express the results of study of the actual organism.
• We could use the model to run a simulation of the activity of the nervous system of the organism. This would allow us to compare our understanding of how the nervous system works with actual measurements.

A logic network model could also serve many purposes:

• We could use the model (augmented with 2-d coordinates) to display the logical layout of the circuit for a computer-aided-design system.
• We could use the model (augmented with 2-d coordinates) to display the physical layout of the circuit for an integrated-circuit-layout system.
• We use the model to simulate the behavior of the circuit in order to test the logical design. This, of course, could be integrated with any of the above.

The model we have been referring to as an epidemic model is actually just a model of a community.

• It could be used to generate traffic for a vehicle simulation as it models people traveling homes, schools, workplaces and stores.

• The same model could be used to simulate the spread of news through a pre-modern community that did not have newspapers or the Internet.

• The same model could be used to model the spread of an epidemic through a community, the particular goal that motivated our discussion of this model.

Our goal for the projects in this course is to build simulations, and the time has come to discuss this in more detail.

Simulation

An old bumper sticker I picked up at a simulation conference said: "Simulationists do it continuously and discretely." The sticker was a joke because, while members of the general public reading the sticker might guess one subject (sex), the actual statement is entirely true when you read "it" as a reference to computer simulation. There are two broad categories of simulation:

• Continuous simulation typically involves use of differential equations to describe the behavior of a physical system, and the simulation itself involves the numerical solution of those equations.

  Continuous simulation models are common in fields as distinct as analog electonics, weather forecasting and macro economics. Here at the University of Iowa, the Hydraulics Institute is largely devoted to continuous simulation of fluid flow. Their building was built in the era when their research was largely done using actual tanks of water and even water from the Iowa River, but today, much of their work is done on computers, simulating not only water flow but also atmospheric flow.

• Discrete-event simulation typically involves systems where the behavior of the system can be described in terms of a sequence of events. Between events, nothing happens, but at each event, parts of the model interact to make changes in the simulation variables.

  Discrete event simulations are common in fields as distinct as logistics and digital logic simulation.

Almost all simulation models are based on simplifying assumptions. Most physics models assume that air is a vacuum and that the earth is flat. You can build bridges with these assumptions, although for medium and large bridges, it is worth checking how the bridge responds to the wind. (The Tacoma Narrows Bridge disaster of 1940 shows what can happen if you forget the wind when you design a large bridge -- it's in Wikipedia, watch the film clip.)

Our distinction between continuous and discrete models is also oversimplified. There are mixed models, for example, where the set of differential equations that describe the behavior of a system changes at discrete events. At each such event, you need to do continuous simulations to predict the times of the next events.

A Highway Network Model

In a highway network, the events we are concerned with are:

• Arrival at intersection: When a car arrives at an intersection, it may have to wait for other cars in the intersection to clear the intersection (we can assume that cars spend only a fixed short time in the intersection), and it may have to wait for the light to change. If it does not have to wait, it immediately schedules a road-entry event for the road its navigation algorithm selects.

• Entry to road: When a car leaves the intersection, it unblocks the intersection, allowing another car to pass through (if one was waiting); it does this by scheduling an entry-to-road event for that car. It also schedules the arrival at intersection event for this car at the other end of the road this car is entering.

• Stoplight change: Each intersection with a stoplight has a light that changes periodically by scheduling stoplight changes at that intersection. With each change, a different road (and the queue of waiting cars for that road) is unblocked, allowing cars in that queue to proceed by scheduling a road-entry event for the first waiting car.

Of course, the model can vary considerably in complexity. A simple model might have a fixed travel time on each road segment, while a more complex simulation might model congestion by having the travel time get longer if the population of a road segment exceeds some threshold, where that threshold may itself depend on the unpopulated travel time.

In a crude model, each car might make random navigation decisions at each intersection. A more complex model might have each car follow a fixed route through the road network, while a really complex model might include cars with adaptive navigation algorithms so that drivers can take alternate routes when congestion slows them on the path they originally planned.

A Digital Logic Model

In a network of digital logic gates connected by wires, we might have the following events:

• Input to a logic gate changes: The values of the other inputs to that gate are checked, using the rule for that gate to compute the new output. For and gates, for example, if both inputs are true, the output is true, otherwise it is false. But, gates do not compute their outputs instantaneously, so if the new output differs from the previous output, a change of the output of the gate is scheduled to happen at a future time based on the time delay of this gate. In real logic used in modern microprocessors, this delay is a fraction of a nanosecond.

• Output of a logic gate changes: When the output of a logic gate changes, this change is transmitted to all of the inputs connected to that output by wires, but wires do not transmit this change instantly. In a perfectly straight wire in a vacuum with no conductors nearby, the change travels down the wire at the speed of light. In real wires, the top speed is closer to half that, and if the signal is being propagated through the body of an integrated circuit chip, it can be 1/10 the speed of light. So, the input changes at the far end of a wire are scheduled to take place some delay after the output of a gate injects a change into a wire.

The key element in the above that needs extra discussion is that if the output of a gate is changed and then changed back very quickly, no output change actually occurs. That is, there is a shortest pulse that the gate can generate on its output.

A Neural Model

In a neural network model, with neurons connected by syapses, we might have the following events:

• Primary synapse fires: The voltage of the neuron at this time is computed, and then it is incremented by the strength of this synapse. If the resulting voltage is above the neuron's threshold, the neuron fires. This schedules synapse-fire events at each outgoing synapse to occur at future times determined by the delay of each synapse.

• Secondary synapse fires: The strength of the destination primary synapse is adjusted by the strength of this secondary synapse. Secondary synapses that operate this way can serve as the basis of long-term memory because they can turn primary synapses on and off.

The key element in the above that needs extra discussion is how the voltage on a neuron changes with time. Between events, the voltage on a neuron decays exponentially, slowly leaking away unless it is pumped up by a synapse firing. So, for each neuron, we record:

• The voltage, v
• The time at which that voltage was known t

Now, if we want to know the voltage at a later time t', we use this formula:

Of course, once you compute the voltage at the new time t', you record the new voltage and the new time so you can work forward from that the next time you need to do this computation. The constant k determines the decay rate of the neuron (it must be negative).

In a simple model, all the neurons might have the same threshold and the same decay rate, and all synapses might have the same strength. More complex models allow these to be varied.

In simpler models, the voltage on a neuron goes to zero when that neuron fires. In more complex models, the threshold has its own decay rate and the threshold goes up when the neuron fires.

In complex models, the strength of a synapse weakens each time it fires because the chemical neurotransmitter is used up by firing. During the time before the next firing, the neurotransmitter can build up toward its resting level. This allows the neural network to get tired if it is too active. (You can actually see this effect at work in your visual pathways. Look at a bright light and then look away, and you will see a negative afterimage that fades as the synapses that were overexcited recharge.)

An Epidemic Model

In an epidemic model, where people move between places and their contact patterns spread some disease, we have the following basic events:

• Person arrive at place: The person is added to the list of people who might be infected because they share that place, and if the person is contageous, the count of infected people in that place is incremented. Also, the person's departure from that place is scheduled at some future time depending on the place, type of person and their infection state. People who are symptomatic (bedridden) do not leave home until recovered. For example, an employee arriving at a workplace schedules a return home at the end of the workday, and a student arriving at a school schedules a return home at the end of the school day. People arriving home may schedule trips to stores in the late afternoon and evening. Students arriving home, if they do not schedule another trip to a store, schedule the next morning's trip to school, and employees arriving home, if they do not schedule a trip to a store, schedule their next trip to work the next morning.

• Person departing from place: The person is removed from the list of people who might be infected because they share that place, and if the person is contageous, the count of infected people at that place is decremented. The distribution of travel times is used to schedule the person's arrival at their destination.

  Any time a person arrives or departs from a place, the number of infected people in that place during the interval since the last arrival or departure, times the length of that interval, times the probability of infection per unit time for the place, gives the probablity that any person in that place during that time will be infected.

• Person is infected: Schedule the time the person's state will change from latent to infectious.

• Person becomes infectious: Schedule the time the person's state will change from infectious to bedridden.

• Person becomes bedridden: Schedule the time the person's state will change from bedridden to either recovered or dead.

• Person becomes recovered: If at home, schedule their next departure for work, school or shopping. If anywhere else, their next trip is already scheduled.

• Person becomes dead: Wherever that person is, remove them from that place, including from the list of people at that place and the count of infectious people there, and do not schedule any further activity for that person.

The above outline is appropriate for a "slow disease," that is, one where, when you begin to feel bad while at work or school, you stay until the end of the day before going home and staying home until better. The model would need modification for a "fast disease," whith a sudden onset that sends people home in midday.

The Event Set

The key to discrete-event simulation is a data structure called the pending-event set. This holds the set of all future events, events that have not yet been simulated, but that have been predicted to occur at future times as a result of events that have already been simulated.

The simulator operates by picking events out of the pending event set in chronological order. Simulating any particular event can involve any combination of changing variables that represent the state of the simulation model, on the one hand, and scheduling future events by adding them to the pending event set.

Some events may only change state variables without scheduling any future events. Other events may schedule future events without making any change to state variables. We can summarize the basic discrete-event simulation algorithm with the following pseudo-code:

// Initialize
PendingEventSet eventSet = new PendingEventSet();
for each event e that known in advance to happen at future time e.time {
    eventSet.add( e );
{

// Simulate
while (!eventSet.isEmpty()) {
    Event e = eventSet.remove(); // e is the new current event at e.time
    simulate e by {
        // update simulation state variables to account for e at e.time
        for each event f that is a consequence of e {
            eventSet.add( f );
        }
        // if event requires, force simulation to terminate
    }
    // simulation terminated because we ran out of events
}

Note that the above code gives two different ways to terminate the simulation, one by running out of events, and the other involving a specific event, perhaps scheduled as part of the initialization to mark the end of time.

Either approach to termination works equally well. If we have a global constant endOfTime, we can make the event set become empty at that time by checking, as each event is scheduled, to see if it happens before or after the end of time. If it happens before, schedule it. If it happens after, simply discard the event notice.

So what is the event set? It is a priority queue sorted by the time that events are scheduled to happen. The order in which events are scheduled to happen has little to do with the times at which it can be accurately predicted that they will happen, so this is not a first-in, first-out queue.

Java provides several classes that can be made to do this job, not all of which are even implementations of the Queue interface. This means that the different candidates for the pending event set in the Java library aren't interchangable. Java's class PriorityQueue is based on a heap implementation. ConcurrentSkipListSet and TreeSet are other options. Sadly, the Java library has PriorityQueue as a final class and not an interface, since there are actually a wide range of algorithms for priority queues that all accomplish the same thing with different tradeoffs between different performance parameters.

We'll luse PriorityQueue here. One of the big differences between the Java alternatives that may concern you is whether the ordering is stable. Stable priority queues guarantee that two elements inserted with equal priority will be dequeued in the order in which they were enqueued. For well formulated simulation models, stable priority queues are not necessary because the order in which simultaneous events are simulated should not matter. In the real world, if two outcomes are possible from simultaneous events, it is highly likely that either outcome is correct. Stability may be useful to endure repeatability and easy debugging, but it may also be misleading in cases where the real world behavior is not repeatable because both outcomes are possible.

If you look up class PriorityQueue in the on-line Java reference material, you'll find that the most elementary constructor for priority queues sorts the elements by their natural ordering. That is, all of the items in the queue must be comparable. There is, however, a constructor described as follows:

PriorityQueue(int initialCapacity, Comparator comparator)

There are two obvious questions here: What is the initial capacity and what is a comparator. Java will automatically grow the queue to any size (until your computer runs out of memory), but by specifying the initial capacity of the queue, you can avoid unnecessary startup costs. For example, if a simulation of a population guarantees that all people will have one pending event (the time of their next scheduled move from here to there), you can use this to give the queue a hint that it will need to be at least that big.

The comparator is a more interesting problem. This must be an object with a compare method that can be used to compare two elements of the queue. Every time the queue implementation needs to compare two items, it will use the compare method of the object you pass. As we have seen, λ notation is a way to pass an object that carries a function, and this works for passing a comparator. Therefore, we can create the priority queue we need like this, assuming that event times of type double:

PriorityQueue eventSet = new PriorityQueue<>(
    initialCapacity,
    (Event e1, Event e2)-> Double.compare( e1.time, e2.time )
);

/** Framework for discrete event simulation.
 */
public class Simulator {

    public static abstract class Event {
        public double time;       // the time of this event
        abstract void trigger(); // the action to take
    }

    private static PriorityQueue<Event> eventSet
        = new PriorityQueue<Event> (
            (Event e1, Event e2) -> Double.compare( e1.time, e2.time )
        );

    /** Call schedule(e) to make e happen at its (future) time
     */
    static void schedule( Event e ) {
        eventSet.add( e );
    }

    /** Call run() after scheduling some initial events
     *  to run the simulation.
     */
    static void run() {
        while (!eventSet.isEmpty()) {
            Event e = eventSet.remove();
	    e.trigger();
        }
    }
}

/** Framework for discrete event simulation.
 */
public class Simulator {

    public interface Action {
        // actions contain the specific code of each event
        void trigger( double time );
    }

    private static class Event {
        public double time; // the time of this event
        public Action act; // what to do at that time
    }

    private static PriorityQueue<Event> eventSet
        = new PriorityQueue<Event> (
            (Event e1, Event e2)-> Double.compare( e1.time, e2.time )
        );

    /** Call schedule to make act happen at time.
     *  Users typically pass the action as a lambda expression:
     *  <PRE>
     *  Simulator.schedule( t, ( double time )-> method( ... time ... ) )
     *  </PRE>
     */
    static void schedule( double time, Action act ) {
        Event e = new Event();
        e.time = time;
        e.act = act;
        eventSet.add( e );
    }

    /** run the simulation.
     *  Call run() after scheduling some initial events to run the simulation.
     */
    static void run() {
        while (!eventSet.isEmpty()) {
            Event e = eventSet.remove();
            e.act.trigger( e.time );
        }
    }
}