#
#   filas.R - A set of R functions for simulating M/M/k queues
#            
#Copyright (C) 2004 Fernando Henrique F. P. Rosa  <feferraz at ime.usp.br>
#                   Vagner Aparecido Pedro Junior <vapjunior at hotmail.com>
#   
# 
#This program is free software; you can redistribute it and/or
#modify it under the terms of version 2 of the GNU General Public License
#as published by the Free Software Foundation.  A copy of this license should
#be included in the file COPYING.
#
#This program is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with this program; if not, write to the Free Software
#Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
#
#
#
# Última data de edição: 20040630
# Descrição: esse conjunto de funções simula diversos processos de fila
# M/M/k. Para descrição da motivação e algum desenvolvimento teórico veja:
# http://www.feferraz.net/files/lista/mae312-trabalho.pdf


filanat <- function(lambda,mi,nat,n,verbose=FALSE,file="") {
        sistema <- numeric(n)
        despertadores <- rep(-1,nat)
	clientes.chegaram <- 0
	clientes.atendidos <- 0
        for (i in 2:n) {
                tempo.atual <- i * 1/10
                chegada = rpois(1,lambda * 1/10)
    		clientes.chegaram <-  clientes.chegaram + chegada 
                sistema[i] <- sistema[i-1]+chegada
                if (sistema[i-1] < nat & chegada > 0) { # verifica se tinha alguem livre em i-1 e aloca chegadas se sim
                        vagas <- nat - sistema[i-1]
                        if (vagas != sum(despertadores == -1)) { stop('non conforming supositions!') }
                        for (j in 1:min(vagas,chegada)) {
                                relogio <- rexp(1,mi)
                                despertadores[which(despertadores == -1)[1]] <- tempo.atual + relogio
				clientes.atendidos <- clientes.atendidos + 1
                        }

                }
                if (sum(tempo.atual >= despertadores & despertadores != -1) > 0) { # verifica se tem algum atendente que liberou alguem de i-1 pra i, se sim, libera o caixa e aloca se tiver alguem na fila
                        for(j in which(tempo.atual >= despertadores & despertadores != -1)) {
                                sistema[i] <- sistema[i]-1
                                despertadores[j] <- -1
                        }
                        if (sum(despertadores == -1) - max(nat-sistema[i],0)  > 0) { # se essa condicao é verdadeira, há despertadores que precisam ser disparados ainda
                                for(j in 1:(sum(despertadores == -1) - max(nat-sistema[i],0))) {
                                        relogio <- rexp(1,mi)
                                        despertadores[which(despertadores == -1)[1]] <- tempo.atual + relogio
	                                clientes.atendidos <- clientes.atendidos + 1
                                }
                        }
                }
                if (verbose) { cat(tempo.atual,despertadores,sistema[i],'\n',file=file,append=TRUE) }
        } 
        list(sistema=sistema,clientes.atendidos=clientes.atendidos,clientes.chegaram=clientes.chegaram)
} 

estimaL <- function(sistema) {
        sum(sapply(unique(sistema),estima.Pt,t=10000,sistema=sistema,inicial=0) * unique(sistema))
}

# exemplos de uso

exemplos <- function() {
	# cria uma simulacao para 60000 minutos de uma fila com 4 atendentes, lambda=1/10 min^-1 
	# mi = 1/20 min^-1 
	fila.quad <- filanat(1/10,1/20,4,600000)
	# cria uma simulacao para 60000 minutos de uma fila com 1 atendente, lambda=1/10 min^-1
	# mi' = 4mi = 4/20 min^-1
	fila.unic <- filanat(1/10,4/20,1,600000)
	
	fila.unic <- filanat(1/10,8/20,1,600000)
        fila.oct <- filanat(1/10,1/20,8,600000)

	# estima E(Pi) - esperança da distribuicao estacionaria
	estimaL(fila.quad)
	#> 2.270097
	estimaL(fila.unic)
        #> 1.082319
	
	# estima o tamanho medio da fila
	mean(fila.quad[which((fila.quad-4)>=0)]-4)
	#[1] 0.9450508
	#> mean(fila.unic)
	#[1] 1.073482

	# distribuicao estacionaria
	#> sapply(sort(unique(fila.unic)),estima.Pt,t=10000,sistema=fila.unic,inicial=0)  
	# [1] 4.790111e-01 2.461459e-01 1.273909e-01 7.257092e-02 4.075567e-02
	# [6] 1.965919e-02 8.018774e-03 3.369915e-03 1.149019e-03 5.481187e-04
	#[11] 2.517286e-04 5.156376e-04 1.096237e-04 1.096237e-04 3.045104e-04
	#[16] 8.932305e-05
	#> sapply(sort(unique(fila.quad)),estima.Pt,t=10000,sistema=fila.quad,inicial=0) 
	# [1] 0.126052460 0.250845761 0.256064445 0.166239362 0.095438203 0.055964319
	# [7] 0.025016308 0.010316003 0.007115008 0.003777478 0.001683936 0.001077112
	#[13] 0.000409606


	# analise por series temporais
	par(mfrow=c(2,1))
	t <- (1:100000)/10
	plot(t,fila.unic[1:100000],type="l",ylim=range(fila.unic[1:100000],fila.quad[1:100000]),main="fila unica")
        plot(t,fila.quad[1:100000],type="l",ylim=range(fila.unic[1:100000],fila.quad[1:100000]),main="fila quad")

	# gerando uma fila vendo o output:
	filanat(1/2,1/5,2,200,verbose=T)

	# gerando o output e salvando num arquivo:
	filanat(1/2,1/5,2,200,verbose=T,file="relatorio.txt")
	
	#lendo o output no R	
	relat <- read.table('relatorio.txt',col.names=c("tempo","despertador 1","despertador 2","estado do sistema"))

}