# This file was formerly stored as `tests/gsDesign_independent_code.R` and is
# now moved to `tests/testthat/helper.R` to follow testthat best practices.

# This script contains independently programmed functions for validating some of 
# the functions of the gsDesign package.
#-------------------------------------------------------------------------------
# gsPP : averages conditional power across a posterior distribution to compute
#       predictive power.
#-------------------------------------------------------------------------------
# validation code Author : Apurva Bhingare
# x:     design object size.
# i:     look position
# zi:    interim test statistic at ith look
# theta: a vector with theta value(s) at which conditional power is to be computed
# wgts:  Weights to be used with grid points in theta.
# r:     Integer value controlling grid for numerical integration
# total: The default of total=TRUE produces the combined probability for all
#        planned analyses after the interim analysis specified in i
#-------------------------------------------------------------------------------

validate_gsPP <- function(x, i, zi, theta, wgts, r, total = total) {
  k0 <- x$k - i
  I0 <- x$n.I[(i + 1):x$k] - x$n.I[i]

  if (x$test.type == 1) {
    a0 <- rep(-3, k0)
  } else {
    a0 <- (x$lower$bound[(i + 1):x$k] - zi * sqrt(x$n.I[i] / x$n.I[(i + 1):x$k])) / 
          sqrt(I0 / x$n.I[(i + 1):x$k])
  }

  b0 <- (x$upper$bound[(i + 1):x$k] - zi * sqrt(x$n.I[i] / x$n.I[(i + 1):x$k])) /
         sqrt(I0 / x$n.I[(i + 1):x$k])

  cp <- gsProbability(k = k0, theta = theta, n.I = I0, a = a0, b = b0,
                      r = r, overrun = 0)

  gsDen <- dnorm(zi, mean = sqrt(x$n.I[i]) * theta) * wgts

  pp <- cp$upper$prob %*% gsDen / sum(gsDen)

  if (total == TRUE) {
    total.pp <- sum(pp)
    return(total.pp)
  }
  else {
    total.pp <- pp
    return(total.pp)
  }
}
#-------------------------------------------------------------------------------
## gsZ: its computes density for interim test statistic value.
#-------------------------------------------------------------------------------
# x:     design object.size.
# theta: a vector with theta value(s) at which conditional power is to be computed
# i:     look position.
# zi:    interim test statistic at ith look
#-------------------------------------------------------------------------------
validate_gsZ <- function(x, theta, i, zi) {
  nIi <- x$n.I[i]
  mu <- theta * sqrt(nIi)
  fz <- matrix(0, nrow = length(zi), ncol = length(mu))

  for (nc in 1:length(mu)) {
    for (nr in 1:length(zi)) {
      fz[nr, nc] <- dnorm(zi[nr], mu[nc])
    }
  }
  return(fz)
}
#-------------------------------------------------------------------------------
# gsPI: computes Bayesian prediction intervals for future analyses corresponding
#      to results produced by gsPP()
#-------------------------------------------------------------------------------
# x: design object size.
# i: look position
# j: specific analysis for which prediction is being made; must be >i and no more than x$k
# zi: interim test statistic at ith look
#-------------------------------------------------------------------------------
validate_gsPI <- function(x, i, j, k, zi, alpha, mu, sigma1sq) {
  n <- x$n.I
  ti <- x$timing
  postmean <- (mu / sigma1sq + zi[i] * sqrt(n[i])) / (1 / sigma1sq + n[i])
  postvar <- 1 / (1 / sigma1sq + n[i])

  b <- zi[i] * sqrt(ti[i])
  pimean <- b + postmean * (ti[j] - ti[i]) * sqrt(x$n.I[k])
  pivar <- (postvar * (x$n.I[j] - x$n.I[i]) + 1) * (ti[j] - ti[i])
  bpi <- pimean + qnorm(c(alpha / 2, 1 - (alpha / 2))) * sqrt(pivar)
  zpi <- bpi / sqrt(ti[j])
  return(zpi) 
}

