Code optimization (first round) #528

yihui · 2025-04-14T18:21:27Z

Close #527

…ata.table; in fact, using data.table will make the code slower)

jdblischak

Do you have some example code for benchmarking that is known to be slow in the app? When I just use gs_design_ahr() |> summary() |> as_gt(), I don't observe much of a difference.

main

git2r::commits(n = 1)[[1]]
## [6b8d5ac] 2025-04-14: Merge pull request #529 from Merck/gt-latex-table-position

library("gsDesign2")
library("microbenchmark")
library("profvis")

system.time(
  gs_design_ahr() |> summary() |> as_gt()
)
## user  system elapsed
## 0.05    0.00    0.08

profvis(gs_design_ahr() |> summary() |> as_gt())

profvis(for (i in 1:10) gs_design_ahr() |> summary() |> as_gt())

microbenchmark(gs_design_ahr() |> summary() |> as_gt())
## Unit: milliseconds
##                             expr     min      lq     mean   median       uq      max neval
##  as_gt(summary(gs_design_ahr())) 56.9447 59.9736 64.96923 62.94655 67.02665 111.5459   100

PR

git2r::commits(n = 1)[[1]]
## [1d3675e] 2025-04-14: Merge branch 'main' into optimize-pw_info

library("gsDesign2")
library("microbenchmark")
library("profvis")

system.time(
  gs_design_ahr() |> summary() |> as_gt()
)
## user  system elapsed
## 0.08    0.00    0.06

profvis(gs_design_ahr() |> summary() |> as_gt())

profvis(for (i in 1:10) gs_design_ahr() |> summary() |> as_gt())

microbenchmark(gs_design_ahr() |> summary() |> as_gt())
## Unit: milliseconds
##                             expr     min      lq     mean   median       uq      max neval
##  as_gt(summary(gs_design_ahr())) 52.4006 58.1329 62.64138 60.54605 63.06605 120.6613   100

yihui · 2025-04-15T16:20:11Z

Yes, to_integer() is slow: #527 (comment)

To make uniroot() faster, I need to make event_diff() faster, which led me to pw_info() as the first target to be optimized. I think event_diff() itself can be optimized (there appears to be extra computation being wasted), and I'll get back to it later.

For now, pw_info() should be faster with this PR, although I still haven't finished it, and there is still room for improvement. The codebase is too complicated to me. I'm still slowly digesting it.

… unnecessary since the cumsum() result is guaranteed to be sorted and unique as long as `duration` is positive should we use `<` instead of `<=` here?

…of filtering `start_time_fr` with `td` again the rewrite here doesn't guarantee equivalence: I actually changed `<=` to `<`, which I feel makes more sense since `start_time_fr` is a simple c() expression, let's just use the expression directly to save a variable name (since naming is always hard)

…mulate the data frames in a list and rbind them all at the end instead

…ll actually be a lot faster a simple benchmark: y = head(penguins) system.time(for (i in 1:1e6) { x = vector('list', 10) for (j in seq_along(x)) x[[j]] = y }) y = head(penguins) system.time(for (i in 1:1e6) { x = list() for (j in seq_along(x)) x[[j]] = y })

…table only for a simple column operation (i.e., adding a `treatment` column)

…l group the computation by `t` next (which will implicitly sort by `t`), and 2) the data is already ordered by `treatment` (via `rbind(control, experimental)` on the previous line, and `control` happens to be "smaller" than `experimental`)

jdblischak · 2025-04-16T15:01:17Z

Yes, to_integer() is slow: #527 (comment)

My apologies. I had seen your comment earlier but then forgot about it when I was performing the benchmarking yesterday.

For now, pw_info() should be faster with this PR

Confirmed! This PR improves the performance, especially for the worst-case scenario (max).

# main
git2r::commits(n = 1)[[1]]
## [6b8d5ac] 2025-04-14: Merge pull request #529 from Merck/gt-latex-table-position
library("gsDesign2")
library("microbenchmark")

benchmark <- function() gs_design_ahr() |> to_integer() |> summary() |> as_gt()
microbenchmark(benchmark())
## Unit: milliseconds
##         expr      min       lq     mean   median       uq      max neval
##  benchmark() 421.8117 447.9604 598.8056 463.4538 482.0702 2983.017   100

