#!/usr/bin/perl
#
# This code is copyright of Matthew Toia
# It was produced to satisfy a programming project
# and should not be used for anything serious
#

use strict;			# I am forcing some sanity onto the code
use Math::Rand48;		# make sure this is installed

####################
# Helper Subs
####################
sub rand_exp { return -log(1-drand48())/$_[0] }

my $last_step=-1;
$|++;
sub progress {
    my ($curtime, $endtime) = @_;
    my $complete = $curtime / $endtime; 
    my $step = int(50 * $complete);
    return if $step == $last_step;
    $last_step=$step;
    printf("\r[%-50s] ", '=' x $step . '-');
}

sub help {
    print <<HERE;
Usage: $0 [options]
Options:
    -mu1 n           set service rate for first queue to n

    -c1 n            set the number of servers in first queue to n

    -mu2 n           set service rate for second queue to n

    -c2 n            set the number of servers in the second queue to n

    -l, -lambda n    set the input rate to n

    -w, -waitsize n  set the maximum number of requests that can be waiting
                     at anyone queue to n

    -S, -simtime n   sets the length of the simulation to n ticks

    -s, -settle n    sets the number of ticks before accounting begins to n

    -h, -help        displays this message
HERE
}
###################
# Helper objects
###################
package request;
my $next_id = 0;

sub new {
    my ($package, $starttime) = @_;
    my $self = {};
    $self->{id} = $next_id++;
    $self->{starttime} = $starttime;
    $self->{endtime} = $starttime;
    bless $self, $package;
    return $self;
}

# Accessors
sub id {
    my $self = shift;
    return $self->{id} unless @_;
    return $self->{id} = shift;
}
sub starttime {
    my $self = shift;
    return $self->{starttime} unless @_;
    return $self->{starttime} = shift;
}
sub endtime {
    my $self = shift;
    return $self->{endtime} unless @_;
    return $self->{endtime} = shift;
}

################
package pqueue;

sub new {
    my $package = shift;
    my $self = [0]; #store logical length at index 0
    bless $self, $package;
}

sub l_child { return $_[0] * 2 }
sub r_child { return $_[0] * 2 + 1 }
sub parent { return int($_[0] /2) }
sub len { return $_[0]->[0] }
sub peek { return $_[0]->len ? $_[0]->[1] : undef }

sub enqueue {
    my $self = shift;
    $self->[0]++;
    $self->[$self->len] = shift;
    $self->bubble_up($self->len);
}

sub dequeue {
    my $self = shift;
    return undef unless $self->len;

    my $ret = $self->[1];
    $self->swap(1, $self->len);
    $self->[0]--;
    $self->bubble_down(1) if $self->len; #only bubble down if there are elements
    return $ret;
}

sub swap { #for bubbling up and down
    my $temp = $_[0]->[$_[1]];
    $_[0]->[$_[1]] = $_[0]->[$_[2]];
    $_[0]->[$_[2]] = $temp;
}

sub bubble_up {
    my $self = shift;
    my $pos = shift;

    while($pos > 1) {
	if($self->[$pos]->endtime <
	   $self->[parent($pos)]->endtime) {
	    $self->swap($pos, parent($pos));
	    $pos = parent($pos);
	} else {
	    last;
	}
    }
}

sub bubble_down {
    my $self = shift;
    my $pos = 1;
    
    while($pos <= $self->len) {
	if(l_child($pos) <= $self->len) {
	    if(r_child($pos) <= $self->len &&      # right child
	       $self->[r_child($pos)]->endtime <   # has priority
	       $self->[l_child($pos)]->endtime &&  #
	       $self->[r_child($pos)]->endtime <   #
	       $self->[$pos]->endtime              #
	       ) {                                 #
		$self->swap($pos, r_child($pos));
		$pos = r_child($pos);
	    }
	    elsif($self->[l_child($pos)]->endtime < # left child has
		  $self->[$pos]->endtime            # priority
		  ) {                               #
		$self->swap($pos, l_child($pos));
		$pos = l_child($pos);
	    }
	    else {
		last;
	    }
	} else {                                    # current pos
	    last;                                   # has priority
	}
    }
}

################
package MMC;
sub new {
    my $package = shift;
    my $self = {};
    $self->{mu} = shift;          # service constant
    $self->{c} = shift;           # no. of servers
    $self->{waitsize} = shift;    # maxmimum wait length before error
    $self->{wait} = pqueue->new;  # wait queue (priority by starttime)
    $self->{proc} = pqueue->new;  # process queue (by endtime)
    bless $self, $package;
    return $self;
}

