accelerate compute_state()

[hoxo-m]

Make compute_state() 2x faster.

actions <- c( 1,  1,  1,  1,  2,  2,  3,  3,  3,  4,  4,   5,  5,   5)
resps   <- c(.2, .1, .0, .3, .2, .4, .1, .6, .8, .5, .8, 1.1, .9, 1.6)

sample_standard_deviation <- RLoptimal:::sample_standard_deviation
compute_state <- RLoptimal:::compute_state

compute_state2 <- function(actions, resps, N_total) {
  resps_per_action <- split(resps, actions)
  count_per_action <- vapply(resps_per_action, length, integer(1L), USE.NAMES = FALSE)
  stopifnot("the number of allocated subjects at each dose should be >= 2" = count_per_action >= 2L)
  
  mean_resps <- vapply(resps_per_action, mean, double(1L), USE.NAMES = FALSE)
  shifted_mean_resps <- mean_resps[-1L] - mean_resps[1L]
  sd_resps <- vapply(resps_per_action, sample_standard_deviation, double(1L), USE.NAMES = FALSE)
  proportion_per_action <- count_per_action / N_total
  
  state <- as.array(c(shifted_mean_resps, sd_resps, proportion_per_action))
  state
}

identical(
  compute_state(actions, resps, 150),
  compute_state2(actions, resps, 150)
)
#> [1] TRUE

library(microbenchmark)

microbenchmark(
  compute_state = compute_state(actions, resps, 150),
  compute_state2 = compute_state2(actions, resps, 150)
)
#> Unit: microseconds
#>           expr   min     lq    mean median     uq   max neval
#>  compute_state 344.8 368.65 412.810 383.15 427.30 831.4   100
#> compute_state2 147.5 160.45 172.628 167.75 178.35 332.1   100