Three-Stage Designs for Phase 1&2 Clinical Trials

#' @keywords internal
convolve_pmf <- function(lhs, rhs) {
  out <- numeric(length(lhs) + length(rhs) - 1L)
  rhs_len <- length(rhs)
  nz <- which(lhs != 0)
  for(i in nz) {
    out[i:(i + rhs_len - 1L)] <- out[i:(i + rhs_len - 1L)] + lhs[i] * rhs
  }
  out
}

#' @keywords internal
stage_total_pmf_from_lookup <- function(t1, stage1_pmf, stage2_pmf) {
  start <- t1[1L]
  stage2_len <- length(stage2_pmf)
  out <- numeric(t1[length(t1)] - start + stage2_len)
  weights <- stage1_pmf[t1 + 1L]
  for(i in seq_along(t1)) {
    offset <- t1[i] - start + 1L
    out[offset:(offset + stage2_len - 1L)] <- out[offset:(offset + stage2_len - 1L)] + weights[i] * stage2_pmf
  }
  list(start = start, pmf = out, total = sum(weights))
}

#' @keywords internal
stage_total_pmf <- function(t1, n1, prob1, n2, prob2) {
  pmf1 <- numeric(n1 + 1L)
  pmf1[t1 + 1L] <- dbinom(t1, n1, prob1)
  list(start = 0L, pmf = convolve_pmf(pmf1, dbinom(0:n2, n2, prob2)))
}

#' @keywords internal
select_total_pmf <- function(pmf, totals, start = 0L) {
  out <- numeric(max(totals) - min(totals) + 1L)
  idx <- totals - start + 1L
  valid <- idx >= 1L & idx <= length(pmf)
  if(any(valid)) {
    pos <- totals[valid] - min(totals) + 1L
    out[pos] <- out[pos] + pmf[idx[valid]]
  }
  list(start = min(totals), pmf = out, total = sum(out))
}

#' @keywords internal
boundary_cdf_from_pmf <- function(pmf, start, boundaries) {
  cdf <- cumsum(pmf)
  idx <- boundaries - start + 1L
  out <- numeric(length(boundaries))
  valid <- idx >= 1L
  if(any(valid)) {
    out[valid] <- cdf[pmin(idx[valid], length(cdf))]
  }
  out
}

#' @keywords internal
boundary_survival_from_pmf <- function(pmf, start, boundaries, total = sum(pmf)) {
  total - boundary_cdf_from_pmf(pmf, start, boundaries)
}

