set.seed(37);
simul = function(N=100, arrival_rate=.2, processing_rate=1, debug=FALSE) {
    Arrival = cumsum(rexp(n=N, rate=arrival_rate));  
    Service = rexp(n=N, rate =processing_rate); # rep(N,x=1); rgamma(N, shape=.1,rate=.1)
    Remaining = rep(N, x=NA);
    Completion = rep(N, x=NA);
    t = 0;
    Working = 0;
    CurrentTask = NA;
    NextArrival = 1;
    while (TRUE) {
        if(debug) print(t)
        if(debug) print(Arrival);
        if(debug) print(Service);
        if(debug) print(Remaining);
        dtA = NA;
        dtC = NA;
        if(length(Arrival[Arrival>t])>0) {
NA
        }
        if(!is.na(CurrentTask)) {
NA
        }
        if(is.na(dtA) & is.na(dtC)) {
            break;
        } 
        dt = min(dtA,dtC,na.rm=T)
        
NA
        #   la date
        t = t + dt;
NA
        if((NextArrival <=N) & (Arrival[NextArrival] <= t)) { ## je met un <= et pas un == car 3-2.9!=0.1 ...
            Remaining[NextArrival] = Service[NextArrival];
            NextArrival = NextArrival + 1;
        }
        #   le remaining 
        if(!is.na(CurrentTask)) {
            Remaining[CurrentTask] = Remaining[CurrentTask] - dt ;
            if(Remaining[CurrentTask] <= 0) {
                Completion[CurrentTask] = t;
                Remaining[CurrentTask] = NA;
            }
            CurrentTask = NA;
            Working = Working + dt

        }
        #   et currentTask (politique d'ordonnancement: FIFO)
        WaitingList=(1:N)[!is.na(Remaining)];
        if(length(WaitingList)>0) {
            CurrentTask = head(WaitingList,n=1);
        }
    }
    if(debug) print(Completion)
    
    return(data.frame(
        Arrival_rate=arrival_rate,
        Processing_rate=processing_rate,
        Arrival = Arrival,
        Service = Service,
        Completion = Completion,
        Notworking = tail(Completion,1)-Working
        ))
        
}

#Question1.Comparaisondelaperformancede S1 etde S2 NA

NA
#+begin_src R :results output :session *R* :exports both
dfs1 = data.frame()
for (r in  seq(from=.1, to=.9, by=.1)) {
    dfs1 = rbind(dfs1,simul(N=1000, arrival_rate=r, processing_rate=1, debug=FALSE));
}
#+end_src

#+begin_src R :results output :session *R* :exports both
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
dfs12 = dfs1 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(Response = mean(Completion - Arrival), 
              Response_se = sd(Completion - Arrival)/sqrt(n()))
#+end_src

dfs2 = data.frame()
for (r in  seq(from=.1, to=.9, by=.1)) {
    dfs2 = rbind(dfs2,simul(N=1000, arrival_rate=r, processing_rate=2, debug=FALSE));
}
#+end_src

#+begin_src R :results output :session *R* :exports both
dfs22 = dfs2 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(Response = mean(Completion - Arrival), 
              Response_se = sd(Completion - Arrival)/sqrt(n()))
#+end_src

#+begin_src R :results output graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R* 
library(ggplot2);
ggplot(data = dfs12, aes(x=Arrival_rate, y=Response)) + geom_line() +
    geom_errorbar(data = dfs12, aes(ymin = Response - 2*Response_se,
                      ymax = Response + 2*Response_se)) +
    geom_errorbar(data = dfs22, aes(ymin = Response - 2*Response_se,
                      ymax = Response + 2*Response_se))

#+end_src

NA

#+begin_src R :results output :session *R* :exports both
dfs1 = data.frame()
for (r in  seq(from=.1, to=.9, by=.1)) {
    dfs1 = rbind(dfs1,simul(N=1000, arrival_rate = r, processing_rate=1, debug=FALSE));
}




#+begin_src R :results output :session *R* :exports both
dfs12 = dfs1 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(Response = mean(Service+ Notworking * 0.5), 
              Response_se = sd(Service+Notworking * 0.5)/sqrt(n()))
#+end_src

dfs2 = data.frame()
for (r in  seq(from=.1, to=.9, by=.1)) {
    dfs2 = rbind(dfs2,simul(N=1000, arrival_rate=r, processing_rate=2, debug=FALSE));
}
#+end_src

