########### Inventory ##############


sS <- function(x){
  littleS <- ceiling(x/40)
  bigS <- littleS + x - (littleS - 1)*40
  list(littleS = littleS, bigS = bigS, Diff = bigS - littleS)
}

MySim <- function(x, n=1, RandomSeed=-1){
  # simulates the (s,S) inventory example of Koenig & Law
  # x in {1,2,...,1600} is the system index
  # n = number of replications
  # RandomSeed sets the initial seed
  # output is average cost for 30 periods
  
  littleS <- ceiling(x/40)
  bigS <- littleS + x - (littleS - 1)*40
  
  if (RandomSeed > 0){set.seed(RandomSeed)}
  
  Y <- rep(0, n)
  for (j in 1:n){
    InvtPos <- bigS
    Cost <- 0
    for (period in 1:30){
      Demand = rpois(1,25)
      if (InvtPos < littleS){
        INext <- bigS
        Cost <- Cost + 32 +3*(bigS - InvtPos)
      }
      else{
        INext <- InvtPos
      }
      if (INext - Demand >= 0){
        Cost <- Cost + INext - Demand
      }
      else{
        Cost <- Cost + 5*(Demand - INext) 
      }
      InvtPos <- INext - Demand
    }
    Y[j] <- Cost/30
  }
  -Y
}


########## SAN ################

MySim <- function(x, n=1, RandomSeed=-1){
  # Simulation of stochastic activity network
  # with given Means for the exponential activity times
  # x in {1,2,3,4,5} is the system index
  # n = number of replications
  # RandomSeed sets the initial seed
  # output is -time to complete network
  
  Means <- matrix(c(0.5,1,1,1,1, 1,0.5,1,1,1, 1,1,0.5,1,1, 0.3,1,1,0.4,1,
                    1,1,1,1,0.5),ncol=5, byrow=F)
  if (RandomSeed > 0){set.seed(RandomSeed)}
  
  Y <- rep(0, n)
  for (j in 1:n){
    A <- rexp(5, rate=1/Means[,x])
    Y[j] <- max(A[1]+A[4], A[1]+A[3]+A[5], A[2]+A[5])
  }
  -Y
}

############## Normal ###################

MySim <- function(x, n=1, RandomSeed=-1){
  # normally distributed data
  # x in {1,2,...,11} is the system index
  # n = number of replications
  # RandomSeed sets the initial seed
  mu <- c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0)
  #mu <- rep(0, 11)
  sigma <- 2
  if (RandomSeed > 0){set.seed(RandomSeed)}
  rnorm(n, mean = mu[x], sd = sigma)
}

############# M/M/1 #####################

MySim <- function(x, n=1, RandomSeed=-1){
  # Simulation of M/M/1 queue with costs for
  # service rate and waiting; arrival rate is 1
  # x in {1,2,...,100} is the system index
  # n = number of replications
  # RandomSeed sets the initial seed
  # output is negative of cost
  mu <- (1:100)/5 + 1
  if (RandomSeed > 0){set.seed(RandomSeed)}
  
  Y <- rep(0, n)
  Wq <- 1/(mu[x]*(mu[x] - 1))
  for (j in 1:n){
    W <- rep(0, 1000)
    D <- Wq
    Snext <- rexp(1, rate=mu[x])
    S <- rexp(1, rate=mu[x])
    for (i in 1:1000){
      D = max(0, D + S - rexp(1, 1))
      W[i] = D + Snext
      S <- Snext
      Snext <- rexp(1, rate=mu[x])
    }
    Y[j] <- mu[x] + 36*mean(W)
  }
  -Y
}