#' @keywords internal
three.opt <- function(alpha1, alpha2, pt, n, sf.param, pe.par, ...){
  # initialization
  nc <- cumsum(n)
  nt <- nc[3]
  as_left <- HSD(alpha1, nc/nt, sf.param)
  as_right <- HSD(alpha2, nc/nt, sf.param)
  # boundary of r1: [0, n1]
  r1_bdry <- 0:nc[1]
  out_r1 <- unique(pbinom(r1_bdry, n[1], pt[1]))
  ind_r1 <- which(out_r1 <= as_left[1])
  # check if conditions hold
  if(length(ind_r1)==0) stop("No optimal design (left side)")
  r1 <- r1_bdry[max(ind_r1)]
  out_r1 <- out_r1[max(ind_r1)]
  # boundary of s1 (r1, n1-1]
  s1_bdry <- r1:(n[1]-1)
  out_s1 <- unique(1 - pbinom(s1_bdry, n[1], pt[2]))
  ind_s1 <- which(out_s1 <= as_right[1])
  # check if condition holds 
  if(length(ind_s1) == 0) stop("No optimal design (right side)")
  s1 <- s1_bdry[min(ind_s1)]
  out_s1 <- out_s1[min(ind_s1)]
  # boundary of r1: [r1, n1+n2-1]
  r2_bdry <- r1:(nc[2]-1)
  t1 <- (r1+1):s1
  left12 <- stage_total_pmf(t1, n[1], pt[1], n[2], pt[1])
  right12 <- stage_total_pmf(t1, n[1], pt[2], n[2], pt[2])
  out_r2 <- boundary_cdf_from_pmf(left12$pmf, left12$start, r2_bdry) + out_r1
  out_r2 <- unique(out_r2)
  ind_r2 <- which(out_r2 <= as_left[2])
  # check if conditions hold
  if(length(ind_r2) == 0) stop("No optimal design (left side)")
  r2 <- r2_bdry[max(ind_r2)]
  out_r2 <- out_r2[max(ind_r2)]
  # boundary of s2
  s2_bdry <- r2:nc[2]
  out_s2 <- boundary_survival_from_pmf(right12$pmf, right12$start, s2_bdry) + out_s1
  out_s2 <- unique(out_s2)
  ind_s2 <- which(out_s2 <= as_right[2])
  # check if conditions hold
  if(length(ind_s2) == 0) stop("No optimal design (right side)")
  s2 <- s2_bdry[min(ind_s2)]
  out_s2 <- out_s2[min(ind_s2)]
  # boundary of r3: [r2, n1+n2+n3=n]
  r3_bdry <- r2:nc[3]
  left23 <- select_total_pmf(left12$pmf, seq.int(r2 + 1L, s2))
  out_r3 <- boundary_cdf_from_pmf(convolve_pmf(left23$pmf, dbinom(0:n[3], n[3], pt[1])), left23$start, r3_bdry) + out_r2
  out_r3 <-unique(out_r3)
  ind_r3 <- which(out_r3 <= as_left[3])
  # check if conditions hold
  if(length(ind_r3) == 0) stop("No optimal design (left side)")
  r3 <- r3_bdry[max(ind_r3)]
  out_r3 <- out_r3[max(ind_r3)]
  # boundary of s3
  s3_bdry <- r3:(nc[3]-1)
  right23 <- select_total_pmf(right12$pmf, seq.int(r2 + 1L, s2))
  out_s3 <- boundary_survival_from_pmf(convolve_pmf(right23$pmf, dbinom(0:n[3], n[3], pt[2])), right23$start, s3_bdry) + out_s2
  out_s3 <- unique(out_s3)
  ind_s3 <-  which(out_s3 <= as_right[3])
  # check if conditions hold
  if(length(ind_s3) == 0) stop("No optimal design (right side)")
  s3 <- s3_bdry[min(ind_s3)]
  out_s3 <- out_s3[min(ind_s3)]
  # save feasible designs & errors
  bdry <- c(r1, r2, r3, s1, s2, s3)
  err <- c(out_r1, out_r2, out_r3, out_s1, out_s2, out_s3)
  pe <- pt[2] + pe.par
  power12 <- stage_total_pmf(t1, n[1], pe, n[2], pe)
  power23 <- select_total_pmf(power12$pmf, seq.int(r2 + 1L, s2))
  emp_power <- 1 - pbinom(s1, n[1], pe) +
    boundary_survival_from_pmf(power12$pmf, power12$start, s2) +
    boundary_survival_from_pmf(convolve_pmf(power23$pmf, dbinom(0:n[3], n[3], pe)), power23$start, s3)
  # merge results
  names(bdry) <- c("r1", "r2", "r3", "s1", "s2", "s3")
  names(err) <- c("alpha11", "alpha12", "alpha13", "alpha21", "alpha22", "alpha23")
  out <- list(bdry = bdry, error = err, pt = pt, n = n, alpha = c(alpha1, alpha2), beta = 1 - emp_power, sf.param = sf.param)
  class(out) <- "2opt"
  return(out)
}