# PR
git2r::commits(n = 1)[[1]]
## [1d3675e] 2025-04-14: Merge branch 'main' into optimize-pw_info
library("gsDesign2")
library("microbenchmark")

benchmark <- function() gs_design_ahr() |> to_integer() |> summary() |> as_gt()
microbenchmark(benchmark())
## Unit: milliseconds
##         expr      min       lq     mean   median       uq      max neval
##  benchmark() 347.3618 374.5193 397.7408 386.4939 413.2976 581.9334   100

nanxstats · 2025-04-16T17:18:45Z

Nice! @jdblischak for benchmarking, maybe use a more demanding example like the app's default design. Code can be copied from the report tab, need to run with the latest development version of gsDesign2 on GitHub because of the recent changes on spending function syntax.

… the data by `time` and `stratum`, which we have already done; the columns seem to be already this order, so no need to reorder them

…ogether

…atum, t)`)

….frame instead (although it doesn't really make much difference)

…et by stratum first, then loop on total_duration, otherwise we may be repeatedly subsetting the data for each total_duration

…clearer

LittleBeannie

Hi @yihui, thank you for your code optimization! Rather than referring to it as a review, I would describe it as a learning opportunity for me to pick up all the tricks! I sat beside the desk for a couple of hours to digest everything and learn. Below, please find my minor comments for your consideration:

R/check_arg.R

LittleBeannie · 2025-04-28T20:22:16Z

R/check_arg.R

-  if (theta[n_analysis] < 0) {
-    stop("gs_design_npe() or gs_power_npe(): final effect size must be > 0!")
-  }
+check_theta <- function(theta, K) {


Please lower case the K, as all gsDesign2 source code is in snake case.

I can definitely make that change, but I feel we should occasionally break the snake case rule. In this particular case, I think it makes more sense to match the mathematical notation. In our math notation, we use K to denote the number of analyses, and lower case k as the subscript. What do you think?

BTW, this is an internal function not for end users, so the coding style consistency matters less.

R/check_arg.R

R/expected_accural.R

LittleBeannie · 2025-04-28T20:54:14Z

R/expected_event.R

    ),
    total_duration = 25,
    simple = TRUE) {
-  # Check input values ----


After the re-coding of expected_event, how much running time can we expect to save? I am okay with the changes in the input checking, but I am cautious about the rest. The computation of the expected event is very complex, and I am concerned that a small change may impact areas we have not previously exposed.

If we decide to proceed with these edits, we should also update this vignette (https://merck.github.io/gsDesign2/articles/story-compute-expected-events.html#organizing-calculations-under-the-piecewise-model) to keep everything in sync. I would prefer not to postpone the merge of this PR. Would it be acceptable to leave the edits for expected_event for now and potentially create another PR if needed?

The code has been restored.

R/pw_info.R

LittleBeannie · 2025-04-28T21:59:37Z

R/pwe.R

-    H = cumsum(h * duration), # cumulative hazard
-    survival = exp(-H) # survival
-  )
+  H <- cumulative_rate(x, duration, rate, last_(rate)) # cumulative hazard


Please avoid the upper case variable name H.

That's what we were using previously, although H was a column name instead of a variable name. I think mathematically it makes sense to denote hazard by h and cumulative hazard by H.

LittleBeannie · 2025-04-28T22:07:05Z

R/pwe.R

+  if (survival[1] > 1) stop("`survival` must not be greater than 1")
+  if (last_(survival) >= 1) stop("`survival` must have at least one value < 1")

-  ans <- tibble(Times = times, Survival = survival) %>%


I like the old way in s2pwe, which is more statistically clear. Since s2pwe is commonly a stand-alone function, it should be okay to keep the old programming?

I understand that you may prefer the old tibble() %>% mutate() %>% select() structure, but let me explain the redundancy in the code.

First, the original code:

tibble(Times = times, Survival = survival) %>% mutate( duration = Times - fastlag(Times, first = 0), H = -log(Survival), rate = (H - fastlag(H, first = 0)) / duration ) %>% select(duration, rate)

Capitalizing names to Times and Survival seems to be unnecessary, so we simplify it to:

tibble( duration = times - fastlag(times, first = 0), H = -log(survival), rate = (H - fastlag(H, first = 0)) / duration ) %>% select(duration, rate)

I have created a simple helper function diff_one() (a shorthand of diff(c(0, x))) for the pattern x - fastlag(x, first = 0) (which appeared in this package several times), so we call diff_one() instead:

tibble( duration = diff_one(times), H = -log(survival), rate = diff_one(H) / duration ) %>% select(duration, rate)

At this point, you realize the only reason we need select() is the intermediate column H, so why not just compute H and create the tibble() directly? There is no need to create a column in tibble() and remove it next.

H <- -log(survival) tibble(duration = diff_one(times), rate = diff_one(H) / duration)

As soon as you know what diff_one() means (as I said, this first-order diff operation is quite common throughout the package), the 2 lines of new code should clearer than the 7 lines of old code.

…` [ci skip]

