%==========================================================================
%SIMULACI�N DE MONTE CARLO: COLA EN SALA DE ESPERA
%Se cuenta con una guardia medica a la que acuden pacierntes donde:
%  1. Hay N consultorios
%  2. Comienza 8am y se simula durante 24hs
%  3. Entre 8 y 22: tiempo de arribo 2+-1 minuto.
%  4. Entre 22 y 8: tiempo de arribo 10+-5 minutos.
%  5. Cada paciente tarda 5 minutos en entrar y 2 en salir.
%  6. El tiempo de atención es 25+-15 minutos.

%Comenzar con N=5 y estimar:
%  a. Tiempo medio de espera de los pacientes.

%NOTA: La unidad de tiempo es el MINUTO
%==========================================================================
clear; close all; clc;

%INICIALIZACION DE VARIABLES-----------------------------------------------
%Cantidad de Consultorios máxima
L=30;
TE_Total=zeros(2,L);
TCI_Total=zeros(2,L);
for N = 1:L

    %Cantidad máxima de pacientes
    M=1000;
    %Generación de tiempos de llegada. Distribución NORMAL con medias y desvios
    %de acuerdo a lo especificado en 3 y 4. Se obtiene un vector que indica los
    %minutos que pasaron hasta hasta la llegada de cada alumno.

    M1 = 2;
    DS1 = 1;
    M2 = 10;
    DS2 = 5;
    TpoLlegada_aux=zeros(M);
    m=1;
    while sum(TpoLlegada_aux)<(24*60 +1)
        if (sum(TpoLlegada_aux)+(8*60))<(22*60)
            TpoLlegada_aux(m) = normrnd(M1,DS1);
        else
            TpoLlegada_aux(m) = normrnd(M2,DS2);
        end
        m=m+1;
    end
    TpoLlegada_aux(TpoLlegada_aux<0)=0;
    TpoLlegada = TpoLlegada_aux(1:m);
    %m es la cantidad de personas que llegaron en 1 día a la guardia.

    %Tiempo de atención. Se obtiene un vector que indica el tiempo que tardara
    %cada paciente en ser atendido una vez dentro del consultorio.
    M3 = 25;
    DS3 = 15;
    TpoAtencion= normrnd(M3,DS3,[1,m]);
    TpoAtencion(TpoAtencion<0)=0;

    %Almacenamiento de datos
    %Hora de llegada. Se considera hora de inicio de simulacion en t=8hs
    HoraLLegada=zeros(1,m);
    %Hora de Inicio del servicio (llegada + espera)
    HoraIniServ=zeros(1,m);
    %Hora de finalizacion (llegada + espera + tiempo de atencion)
    FinServicio=zeros(1,m);
    %Tiempo de espera (Inicio de servicio - Llegada)
    TpoEspera=zeros(1,m);
    %Tiempo de espera de consultorio (Inicio de servicio - Llegada)
    TpoConsultIdle=zeros(1,m);
    %Cantidad de pacientes en cola ante la llegada de uno nuevo
    EnCola=zeros(1,m);

    %PROCESO DE SIMULACION-----------------------------------------------------
    %Primeros N pacientes: Los consultorios se encuentran vacios, por lo que el
    %tiempo de espera es 0. La hora de llegada coincide con el tiempo de
    %llegada. La hora de entrada al consultorio con la hora de llegada + 5min
    %que le lleva entrar. La hora en que el consultorio queda libre es el
    %tiempo de llegada + 5min + tiempo de atencion +2min que le lleva salir.
    if N==5
        fprintf('Persona &Hora LLegada (hr) &Hora de inicio de servicio (hr) &Fin de servicio (hr) &Tiempo de espera (min) &Cola\\\\ \n')
    end
    HoraLlegada = cumsum(TpoLlegada);
    for i = 1:N
        TpoConsultIdle(i)=HoraLlegada(i);
        HoraIniServ(i)=HoraLlegada(i);
        FinServicio(i)=HoraIniServ(i)+5+TpoAtencion(i)+2;
        if N==5
           fprintf('%d & \t\t%6.2f &  \t\t%6.2f & \t\t\t%6.2f & \t\t%6.2f & \t\t%d \\\\ \n',...
                    i,HoraLlegada(i)/60,HoraIniServ(i)/60,FinServicio(i)/60,TpoEspera(i),EnCola(i));
        end
    end
    HoraLiberacionCons = FinServicio(1:N);

    for k=N+1:m,
       %Hora de inicio del servicio. Si el consultorio esta libre (la hora de
       %llegada es mayor que la hora de fin de alguno de los servicios) se le
       %asigna la hora de llegada al inicio del servicio. En caso contrario,
       %solo podra comenzar luego del fin del servicio anterior.
       [c,i]=min(HoraLiberacionCons);
       if HoraLlegada(k)>=c
           TpoConsultIdle(k)=HoraLlegada(k)-c;
           HoraIniServ(k)=HoraLlegada(k);
       else
           HoraIniServ(k)=c;
           TpoEspera(k)=c-HoraLlegada(k);
       end
       FinServicio(k)=HoraIniServ(k)+5+TpoAtencion(k)+2;
       HoraLiberacionCons(i)=FinServicio(k);

       %Calculo de gente en cola al momento de la llegada
       EnCola(k)=0;
       %Para calcular la gente que esta en la cola hasta el momento se recorre
       %la lista. Si en cada posicion la hora de llegada es inferior a la hora 
       %de fin del servicio anterior, significa que hay alumnos en espera. 
       for j=1:k-1
           if HoraIniServ(j)>HoraLlegada(k)
               EnCola(k)=EnCola(k)+1;
           end
       end
       %Visualizacion de datos
       if N==5
          fprintf('%d & \t\t%6.2f &  \t\t%6.2f & \t\t\t%6.2f & \t\t%6.2f & \t\t%d \\\\ \n',...
                   k,HoraLlegada(k)/60,HoraIniServ(k)/60,FinServicio(k)/60,TpoEspera(k),EnCola(k));
       end
    end

    %EVALUACION DE LA SIMULACION-----------------------------------------------
    %En todos los item se evaluan media y desvio 
    if N==5
        fprintf('------------------------------------------------------------------------------------------------ \n'); 
        fprintf('Tiempo Medio de Llegada & \tTiempo Medio de Espera & \tTiempo Medio de Inactividad & \t\t Cola Maxima \\\\ \n'); 
        fprintf('   %2.1f+/-%2.1f  &   \t\t%2.1f+/-%2.1f  &  \t\t%2.1f+/-%2.1f  & \t\t  %2.1f \\\\ \n',...
                 mean(TpoLlegada),std(TpoLlegada), mean(TpoEspera),std(TpoEspera),mean(TpoConsultIdle),std(TpoConsultIdle),max(EnCola));
    end
    TE_Total(1,N)=mean(TpoEspera);
    TE_Total(2,N)=std(TpoEspera);
    TCI_Total(1,N)=mean(TpoConsultIdle);
    TCI_Total(2,N)=std(TpoConsultIdle);
    if N==5
        figure()
        plot(m,TpoEspera,'g')
        title(['Tiempo de espera para 5 consultorios, con ' num2str(m) ' pacientes'])
        xlabel('Paciente');
        ylabel('Espera (min)');
    end