#-------------------------------------------------------------------------------
# gsBoundCP: (function description)
#-------------------------------------------------------------------------------
# x: design object size.
#-------------------------------------------------------------------------------
validate_gsBoundCP <- function(x) {
  i0 <- x$k - 1
  theta <- x$delta
  thetahat_hi <- thetahat_hi <- rep(0, i0)
  CP_lo <- CP_hi <- rep(0, i0)

  if (x$test.type == 1) {
    thetahat_hi <- thetahat_low <- if (theta != "thetahat") {
      rep(theta, i0)
    } else {
      x$upper$bound[1:i0] / sqrt(x$n.I[1:i0])
    }
  } else {
    thetahat_hi <- if (theta != "thetahat") {
      rep(theta, i0)
    } else {
      x$upper$bound[1:i0] / sqrt(x$n.I[1:i0])
    }

    thetahat_low <- if (theta != "thetahat") {
      rep(theta, i0)
    } else {
      x$lower$bound[1:i0] / sqrt(x$n.I[1:i0])
    }
  }

  for (i in 1:i0) {
    if (x$test.type > 1) {
      xhi <- gsCP(x, thetahat_hi[i], i, x$upper$bound[i])
      CP_hi[i] <- sum(xhi$upper$prob)

      xlo <- gsCP(x, thetahat_low[i], i, x$lower$bound[i])
      CP_lo[i] <- sum(xhi$lower$prob)

      Bounds <- data.frame(CP_lo, CP_hi)
    } else {
      xhi <- gsCP(x, thetahat_hi[i], i, x$upper$bound[i])
      CP_hi[i] <- sum(xhi$upper$prob)
      Bounds <- CP_hi
    }
  }
  return(Bounds)
}

#-------------------------------------------------------------------------------
# gsCPOS:  gsCPOS() assumes no boundary has been crossed before and including an
#        interim analysis of interest, and computes the probability of success based on this event
#-------------------------------------------------------------------------------
# x:     design object size.
# theta: a vector with theta value(s) at which conditional power is to be computed
# wgts:  Weights to be used with grid points in theta
# i0:    look position
#-------------------------------------------------------------------------------

validate_gsCPOS <- function(x, theta, wgts, i0 = i) {
  one <- rep(0, x$k)
  one[1:i0] <- rep(1, i0)
  zero <- 1 - one

  y <- gsProbability(theta = theta, d = x)

  pAi <- 1 - one %*% (y$upper$prob + y$lower$prob) %*% wgts
  pABi <- zero %*% (y$upper$prob) %*% wgts

  CPOS <- pABi / pAi
  return(CPOS)
}

#-------------------------------------------------------------------------------
# gsPosterior: gsPosterior() computes the posterior density for the group sequential
#             design parameter of interest given a prior density and an interim
#             outcome that is exact or in an interval
#-------------------------------------------------------------------------------
# x:     design object size.
# theta: a vector with theta value(s) at which conditional power is to be computed;
# i:    look position
# zi:   interim test statistic at ith look
#-------------------------------------------------------------------------------

validate_gsPosterior <- function(x, theta, density, gridwgts, wgts, i, zi) {
  if (!is.null(gridwgts)) {
    nIi <- x$n.I[i]
    mu <- theta * sqrt(nIi)
    fz <- matrix(0, nrow = length(zi), ncol = length(mu))


    for (nc in 1:length(mu)) {
      for (nr in 1:length(zi)) {
        fz[nr, nc] <- dnorm(zi[nr], mu[nc])
      }
    }

    p <- (fz * density)
    marg <- sum(fz * density * gridwgts)
    posterior <- p / marg
  } else {
    gridwgts <- rep(1, length(theta))
    nIi <- x$n.I[i] ## ss at the ith look
    mu <- theta * sqrt(nIi) ## theta is standardized delta
    fz <- matrix(0, nrow = length(zi), ncol = length(mu))


    for (nc in 1:length(mu)) {
      for (nr in 1:length(zi)) {
        fz[nr, nc] <- dnorm(zi[nr], mu[nc])
      }
    }

    p <- (fz * density * gridwgts)
    marg <- sum(fz * density * gridwgts)
    posterior <- p / marg
  }

  return(posterior)
}
#-------------------------------------------------------------------------------
# sfPoints: The function sfPoints implements a spending function with values
#           specified for an arbitrary set of specified points.
#-------------------------------------------------------------------------------
# Independent code Author : Apurva Bhingare
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# t :    A vector of time points (information fraction) with increasing
#        values from >0 and <=1.
# param: A real vector of the same length as t specifying the cumulative proportion (between 0 and 1)
#        of spending to corresponding to each point in t
#-------------------------------------------------------------------------------
validate_sfPoints <- function(alpha, t, param) {
  t[t > 1] <- 1
  k <- length(t)
  j <- length(param)

  if (j == k - 1) {
    param <- c(param, 1)
    j <- k
  }

  spend <- alpha * param
  return(spend)
}