#' @keywords internal
two.opt <- function(alpha1, alpha2, pt, n, sf.param, pe.par, ...){
  # initialization
  nc <- cumsum(n)
  nt <- nc[2]
  as_left <- HSD(alpha1, nc/nt, sf.param)
  as_right <- HSD(alpha2, nc/nt, sf.param)
  comb <- NULL
  err <- NULL
  # boundary of r1: [0, n1]
  r1_bdry <- 0:nc[1]
  out_r1 <- unique(pbinom(r1_bdry, n[1], pt[1]))
  ind_r1 <- which(out_r1 <= as_left[1])
  # check if conditions hold
  if(length(ind_r1)==0) stop("No optimal design (left side)")
  r1 <- r1_bdry[max(ind_r1)]
  out_r1 <- out_r1[max(ind_r1)]
  # boundary of s1 (r1, n1-1]
  s1_bdry <- r1:(n[1]-1)
  out_s1 <- unique(1 - pbinom(s1_bdry, n[1], pt[2]))
  ind_s1 <- which(out_s1 <= as_right[1])
  # check if condition holds 
  if(length(ind_s1) == 0) stop("No optimal design (right side)")
  s1 <- s1_bdry[min(ind_s1)]
  out_s1 <- out_s1[min(ind_s1)]
  # boundary of r1: [r1, n1+n2-1]
  r2_bdry <- r1:(nc[2]-1)
  t1 <- (r1+1):s1
  left12 <- stage_total_pmf(t1, n[1], pt[1], n[2], pt[1])
  right12 <- stage_total_pmf(t1, n[1], pt[2], n[2], pt[2])
  out_r2 <- boundary_cdf_from_pmf(left12$pmf, left12$start, r2_bdry) + out_r1
  out_r2 <- unique(out_r2)
  ind_r2 <- which(out_r2 <= as_left[2])
  # check if conditions hold
  if(length(ind_r2) == 0) stop("No optimal design (left side)")
  r2 <- r2_bdry[max(ind_r2)]
  out_r2 <- out_r2[max(ind_r2)]
  # boundary of s2
  s2_bdry <- r2:nc[2]
  out_s2 <- boundary_survival_from_pmf(right12$pmf, right12$start, s2_bdry) + out_s1
  out_s2 <- unique(out_s2)
  ind_s2 <- which(out_s2 <= as_right[2])
  # check if conditions hold
  if(length(ind_s2) == 0) stop("No optimal design (right side)")
  s2 <- s2_bdry[min(ind_s2)]
  out_s2 <- out_s2[min(ind_s2)]
  bdry <- c(r1, r2, s1, s2)
  # calculate type-2 error with pt = pt + 0.2
  pe <- pt[2] + pe.par
  err <- c(out_r1, out_r2, out_s1, out_s2)
  power12 <- stage_total_pmf(t1, n[1], pe, n[2], pe)
  emp_power <- 1 - pbinom(s1, n[1], pe) + boundary_survival_from_pmf(power12$pmf, power12$start, s2)
  # merge results
  names(bdry) <- c("r1", "r2", "s1", "s2")
  names(err) <- c("alpha11", "alpha12", "alpha21", "alpha22")
  out <- list(bdry= bdry, error = err, pt = pt, n = n, alpha = c(alpha1, alpha2), beta = 1 - emp_power, sf.param = sf.param)
  class(out) <- "2opt"
  return(out)
}

#' @keywords internal
right.two.opt <- function(alpha, pt, n, sf.param, ...){
  # initialization
  nc <- cumsum(n)
  nt <- nc[2]
  as_right <- HSD(alpha, nc/nt, sf.param)
  # boundary of s1 (0, n1-1]
  s1_bdry <- 0:(n[1]-1)
  out_s1 <- 1-pbinom(s1_bdry, n[1], pt)
  ind_s1 <- which(out_s1 <= as_right[1])
  # check if condition holds 
  if(length(ind_s1)==0) stop("No optimal design (right side)")
  s1 <- s1_bdry[min(ind_s1)]
  out_s1 <- out_s1[min(ind_s1)]
  t1 <- 0:s1
  # boundary of s2
  s2_bdry <- s1:(nc[2]-1)
  right12 <- stage_total_pmf(t1, n[1], pt, n[2], pt)
  out_s2 <- boundary_survival_from_pmf(right12$pmf, right12$start, s2_bdry) + out_s1
  out_s2 <- unique(out_s2)
  ind_s2 <- which(out_s2 <= as_right[2])
  # check if conditions hold
  if(length(ind_s2) == 0) stop("No optimal design (right side)")
  s2 <- s2_bdry[min(ind_s2)]
  out_s2 <- out_s2[min(ind_s2)]
  bdry <- c(s1, s2)
  err <- c(out_s1, out_s2)
  # merge results
  names(bdry) <- c("s1", "s2")
  names(err) <- c("alpha11", "alpha12")
  out <- list(bdry = bdry, error = err, pt = pt, n = n, sf.param = sf.param, alpha = alpha)
  class(out) <- "1opt"
  return(out)
}

