#
#   CTMC.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: 20040831
# 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


CTMC <- function(lambda,mi,n) {
	sistema <- numeric(n)
	for (i in 2:n) {
		tempo.atual <- i * 1/10
                chegada = rpois(1,lambda * 1/10)
		if (sistema[i-1] == 0) {
			if (chegada > 0) {
				sistema[i] = 1
				relogio <- rexp(1,mi)
				despertador <- tempo.atual + relogio
			}
		}
		else {
			if (tempo.atual >= despertador) {
				sistema[i] = 0
			}
			else {
				sistema[i] = 1
			}
		}
	}
	sistema
}

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

estima.Pt <- function(t,sistema,inicial,final) {
        sucessos <- 0
        L <- length(sistema)
        if (t > L/10) { stop('unable to estimate t over larger interval than sampled.') }
        tamanho.janela <- t*10
        N <- L-tamanho.janela+1
        for (i in 1:N) {
                a <- sistema[i]
                b <- sistema[i+(tamanho.janela)-1]
		if (a != inicial) { N <- N -1 }
		else if (b == final) { sucessos <- sucessos + 1 }
        }
        sucessos/N
}

grafico <- function(x,y,inic=0,fim=0,est=0.6,lx=36,ly=0.4) {
	z <- smooth.spline(x,y)
        plot(x,y,xlim=c(0,50),ylim=c(0,1),xlab='t',ylab='f(t)',main=paste('Simulação de Probabilidades de Transição P',inic,fim,'(t)',sep=''))
        curve(predict(z,x)$y,add=T,col="blue")
        curve(Pt.teorica(x,inic,fim),add=T,col="red")
        abline(h=est,lty=2)
	legend(lx,ly,legend=c('simulada',expression('teórica'),bquote(pi * .(fim))),lty=c(1,1,2),col=c('blue','red','black'))
}

Pt.teorica <- function(x,inic,fim) {
 	if (inic == 0 & fim == 0) { mi/(mi+lambda) + lambda/(mi+lambda)*exp(-(lambda+mi)*x) }
	else if (inic == 0 & fim == 1) {1 - ( mi/(mi+lambda) + lambda/(mi+lambda)*exp(-(lambda+mi)*x)) }
        else if (inic == 1 & fim == 1) {( lambda/(mi+lambda) + mi/(mi+lambda)*exp(-(lambda+mi)*x)) }
	else if (inic == 1 & fim == 0) {1-( lambda/(mi+lambda) + mi/(mi+lambda)*exp(-(lambda+mi)*x)) }
	else { stop('not implemented') }
}