#-------------------------------------------------------------------------------
#### sfLinear : The function sfLinear() allows specification of a piecewise
#               linear spending function.
#-------------------------------------------------------------------------------
# Independent code Author : Apurva Bhingare
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# t :    A vector of points with increasing values from 0 to 1, inclusive.
# param: A vector with a positive, even length. Values must range from 0 to 1, inclusive.
#-------------------------------------------------------------------------------
validate_sfLinear <- function(alpha, t, param) {
  j <- length(param)
  k <- j / 2

  s <- t
  s[t <= 0] <- 0
  s[t >= 1] <- 1

  index <- (0 < t) & (t <= param[1])
  s[index] <- param[k + 1] * t[index] / param[1]

  index <- (1 > t) & (t >= param[k])
  s[index] <- param[j] + (t[index] - param[k]) / (1 - param[k]) * (1 - param[j])


  if (k > 1) {
    for (i in 2:k)
    {
      index <- (param[i - 1] < t) & (t <= param[i])
      s[index] <- param[k + i - 1] + (t[index] - param[i - 1]) /
        (param[i] - param[i - 1]) *
        (param[k + i] - param[k + i - 1])
    }
  }
  spend <- alpha * s
  return(spend)
}

#-------------------------------------------------------------------------------
#### sfStep : The function sfStep() specifies a step function spending function
#-------------------------------------------------------------------------------
# Independent code Author : Apurva Bhingare (modified by K Anderson, 5/26/2022)
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# t :    A vector of time points (information fraction) with increasing
#        values from >0 and <=1.
# param: A vector with a positive, even length. Values must range from 0 to 1, inclusive.
#-------------------------------------------------------------------------------
validate_sfStep <- function(alpha, t, param) {
  j <- length(param)

  k <- j / 2

  s <- t

  s[t < param[1]] <- 0
  s[t >= param[k]] <- param[j]
  
  if (k > 1){
    for (i in 1:(k-1)) s[(param[i] <= t) & (t < param[i + 1])] <- param[k + i]
  }
  s[t >= 1] <- 1
  
  spend <- alpha * s
  return(spend)
}

#--------------------------------------------------------------------------------
## sfTDist : The function sfTDist() provides perhaps the maximum flexibility
#           among spending functions provided in the gsDesign package.
#-------------------------------------------------------------------------------
# Independent code Author : Apurva Bhingare
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# t :    A vector of time points (information fraction) with increasing
#        values from >0 and <=1.
# param: A parameter vector of length 3 (5) for specifying t-ditribution, where third (fifth)
#        parameter gives the df
#-------------------------------------------------------------------------------
# Test sfTDist for param of length 3.
validate_sfTDist <- function(alpha, t, param) {
  if (length(param) == 3) {
    t[t > 1] <- 1

    a <- param[1]
    b <- param[2]
    df <- param[3]

    x <- qt(t, df)
    y <- pt(a + b * x, df)
    spend <- alpha * y
  } else {
    if (length(param) == 5) {
      t[t > 1] <- 1
      t0 <- param[1:2]
      p0 <- param[3:4]
      df <- param[5]

      x <- qt(t0, df)
      y <- qt(p0, df)
      b <- (y[2] - y[1]) / (x[2] - x[1])
      a <- y[2] - b * x[2]

      x <- qt(t, df)
      y <- pt(a + b * x, df)
      spend <- alpha * y
    }
    else {
      stop("Check parameter specification")
    }
  }

  return(spend)
}

#-------------------------------------------------------------------------------
## sfTruncated : The functions sfTruncated() and sfTrimmed apply any other
#               spending function over a restricted range. This allows eliminating spending for
#               early interim analyses when you desire not to stop for the bound being specified;
#-------------------------------------------------------------------------------
# Independent code Author : Apurva Bhingare
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# tx :    A vector of time points (information fraction) with increasing
#        values from >0 and <=1.
# param: a list containing the elements sf (a spendfn object such as sfHSD),
#-------------------------------------------------------------------------------
validate_sfTruncated <- function(alpha, tx, param.list) {
  trange <- param.list$trange
  param <- param.list$param
  sf <- param.list$sf

  spend <- as.vector(rep(0, length(tx)))
  spend[tx >= trange[2]] <- alpha
  indx <- trange[1] < tx & tx < trange[2]
  ttmp <- (tx[indx] - trange[1]) / (trange[2] - trange[1])

  if (max(indx)) {
    stmp1 <- sf(alpha = alpha, t = ttmp, param)
    spend[indx] <- stmp1$spend
  }

  spend.truncated <- spend
  return(spend.truncated)
}

