// Author:  Douglas Wilhelm Harder
// Copyright (c) 2012 by Douglas Wilhelm Harder.  All rights reserved.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <assert.h>
#include "Queue.h"

struct data {
	double arrival_time;
	double start_time;
	double finish_time;
};




// Mean arrival rate
#define LAMBDA (85.0/60.0)

// Mean service rate
#define MU     (1.0/0.666)

#define RHO    ((LAMBDA)/(MU))

// Maximum number of processes generated
#define N      100000

int main() {
	assert( LAMBDA < MU );

	struct data request[N];

	srand48( time(NULL) );

	int i;

	/////////////////////////////////////////////////////////////
	// Generate the random arrivals from a Poisson Distrubtion //
	/////////////////////////////////////////////////////////////

	request[0].arrival_time = 0.0;

	for ( i = 1; i < N; ++i ) {
		// drand48() returns [0, 1);
		// we need (0, 1] so use 1 - drand48()

		request[i].arrival_time = request[i - 1].arrival_time
			- log(1.0 - drand48())/LAMBDA;
	}

	/////////////////////////////////////////////////////
	// Validate that the data is approximately Poisson //
	/////////////////////////////////////////////////////

	// Go out six standard deviations
	int max_poisson = LAMBDA + (int) ceil( 6*sqrt( LAMBDA ) );

	double poisson[max_poisson];
	for ( i = 0; i < max_poisson; ++i ) {
		poisson[i] = 0.0;
	}

	int j = 1;
	i = 0;

	while( 1 ) {
		int count = 0;

		while( request[i].arrival_time < j ) {
			++count;
			++i;

			if ( i == N ) {
				goto finish;
			}
		}

		if ( count < max_poisson ) {
			++( poisson[count] );
		}

		++j;
	}

	finish: ;

	// These two values should be roughly equal

	double factorial = 1.0;

	for ( i = 0; i  < max_poisson; ++i ) {
		printf( "%d\t%f %f\n", i, poisson[i]/j, pow( LAMBDA, i )/factorial*exp( -LAMBDA ) );

		factorial *= i + 1.0;
	}

	/////////////////////////////////
	// Now perform the simulation. //
	/////////////////////////////////

	struct Queue ready_list;
	init( &ready_list, 100 );

	double curr_time = 0.0;
	int to_enqueue = 0;

	for ( i = 0; i < N; ++i ) {
		if ( size( &ready_list ) == 0 ) {
			push( &ready_list, &(request[to_enqueue]) );
			++to_enqueue;
		}

		struct data *next_event = (struct data *)front( &ready_list );
		pop( &ready_list );

		if ( next_event->arrival_time > curr_time ) {
			curr_time = next_event->arrival_time;

			while ( ( to_enqueue < N ) && ( request[to_enqueue].arrival_time < curr_time ) ) {
				push( &ready_list, &(request[to_enqueue]) );
				++to_enqueue;
			}
		}

		next_event->start_time = curr_time;
		curr_time += 1/MU;
		next_event->finish_time = curr_time;

		while ( ( to_enqueue < N ) && ( request[to_enqueue].arrival_time < curr_time ) ) {
			push( &ready_list, &(request[to_enqueue]) );
			++to_enqueue;
		}
	}

	deinit( &ready_list );

	//////////////////////////////////////
	// Collect and print the statistics //
	//////////////////////////////////////

	double average_time_in_the_system = 0.0;

	for ( i = 0; i < N; ++i ) {
		average_time_in_the_system += request[i].finish_time - request[i].arrival_time;
	}

	average_time_in_the_system /= N;

	printf( "\nExpected mean time in the system:  %f\n", (2.0 - RHO)/2.0/MU/(1 - RHO) );
	printf(   "Actual average time in the system: %f\n", average_time_in_the_system );

	return 0;
}