/* Discrete-event simulation of Banff park entry queue   */
/*                                                       */
/* Written by Carey Williamson, University of Calgary    */
/*                                                       */
/* Usage: cc -o parksim parksim.c -lm                    */
/*        ./parksim                                      */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h> /* use man 3 log for the math functions */

#define NUM_CARS 10000   /* Number of cars to simulate */
#define END_OF_TIME 999999.9 /* Time far in the future */

/* Debugging output */
/* #define DEBUG 1 */

/* Queue stats */
/* #define STATS 1 */
#define REPORT_INTERVAL 1000.0
/* #define BATCHES 1 */
#define WARMUP_INTERVAL 0.0

/***********************************************************************/
/*                 RANDOM NUMBER GENERATION STUFF                      */
/***********************************************************************/

/* Parameters for random number generation. */
#define MAX_INT 2147483647       /* Maximum positive integer 2^31 - 1 */

/* Generate a random floating point value uniformly distributed in [0,1] */
float Uniform01()
  {
    float randnum;

    /* get a random positive integer from random() */
    randnum = (float) 1.0 * random();

    /* divide by max int to get something in 0..1  */
    randnum = (float) randnum / (1.0 * MAX_INT); 

    return( randnum );
  }

/* Generate a random floating point number from an exponential    */
/* distribution with mean mu.                                     */

float Exponential(mu)
    float mu;
  {
    float randnum, ans;

    randnum = Uniform01();
    ans = -(mu) * log(randnum);

    return( ans );
  }

/***********************************************************************/
/*                 MAIN PROGRAM                                        */
/***********************************************************************/