#-------------------------------------------------------------------------------
# sfTrimmed : sfTrimmed simply computes the value of the input spending function
#            and parameters in the sub-range of [0,1]
#-------------------------------------------------------------------------------
# Independent code Author : Apurva Bhingare
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# tx :   A vector of time points (information fraction) with increasing
#        values from >0 and <=1.
# param: a list containing the elements sf (a spendfn object such as sfHSD),
#-------------------------------------------------------------------------------
validate_sfTrimmed <- function(alpha, tx, param.list) {
  trange <- param.list$trange
  param <- param.list$param
  sf <- param.list$sf

  spend <- as.vector(rep(0, length(tx)))
  spend[tx >= trange[2]] <- alpha
  indx <- trange[1] < tx & tx < trange[2]

  if (max(indx)) {
    stmp2 <- sf(alpha = alpha, t = tx[indx], param)
    spend[indx] <- stmp2$spend
  }

  spend.trimmed <- spend
  return(spend.trimmed)
}

#-------------------------------------------------------------------------------
## sfGapped : sfGapped() allows elimination of analyses after some time point in the trial
#-------------------------------------------------------------------------------
# Independent code Author : Apurva Bhingare
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# tx :   A vector of time points (information fraction) with increasing
#        values from >0 and <=1.
# param: a list containing the elements sf (a spendfn object such as sfHSD),
#-------------------------------------------------------------------------------
validate_sfGapped <- function(alpha, tx, param.list) {
  trange <- param.list$trange
  param <- param.list$param
  sf <- param.list$sf

  spend <- as.vector(rep(0, length(tx)))
  spend[tx >= trange[2]] <- alpha

  indx <- trange[1] > tx

  if (max(indx)) {
    s <- sf(alpha = alpha, t = tx[indx], param)
    spend[indx] <- s$spend
  }
  indx <- (trange[1] <= tx & trange[2] > tx)
  if (max(indx)) {
    spend[indx] <- sf(alpha = alpha, t = trange[1], param)$spend
  }
  spend.gapped <- spend
  return(spend.gapped)
}

#-------------------------------------------------------------------------------
## spendingFunction: spendingFunction functions in general produce an object of type spendfn.
#-------------------------------------------------------------------------------
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# tx :   A vector of time points (information fraction) with increasing
#        values from >0 and <=1.
# param: A single real value or a vector of real values specifying the spending function parameter(s).
#-------------------------------------------------------------------------------
validate_spendingFunction <- function(alpha, tx, param) {
  spend <- alpha * tx
  return(spend)
}

#-------------------------------------------------------------------------------
############
#sfLogistic
############
# Independent code Author : Apurva Bhingare
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# t :    A vector of time points (information fraction) with increasing
#        values from >0 and <=1.
# param: A vector with a positive, even length. Values must range from 0 to 1, inclusive.
#-------------------------------------------------------------------------------
validate_sfLogistic <- function(alpha, t, param)
{
  if (length(param) == 2 )
  {
    a<-param[1]
    b<-param[2]
    
    xv<-qlogis(t)
    sp<-alpha*plogis(a+b*xv)
    
  }
  else
  { 
    if (length(param) == 4){
      t0<-param[1:2]
      p0<-param[3:4]
      
      xv<-qlogis(t0)
      y<-qlogis(p0)
      
      b<-(y[2]-y[1])/(xv[2]-xv[1])
      a<-y[2]-b*xv[2]
      
      xv<-qlogis(t)
      sp<-alpha*plogis(a+b*xv)
    }
    else
    {
      stop("Check parameter specification")
    }
  }
  return(sp)
}

#-------------------------------------------------------------------------------
###########
#sfCauchy
###########
# Independent code Author : Apurva Bhingare
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# t :    A vector of time points (information fraction) with increasing
#        values from >0 and <=1.
# param: A vector with a positive, even length. Values must range from 0 to 1, inclusive.
#-------------------------------------------------------------------------------
validate_sfCauchy <- function(alpha, t, param) {
  if (length(param) == 2) {
    a <- param[1]
    b <- param[2]
    
    xv <- qcauchy(t)
    sp <- alpha * pcauchy(a + b * xv)
  }
  else {
    if (length(param) == 4) {
      t0 <- param[1:2]
      p0 <- param[3:4]
      
      xv <- qcauchy(t0)
      y <- qcauchy(p0)
      
      b <- (y[2] - y[1]) / (xv[2] - xv[1])
      a <- y[2] - b * xv[2]
      
      xv <- qcauchy(t)
      sp <- alpha * pcauchy(a + b * xv)
    }
    else {
      stop("Check parameter specification")
    }
  }
  return(sp)
}

