ex84x.c

/* ------------------------------------------------------------------------- 
 * This program is a next-event simulation of a single-server FIFO service
 * node with immediate feedback using Exponentially distributed interarrival
 * times and service times (i.e., a M/M/1 queue).  
 *
 * The program provides steady-state statistics using the method of
 * batch means.
 *
 * Name            : ex84x.c
 * Author          : Steve Park & Dave Geyer
 * Additional Impl.: Rob Jones
 * Language        : ANSI C
 *
 * Change History:
 *   Date      Reason
 *   --------  -------------------------------------------------------------
 *   02-18-97  Revised to use method of batch means  
 *
 * Latest Revision : 02-18-97
 * ------------------------------------------------------------------------- 
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "rngs.h"              
#include "rvgs.h"
#include "rvms.h"
#include "sort.h"

#define PLOT
/* #define OUTPUT_FOR_ESTIMATE */

#define BATCHSIZE      16             /* batch size           */
#define NUMBATCHES     1024           /* number of batches    */
#define LOC            0.95           /* level of confidence  */
#define NEVER       1e300          /* a large number       */

static constexpr double lambda = 1.0;   /* job arrival rate     */
static constexpr double nu     = 1.25;  /* job service rate     */
static constexpr double BETA   = 0.0;   /* feedback probability */

int rns;

/* ============================== */
   static double Min(double a, double c)
/* ============================== */
{ 
  if (a < c)
    return (a);
  else
    return (c);
} 

/* ======================== */
   static double GetInterarrival()              
/* ======================== */          
{
  SelectStream(rns); 
  return (Exponential(1.0/lambda));
} 

/* =================== */
   static double GetService()                 
/* =================== */             
{
  SelectStream(rns);
  return (Exponential(1.0/nu));
} 

/* ================= */
   static int GetFeedback()
/* ================= */
{
  SelectStream(2);
  if (Random() < BETA)
    return (1);
  else
    return (0);
}
  
/* ============== */
   static int sssim(void)
/* ============== */
{
  struct {
    double arrival;         /* next arrival time                        */
    double completion;      /* next completion time                     */
    double current;         /* current time                             */
    double next;            /* next (most imminent) event time          */
    double interval;        /* time width of batch                      */
  } t;
  long   i;                 /* batch index                              */
  long   j;                 /* in-batch data index                      */
  long   number = 0;        /* number in the node                       */
  long   aindex;            /* used to count arriving jobs              */
  long   index;             /* used to count departed jobs              */
  long   seed = 123456789;  /* random number generator seed             */
  double *wait;             /* to gather sample wait times              */
  long   capacity;          /* wait[1..capacity] is allocated           */
  double aw;                /* batch mean of wait times                 */
  double avg_w, var_w;      /* mean, var of wait times from batch means */
  double delta, u, tstar = 0.0; 
  /* 
  *  initialize 
  */
  PlantSeeds(seed);
  /* wait[1..capacity], indexed by job number.  Arrivals run ahead of
     departures, so the array grows on demand.  (2026: was a fixed array
     of BATCHSIZE*NUMBATCHES, which the final arrivals overran.) */
  capacity = BATCHSIZE*NUMBATCHES;
  wait = (double *)calloc(capacity + 1, sizeof(double));
  if (wait == NULL) {
    fprintf(stderr, "Error allocating memory.\n");
    exit(EXIT_FAILURE);
  }
  t.current    = 0.0;                           /* set the clock             */
  t.arrival    = t.current + GetInterarrival(); /* schedule first arrival    */
  t.completion = NEVER;                      /* first event != completion */
  aindex = 0; 
  index  = 0; 
  avg_w  = 0.0;
  var_w  = 0.0;
  /* 
  *  constraint checks 
  */
  if (NUMBATCHES < 2) {
    printf("ERROR: must have at least 2 batch means.\n");
    return (0);
  }

  for (i = 1; i <= NUMBATCHES; i++) {
    /* 
    *  re-initialize the batch statistics
    */
    aw = 0.0;
    j  = 0;

    while (j < BATCHSIZE) {
      t.next = Min(t.arrival, t.completion);  /* next event time    */
      t.current = t.next;                     /* advance the clock  */

      if (t.current == t.arrival)  {          /* process an arrival */
        aindex++;
        if (aindex > capacity) {              /* make room for this job */
          double *grown = (double *)realloc(wait, (2*capacity + 1)*sizeof(double));
          if (grown == NULL) {
            fprintf(stderr, "Error allocating memory.\n");
            exit(EXIT_FAILURE);
          }
          wait      = grown;
          capacity *= 2;
        }
        /* store arrival time for this job */
        wait[aindex] = t.arrival;
        number++;
        t.arrival = t.current + GetInterarrival();
        if (number == 1)
          t.completion = t.current + GetService();
        } 
      else {                                  /* process a completion */
        if (GetFeedback())
          t.completion = t.current + GetService();
        else {
          index++;
          number--;
          /* store wait time for this job */
          wait[index] = t.completion - wait[index];

          /* Welford's algorithm to compute batch mean of wait times */
          aw += (wait[index] - aw)/(++j);

          if (number > 0)
            t.completion = t.current + GetService();
          else
            t.completion = NEVER;
         }
      }
    } /*end while*/ 
  /* 
  *  Welford's algorithm to compute mean/variance 
  *  of wait times from batch means...
  */
  delta  = aw - avg_w; 
  var_w += delta*delta*(i-1)/i; 
  avg_w += delta/i;

#ifdef OUTPUT_FOR_ESTIMATE
  printf("%.6f\n",aw); 
#endif

  }/* end for */

  i--;                                   /* i = the number of data points   */

  if (i > 1) {
    u     = 1.0 - 0.5*(1.0 - LOC);       /* interval parameter              */
    tstar = idfStudent(i-1,u);           /* critical value of t             */
  } 
  delta   = tstar*sqrt(var_w/(i*(i-1))); /* half-width of interval estimate */

#ifdef PLOT
  printf("%d %6.4f %6.4f\n", rns, avg_w, delta);
#endif
#ifndef PLOT 
#ifndef OUTPUT_FOR_ESTIMATE 
  printf("\n");
  printf("number of data points =  %ld\n", i);
  printf("batch size........... =  %ld\n", BATCHSIZE); 
  printf("confidence level..... =  %d%%\n", (int)(100.0*LOC+0.5)); 
  printf("arrival rate......... =  %g\n", lambda); 
  printf("service rate......... =  %g\n", nu); 
  printf("feedback prob........ =  %g\n", BETA); 
  printf("\n");
  printf("    AVERAGE WAIT..... = %7.4f +/- %6.4f\n", avg_w, delta);
  printf("\n");
#endif
#endif

  free(wait);

  return (0);
}


/* ============== */
   int main(void)
/* ============== */
{
  for (rns = 0; rns < 10; rns++) {
    sssim();
  }
  return (0);
}