Optimize gs_spending_bound()

[nanxstats]

In gsDesign, the C function hupdate was not directly called in R, but indirectly called via a C function gsbound, which does the memory allocation and management effectively.

In gsDesign2, the corresponding part of the logic is currently implemented in R, which could be a bottleneck to slow things down, to the extent that a prospect web interface's responsiveness will be affected.

So, at some point, we will want to rewrite gsbound (and gsbound1) in C++.

[LittleBeannie]

Hi @nanxstats , is there any ideal person we can assign this ticket to? Let us do profiling first before the C++ rewrite.

[nanxstats]

@LittleBeannie No one in particular came to mind, do you have a suggestion?

[LittleBeannie]

I guess you are the only expert in this area in BAAMR 👍

@LittleBeannie No one in particular came to mind, do you have a suggestion?

[elong0527]

Per team discussion, @nanxstats will drive the implementation.

[LittleBeannie]

On hold for now

[jdblischak]

I did some benchmarking and profiling. I don't think the R code in gs_spending_bound() is a major bottleneck, at least in the examples I tried. It would be helpful for me to have some extra examples of code that would typically be executed via the web interface.

First I benchmarked and profiled the example in ?gs_spending_bound, which calls gs_power_ahr(). There are many bottlenecks, but gs_spending_bound() and hupdate() barely register. The biggest cumulative bottlenecks are expected_time(), expected_event(), gs_info_ahr(), pw_info(), and ahr() (note that these are likely nested, so I'd need to dig deeper to find the best place to optimize first)

devtools::load_all()
library("microbenchmark")
library("profvis")

microbenchmark(
  gs_power_ahr = gs_power_ahr(
    analysis_time = c(12, 24, 36),
    event = c(30, 40, 50),
    binding = TRUE,
    upper = gs_spending_bound,
    upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
    lower = gs_spending_bound,
    lpar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL)
  ),
  times = 50
)
## Unit: milliseconds
## expr     min       lq     mean   median       uq      max neval
## gs_power_ahr 863.136 1004.379 1144.783 1104.985 1238.209 1835.487    50

profvis(
  gs_power_ahr(
    analysis_time = c(12, 24, 36),
    event = c(30, 40, 50),
    binding = TRUE,
    upper = gs_spending_bound,
    upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
    lower = gs_spending_bound,
    lpar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL)
  )
)

profvis({
  for (i in 1:10)
  gs_power_ahr(
    analysis_time = c(12, 24, 36),
    event = c(30, 40, 50),
    binding = TRUE,
    upper = gs_spending_bound,
    upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
    lower = gs_spending_bound,
    lpar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL)
  )
})

Second, I focused in on gs_spending_bound() by using the old examples that were removed in #261. This identified only two bottlenecks:

The biggest bottleneck was calling out to gsDesign::sfLDOF()

gsDesign2/R/gs_spending_bound.R

Line 134 in 50067c4

spend <- par$sf(alpha = par$total_spend, t = timing, param = par$param)$spend

The second bottleneck, which was only a fraction of the above, was calling out to hupdate(), which has already been converted to Rcpp

gsDesign2/R/RcppExports.R

Line 13 in 50067c4

.Call(`_gsDesign2_hupdate_rcpp`, r, theta, I, a, b, thetam1, Im1, gm1)

old_examples <- function() {
  info <- (1:3) * 10
  info_frac <- info / max(info)
  k <- length(info_frac)
  
  # 1st analysis
  a1 <- gs_spending_bound(
    k = 1, efficacy = FALSE, theta = 0,
    par = list(sf = gsDesign::sfLDOF, total_spend = 0.025, timing = info_frac, param = NULL),
    hgm1 = NULL
  )
  b1 <- gs_spending_bound(
    k = 1, efficacy = TRUE, theta = 0,
    par = list(sf = gsDesign::sfLDOF, total_spend = 0.025, timing = info_frac, param = NULL),
    hgm1 = NULL
  )
  cat("The (lower, upper) boundary at the 1st analysis is (", a1, ", ", b1, ").\n")
  
  # 2nd analysis
  a2 <- gs_spending_bound(
    k = 2, efficacy = FALSE, theta = 0,
    par = list(sf = gsDesign::sfLDOF, total_spend = 0.025, timing = info_frac, param = NULL),
    hgm1 = h1(r = 18, theta = 0, info = info[1], a = a1, b = b1)
  )
  b2 <- gs_spending_bound(
    k = 2, efficacy = TRUE, theta = 0,
    par = list(sf = gsDesign::sfLDOF, total_spend = 0.025, timing = info_frac, param = NULL),
    hgm1 = h1(r = 18, theta = 0, info = info[1], a = a1, b = b1)
  )
  cat("The upper boundary at the 2nd analysis is (", a2, ", ", b2, ").\n")
}

profvis(for (i in 1:10) old_examples())

In summary, at least for the example call to gs_power_ahr(), converting more of gs_spending_bound() to C++ is unlikely to make a noticeable difference.

[keaven]

I think you’re figuring out where investment may be most beneficial

[jdblischak]

I think you’re figuring out where investment may be most beneficial

After looking at the profiling results more carefully, I plan to start by optimizing ahr() and pw_info()

[nanxstats]

The initial data.table rewrite seems promising to cut the running time substantially, which is a major source of bottlenecks.

As there is not a current goal to have a complete C++ rewrite of these functions (which would be a much bigger effort), let's close this for now as completed as the primary goal is achieved and we have a clear direction to go.

We can open a new one if and when additional optimization is required.

Thanks @jdblischak for the investigations!