#-------------------------------------------------------------------------------
###############
#sfExtremeValue
###############
# Independent code Author : Apurva Bhingare
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# t :    A vector of time points (information fraction) with increasing
#        values from >0 and <=1.
# param: A vector with a positive, even length. Values must range from 0 to 1, inclusive.
#-------------------------------------------------------------------------------
validate_sfExtremeValue <- function(alpha, t, param)
{
  if(length(param) == 2){
    a<-param[1]
    b<-param[2]
    
    x<- -log(-log(t))
    sp<-alpha*exp(-exp(-(a+b*x)))
  }
  else {
    if(length(param) == 4){
      
      t0<-param[1:2]
      p0<-param[3:4]
      
      x<- -log(-log(t0))
      y<- -log(-log(p0))
      
      b<-(y[2]-y[1])/(x[2]-x[1])
      a<-y[2]-b*x[2]
      
      x<- -log(-log(t))
      sp<-alpha*exp(-exp(-(a+b*x)))
    }
    else{
      stop("Check parameter specification")
    }
  }
  return(sp)
}

#-------------------------------------------------------------------------------
#################
#sfExtremeValue2
#################
# Independent code Author : Apurva Bhingare
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# t :    A vector of time points (information fraction) with increasing
#        values from >0 and <=1.
# param: A vector with a positive, even length. Values must range from 0 to 1, inclusive.
#-------------------------------------------------------------------------------
validate_sfExtremeValue2 <- function(alpha, t, param)
{
  if(length(param) == 2){
    a<-param[1]
    b<-param[2]
    
    xv<- log(-log(1-t))
    sp<-alpha*(1-exp(-exp((a+b*xv))))
  }
  else {
    if(length(param) == 4){
      t0<-param[1:2]
      p0<-param[3:4]
      
      xv<- log(-log(1-t0))
      y<- log(-log(1-p0))
      
      b<-(y[2]-y[1])/(xv[2]-xv[1])
      a<-y[2]-b*xv[2]
      
      xv<- log(-log(1-t))
      sp<-alpha*(1-exp(-exp((a+b*xv))))
    }
    else{
      stop("Check parameter specification")
    }
  }
  return(sp)
}

#-------------------------------------------------------------------------------
###########
#sfBetaDist
###########
# Independent code Author : Apurva Bhingare
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# t :    A vector of time points (information fraction) with increasing
#        values from >0 and <=1.
# param: A vector with a positive, even length. Values must range from 0 to 1, inclusive.
#-------------------------------------------------------------------------------
validate_sfBetaDist <- function(alpha, t, param){
  t[t > 1] <- 1
  sp <- alpha * stats::pbeta(t, param[1], param[2])
  return(sp)
}

#-------------------------------------------------------------------------------
#################
#sfNormal
#################
# Independent code Author : Apurva Bhingare
# alpha: Type I error (or Type II error) specification takes values between 0 and 1.
# t :    A vector of time points (information fraction) with increasing
#        values from >0 and <=1.
# param: A vector with a positive, even length. Values must range from 0 to 1, inclusive.
#-------------------------------------------------------------------------------
validate_sfNormal <- function(alpha, t, param)
{
  if(length(param) == 2){
    a<-param[1]
    b<-param[2]
    
    xv <- qnorm(1 * (!is.element(t, 1)) * t)
    y <- pnorm(a + b * xv)
    sp <- alpha * (1 * (!is.element(t, 1)) * y + 1 * is.element(t, 1))
  }
  else {
    if(length(param) == 4){
      t0<-param[1:2]
      p0<-param[3:4]
      
      xv<-qnorm(t0)
      y<-qnorm(p0)
      
      b<-(y[2]-y[1])/(xv[2]-xv[1])
      a<-y[2]-b*xv[2]
      
      xv <- qnorm(1 * (!is.element(t, 1)) * t)
      y <- pnorm(a + b * xv)
      sp <- alpha * (1 * (!is.element(t, 1)) * y + 1 * is.element(t, 1))
    }
    else{
      stop("Check parameter specification")
    }
  }
  return(sp)
}