int main()
  {
    int i, arrivals, departures, qsize, frontcar, maxq;
    float now, nextarrival, nextdeparture, fdone, gdone, rand;
    float arrivaltimes[NUM_CARS], servicetimes[NUM_CARS];
    float waittimes[NUM_CARS], sojourntimes[NUM_CARS];
    float meanwait, varwait, stdwait;
    float last, meanqueue, varqueue;
    float nextreport, lastreport;
    int warmupdone;

    /* Seed the PRNG */
    srandom(1234567);       /* for a fixed seed every time you run it */
    srandom(time(NULL));    /* for different seed every time you run it */

    /* Make huge arrays with all the timing information */
    for( i = 0; i < NUM_CARS; i++ )
      {
	if( i == 0 )
	  arrivaltimes[i] = Exponential(2.0);
	else arrivaltimes[i] = arrivaltimes[i-1] + Exponential(2.0/1.0);
#ifdef DETERMINISTIC
	servicetimes[i] = 1.5;
#endif
#ifdef HYPEREXPONENTIAL
	if( Uniform01() < 0.5 )
	  servicetimes[i] = Exponential(1.0);
	else servicetimes[i] = Exponential(2.0);
#endif
	servicetimes[i] = Exponential(1.5);
#ifdef DEBUG
	printf("Car %d arrival time %8.6f service time %8.6f\n",
	       i, arrivaltimes[i], servicetimes[i]);
#endif
      }

    /* Initialization */
    arrivals = 0;
    departures = 0;
    qsize = 0;
    meanqueue = 0.0;
    varqueue = 0.0;
    last = 0.0;
    maxq = 0;
    now = 0.0;
    frontcar = 0;
    nextarrival = arrivaltimes[frontcar];
    nextdeparture = END_OF_TIME;
    nextreport = REPORT_INTERVAL;
    lastreport = 0.0;
    warmupdone = 0;

#ifdef DEBUG
    printf("Beginning trace-driven simulation now...\n");
#endif

    while( departures < NUM_CARS )
      {
	/* determine next event */
	if( nextarrival < nextdeparture )
	  {
	    now = nextarrival;
#ifdef DEBUG
	    printf("Car %d arrives at time %8.6f\n", arrivals, now);
#endif

	    /* Update and/or report the current queue size */
	    meanqueue += (now - last)*qsize;
	    varqueue += (now - last)*qsize*qsize;
	    last = now;
	    qsize++;
#ifdef DEBUG
	    printf("Lineup at booth now has %d cars!\n", qsize);
#endif
#ifdef STATS
	    printf("%8.6f %d\n", now, qsize);
#endif
	    if( qsize > maxq )
	      maxq = qsize;

	    /* schedule departure event, if queue was empty before */
	    if( qsize == 1 )
	      {
		frontcar = arrivals;
		nextdeparture = now + servicetimes[frontcar];
#ifdef DEBUG
		printf("Car %d will be completing service at time %8.6f\n",
		       frontcar, nextdeparture);
#endif
	      }

	    /* schedule next arrival, if any */
	    arrivals++;
	    if( arrivals < NUM_CARS )
	      nextarrival = arrivaltimes[arrivals];
	    else nextarrival = END_OF_TIME;
	  }
	else   /* departure */
	  {
#ifdef DEBUG
	    printf("Next event is service completion for car %d at time %8.6f\n",
		   frontcar, nextdeparture);
#endif
	    now = nextdeparture;
	    sojourntimes[frontcar] = now - arrivaltimes[frontcar];
	    waittimes[frontcar] = sojourntimes[frontcar] - servicetimes[frontcar];
	    departures++;
	    meanqueue += (now - last)*qsize;
	    varqueue += (now - last)*qsize*qsize;
	    last = now;
	    qsize--;
#ifdef DEBUG
	    printf("Lineup at booth now has %d cars\n", qsize);
#endif
#ifdef STATS
	    printf("%8.6f %d\n", now, qsize);
#endif
	    if( qsize > 0 )
	      {
		frontcar = departures;
		nextdeparture = now + servicetimes[frontcar];
#ifdef DEBUG
		printf("Front car %d has service time %8.6f departure %8.6f\n",
		       frontcar, servicetimes[frontcar], nextdeparture);
#endif
	      }
	    else
	      {
		nextdeparture = END_OF_TIME;
	      }
	  }
	/* check to see if data dump is needed */
	if( now > nextreport )
	  {
#ifdef BATCHES
	    printf("Batch Time: %8.6f  Mean queue size: %5.3f\n",
		   now, meanqueue/(now - lastreport));
	    meanqueue = 0.0;
	    varqueue = 0.0;
	    lastreport = now;
#else
	    if( now > WARMUP_INTERVAL + REPORT_INTERVAL )
	      printf("Cumulative Time: %8.6f  Mean queue size: %5.3f\n",
		     now, meanqueue/(now - WARMUP_INTERVAL));
	    else printf("Cumulative Time: %8.6f  Mean queue size: %5.3f\n",
		   now, meanqueue/now);
#endif
	    nextreport += REPORT_INTERVAL;
	  }

	/* see if warmup period is done or not */
	if( (warmupdone == 0) && (now > WARMUP_INTERVAL) )
	  {
	    meanqueue = 0.0;
	    varqueue = 0.0;
	    warmupdone++;
	  }

      }

    printf("\n");
    printf("Banff Park Entry Simulation Model:\n");
    printf("Number of cars: %d\n", NUM_CARS);
    printf("Time: %8.6f  Arrivals: %d  Departures: %d  Queue: %d\n",
	   now, arrivals, departures, qsize);
    meanqueue = meanqueue/(now-WARMUP_INTERVAL);
    varqueue = varqueue/(now-WARMUP_INTERVAL);
    varqueue -= meanqueue*meanqueue;
    printf("Mean queue size: %5.3f\n", meanqueue);
    printf("Variance of queue size: %5.3f\n", varqueue);
    printf("Maximum queue size observed: %d\n", maxq);

#ifdef DEBUG
    /* print detailed report */
    printf(" Car    Arrival     Wait        Service     Sojourn\n");
    for( i = 0; i < NUM_CARS; i++ )
      {
	printf(" %3d  %10.6f  %10.6f  %10.6f  %10.6f\n",
	       i, arrivaltimes[i], waittimes[i],
	       servicetimes[i], sojourntimes[i]);
      }
#endif
    
    /* calculate and report statistics */
    meanwait = 0.0;
    for( i = 0; i < NUM_CARS; i++ )
      meanwait += waittimes[i];
    meanwait /= NUM_CARS;
    varwait = 0.0;
    for( i = 0; i < NUM_CARS; i++ )
      varwait += (waittimes[i] - meanwait)*(waittimes[i] - meanwait);
    varwait /= NUM_CARS;
    stdwait = sqrt(varwait);
    printf("Mean wait time: %8.6f\n", meanwait);
    printf("Standard deviation of wait time: %8.6f\n", stdwait);
  }