[ci skip]

yihui

I've restored expected_event() as requested, and will create a separate PR (unless we are going to call gsDesign::eEvents() in it).

To me, the only minor issue left is the variable name K. I have renamed it to k but I feel we should break the snake case rule here and use capital K to match the math notation: #528 (comment)

@LittleBeannie If don't have objections, I'll revert the commit 7c18ad1 and this PR will be ready to merge. If you have a strong preference on the lowercase k, the PR can be merged right now.

yihui · 2025-05-02T18:51:47Z

R/expected_event.R

    ),
    total_duration = 25,
    simple = TRUE) {
-  # Check input values ----


The code has been restored.

LittleBeannie

Thank you for explaining your thoughts to me, @yihui! I now have a better understanding of the simplified code. I appreciate the significant reduction in code complexity, and I’m grateful for your efforts in completing this extensive task!

jdblischak · 2025-05-09T18:39:14Z

Since some updates to expected_event() were removed, I re-ran my benchmarking code from #528 (comment). The results are still similar. A ~2x speedup from a cold run and a ~10x speedup from a cached run.

# main
git2r::commits(n = 1)[[1]]
## [a22f40d] 2025-05-08: Merge pull request #543 from Merck/533-the-design-of-the-weight-argument-in-gs_xxx_wlr

system.time(benchmark())
## user  system elapsed
## 1.91    0.06    1.96
microbenchmark(benchmark())
## Unit: seconds
##         expr      min       lq     mean   median     uq      max neval
##  benchmark() 1.819074 1.949909 3.403863 2.061791 2.6421 12.36416   100

# PR
git2r::commits(n = 1)[[1]]
## [380f6c4] 2025-05-07: revert 7c18ad1: k -> K

# cold
clrhash(gsDesign2:::fun_hash)
system.time(benchmark())
## user  system elapsed
## 1.02    0.03    1.17

# cached
microbenchmark(benchmark())
## Unit: milliseconds
##         expr      min       lq     mean   median       uq      max neval
##  benchmark() 167.3944 173.5312 180.8217 176.8842 181.3281 276.1453   100

jdblischak

Amazing work @yihui!

jdblischak · 2025-05-09T18:45:37Z

R/gsDesign2-package.R

+#' @importFrom stats pnorm qnorm setNames uniroot
 #' @importFrom tibble tibble
-#' @importFrom utils tail
+#' @import utils


This would be a good discussion topic for a future group meeting. I prefer to explicitly list each imported function like we currently do for {data.table} and {dplyr}.

For base packages, usually I fully import them. For 3rd party packages, usually I import functions explicitly. I'm okay with selective imports for base packages, too. In fact, the stats package has to be imported selectively to avoid conflicts between stats::filter and dplyr::filter.

yihui added 9 commits April 14, 2025 13:12

don't convert enroll_rate/fail_rate to data.table or add extra classes

24e28d5

use data frames (these objects are fairly small and not worth using d…

ff39138

…ata.table; in fact, using data.table will make the code slower)

1/x is equivalent to x^-1 but faster

8fe60db

Merge branch 'gt-latex-table-position' into optimize-pw_info

9ca055c

can we use the column name "t" directly?

885f2f2

assign row names to NULL to restore to the natural names 1:n

b584d46

factor out the common line in the if-else branches

1e6c5b3

Merge branch 'main' into optimize-pw_info

1d3675e

rewrite the if-else statement

3a6e8f8

jdblischak reviewed Apr 15, 2025

View reviewed changes

yihui added 10 commits April 15, 2025 16:00

call cumsum() only once instead of 3 times; sort()/unique() should be…

d228509

… unnecessary since the cumsum() result is guaranteed to be sorted and unique as long as `duration` is positive should we use `<` instead of `<=` here?

let the arguments time and stratum automatically recycle the scalars

5121301

prepare start_time_fr earlier

22507a5

growing tbl_n in every step of the loop is normally very slow; accu…

ba7d342

…mulate the data frames in a list and rbind them all at the end instead

no need to create the fail object; use fail_c

e417812

it's probably an overkill to convert event_c and event_e to data.…

a1d7da8

…table only for a simple column operation (i.e., adding a `treatment` column)

make sure tbl_n is data.table

3519b8b

yihui added 7 commits April 16, 2025 15:33

the by operation doesn't mutate .SD, so it is equivalent to sorting…

801a05e

… the data by `time` and `stratum`, which we have already done; the columns seem to be already this order, so no need to reorder them

make one fewer copy of ans by applying the row and column filters t…

117c506

…ogether

ln_hr is unused

e2648be

move the recompute outside the inner loop (by = .(t) -> `by = .(str…

6999e11

…atum, t)`)

then move the recompute further outside the outer loop

38d28e8

a single column operation may not worth using data.table, so use data…

e77b4b6

….frame instead (although it doesn't really make much difference)

swap the inner and outer loops to save time on subsetting, i.e., subs…

36a9e82

…et by stratum first, then loop on total_duration, otherwise we may be repeatedly subsetting the data for each total_duration

show argument values in the error messages to make messages a little …

e8abc39

…clearer

yihui force-pushed the optimize-pw_info branch from 512eba6 to e8abc39 Compare April 28, 2025 13:53

LittleBeannie reviewed Apr 28, 2025

View reviewed changes

yihui added 8 commits May 2, 2025 08:45

show different stratum values in the error message

82f3fac

clarify error message in check_theta() [ci skip]

adebd58

use tail(theta, 1) to retrieve the last theta, instead of `theta[K]…

a4cc416

…` [ci skip]

lower case k

7c18ad1

info -> number of analyses

3b8e67e

rates -> rate_group

9a4ac32

k -> i

2defb5e

diff2 -> diff_one

2353c96

[ci skip]

yihui force-pushed the optimize-pw_info branch 2 times, most recently from 35bd79b to 7fffe00 Compare May 2, 2025 18:00

yihui commented May 2, 2025

View reviewed changes

yihui force-pushed the optimize-pw_info branch from b1ec95e to 6a6ad07 Compare May 2, 2025 19:27

yihui added 4 commits May 2, 2025 16:05

temporarily restore expected_event()

7c6e0d9

add news and bump version [ci skip]

b318c16

let check_theta() return the possibly recycled theta values

3a0e19d

similarly, simplify checks for test_upper and test_lower

2c3cd4d

yihui force-pushed the optimize-pw_info branch from 08230d1 to 2c3cd4d Compare May 2, 2025 21:05

yihui added 2 commits May 5, 2025 10:50

Merge branch 'main' into optimize-pw_info

fcaedc9

show event and analysis_time values in the error message

0fdd8e7

LittleBeannie approved these changes May 7, 2025

View reviewed changes

revert 7c18ad1: k -> K

380f6c4

jdblischak approved these changes May 9, 2025

View reviewed changes

LittleBeannie merged commit 3430b56 into main May 9, 2025
15 of 17 checks passed

LittleBeannie deleted the optimize-pw_info branch May 9, 2025 19:00

This was referenced May 12, 2025

Tests failing with dev version of {gsDesign2} Merck/simtrial#328

Closed

Relax test tolerance for summary.simtrial_gs_wlr() Merck/simtrial#329

Merged

Code optimization (first round) #528

Code optimization (first round) #528

Uh oh!

Conversation

yihui commented Apr 14, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jdblischak left a comment

Choose a reason for hiding this comment

Uh oh!

yihui commented Apr 15, 2025

Uh oh!

jdblischak commented Apr 16, 2025

Uh oh!

nanxstats commented Apr 16, 2025

Uh oh!

LittleBeannie left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

yihui May 2, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

yihui left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

LittleBeannie left a comment

Choose a reason for hiding this comment

Uh oh!

jdblischak commented May 9, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jdblischak left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

yihui commented Apr 14, 2025 •

edited

Loading

yihui May 2, 2025 •

edited

Loading

jdblischak commented May 9, 2025 •

edited

Loading