#' @keywords internal
right.three.opt <- function(alpha, pt, n, sf.param, ...){
  # initialization
  nc <- cumsum(n)
  nt <- nc[3]
  as_right <- HSD(alpha, nc/nt, sf.param)
  # boundary of s1 (0, n1-1]
  s1_bdry <- 0:(n[1]-1)
  out_s1 <- 1-pbinom(s1_bdry, n[1], pt)
  ind_s1 <- which(out_s1 <= as_right[1])
  # check if condition holds 
  if(length(ind_s1)==0) stop("No optimal design (right side)")
  s1 <- s1_bdry[min(ind_s1)]
  out_s1 <- out_s1[min(ind_s1)]
  # loop over all s1 
  t1 <- 0:s1
  # boundary of s2
  s2_bdry <- s1:(nc[2]-1)
  right12 <- stage_total_pmf(t1, n[1], pt, n[2], pt)
  out_s2 <- boundary_survival_from_pmf(right12$pmf, right12$start, s2_bdry) + out_s1
  out_s2 <- unique(out_s2)
  ind_s2 <- which(out_s2 <= as_right[2])
  # check if conditions hold
  if(length(ind_s2) == 0) stop("No optimal design (right side)")
  s2 <- s2_bdry[min(ind_s2)]
  out_s2 <- out_s2[min(ind_s2)]
  # boundary of s3
  s3_bdry <- s2:(nc[3]-1)
  right23 <- select_total_pmf(right12$pmf, seq.int(0L, s2))
  out_s3 <- boundary_survival_from_pmf(convolve_pmf(right23$pmf, dbinom(0:n[3], n[3], pt)), right23$start, s3_bdry) + out_s2
  out_s3 <- unique(out_s3)
  ind_s3 <-  which(out_s3 <= as_right[3])
  # check if conditions hold
  if(length(ind_s3) == 0) stop("No optimal design (right side)")
  s3 <- s3_bdry[min(ind_s3)]
  out_s3 <- out_s3[min(ind_s3)]
  # save feasible designs & errors
  bdry <- c(s1, s2, s3)
  err <- c(out_s1, out_s2, out_s3)
  # merge results
  names(bdry) <- c("s1", "s2", "s3")
  names(err) <- c("alpha11", "alpha12", "alpha13")
  out <- list(bdry = bdry, error = err, pt = pt, n = n, sf.param = sf.param, alpha = alpha)
  class(out) <- "1opt"
  return(out)
}

