#
#   funilaria.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: 20040701
# 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


estacionaria.Teorica <- function(lambda,mi1,mi2,estado) {
	# vamos montar o sistema com as eqs de balanceamento e
	# resove-lo. (pag 428 Ross/notas de aula)
	# aqui eu só peguei o sistema e escrevi em equacao matricial pra resolve-lo, mas é legal detalhar isso no trabalho tambem, mostrar as equacoes na forma normal e dai a forma matricial
	A <- matrix(c(lambda,-lambda,0,0,1,-mi2,0,lambda+mi2,-lambda,1,0,mi1,-mi1,0,1,0,-mi2,0,mi1+mi2,1,0,0,-mi2,0,1),ncol=5)
	colnames(A) <- c("P00","P01","P10","P11","PE1")
	b  <- matrix(c(0,0,0,0,1),ncol=1)
	colnames(b) <- "P estacionaria"
	resultado <- solve(A,b)
	if (missing(estado)) { resultado }
	else { resultado[paste("P",estado,sep=''),]  }
}


funilaria <- function(lambda,mi1,mi2,n,verbose=FALSE,file="") {
        clientes.perdidos <- clientes.atendidos <- 0
	sistema <- rep("00",n)
	desp.mi1 <- desp.mi2 <- -1
	
	for(i in 2:n) {
		tempo.atual <- i * 1/10
		chegada = rpois(1,lambda * 1/10)
		sistema[i] <- sistema[i-1]
		if (tempo.atual >=  desp.mi1 & tempo.atual >=  desp.mi2 & desp.mi1 != -1 & desp.mi2 != -1) {  # tocam mi1 e mi2
			sistema[i] <- "01"
			desp.mi1 <- -1
			desp.mi2 <- tempo.atual+rexp(1,mi2)
		}
		else if (tempo.atual >=  desp.mi2 & desp.mi2 != -1) { # toca só mi2
			if (sistema[i-1] == "E1") { sistema[i] <- "01";  desp.mi2 <- tempo.atual+rexp(1,mi2) }
			else if (sistema[i-1] == "11") { sistema[i] <- "10";  desp.mi2 <- -1}
			else { sistema[i] <- "00";  desp.mi2 <- -1 }
		}
                else if (tempo.atual >=  desp.mi1 & desp.mi1 != -1) { # toca só mi1
			if (sistema[i-1] == "11") { sistema[i] <- "E1" }
			else { sistema[i] <- "01"; desp.mi2 <- tempo.atual+rexp(1,mi2) }
			desp.mi1 <- -1
		}
		if (tempo.atual < desp.mi1) { # nao toca mi1 (mi1 ligado)
			substr(sistema[i],1,1) <- "1"
		}
		if (tempo.atual < desp.mi2) { # nao toca mi2 (mi2 ligado)
                        substr(sistema[i],2,2) <- "1"

		}

		# qq modificacao q ocorreu no sistema por despertadores tocando	
		# acabou. agora resta ver se chegou alguem
		if (chegada > 0 & substr(sistema[i],1,1) == "0") {
			clientes.atendidos <- clientes.atendidos + 1
			substr(sistema[i],1,1) <- "1"
			relogio <- rexp(1,mi1)
			desp.mi1 <- tempo.atual + relogio
		}
		else if (chegada > 0 & substr(sistema[i],1,1) == "1") {
			clientes.perdidos <- clientes.perdidos + chegada
		}
                if (verbose) { cat(tempo.atual,desp.mi1,desp.mi2,sistema[i],'\n',file=file,append=TRUE) }

	}
	list(sistema=sistema,clientes.perdidos=clientes.perdidos,clientes.atendidos=clientes.atendidos,lambda=lambda,mi1=mi1,mi2=mi2,n=n)
}

analise.desempenho <- function(funil.list,custo.amass,custo.pint,preco) {
        # custo por minuto para desamassar: custo.amass
        # custo por minuto para pintar: custo.pint
        # preco por cliente atendido: preco
        tempo.amass <- sum(funil.list$sistema == "10" | funil.list$sistema == "11")/10
        tempo.pint <- sum(funil.list$sistema == "E1" | funil.list$sistema == "11" | funil.list$sistema == "01")/10

        custo.total <- custo.amass * tempo.amass + custo.pint * tempo.pint
        lucro.total <- funil.list$clientes.atendidos*preco - custo.total
        lucro.medio <- (10*lucro.total/funil.list$n)

	x <- list(sistema=funil.list$sistema,mi1=funil.list$mi1,mi2=funil.list$mi2,clientes.perdidos=funil.list$clientes.perdidos,clientes.atendidos=funil.list$clientes.atendidos,lambda=funil.list$lambda,n=funil.list$n,tempo.amass=tempo.amass,tempo.pint=tempo.pint,custo.total=custo.total,lucro.total=lucro.total,lucro.medio=lucro.medio,custo.amass=custo.amass,custo.pint=custo.pint,preco=preco)
	class(x) <- c("analise.desempenho", "list")
	return(x)
}

print.analise.desempenho <- function(x) {
	cat("\n      Analise de Desempenho do Sistema - Funilaria e Pintura\n\n")
	cat("  Modelo teórico \n")
	cat("   Chegadas: PPP(lambda = ",x$lambda,")\n",sep='')
	cat("   Tempo para funilaria: exp(mi1=",x$mi1,")\n",sep='')
        cat("   Tempo para pintura: exp(mi2=",x$mi2,")\n",sep='')
	cat("   Minutos simulados: ",x$n/10,"\n\n")
	cat("  Simulação\n")
	cat("   Tempo médio \t Funilaria \t Pintura\n")
	cat("   esperado \t ",1/x$mi1,"\t\t",1/x$mi2,"\n")
        cat("   observado \t ",x$tempo.amass/x$clientes.atendidos,"\t",x$tempo.pint/x$clientes.atendidos,"\n\n")
	cat("  Balanço financeiro\n")
	cat("   Custo para desamassar: \tR$ ",format(x$custo.amass,nsmall=2),"por min.\n")
        cat("   Custo para pintar: \tR$",format(x$custo.pint,nsmall=2),"por min.\n")
	cat("   Preco por cliente: \tR$",format(x$preco,nsmall=2),"\n")
	cat("   Clientes atendidos:",x$clientes.atendidos," Clientes perdidos:",x$clientes.perdidos,"\n")
	cat("   Custo total: R$",format(x$custo.total,nsmall=2),"Recebimento total: R$ ",format(x$clientes.atendidos*x$preco,nsmall=2),"\n")
	cat("   Lucro total: R$",format(x$lucro.total,nsmall=2),"\n") 
	cat("   Lucro médio por cliente: \tR$",format(x$lucro.total/x$clientes.atendidos),"\n")
	cat("   Lucro médio por minuto: \tR$",format(x$lucro.medio),"\n")
	return(invisible(x))
}