/*
 * $Id: mm1.c,v 1.3 2009/02/10 02:14:44 dxt22 Exp $
 *
 * Date: February 9, 2009
 * Modified by: Dmitri Tikhonov <dxt22@njit.edu>
 *
 * CS661, Homework 2.
 *
 * DESCRIPTION:
 *
 * The modifications transformed the original program from a simulation of
 * a single server with a single queue to a simulation of a series of servers,
 * each with a single queue:
 *
 *  +----------+        +----------+         +----------+
 *  | SERVER 1 | --->>  | SERVER 2 |  >...>  | SERVER N |
 *  +----------+        +----------+         +----------+
 *
 * Parts enter the system by arriving at server 1 queue and, after being
 * serviced, go to the next server and so on until the last server.
 *
 * PARTICULARS:
 *
 * - Constant NSERVERS defines the number of servers for the system.  The
 *   program must be recompiled for the changes to take effect.  Supported
 *   are values of 1 and higher.  NSERVERS = 1 reduces the program to the
 *   original.
 *
 * - To get a bunch of debugging output, turn on TRACE constant.
 *
 * - In order to get the same results for SERVER 1 in a run with any number
 *   of servers (as long as the same parameters are provided, of course),
 *   function expon() has been modified to allocate a stream of random values
 *   per server.
 *
 * - arrive() function has been modified to only schedule EVENT_NEW_ARRIVAL
 *   for the first server.
 *
 * The latest version of this program is available on the web:
 *   http://web.njit.edu/~dxt22/cs661/
 */

/* External definitions for single-server queueing system, fixed run length. */

#include <assert.h>

#include <stdio.h>

#include <stdlib.h>

#include <math.h>


#define Q_LIMIT 100  /* Limit on queue length. */

#define BUSY      1  /* Mnemonics for server's being busy */

#define IDLE      0  /* and idle. */



/* The following 3 declarations are for use of the random-number generator
   lcgrand and the associated functions lcgrandst and lcgrandgt for seed
   management.  This file (named lcgrand.h) should be included in any program
   using these functions by executing
       #include "lcgrand.h"
   before referencing the functions. */

float lcgrand(int stream);
void  lcgrandst(long zset, int stream);
long  lcgrandgt(int stream);



/* Number of queue/server pairs in the system: server 1 gets index 0,
 * server 2 gets index 1, etc.
 */
#define NSERVERS 2


#define NEVENTS (2 + NSERVERS)


#define EVENT_EMPTY_QUEUE 0

/* New arrivals are to server 1 only */
#define EVENT_NEW_ARRIVAL 1

/* Generate any number of departures in between */
#define EVENT_DEPART_SERVER(_i_) ((_i_) + EVENT_NEW_ARRIVAL + 1)

#define CONVERT_EVENT_TO_SERVER_ID(_event_) ((_event_) - EVENT_NEW_ARRIVAL - 1)

#define EVENT_LEAVE_SYSTEM (NEVENTS - 1)

#define EVENT_END          NEVENTS


static char _event_name_str[0x20];
#define EVENT_TO_STRING(_e_) ((const char *) (                          \

    EVENT_EMPTY_QUEUE  == (_e_) ? "EVENT_EMPTY_QUEUE"   :               \
    EVENT_NEW_ARRIVAL  == (_e_) ? "EVENT_NEW_ARRIVAL"   :               \
    EVENT_LEAVE_SYSTEM == (_e_) ? "EVENT_LEAVE_SYSTEM"  :               \
    EVENT_END          == (_e_) ? "EVENT_END"           :               \
    (snprintf(_event_name_str, sizeof(_event_name_str),                 \
        "EVENT_DEPART_SERVER_%d", CONVERT_EVENT_TO_SERVER_ID(_e_) + 1), \
        _event_name_str)))

int   num_custs_delayed[NSERVERS] = {}
    , num_in_q[NSERVERS] = {}
    , server_status[NSERVERS] = {}
    ;

float area_num_in_q[NSERVERS] = {}
    , area_server_status[NSERVERS] = {}
    , mean_service[NSERVERS] = {}
    , mean_interarrival // Only for server #1, others are dictated by
                        // preceding queues.
    , sim_time
    , time_arrival[NSERVERS][Q_LIMIT + 1] = {}
    , time_end
    , time_last_event
    , time_next_event[1 + NEVENTS] = {}
    , total_of_delays[NSERVERS] = {}
    ;

FILE  *infile, *outfile;

void  initialize(void);

/* Returns next event type */
int   timing(void);

/* Arrive to server id `id'.  If first server (id = 0), new arrival event
 * is scheduled.
 */
void  arrive(int id);

/* Depart from server `id'.
 */
void  depart(int id);

void  report(void);
void  update_time_avg_stats(void);