#-------------------------------------------------------------------------------
#varBinomial
#-----------
# ratio: sample size ratio for group 2 divided by group 1
# x: Number of "successes" in the combined control and experimental groups
# n: Number of observations in the combined control and experimental groups
# delta:
# (1) scale="RR", delta0 is the logarithm of the relative risk of event rates (p10 / p20) 
# (2) scale="OR", delta0 is the difference in natural logarithm of the odds-ratio 



validate_varBinom_rr<-function(x, n, delta , ratio, scale = "rr"){
  
  a <- ratio+1
  RR <- exp(delta)
  p <- x/n
  
  if(delta==0){
    
    v <- (((1 - p) / p) * a^2/(a-1))/n
    
  }else{
    
    p1 <- p * a / (ratio * RR + 1)
    p2 <- RR * p1
    
    
    t1 <- a * RR
    t2 <- -(RR * ratio + p2 * ratio + 1 + p1 * RR)
    t3  <- ratio * p2 + p1
    
    p10 <- ((-t2 - sqrt(t2^2 - 4 * t1 * t3)) / 2)/t1
    p20 <- RR * p10
    
    v <- (a * ((1 - p10) / p10 + (1 - p20) / ratio / p20))/n
    
  }
  return(v)
  
}

validate_varBinom_or<-function(x, n, delta , ratio, scale = "OR"){
  
  v <- rep(0, max(length(delta), length(x), length(ratio), length(n)))
  
  p <- x/n
  OR <- exp(delta)
  
  a <- (ratio + 1)
  b<- OR - 1
  c <- -p * a
  d <- 1 + (a-1) * (b+1) + (b) * c
  
  p10 <- (-d + sqrt(d^2 - 4 * b * c)) / 2 / b
  p20 <- (b+1) * p10 / (1 + p10 * b)
  
  v <- (a* (1 / p10 / (1 - p10) + 1 / p20 / (1 - p20) / (a-1)))/n
  v[delta == 0] <- (1 / p / (1 - p) * (1 + 1 / (a-1)) * a)/n
  
  return(v)
}


validate_varBinom_Diff<-function(x, n, delta , ratio){
  
  a <- ratio+1
  p <- x/n
  
  if(delta==0){
    
    x0 <- n*p/a
    x1 <- x0*(a-1)
    
    n0 <- n/(1+ratio)
    n1 <- n-n0
    
    R1<-x1/n1
    R0<-x0/n0
    
    v <- (R1*(1-R1)/n1 + R0*(1-R0)/n0)
  }
  
  return(v)
}
#-------------------------------------------------------------------------------
##eEvents
#--------
## These independently coded functions have been implemented from the
## Lachin and Foulkes (1986) paper


## compute expected value of status (event/censor)
expect_status <- function(lambda, drprate, maxstudy, accdur) {
  rtsum <- lambda + drprate
  t1 <- lambda / rtsum
  t2n <- exp(-rtsum * (maxstudy - accdur)) - exp(-rtsum * maxstudy)
  t2d <- rtsum * accdur
  t2 <- 1 - (t2n / t2d)
  op <- t1 * t2
  return(op)
}

