Simulation with Cofunctions
===========================

In this example, I'm going to develop a process-oriented discrete-event
simulation kernel.

This will be very similar to my earlier scheduler example. The main difference
is that we will need to provide a way for a process to suspend itself until
a given amount of simulation time has elapsed.

The Kernel
----------

First we will want a variable to keep track of the current simulation time.
While we're at it, we'll define a public function for accessing it as well.

  current_time = 0.0
  
  def now():
    return current_time

Next we will need a way to represent the queue of scheduled processes,
together with the times at which they are due to run. For this, we will use
a heap queue of tuples (t, p), where t is a time and p is a process.

  from heapq import heappush, heappop
  event_queue = []

Scheduling a process at a given time consist of inserting a tuple into the
event queue.

  def schedule(process, time):
    heappush(event_queue, (time, process))

Now we're ready to write the main loop. We repeatedly retrieve the process
due to run next, advance the simulation time to its scheduled time, and give
it a chance to run. We will also store the current process in a global,
because it will be useful to know later.

  current_process = None

  def run():
    global current_time, current_process
    while event_queue:
      current_time, current_process = heappop(event_queue)
      try:
        next(current_process)
      except StopIteration:
        pass

When there are no processes left in the event queue, there is nothing more to
do, so run() simply returns.

Here we're making an assumption that, each time it is run, a process will do
something to either re-schedule itself at a later time, or suspend itself on a
queue somewhere so that it can be re-scheduled later. If it simply yields
without doing either of those things, it will vanish from the system. This
could be fixed in a variety of ways, but in the interests of simplicity we will
leave things as they are.

At this point, we have everything we absolutely need for discrete-event
simulation, but things aren't quite as convenient to use as they could be,
so let's define a few higher-level interfaces using the primitives we have
so far.

One very common thing to do is suspend the current process for a specified
length of time, an operation traditionally known as "holding". This will be
a cofunction, because it requires the current process to suspend itself.

  codef hold(delay):
    schedule(current_process, now() + delay)
    yield

Another thing that it will be common to do is schedule a process to run
at the current time. This is an ordinary function, because it applies to
another process and doesn't need to suspend the calling process.

  def wakeup(process):
    schedule(process, now())

We would also like a convenient way of getting a process started.
The schedule() function as we have defined it expects an already-started
generator, whereas we would like to be able to pass a cofunction.

  def start(cofunction, *args, **kwds):
    wakeup(costart(cofunction, *args, **kwds))

Restaurant 1
------------

Now we can do some simulation. Let's model a restaurant. We'll start by
writing a process to generate customers at random intervals.

  import random
  random.seed(12345)

  codef generate_customers(howmany):
    for i in range(howmany):
      print("Generating a customer at", now())
      start(customer, i)
      hold(random.expovariate(1/5))

To begin with, our customers will just start eating straight away. We'll
make it more realistic later.

  codef customer(i):
    print("Customer", i, "starting to eat spam")
    hold(random.normalvariate(10, 5))
    print("Customer", i, "finished eating at", now())

To get things going, we start the customer generator, then run the
scheduler.

  start(generate_customers, 5)
  run()

Typical output:

  Generating a customer at 0.0
  Customer 0 starting to eat spam
  Generating a customer at 4.37790524776
  Customer 1 starting to eat spam
  Customer 0 finished eating at 7.26531575703
  Generating a customer at 12.5861265916
  Customer 2 starting to eat spam
  Customer 1 finished eating at 15.0533030872
  Customer 2 finished eating at 21.2725385233
  Generating a customer at 23.0127452598
  Customer 3 starting to eat spam
  Generating a customer at 31.7463824577
  Customer 4 starting to eat spam
  Customer 3 finished eating at 33.7204064476
  Customer 4 finished eating at 40.2268060542

Resources
---------

Our restaurant model is a bit optomistic at the moment, because customers
get served immediately and there is an infinite number of tables. To fix
this, we need to be able to model finite resources that processes can wait
for.

To do this, it will be convenient to add another primitive to our scheduler.
It takes a list representing a queue and adds the current process to the
tail of the queue.

  def enqueue(queue):
    queue.append(current_process)

Now we can create a class for modelling a resource. It will have a number
representing the amount available, and a queue of processes waiting for the
resource.

  class Resource:
  
    def __init__(self, capacity):
      self.available = capacity
      self.queue = []
    
The next method is a little subtle. When a process tries to acquire some of the
resource, it adds itself to the queue and then waits until it reaches the head.
Then, it needs to stay at the head of the queue and wait until the requested
amount of the resource becomes available. Once that happens, it can remove
itself from the queue and take the resource.

    codef acquire(self, amount):
      enqueue(self.queue)
      if len(self.queue) > 1:
        yield
      while amount > self.available:
        yield
      self.queue.pop(0)
      self.available -= amount
    
The following method cooperates with the previous one. Whenever a process
releases some of the resource, and there is another process at the head of
the queue, it wakes the other process up so that it can check whether there
is enough of the resource available for it.

    def release(self, amount):
      self.available += amount
      if self.queue:
        wakeup(self.queue[0])

Restaurant 2
------------

Now we can model our restaurant more accurately. It will have a Resource
representing tables, and another one representing waiters. A customer will
first acquire a table, and then a waiter. When he gets served, he will
release the waiter and eat. Finally, he will release the table and leave.

  tables = Resource(8)
  waiters = Resource(3)
  
  codef customer(i):
    print("Customer", i, "arriving at", now())
    tables.acquire(1)
    print("Customer", i, "sits down at a table at", now())
    waiters.acquire(1)
    print("Customer", i, "orders spam at", now())
    hold(random.normalvariate(20, 2))
    waiters.release(1)
    print("Customer", i, "gets served spam at", now())
    hold(random.normalvariate(10, 5))
    print("Customer", i, "finished eating at", now())
    tables.release(1)