#' @keywords internal
zhong.three <- function(alpha1, alpha2, beta, pc, pe, frac_n1 = c(0.1, 0.3), frac_n2 = c(0.2,0.4), sf.param = 1, stop.eff = FALSE, show = TRUE, nmax = 100, ...) {
  if(length(pc) == 1) {
    pc <- rep(pc, 2)
  }
  if(length(alpha1) == 1) {
    alpha1_ini <- c(1, 1, alpha1)
    alpha2_ini <- c(1, 1, alpha2)
  } else {
    stop("alpha1 and alpha2 should be a single value")
  }
  n_count <- 0
  for(nt in 3:nmax){
    results <- list()
    idx <- 0L
    n1_bdry <- floor(nt*frac_n1[1]):ceiling(nt*frac_n1[2])
    n2_bdry <- floor(nt*frac_n2[1]):ceiling(nt*frac_n2[2]) 
    # remove n1=0, n2=0, from boundary 
    n1_bdry <- n1_bdry[n1_bdry!=0]
    n2_bdry <- n2_bdry[n2_bdry!=0]
    n1_n2 <- expand.grid(n1_bdry, n2_bdry)
    # remove n3=0
    sum_n1_n2 <- rowSums(n1_n2)
    id_rm <- sum_n1_n2 < nt 
    n3_bdry <- nt-sum_n1_n2[id_rm]
    n <- as.matrix(cbind(n1_n2[id_rm,], n3_bdry))
    colnames(n) <- c("n1", "n2", "n3")
    for(i in 1:nrow(n)){
      n1 <- n[i, 1]
      n2 <- n[i, 2]
      n3 <- n[i, 3]
      if(is.null(sf.param)) {
        alpha1 <- alpha1_ini
        alpha2 <- alpha2_ini
      }	else {
        alpha1 <- HSD(alpha1_ini[3], cumsum(n[i, ])/nt, sf.param)
        alpha2 <- HSD(alpha2_ini[3], cumsum(n[i, ])/nt, sf.param)
      }
      for(r1 in (n1-1):0) {
        if(stop.eff) s1_lb <- r1 + 1
        else s1_lb <- n1
        # check L1 
        L1 <- pbinom(r1, n1, pc[1])
        if(L1 <= alpha1[1] & L1 <= alpha1[3]){
          for(s1 in s1_lb: n1){
            R1 <- 1 - pbinom(s1, n1, pc[2])
            if(R1 <= alpha2[1] & R1 <= alpha2[3]) {
              t1 <- (r1+1) : s1
              left12 <- stage_total_pmf(t1, n1, pc[1], n2, pc[1])
              right12 <- stage_total_pmf(t1, n1, pc[1], n2, pc[2])
              power12 <- stage_total_pmf(t1, n1, pe, n2, pe)
              pmf3_left <- dbinom(0:n3, n3, pc[1])
              pmf3_right <- dbinom(0:n3, n3, pc[2])
              pmf3_power <- dbinom(0:n3, n3, pe)
              stage1_power <- 1 - pbinom(s1, n1, pe)
              for(r2 in (s1+n2) : r1){
                L2 <- L1 + boundary_cdf_from_pmf(left12$pmf, left12$start, r2)
                if(L2 <= alpha1[2] & L2 <= alpha1[3]) {
                  if(stop.eff) {
                    s2_lb <- max(r2, s1)
                  } 
                    else {
                      s2_lb <- s1+n2
                    } 
                  for(s2 in s2_lb : (s1+n2)) {
                    R2 <- R1 + boundary_survival_from_pmf(right12$pmf, right12$start, s2)
                    if(R2 <= alpha2[2] & R2 <= alpha2[3]) {
                      left23 <- select_total_pmf(left12$pmf, seq.int(r2 + 1L, s2))
                      power23 <- select_total_pmf(power12$pmf, seq.int(r2 + 1L, s2))
                      r3_bdry <- (n3 + s2) : r2
                      l3_vals <- L2 + boundary_cdf_from_pmf(convolve_pmf(left23$pmf, pmf3_left), left23$start, r3_bdry)
                      s3_all <- r2:(s2 + n3)
                      right3_vals <- R2 + boundary_survival_from_pmf(convolve_pmf(left23$pmf, pmf3_right), left23$start, s3_all)
                      power3_vals <- stage1_power +
                        boundary_survival_from_pmf(power12$pmf, power12$start, s2) +
                        boundary_survival_from_pmf(convolve_pmf(power23$pmf, pmf3_power), power23$start, s3_all)
                      for(rr in seq_along(r3_bdry)){
                        r3 <- r3_bdry[rr]
                        L3 <- l3_vals[rr]
                        if(L3 <= alpha1[3]){
                          if(stop.eff) s3_lb <- max(r3, s2)
                          else s3_lb <- r3
                          start_idx <- s3_lb - r2 + 1L
                          for(ss in start_idx:length(s3_all)){
                            s3 <- s3_all[ss]
                            R3 <- right3_vals[ss]
                            if(R3 <= alpha2[3]){
                              Rpe <- power3_vals[ss]
                            } else next
                            if(Rpe >= 1-beta){
                              idx <- idx + 1L
                              results[[idx]] <- c(L1, L2, L3, R1, R2, R3, 1 - Rpe, r1, r2, r3, s1, s2, s3, n1, n2, n3)
                            } else break
                          }
                        } else next
                      } 
                    } else next
                  }
                } else next
              } 
            } else next
          }  
        } else next
      } 
    }
    if(idx > 0L) break
    # if(n_count > n.ratio * n) break
    # if(n_count >= 1) break
    if(show) print(paste("current sample size is", nt))
  }
  out <- do.call(rbind, rev(results))
  colnames(out) <- c("alpha11", "alpha12", "alpha13", "alpha21", "alpha22", "alpha23", "beta", "r1", "r2", "r3", "s1", "s2", "s3", "n1", "n2", "n3")
  # out <- out[order(-out[,3], -out[,6], out[,7]), ]
  out <- out[order(-out[,1], -out[,2], -out[,3], -out[,4], -out[,5], -out[,6], out[,7]), ]
  bdry <- out[1, ][8:13]
  error <- out[1, ][1:7]
  nn <- out[1, ][14:16]
  names(bdry) <- c("r1", "r2", "r3", "s1", "s2", "s3")
  names(nn) <- c("n1", "n2", "n3")
  names(error) <- c("alpha11", "alpha12", "alpha13", "alpha21", "alpha22", "alpha23", "beta")
  return(list(bdry = bdry, error = error, n = nn, complete = out))
}