## compute # of events per arm
## Inputs: lambda - the hazard rate on a single arm
##         drprate - dropout rate
##         maxstudy - max. study duration
##         accdur - accrual duration (in the same units as maxstudy)
##         totalSS - total sample size
##         rndprop - randomisation proportion
##
## Output: events - the expected number of events
expect_ev_arm <- function(lambda, drprate, maxstudy, accdur,
                          totalSS, rndprop) {
  expst <- expect_status(lambda, drprate, maxstudy, accdur)
  events <- totalSS * rndprop * expst
  return(events)
}
#------------------------------------------------------------------------
#gsNormalGrid
#-------------
#' Grid points for group sequential design numerical integration
#'
#' Points and weights for Simpson's rule numerical integration from
#' p 349 - 350 of Jennison and Turnbull book.
#' This is not used for arbitrary integration, but for the canonical form of Jennison and Turnbull.
#' mu is computed elsewhere as drift parameter times sqrt of information.
#' Since this is a lower-level routine, no checking of input is done; calling routines should
#' ensure that input is correct.
#' Lower limit of integration can be \code{-Inf} and upper limit of integration can be \code{Inf}
#'
#' @details
#' Jennison and Turnbull (p 350) claim accuracy of \code{10E-6} with \code{r=16}.
#' The numerical integration grid spreads out at the tail to enable accurate tail probability calcuations.
#'
#'
#' @param r Integer, at least 2; default of 18 recommended by Jennison and Turnbull
#' @param mu Mean of normal distribution (scalar) under consideration
#' @param a lower limit of integration (scalar)
#' @param b upper limit of integration (scalar \code{> a})
#'
#' @return A \code{tibble} with grid points in \code{z} and numerical integration weights in \code{w}
#' @export
#'
#' @examples
#' library(dplyr)
#' # approximate variance of standard normal (i.e., 1)
#' gridpts() %>% summarize(var = sum(z^2 * w * dnorm(z)))
#'
#' # approximate probability above .95 quantile (i.e., .05)
#' gridpts(a = qnorm(.95), b = Inf) %>% summarize(p05 = sum(w * dnorm(z)))
gridpts <- function(r = 18, mu = 0, a = -Inf, b = Inf){
  # Define odd numbered grid points for real line
  x <- c(mu - 3 - 4 * log(r / (1:(r - 1))),
         mu - 3 + 3 * (0:(4 * r)) / 2 / r,
         mu + 3 + 4 * log(r / (r - 1):1)
  )
  # Trim points outside of [a, b] and include those points
  if (min(x) < a) x <- c(a, x[x > a])
  if (max(x) > b) x <- c(x[x < b], b)
  
  # Define even numbered grid points between the odd ones
  m <- length(x)
  y <- (x[2:m] + x[1:(m-1)]) / 2
  
  # Compute weights for odd numbered grid points
  i <- 2:(m-1)
  wodd <- c(x[2] - x[1],
            (x[i + 1] - x[i - 1]),
            x[m] - x[m - 1]) / 6
  
  weven <- 4 * (x[2:m] - x[1:(m-1)]) / 6
  
  # Now combine odd- and even-numbered grid points with their
  # corresponding weights
  z <- rep(0, 2*m - 1)
  z[2 * (1:m) - 1] <- x
  z[2 * (1:(m-1))] <- y
  w <- z
  w[2 * (1:m) - 1] <- wodd
  w[2 * (1:(m-1))] <- weven
  
  return(tibble::tibble(z=z, w=w))
}
#-------------------------------------------------------------------------------

## This script contains independently programmed functions for validating some of the functions of the gsDesign package.

#########################################################################################
# This Function validates z2Z.
# Independent Code Author: Apurva
# Independent Code Date: 11/11/2020
#########################################################################################

## x - Base case Design
## z10 - Stage 1 statistics
## n20 - stage incremental sample size

validate_z2z <- function(x, z10, n20) {
  n10 <- (x$n.I[1])
  w1 <- sqrt(n10 / (n10 + n20))
  w2 <- sqrt(n20 / (n10 + n20))
  z20 <- x$upper$bound[2]
  z2incr <- (z20 - w1 * z10) / w2
  return(z2incr)
}





#########################################################################################
# This function validates z2NC.
# Independent Code Author: Apurva
# Independent Code Date: 11/11/2020
#########################################################################################

## x - Base case Design
## z10 - Stage 1 statistics
## info.frac - information fraction


validate_z2NC <- function(x, z10, info.frac) {
  z20 <- x$upper$bound[2]
  z2incr <- (z20 - sqrt(info.frac[1]) * z10) / sqrt(1 - info.frac[1])
  return(z2incr)
}




#########################################################################################
# This function validates z2Fisher.
# Independent Code Author:Apurva
# Independent Code Date:11/11/2020
#########################################################################################
#########################################################################################
## Using the relation
## P(-2 log (p1p(2)) >= ChiSq_a4)=alpha
## where,
## p1= p.value at stage 1
## p(2)= incremental p.value at stage 2
## ChiSq_a4= Upper alpha quantile of chi-square distribution with df=4
#########################################################################################
## x - Base case Design
## z1 - Stage 1 statistics

validate_z2Fisher <- function(x, z1) {
  alpha <- 1 - pnorm(x$upper$bound[2])
  qalpha <- qchisq(alpha, df = 4, lower.tail = FALSE)
  y <- exp(-0.5 * qalpha - log(pnorm(-z1)))
  zzfsr <- qnorm(y, lower.tail = FALSE)
  return(zzfsr)
}




#########################################################################################
# This Function validates comp_bdry
# Reference: Sequential Analysis, Abraham Wald, 1947
# Independent Code Author: Imran
# Independent Code Date: 11/11/2020
#########################################################################################
## alpha - type 1 error
## beta - type 2 error
## p0 - proportion under null hypothesis
## p1 - proportion under alternate hypothesis
## n- Sample size

