# I - E = c
# R - I ~ Ga(gamma,delta)

simSEIR.Non.Markov.gamma <- function(N, beta, gamma, delta, c) {
   
   # initial number of infectives, susceptibles, exposed;
   I <- 1
   S <- N-1
   E <- 0 
   
   # recording time;
   t <- 0;
   times <- c(t);
   
   # create a vector containing the removal times of all the current infectives.
   k <- rgamma(1, gamma, delta)
   r <- k 
   
   # a vector to store the (potential) infection times of all the current exposed individuals.
   I.t <- NULL
   
   # a vector which records the type of event (0 = exposure, 1=infection, 2=removal)
   type <- c(1);
   
   # a counter for labelling the individuals
   lambda <- 1;
   
   # a vector to store the labels
   labels <- c(1);
   
   no.R <- 0
   
   while (E > 0 || I > 0) {
      
      ############################################
      # simulate times to the next possible events
      ############################################
      
      print(c(E, I, S, no.R))
      
      # time to next exposure
      if (S > 0 && I > 0) {
         T.to.E  <- rexp(1, (beta/N)*I*S)
      }
      else {
         T.to.E <- Inf;
      }
      
      # time to the next infection
      if (E > 0) {
         T.to.I <- min(I.t, na.rm=TRUE)
      }
      else {
         T.to.I <- Inf
      }

      # time to next removal
      if (I > 0) {
         T.to.R <- min(r, na.rm=TRUE);
      }
      else {
         T.to.R <- Inf
      }

      # check which of the three events happens first
      T.to.next.event <- min(t + T.to.E, T.to.I, T.to.R)
      type.next.event <- which.min(c(t + T.to.E, T.to.I, T.to.R))
      
      
      if (type.next.event==1) {
         
         # exposure happens first
         E <- E + 1
         S <- S - 1;
         I.t <- append(I.t, T.to.next.event + c) # record the time to infection
         k <- rgamma(1, gamma, delta)
         r <- append(r, T.to.next.event + c + k) # record the time to removal
         
         lambda <- lambda + 1;
         labels <- append(labels, lambda)
         type <- append(type, 0);
         
         times <- append(times, T.to.next.event);
         t <- T.to.next.event
         
      }
      else if (type.next.event == 2) {
         
         # infection happens first
         E <- E - 1
         I <- I + 1
         
         type <- append(type, 1);
         
         index.min.I <- which(min(I.t, na.rm=TRUE) == T.to.I)
         I.t[index.min.I] <- NA
         labels <- append(labels, index.min.I)
         times <- append(times, T.to.I);
         t <- T.to.I
      }
      else{
         
         no.R <- no.R + 1
         #removal occurs
         I <- I-1
         
         type <- append(type, 2);
         index.min.r <- which(min(r, na.rm=TRUE)==r)
         r[index.min.r] <- NA
         
         labels <- append(labels, index.min.r)
         times <- append(times, T.to.R);
         t <- T.to.R
      }
      
   }
   
   res <- list("t"=times, "type"=type, "labels"=labels);
   res
}