###  
# n is the number of simulations. The code generates n simulations, each one has arms tasks. Each simulation will be 'played' play times. There's no relation between n and play as far as i know.

# answering your question: we need n x play to plot, first by taking the 'row average'
# apply(results$greedy, 1, mean) takes the mean of all plays (all columns)

# average reward by play, each simulation will be played 'play' times
# its very confusing
# lets make n = 1, just one simulation

# a machine has 10 slots, play this machine play times and average the results
# thats what i do, but i dont have yet the reward for each arm, so how do i know the greedy (max) value?



get.task = function(arms = 10, u = 0, sdev = 1){
## assuming each action has a normal distribution 
  rnorm(arms, u, sdev)
}

get.rewards = function(arms = 10, n = 2, eps0 = 0.01, eps1 = 0.1, play = 1000){
# n is the number of actions (different options)
# plays is the number of times you play each action (simulatively)

  # x is the rewards matrix, each row is a simulation of a reward for that play (one row is one play).
  #  each column corresponds to an arm (or an action choice)
  rewards.per.play <- t(replicate(play, get.task(arms=arms)))
  # add colnames for visual aid when debugging
  colnames(rewards.per.play) <- paste0("Arm", ifelse((1:arms)<10, 0, ""), 1:arms)

  ## GREEDY 
  ## At each play, we simply select the arm with the highest running average.

  # calculate the running average
  running.avgs <- apply(x, 2, cumsum) / (1:play)

  # for each play, select the running best
  running.best.arm <- apply(running.avgs, 1, which.max) 
  selection.matrix <- cbind(1:play, running.best.arm)
  gredy.rewards <-  rewards.per.play[selection.matrix]

    # NOTE: Have a look at `running.best.arm` to see what's going on 
    { cat ("The running best for the first 21 plays are are:\n\tPlayNo\tBestArm\n\t\t") ;cat(paste(Play.Number=1:21, Best.Arm=running.best.arm[1:21], sep="\t\t"), sep="\n\t\t")  }



  ## EPSILON_1
  ## At each play, we pull a random number, call it x.
  ##    If x < epsi1, we will select an arm at random for that play.
  ##    Otherwise, we will select the arm with the highest running avg.   
  ## Note that the `rewards.per.play` does not change, but the 
  ##      `running.avgs` does change once we have chosen





  # for each arm, we want the _RUNNING AVERAGE_ of all the plays for that arm  
  
  reward.greedy = matrix(0, nrow = play, ncol = n)
  reward.eps0 = matrix(0, nrow = play, ncol = n)
  reward.eps1 = matrix(0, nrow = play, ncol = n)

#  for (i in 1:n){
    play.greedy = numeric(length = play)
    play.eps0 = numeric(length = play)
    play.eps1 = numeric(length = play)

  
    # setup epsi1 and epsi2
    play.eps0 <- play.eps1 <- play.greedy

    for (j in 1:play){
      play.greedy[j] = max(x.av[j, ])
      temp.prob = runif(1)
      ifelse(temp.prob < )
      if (temp.prob < eps0) play.eps0[j] = sample(x[, i], 1) else play.eps0[j] =  play.greedy[j]
      if (temp.prob < eps1) play.eps1[j] = sample(x[, i], 1) else play.eps1[j] =  play.greedy[j]
    }

    reward.greedy[, i] = play.greedy
    reward.eps0[, i] = play.eps0
    reward.eps1[, i] = play.eps1
  }

 return(list(greedy = reward.greedy, eps0 = reward.eps0, eps1 = reward.eps1))
}

run.simulation = function(eps0 = 0.01, eps1 = 0.1){

  results = get.rewards(arms=3, n=1000, play=500, eps0 = eps0, eps1 = eps1)
  greedy = sort(apply(results$greedy, 2, mean))
  eps0 = sort(apply(results$eps0, 2, mean))
  eps1 = sort(apply(results$eps1, 2, mean))

  plot.ts(greedy, col = 'red', lwd = 2, ylim = range(greedy, eps0, eps1), ylab = 'Average reward')
  lines(eps0, col = 'orange', lwd = 2)
  lines(eps1, col = 'navy', lwd = 2)
  grid(col = 'darkgray')

  legend('topleft', c('greedy', 'eps 0.1', 'eps 0.01'), col = c('red', 'orange', 'navy'),
         lwd = 2, cex = 0.8)
}