/* Since we want to keep the same numbers for the first server no matter
 * which kind of simulation we run: one server, two servers, etc., each
 * server will get its own set of random numbers.  `id' in this case is 0
 * for server 1, 1 for server 2, and so on.
 */
float expon(float mean, int id);

/* Used for debugging */
void  short_report(void);

#define TRACE 0


#if TRACE /* [ */

static void
print_state (void)
{
    printf("TIME: %.3f -- ", sim_time);
    int i;
    for (i = 0; i < NSERVERS; ++i) {
        printf("  SERVER %d: %s, NQ: %d", i + 1,
            (server_status[i] ? "BUSY" : "IDLE"), num_in_q[i]);
    }
    printf("\n");
    return;
}

static void
print_event_list (void)
{
    int i;
    printf("EVENTS: ");
    for (i = 1; i <= NEVENTS; ++i)
        printf("%s: %f; ", EVENT_TO_STRING(i), time_next_event[i]);
    printf("\n");
}
#else

#define print_state()

#define print_event_list()

#endif /* ] TRACE */



int
main(void)
{
    /* Open input and output files. */

    infile  = fopen("mm1.in",  "r");
    outfile = fopen("mm1.out", "w");

    /* Read input parameters. */

    fscanf(infile, "%f", &mean_interarrival);
    int i;
    for (i = 0; i < NSERVERS; ++i)
        fscanf(infile, " %f", &mean_service[i]);
    fscanf(infile, " %f", &time_end);


    /* Write report heading and input parameters. */

    if (1 == NSERVERS) {
        fprintf(outfile, "Single-server queueing system with fixed run");
        fprintf(outfile, " length\n");
    } else {
        fprintf(outfile, "%d-server queueing system (in series), fixed run "
                         "length\n", NSERVERS);
    }

    fprintf(outfile, "Mean interarrival time%11.3f minutes\n",
            mean_interarrival);

    if (1 == NSERVERS) {
        fprintf(outfile, "Mean service time %.3f minutes\n\n", mean_service[0]);
    } else {
        int i;
        for (i = 0; i < NSERVERS; ++i)
            fprintf(outfile, "Mean service time at server #%d: %.3f minutes\n",
                i + 1, mean_service[i]);
    }
    fprintf(outfile, "Length of the simulation: %.3f minutes\n\n", time_end);

    /* Initialize the simulation. */

    initialize();

    /* Run the simulation until it terminates after an end-simulation event
       (EVENT_END) occurs. */

    int next_event_type;
    do {
        print_event_list();
        print_state();

        /* Determine the next event. */

        next_event_type = timing();
        assert(0 <= next_event_type && next_event_type <= NEVENTS);
        if (EVENT_EMPTY_QUEUE == next_event_type) {
            fprintf(outfile, "\nEvent list empty at time %f\n", sim_time);
            break;
        }
        if (TRACE)
            printf("-> Got event %s\n", EVENT_TO_STRING(next_event_type));

        /* Update time-average statistical accumulators. */

        update_time_avg_stats();

        /* Invoke the appropriate event function. */

        switch (next_event_type) {
            case EVENT_NEW_ARRIVAL:
                arrive(0);
                break;
            case EVENT_LEAVE_SYSTEM:
                depart(NSERVERS - 1);
                break;
            case EVENT_END:
                /* Do nothing, will break out of the loop, see while cond */
                break;
            default: {
                assert(next_event_type > EVENT_NEW_ARRIVAL);
                assert(next_event_type < EVENT_LEAVE_SYSTEM);

                /* Go from one server to the next */
                const int id = CONVERT_EVENT_TO_SERVER_ID(next_event_type);
                depart(id);
                arrive(id + 1);
                break;
            }
        }

        if (TRACE)
            short_report();

    /* If the event just executed was not the end-simulation event (type 3),
       continue simulating.  Otherwise, end the simulation. */

    } while (next_event_type != EVENT_END);

    report();

    fclose(infile);
    fclose(outfile);

    return 0;
}


void initialize(void)  /* Initialization function. */
{
    /* Initialize the simulation clock. */
    sim_time = 0.0;

    time_last_event = 0.0;

    int i;
    for (i = 0; i < NSERVERS; ++i) {
        /* Initialize the state variables. */

        server_status[i]   = IDLE;
        num_in_q[i]        = 0;

        /* Initialize the statistical counters. */

        num_custs_delayed[i]  = 0;
        total_of_delays[i]    = 0.0;
        area_num_in_q[i]      = 0.0;
        area_server_status[i] = 0.0;
    }

    /* Initialize event list.  Since no customers are present, the departure
       (service completion) event is eliminated from consideration.  The end-
       simulation event (type 3) is scheduled for time time_end. */

    time_next_event[EVENT_NEW_ARRIVAL] = sim_time + expon(mean_interarrival, 0);
    for (i = EVENT_NEW_ARRIVAL + 1; i < EVENT_END; ++i)
        time_next_event[i] = 1.0e+30;
    time_next_event[EVENT_END] = time_end;
}