validate_comp_bdry <- function(alpha, beta, p0, p1, n) {
  A <- (1 - beta) / alpha
  B <- beta / (1 - alpha)

  tmpratio <- (1 - p1) / (1 - p0)
  bndry_den <- log((p1 / p0) * (1 / tmpratio))
  LB <- (log(B) - (n * log(tmpratio))) / bndry_den
  UB <- (log(A) - (n * log(tmpratio))) / bndry_den
  boundaries <- list("LB" = LB, "UB" = UB)

  return(boundaries)
}


#########################################################################################
## Reference: Sequential Analysis, Abraham Wald, 1947
# This Function validates binomialSPRT.
# Independent Code Author: Imran
# Independent Code Date: 11/11/2020
#########################################################################################
## alpha - type 1 error
## beta - type 2 error
## p0 - proportion under null hypothesis
## p1 - proportion under alternate hypothesis
## nmin - minimum sample size for doing ratio test
## nmax - maximum sample size for doing ratio test


Validate_comp_sprt_bnd <- function(alpha, beta, p0, p1, nmin, nmax) {
  lbnd <- c()
  ubnd <- c()

  for (i in nmin:nmax)
  {
    x <- validate_comp_bdry(alpha, beta, p0, p1, i)
    lbnd <- c(lbnd, x$LB)
    ubnd <- c(ubnd, x$UB)
  }

  df <- data.frame(
    n = nmin:nmax, lbnd = lbnd, ubnd = ubnd,
    lbrnd = floor(lbnd), ubrnd = ceiling(ubnd)
  )
  slope <- as.numeric(coef(lm(ubnd ~ n, data = df))["n"]) ## the slope for the boundary lines
  lintercept <- as.numeric(coef(lm(lbnd ~ n, data = df))["(Intercept)"]) ## y-intercept of lower boundary
  uintercept <- as.numeric(coef(lm(ubnd ~ n, data = df))["(Intercept)"]) ## y-intercept of upper boundary

  result <- list("bounds" = df, "slope" = slope, "lowint" = lintercept, "upint" = uintercept)
  return(result)
}

#-------------------------------------------------------------------------------
#save_gg_plot() : Function to save plots created with ggplot2 package and saved
#                 with ggsave() function.
#
# plot	: Plot to save, defaults to last plot displayed.
# path	:  Path of the directory to save plot to: path and filename are 
#        combined to create the fully qualified file name.
# width : the width of the device (in inches).
# height: the height of the device (in inches).
# dpi	  : Plot resolution.

# Plot size in units ("in", "cm", or "mm")
#-------------------------------------------------------------------------------
save_gg_plot <- function(code, width = 4, height = 4) {
  path <- tempfile(fileext = ".png")
  ggplot2::ggsave(path, plot = code, width = width, height = height, dpi = 72, units = "in")
  path
}

#-------------------------------------------------------------------------------
# save_gr_plot(): Function to save plots created with graphics package and saved
#                 with png() function.
# width: the width of the device (in pixels).
# height: the height of the device (in pixels).
#-------------------------------------------------------------------------------
save_gr_plot <- function(code, width = 400, height = 400) {
  path <- tempfile(fileext = ".png")
  png(path, width = width, height = height)
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))
  code
  invisible(dev.off())
  path
}
#-------------------------------------------------------------------------------
#ssrCP : ssrCP() adapts 2-stage group sequential designs to 2-stage sample size
# re-estimation designs based on interim analysis conditional power.

# x  : design object 
# beta : type 2 error
# delta : to be used for true effect size delta
#-------------------------------------------
validate_ssrCP <- function(x, z1){
  
  beta <- x$beta
  delta <-z1/ sqrt(x$n.I[1])
  
  n1 <- x$n.I[1]
  n2 <- x$n.I[2]
  SE <- delta/z1
  
  
  res1 <- (n1 *SE^2)/ delta^2
  res2 <- (x$upper$bound[2] * sqrt(n2) - z1 * sqrt(n1)) / 
    (sqrt(n2 - n1)) - (qnorm(beta))
  
  final_res <- n1 + (res1 * res2 * res2)
  CP <- condPower(z1 = z1, n2 = n1, z2 = z2NC, theta = NULL, x = x)
  expected_out <- list(final_res, CP)
  return(expected_out)
}
#------------------------------------------------------------------------------