#+begin_src R :results output :session *R* :exports both
dfs22 = dfs2 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(Response = mean(Service*4 + Notworking * 0.25), 
              Response_se = sd(Service*4+ Notworking * 0.25)/sqrt(n()))

#+begin_src R :results output graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R* 
ggplot(dfs12, aes(x=Arrival_rate, y=Response)) + geom_line(data = dfs12) +
    geom_errorbar(data = dfs12, aes(ymin = Response - 2*Response_se,
                      ymax = Response + 2*Response_se)) + geom_line(data = dfs22) +
    geom_errorbar(data = dfs22, aes(ymin = Response - 2*Response_se,
                      ymax = Response + 2*Response_se))

#+end_src

NA

NA

NA NA

NA
simul2 = function(N=100, arrival_rate=0.2) {
    Arrival = cumsum(rexp(n=N, rate=arrival_rate))
    Service1 = rexp(n=N, rate =1)
    Service2 = rexp(n=N, rate =2)
    Remaining = rep(N, x=NA)
    Completion = rep(N, x=NA)
    t = 0
    CurrentTask1 = NA
    CurrentTask2= NA
    NextArrival = 1
    Temp1 = 0
    Temp2 = 0
    while (TRUE) {
        dtA = NA
        dtC = NA
        if(length(Arrival[Arrival>t])>0) {
NA
        }
        if(!is.na(CurrentTask1) || (!is.na(CurrentTask2))) {
NA
        }
        if(is.na(dtA) & is.na(dtC)) {
            break
        }
        dt = min(dtA,dtC,na.rm=T)
       
NA
        #   la date
        t = t + dt
NA
        if((NextArrival <=N) & (Arrival[NextArrival] <= t)) { ## je met un <= et pas un == car 3-2.9!=0.1 ...
NA
            NextArrival = NextArrival + 1
        }
       
        #   le remaining du serveur 1
        if(!is.na(CurrentTask1)) {
            Remaining[CurrentTask1] = Remaining[CurrentTask1] - dt
            if(Remaining[CurrentTask1] <= 0) {
                Completion[CurrentTask1] = t
                Remaining[CurrentTask1] = NA
                CurrentTask1=NA
            }
            Temp1 = Temp1 + dt
        }
        #   le remaining du serveur 2
        if(!is.na(CurrentTask2)) {
            Remaining[CurrentTask2] = Remaining[CurrentTask2] - dt
            if(Remaining[CurrentTask2] <= 0) {
                Completion[CurrentTask2] = t
                Remaining[CurrentTask2] = NA
                CurrentTask2=NA
            }
            Temp2 = Temp2 + dt
        }
       
        #   et currentTask (politique d'ordonnancement: FIFO), sur les deux serveurs
        WaitingList=(1:N)[!is.na(Remaining)]
       
        if(length(WaitingList)>0) {
         
            if (is.na(CurrentTask1)) {
              if(!is.na(CurrentTask2)) {
                CurrentTask1 = WaitingList[2]
              } else {
                CurrentTask1 = head(WaitingList,n=1)
              }
              Remaining[CurrentTask1] = Service1[CurrentTask1]
            }
         
            if (is.na(CurrentTask2)) {
              if(!is.na(CurrentTask1)){
                CurrentTask2 = WaitingList[2]
              } else {
                CurrentTask2 = head(WaitingList,n=1)
              }
              Remaining[CurrentTask2] = Service2[CurrentTask2]
            }
        }
    }
    
    Energy = Temp1*1 + (tail(Completion,1)-Temp1)*.5 + Temp2*4 + (tail(Completion,1)-Temp2)*0.25
    return(data.frame(
        Arrival_rate=arrival_rate,
        Arrival = Arrival,
        Completion = Completion,
        Energy = Energy
        ))
}

dfs3 = data.frame();
for (r in  seq(from=.1, to=.9, by=.1)) {
    dfs3 = rbind(dfs3,simul2(N=1000, arrival_rate=r));
}
#+end_src
#+begin_src R :results output :session *R* :exports both
dfs32 = dfs3 %>% group_by(Arrival_rate) %>%
    summarize(Response = mean(Completion - Arrival), 
              Response_se = sd(Completion - Arrival)/sqrt(n()))
#+end_src
ggplot(data = dfs32, aes(x=Arrival_rate, y=Response)) + geom_line() +
    geom_errorbar(aes(ymin = Response - 2*Response_se,
                      ymax = Response + 2*Response_se))

ggplot(dfs3, aes(x=Arrival_rate, y=Energy))+geom_line();

NA

NA

NA