int
timing (void)
{
    int   i;
    float min_time_next_event = 1.0e+29;

    int next_event_type = 0;

    /* Determine the event type of the next event to occur. */

    for (i = 1; i <= NEVENTS; ++i) {
        /*
        fprintf(stderr, "%s: i: %d, time: %f\n", __FUNCTION__, i,
            time_next_event[i]);
            */
        if (time_next_event[i] < min_time_next_event) {
            min_time_next_event = time_next_event[i];
            next_event_type     = i;
        }
    }

    /* Check to see whether the event list is empty. */

    if (next_event_type == 0)
        return EVENT_EMPTY_QUEUE;

    /* The event list is not empty, so advance the simulation clock. */

    sim_time = min_time_next_event;
    return next_event_type;
}


void
arrive (int id)  /* Arrival event function. */
{
    if (0 == id) {
        /* Only schedule next arrival from outside the system if this is
         * server #1.
         */
        time_next_event[EVENT_NEW_ARRIVAL] =
            sim_time + expon(mean_interarrival, id);
    }

    /* Check to see whether server is busy. */

    if (server_status[id] == BUSY) {

        /* Server is busy, so increment number of customers in queue. */

        ++num_in_q[id];

        /* Check to see whether an overflow condition exists. */

        if (num_in_q[id] > Q_LIMIT) {

            /* The queue has overflowed, so stop the simulation. */

            fprintf(outfile, "\nOverflow of the array time_arrival at");
            fprintf(outfile, " time %f", sim_time);
            exit(2);
        }

        /* There is still room in the queue, so store the time of arrival of the
           arriving customer at the (new) end of time_arrival. */

        time_arrival[id][num_in_q[id]] = sim_time;
    }

    else {

        /* Server is idle, so arriving customer has a delay of zero.  (The
           following two statements are for program clarity and do not affect
           the results of the simulation.) */

        float delay            = 0.0;
        total_of_delays[id] += delay;

        /* Increment the number of customers delayed, and make server busy. */

        ++num_custs_delayed[id];
        // printf("Server %d, delayed %d\n", id, num_custs_delayed[id]);
        server_status[id] = BUSY;

        /* Schedule a departure (service completion). */

        time_next_event[EVENT_DEPART_SERVER(id)] =
            sim_time + expon(mean_service[id], id);
    }
}


void
depart (int id)  /* Departure event function. */
{
    /* Check to see whether the queue is empty. */

    if (num_in_q[id] == 0) {

        /* The queue is empty so make the server idle and eliminate the
           departure (service completion) event from consideration. */

        server_status[id]      = IDLE;
        time_next_event[EVENT_DEPART_SERVER(id)] = 1.0e+30;
    }

    else {

        /* The queue is nonempty, so decrement the number of customers in
           queue. */

        --num_in_q[id];

        /* Compute the delay of the customer who is beginning service and update
           the total delay accumulator. */

        const float delay    = sim_time - time_arrival[id][1];
        total_of_delays[id] += delay;

        /* Increment the number of customers delayed, and schedule departure. */

        ++num_custs_delayed[id];
        time_next_event[EVENT_DEPART_SERVER(id)] =
            sim_time + expon(mean_service[id], id);

        /* Move each customer in queue (if any) up one place. */

        int i;
        for (i = 1; i <= num_in_q[id]; ++i)
            time_arrival[id][i] = time_arrival[id][i + 1];
    }
}


void
report(void)  /* Report generator function. */
{
    /* Compute and write estimates of desired measures of performance. */

    int id;
    for (id = 0; id < NSERVERS; ++id) {
        fprintf(outfile, "Queue/Server %d:\n", id + 1);
        fprintf(outfile, "Average delay in queue: %11.3f minutes\n",
                total_of_delays[id] / num_custs_delayed[id]);
        fprintf(outfile, "Average number in queue: %10.3f\n",
                area_num_in_q[id] / sim_time);
        fprintf(outfile, "Server utilization: %15.3f\n",
                area_server_status[id] / sim_time);
        fprintf(outfile, "Number of delays completed:   %d\n",
                num_custs_delayed[id]);
        fprintf(outfile, "\n\n");
    }
}


/* Compact form of the report: used for debugging */
void
short_report (void)
{
    printf("SHORT REPORT:\n");
    int id;
    for (id = 0; id < NSERVERS; ++id) {
        printf("S #%d: AD %.3f; AQ: %.3f; SU: %.3f; NDC: %d\n",
            id + 1,
            total_of_delays[id] / num_custs_delayed[id],
            area_num_in_q[id] / sim_time,
            area_server_status[id] / sim_time,
            num_custs_delayed[id]
        );
    }
}