end
%%
figure()
x=1:L;
x2=2:L;
[AX,H1,H2] = plotyy(x,TE_Total(1,:),x2,diff(TE_Total(1,:)));
set(get(AX(1),'Ylabel'),'String','Tiempo de espera (min)')
set(get(AX(2),'Ylabel'),'String','Pendiente del tiempo de espera')
set(get(AX(1),'Xlabel'),'String','Cantidad de consultorios');

hold(AX(1),'on');
errorbar(TE_Total(1,:),TE_Total(2,:),'-om','LineWidth',2,'MarkerSize',6,...
    'MarkerEdgeColor','red','MarkerFaceColor','red')
ylim([0 12000])
title('Tiempo de espera promedio segun cantidad de consultorios')

figure()
x=1:L;
x2=2:L;
[AX,H1,H2] = plotyy(x,TE_Total(1,:),x,TCI_Total(1,:));
set(get(AX(1),'Ylabel'),'String','Tiempo de espera (min)')
set(get(AX(2),'Ylabel'),'String','Tiempo Libre en Consultorio (min)')
set(get(AX(1),'Xlabel'),'String','Cantidad de consultorios');

hold(AX(2),'on');
errorbar(AX(2),x,TCI_Total(1,:),TCI_Total(2,:),'-ob','LineWidth',2,'MarkerSize',6,...
    'MarkerEdgeColor','cyan','MarkerFaceColor','cyan')
ylim(AX(2),[0 130]);
hold(AX(2),'off');

hold(AX(1),'on');
errorbar(AX(1),x,TE_Total(1,:),TE_Total(2,:),'-om','LineWidth',2,'MarkerSize',6,...
    'MarkerEdgeColor','red','MarkerFaceColor','red')
ylim(AX(1),[0 12000]);
title('Tiempo de espera promedio segun cantidad de consultorios')
hold(AX(1),'off');