# Accessors
sub c {
    my $self = shift;
    return $self->{c} unless @_;
    return $self->{c} = shift;
}
sub mu {
    my $self = shift;
    return $self->{mu} unless @_;
    return $self->{mu} = shift;
}


sub enter {			# Adds to wait queue, and then fills in
    my $self = shift;		# process queue if there is room
    my $r = shift;
    my $curtime = shift;
    die "Wait size exceeded\n" if $self->{wait}->len() > $self->{waitsize};
    $self->{wait}->enqueue($r);
    $self->fill($curtime);
    return 1;
}

sub fill {			# Moves requests from wait queue into
    my $self = shift;		# process queue
    my $curtime = shift;
    while($self->{wait}->len > 0 &&
	  $self->{proc}->len < $self->c) {
	my $r = $self->{wait}->dequeue;
	$r->endtime($curtime + main::rand_exp($self->mu)); # what the endtime will be
	$self->{proc}->enqueue($r);
    }
}

sub exit {			# removes from processing queue
    my $self = shift;
    return $self->{proc}->dequeue;
}

sub next_event {		# returns the time the next event will end
    my $self = shift;
    return $self->{proc}->peek ? 
	$self->{proc}->peek->endtime :
	undef;
}

################
package main;

###############
# Main Program
#################

# The following scalars set up the experiment
my $waitsize;
my $nexttime;
my $settle = 10;
my $simtime = 200;
my $curtime;
my $lambda = 100;
my $mu1 = 20;
my $mu2 = 5;
my $c1 = 500;
my $c2 = 4500;
my $nextq1;
my $nextq2;
my $sum;
my $ssq;
my $count;

# Proccess Command Line Options
while($_ = shift) {
    if(/^-{1,2}waitsize|^-w/) { $waitsize = shift }
    elsif(/^-{1,2}settle|^-s/) { $settle = shift }
    elsif(/^-{1,2}simtime|^-S/) { $simtime = shift }
    elsif(/^-{1,2}lambda|^-l/) { $lambda = shift }
    elsif(/^-{1,2}mu1/) { $mu1 = shift }
    elsif(/^-{1,2}mu2/) { $mu2 = shift }
    elsif(/^-{1,2}c1/) { $c1 = shift }
    elsif(/^-{1,2}c2/) { $c2 = shift }
    elsif(/^-{1,2}help|^-h/) { help; exit 0; }
    else { warn "Unkown option $_\n" }
}

my $q1 = MMC->new($mu1, $c1, $waitsize);
my $q2 = MMC->new($mu2, $c2, $waitsize);

# Report what we're working with
print "mu1=$mu1 c1=$c1\n";
print "mu2=$mu2 c2=$c2\n";
print "lambda=$lambda simtime=$simtime settle=$settle\n";

while($curtime < $simtime) {
    $nextq1 = $q1->next_event || $simtime+1;
    $nextq2 = $q2->next_event || $simtime+1;

    if($nexttime < $nextq1 && $nexttime < $nextq2) { #new request enters
	$curtime = $nexttime;
	my $r = request->new($curtime);
	$q1->enter($r, $curtime);
	$nexttime += rand_exp($lambda);
    } elsif($nextq1 < $nexttime && $nextq1 < $nextq2) { #request dispatched
	$curtime = $nextq1;
	my $r = $q1->exit;
	$q2->enter($r, $curtime);
    } elsif($nextq2 < $nexttime && $nextq2 < $nexttime) { #search completed
	$curtime = $nextq2;
	my $r = $q2->exit;
	if($curtime>$settle) {	# Save accounting information
	    my $diff = $r->endtime - $r->starttime;
	    $sum += $diff;
	    $ssq += $diff*$diff;
	    $count++;
	}
    }    
    progress($curtime, $simtime); # give user some feedback
}

#report results
my $mean = $sum/$count;
my $stdev = sqrt($ssq/$count - $mean * $mean);

printf("\nCount:\t%d\n", $count);
printf("Sum:\t%d\n", $sum);
printf("SSq:\t%d\n", $ssq); 
printf("Mean:\t%.3f\n", $mean);
printf("STDev:\t%.3f\n", $stdev);