/*----------------------------------------------------------------------------* TIME MANAGEMENT This module manages time through an event scheduler/dispatcher, efficient even for queues containing hundreds of millions of events. This algorithm may be used for scheduling actual events, as in a real-time operating system, or for simulated events. In the former case, it is synchronized with a real-time clock. In the latter, the next event is dispatched immediately and the simulated clock is set to the time of that event. The simulation system thus jumps forward from one event to the next, skipping intervening times and never wasting computer cycles when nothing is happening. It is easy to understand the scheduling algorithm used here with a physical analogy. Suppose we have half a million events to be scheduled over the next five years, and that they appear more-or-less randomly throughout that period. We want to be able to (1) schedule new events, (2) cancel existing events, and (3) notify some dispatcher as the time for each event arrives. Image a series of bins, one per minute for an entire year. The first bin represents the first minute after midnight on New Year's Day, the second bin represent the second minute, and so forth to the last bin, which represents the last minute on December 31st of the year. That is 366 days x 24 hours/day x 60 minutes/hour = 527,040 bins total, each labeled with the month, day, hour, and minute that it represents. With half a million events to be scheduled, that is less than one per bin on average. Also imagine that each event has a ticket with its event number and its scheduled time, represented at least to the nearest second. Then here is the procedure for handling events: 1. Scheduling a new event. a. Go to the bin representing the month, day, hour, and minute for the event, ignoring the year and the second. b. Drop the event's ticket into the bin. c. That took only a single operation. 2. Cancelling an existing event. a. Go to the bin representing the month, day, hour, and minute for the event, ignoring the year and the second. b. Remove all tickets from the bin and flip through them to find the ticket for the event in question. c. Destroy that ticket and place the others back in the bin. d. That took one operation per ticket in the bin. But on average there is only one ticket in the bin. 3. Dispatching the next event. a. Go to the bin representing the current time of day. b. Remove the tickets from the bin and sort them into chronological order. c. Place any tickets for future years back in the bin. d. Hand each ticket for this year to a dispatcher as its scheduled time arrives. e. If any tickets for the current bin arrive while the bin is being dispatched, put them in their proper position among the other tickets for dispatching now or in the future. f. That took one operation per ticket in the bin. But on average there is only one ticket in the bin. The basic algorithm underlying these procedures was discovered in the early days of computer science. Today it would be called a "hash coded search with modulus key." It allows near-instantaneous access to any entry and keeps the entries in order, modulo the number of entries in the list. The first paper clearly describing this kind of search is by Arnold Dumey (1956), entitled "Indexing for rapid random access memory systems," in Computers and Automation 6:6-9. It is interesting that "rapid" in that day meant something like 1/3 second per access! The running time of the above three procedures, remarkably, is not a function of the number of pending events! It is only a function of the ratio between the number of pending events and the number of bins. Therefore, with sufficient memory, the procedures are maximally efficient. In computer parlance, they are called "Order 1", or "O(1)". If there are a hundred million events pending randomly distributed into a hundred million bins, for example, the average number in each bin is 1 and the average in each occupied bin is $1/(1-1/e) = 1.58$. Therefore, essentially no searching/sorting occurs, making the method exceedingly fast. It also has modest space requirements, with only two additional integers needed per event (one for a head-of-list index and another for a forward list index). Therefore, assuming four bytes per integer, managing 100 million events requires less than 3/4 gigabyte of memory, well within the reach even of laptop computers. ($8\times10^8/2^{30}$ GB). Assuming the events occur at random times, then a portion $1/e = 0.36787944$ of the bins will be empty, the same portion will contain a single entry, $1/2e^2 = 0.18393972$ will contain two entries, and so on according to the Poisson distribution. Below is an actual example with one million bins and one million pending events. The expected numbers are Poisson values. The observed numbers are from the end of an SIR disease run, after many millions of events had been added, cancelled, and dispatched. There is no significant clustering to slow things down. Here practice corresponds well with theory. N Observed Expected 0 371047 367879 1 364903 367879 2 181992 183940 3 62032 61313 4 15998 15328 5 3256 3066 6 675 511 7 82 73 8 14 9 9 1 1 As implemented in this module, The bins need not correspond to standard time units such as a minutes. The bin representing the current time is kept sorted so that each event to be dispatched need only be picked from the front of the list. That means that any time a new event is added to the current bin, that bin must be resorted. In this module the three functions that correspond to the above are 'EventSchedule', 'EventCancel,' and 'EventNext'. They work essentially as described above. */ /*----------------------------------------------------------------------------* DATA STRUCTURES The module has two main data structures, a circular series of time bins 'Q[h]', which index one event scheduled for the time slots 'h' represented by each bin, and a series of singly linked lists 'P[i]', which define the earliest event for each individual, linked to their proper time slots. Note that each bin 'Q[h]' represents many related time slots, all equal modulo the width of the series of time bins, 'Qw', described further below. The number of links 'P[i]' must be equal to the number of individuals, and is indexed by individual number. but the number of time bins 'Q[h]' may be smaller or larger than the number of individuals. The size of 'Q' is a matter of optimization. It is typical to have one time bin for each event that could be scheduled, meaning each bin will represent a single event on average. That is optimal if all bin operations are of equal speed. The width 'Qw' of all bins combined is also a matter of optimization. If it is much too large, events will tend to cluster near the bin being dispatched. If it is much too small, events will tend to spread out, with most bins containing events that are for the more distant future, not immediate events. A suitable value for 'Qw' can be found by knowledge of the system being analysed, or by experimental trials to find a good speed of operation. The linked lists of events in 'P' are not maintained in any particular order, for speed of addition, but before a bin is dispatched, it is sorted into order by event times, for speed of dispatching. */ #include "common.h" #define PINIT if(run1==0) EventInit(); #define PEMPTY -1 //Marker for bins containing no linkages. #define TN (INDIV+0) //Maximum number of time bins. #define PN (INDIV+3) //Maximum number of time bin forward indexes. #define TW 20 //Time width of all bins combined (for optimization). dec t; //Current time, last dispatched event. static int run1; //Flag to detect if the routine is being reused. static dec T[PN]; //Time for each scheduled event. static int P[PN]; //Forward indexes within bins, ending with zero. static int Q[TN]; //First index for the bin, with zero for empty bins. static int Qn = TN; //Number of elements in 'Q'. static dec Qw = TW; //Interval of time represented for each cycle in 'Q'. static int Qi = 0; //Index of the immediate time bin. static int Qo = 1; //Flag set if the immediate bin is in order. static int Qe = 0; //Number of events in all bins. static dec Qt0 = 0; //Earliest time representable this cycle in 'Q'. static dec Qt1 = TW; //Earliest time beyond this cycle in 'Q'. /*----------------------------------------------------------------------------* INITIALIZE STATIC DATA STRUCTURES This routine must be called if an entirely new simulation is to be started in the absence of the program being restarted from the beginning. It makes the module serially reusable, and is included to cover a deficiency in MPI whereby 'system' and 'popen' corrupt the system. ENTRY: It is time to start an entirely new run. EXIT: Any prior data have been wiped clean. */ EventInit() { int i; for(i=0; i0' Element 'i' occurs after element 'j'. */ int o1(int i, int j) { dec w = T[i]-T[j]; return w<0? -1: w>0? 1: 0; } /*----------------------------------------------------------------------------* LOCATE NEXT EVENT ENTRY: 'T' contains the time of the next event for every event scheduled. The bin structure is properly initialized. EXIT: 'EventNext' contains the number of the next event. If zero, no events are scheduled. 't' contains the time of the next event, if 'NextEvent' is not zero. WORK: The scheduling data structures are prepared as described above. */ int EventNext() { int j, n; PINIT; //Initialize if necessary. while(Qe>0) { for(; Qi=PN) Error1(734.1, "n=",n); //Check the index and make sure an if(P[n]!=PEMPTY) Error1(735.1, "n=",n); //event is not already scheduled if(te",te); //and is not in the past. T[n] = te; //Record the time of the new event. tr = (te-Qt0)/Qw; tr -= (int)tr; //Convert the time to a bin number i = tr*Qn; if(i==Qi) Qo = 0; //and mark for sorting if needed. P[n] = Q[i]; Q[i] = n; //Add the event to the list for Qe += 1; //that bin and increment the number } //of events. /*----------------------------------------------------------------------------* CANCEL EXISTING EVENT ENTRY: 'n' contains the number (starting with 1) of the event to be cancelled. 'T[n]' contains the scheduled time of the event. EXIT: The event has been removed from the list. WORK: The scheduling data structures are prepared as described above. */ EventCancel(int n) { int i, j, jp; dec tr; PINIT; //Initialize if necessary. if(n<1||n>=PN) Error1(734.2, "n=",n); //Check the index and make sure an if(P[n]==PEMPTY) Error1(736.2, "n=",n); //event is scheduled. tr = (T[n]-Qt0)/Qw; tr -= (int)tr; //Convert the time to a bin number, i = tr*Qn; //modulo the duration of the cycle. if(cancel1(n, i)) return; //Remove it from its normal bin. i = (i-1+Qn) % Qn; //If it is not there, check the if(cancel1(n, i)) return; //bin below (from rounding error). i = (i+2+Qn) % Qn; //Finally, check the bin above. if(cancel1(n, i)) return; Error1(818., "n=",n); //If the specified event was not in } //the list, signal an error. int cancel1(int n, int i) { int j, jp; for(jp=0,j=Q[i]; j>0; jp=j,j=P[j]) //Scan the list of pending events if(j==n) //in this bin and remove the { if(jp>0) P[jp] = P[j]; //specified event. (The average else Q[i] = P[j]; //number events in a non-empty bin P[j] = PEMPTY; //is 1.5) Qe -= 1; if(Qe<0) Error2(819., "n=",n, " bin=",i); return 1; } return 0; } /* Note: It is necessary to check not just the expected bin when an event is to be deleted, but the two adjacent bins as well if the event does not appear. This is because of the vagaries of finite-precision computer arithmetic. Two bins are separated by a knife-edge, and what is calculated one time as 1.0 may be calculated in a slightly different manner as 0.999999999999999, dropping it in an adjacent bin. In simulations testing this condition, it happened only once per ten million times or less, so there is no speed issue here, only a correctness issue. */ /*----------------------------------------------------------------------------* RENUMBER EXISTING EVENT This routine renumbers events. It can be called, for example, to reuse an entry when it becomes available, for example from a simulated death, shifting an existing individual from the end of the array keep the array compact. ENTRY: 'n' contains the new index number, which has no event scheduled. 'm' contains the current index number of the event. There is an event scheduled for 'm'. EXIT: 'n' is the new index number of the individual. The event originally scheduled as 'm' is re-scheduled as 'n'. Event 'm' no longer has an event scheduled and the index is free to be reused. */ EventRenumber(int n, int m) { if(n<1||n>=PN) Error1(734.3, "n=",n); //Check the indexes and make sure if(m<1||m>=PN) Error1(734.4, "n=",m); //they are in range. if(n!=m) { T[n] = T[m]; //Transfer the time. EventCancel(m); //Cancel the old number. EventSchedule(n, T[n]); } //Reschedule as the new number. } /*----------------------------------------------------------------------------* DISPLAY PROFILE This routine displays the distribution of schedule events---how frequently 0, 1, 2, 3, ... events occur in bins, and compares it with a Poisson distribution, $exp(-x) x^n / n!$ ENTRY: 'label' contains an optional label to be displayed with the results, such as "Initial" or "Final". The bins have been initialized. EXIT; 'EventProfile' contains the amount of memory allocated for the main data structures. */ #define PROF 1001 int EventProfile(char *label) { int i, j, n, imax, prof[PROF]; dec nexp, lambda, eml, ln, nf; if(label==0||label[0]==0) label = "Bin"; //Establish a default label. for(i=0; i0; j=P[j],n++) //have no entries, one entry, two if(j<1||j>=PN||n>PN) //entries, etc. Error1(820.2, "j=",j); if(n>PROF-1) n = PROF-1; prof[n] += 1; } for(i=imax=0; i0.5) //randomness of the Poisson. printf("%4d%c%9d %10.0f\n", i, i