void update_time_avg_stats(void)  /* Update area accumulators for time-average
                                     statistics. */
{
    float time_since_last_event;

    /* Compute time since last event, and update last-event-time marker. */

    time_since_last_event = sim_time - time_last_event;
    time_last_event       = sim_time;

    int i;
    for (i = 0; i < NSERVERS; ++i) {
        /* Update area under number-in-queue function. */

        area_num_in_q[i]      += num_in_q[i] * time_since_last_event;

        /* Update area under server-busy indicator function. */

        area_server_status[i] += server_status[i] * time_since_last_event;
    }
}


float expon(float mean, int id)  /* Exponential variate generation function. */
{
    /* Return an exponential random variate with mean "mean". */

    return -mean * log(lcgrand(id + 1));
}





/* Prime modulus multiplicative linear congruential generator
   Z[i] = (630360016 * Z[i-1]) (mod(pow(2,31) - 1)), based on Marse and Roberts'
   portable FORTRAN random-number generator UNIRAN.  Multiple (100) streams are
   supported, with seeds spaced 100,000 apart.  Throughout, input argument
   "stream" must be an int giving the desired stream number.  The header file
   lcgrand.h must be included in the calling program (#include "lcgrand.h")
   before using these functions.

   Usage: (Three functions)

   1. To obtain the next U(0,1) random number from stream "stream," execute
          u = lcgrand(stream);
      where lcgrand is a float function.  The float variable u will contain the
      next random number.

   2. To set the seed for stream "stream" to a desired value zset, execute
          lcgrandst(zset, stream);
      where lcgrandst is a void function and zset must be a long set to the
      desired seed, a number between 1 and 2147483646 (inclusive).  Default
      seeds for all 100 streams are given in the code.

   3. To get the current (most recently used) integer in the sequence being
      generated for stream "stream" into the long variable zget, execute
          zget = lcgrandgt(stream);
      where lcgrandgt is a long function. */

/* Define the constants. */

#define MODLUS 2147483647

#define MULT1       24112

#define MULT2       26143


/* Set the default seeds for all 100 streams. */

static long zrng[] =
{         1,
 1973272912, 281629770,  20006270,1280689831,2096730329,1933576050,
  913566091, 246780520,1363774876, 604901985,1511192140,1259851944,
  824064364, 150493284, 242708531,  75253171,1964472944,1202299975,
  233217322,1911216000, 726370533, 403498145, 993232223,1103205531,
  762430696,1922803170,1385516923,  76271663, 413682397, 726466604,
  336157058,1432650381,1120463904, 595778810, 877722890,1046574445,
   68911991,2088367019, 748545416, 622401386,2122378830, 640690903,
 1774806513,2132545692,2079249579,  78130110, 852776735,1187867272,
 1351423507,1645973084,1997049139, 922510944,2045512870, 898585771,
  243649545,1004818771, 773686062, 403188473, 372279877,1901633463,
  498067494,2087759558, 493157915, 597104727,1530940798,1814496276,
  536444882,1663153658, 855503735,  67784357,1432404475, 619691088,
  119025595, 880802310, 176192644,1116780070, 277854671,1366580350,
 1142483975,2026948561,1053920743, 786262391,1792203830,1494667770,
 1923011392,1433700034,1244184613,1147297105, 539712780,1545929719,
  190641742,1645390429, 264907697, 620389253,1502074852, 927711160,
  364849192,2049576050, 638580085, 547070247 };

/* Generate the next random number. */

float lcgrand(int stream)
{
    long zi, lowprd, hi31;

    zi     = zrng[stream];
    lowprd = (zi & 65535) * MULT1;
    hi31   = (zi >> 16) * MULT1 + (lowprd >> 16);
    zi     = ((lowprd & 65535) - MODLUS) +
             ((hi31 & 32767) << 16) + (hi31 >> 15);
    if (zi < 0) zi += MODLUS;
    lowprd = (zi & 65535) * MULT2;
    hi31   = (zi >> 16) * MULT2 + (lowprd >> 16);
    zi     = ((lowprd & 65535) - MODLUS) +
             ((hi31 & 32767) << 16) + (hi31 >> 15);
    if (zi < 0) zi += MODLUS;
    zrng[stream] = zi;
    return (zi >> 7 | 1) / 16777216.0;
}


void lcgrandst (long zset, int stream) /* Set the current zrng for stream
                                          "stream" to zset. */
{
    zrng[stream] = zset;
}


long lcgrandgt (int stream) /* Return the current zrng for stream "stream". */
{
    return zrng[stream];
}



syntax highlighted by Code2HTML, v. 0.9.1