#' @keywords internal
zhong.two <- function(alpha1, alpha2, beta, pc, pe, stop.eff, sf.param, show, nmax, n.choice, frac_n1,...){
  if(length(pc) == 1) {
    pc <- rep(pc, 2)
  }
  if(length(alpha1) == 1) {
    alpha1_ini <- c(1, alpha1)
    alpha2_ini <- c(1, alpha2)
  } else {
    stop("alpha1 and alpha2 should be a single value")
  }
  
  for(n in 2 : nmax) {
    results <- list()
    idx <- 0L
    n_count <- 0
    n1_bdry <- max(1, floor(n*frac_n1[1])):ceiling(n*frac_n1[2])
    for(n1 in n1_bdry) {
      n2 <- n - n1
      pmf2_left <- dbinom(0:n2, n2, pc[1])
      pmf2_right <- dbinom(0:n2, n2, pc[2])
      pmf2_power <- dbinom(0:n2, n2, pe)
      cdf1_left <- pbinom(0:n1, n1, pc[1])
      tail1_right <- 1 - pbinom(0:n1, n1, pc[2])
      tail1_power <- 1 - pbinom(0:n1, n1, pe)
      db1_left <- dbinom(0:n1, n1, pc[1])
      db1_right <- dbinom(0:n1, n1, pc[2])
      db1_power <- dbinom(0:n1, n1, pe)
      if(!is.null(sf.param)) {
        alpha1 <- HSD(alpha1_ini[2], c(n1, n)/n, sf.param)
        alpha2 <- HSD(alpha2_ini[2], c(n1, n)/n, sf.param)
      } else {
        alpha1 <- alpha1_ini
        alpha2 <- alpha2_ini
      }
      for(r1 in (n1-1):0) {
        if(stop.eff) s1_lb <- r1 + 1
        else  s1_lb <- n1
        # check L1 
        L1 <- cdf1_left[r1 + 1L]
        if(L1 <=  alpha1[1] & L1 <= alpha1[2]){
          for(s1 in s1_lb: n1){
            R1 <- tail1_right[s1 + 1L]
            if(R1 <= alpha2[1] & R1 <= alpha2[2]) {
              t1 <- (r1+1) : s1
              left12 <- stage_total_pmf_from_lookup(t1, db1_left, pmf2_left)
              right12 <- stage_total_pmf_from_lookup(t1, db1_right, pmf2_right)
              power12 <- stage_total_pmf_from_lookup(t1, db1_power, pmf2_power)
              left12_cdf <- cumsum(left12$pmf)
              right12_tail <- right12$total - cumsum(right12$pmf)
              power12_tail <- power12$total - cumsum(power12$pmf)
              stage1_power <- tail1_power[s1 + 1L]
              r2_vals <- (s1 + n2):r1
              r2_idx <- r2_vals - left12$start + 1L
              l2_add <- numeric(length(r2_vals))
              valid_r2 <- r2_idx >= 1L
              l2_add[valid_r2] <- left12_cdf[pmin(r2_idx[valid_r2], length(left12_cdf))]
              L2_vals <- L1 + l2_add
              feasible_r2 <- which(L2_vals <= alpha1[2])
              if(length(feasible_r2) == 0L) {
                next
              }
              boundary_vals <- r1:(s1 + n2)
              boundary_idx_right <- boundary_vals - right12$start + 1L
              boundary_idx_power <- boundary_vals - power12$start + 1L
              r2_add <- numeric(length(boundary_vals))
              rpe_add <- numeric(length(boundary_vals))
              valid_right <- boundary_idx_right >= 1L
              valid_power <- boundary_idx_power >= 1L
              r2_add[valid_right] <- right12_tail[pmin(boundary_idx_right[valid_right], length(right12_tail))]
              rpe_add[valid_power] <- power12_tail[pmin(boundary_idx_power[valid_power], length(power12_tail))]
              R2_all <- R1 + r2_add
              Rpe_all <- stage1_power + rpe_add
              for(rr in feasible_r2){
                r2 <- r2_vals[rr]
                L2 <- L2_vals[rr]
                if(stop.eff) {
                  s2_l <- max(r2, s1)
                } else {
                  s2_l <- r2
                }
                pos <- seq.int(s2_l - r1 + 1L, length(boundary_vals))
                alpha_ok <- which(R2_all[pos] <= alpha2[2])
                if(length(alpha_ok) == 0L) {
                  next
                }
                pos <- pos[alpha_ok[1L]:length(pos)]
                power_ok <- Rpe_all[pos] >= 1 - beta
                keep_len <- match(FALSE, power_ok)
                if(is.na(keep_len)) {
                  keep_pos <- pos
                } else if(keep_len == 1L) {
                  next
                } else {
                  keep_pos <- pos[seq_len(keep_len - 1L)]
                }
                idx <- idx + 1L
                results[[idx]] <- cbind(
                  L1,
                  L2,
                  R1,
                  R2_all[keep_pos],
                  1 - Rpe_all[keep_pos],
                  r1,
                  r2,
                  s1,
                  boundary_vals[keep_pos],
                  n1,
                  n2
                )
              }
            } else next
          }
        } else next
      }
    }
    if(idx > 0L) {
      out <- do.call(rbind, results)
      if(stop.eff){
        en <- apply(out, 1, function(x) x[10] + x[11]*(1-x[3]-x[1]))
      } else {
        en <- apply(out, 1, function(x) x[10] + x[11]*(1-x[1]))
      }
      out <- cbind(out, en)
      out <- round(out, 4)
      out <- out[order(out[,10], -out[, 2], -out[, 4], out[, 5], -out[, 1], -out[, 3], -out[, 8], -out[, 9]), ]
      out <- do.call(rbind, by(out, out[, 10], FUN=function(x) head(x, 1)))
      n_count <- nrow(out)
    }
    if(n_count > n.choice) break
    if(show) print(paste("current sample size is", n))
  }
  out <- as.matrix(out)
  out <- out[order(out[,10], -out[, 2], -out[, 4], -out[, 5]), ]
  colnames(out) <- c("alpha11", "alpha12", "alpha21", "alpha22", "beta", "r1", "r2", "s1", "s2", "n1", "n2", "EN")
  opt <- out[which.min(out[, 12]), ]
  return(list(bdry = opt[6:9], error = opt[1:5], n = opt[10:11], complete = out))
}





#' @keywords internal
HSD <- function (alpha, t, param) {
  t[t > 1] <- 1
  spend <- if (param == 0) t * alpha else alpha * (1 - exp(-t * param))/(1 - exp(-param))
  return(spend)
}