############################################################################

.chktargsint <- function(targs) {

   if (length(targs) > 3L)
      stop("Length of 'targs' argument must be <= 3.", call.=FALSE)

   if (.is.vector(targs)) {
      if (is.null(names(targs))) {
         names(targs) <- c("tau2", "lower", "upper")[seq_along(targs)]
         targs <- as.list(targs)
      } else {
         targs <- list(tau2=unname(targs[startsWith(names(targs), "t")]), lower=unname(targs[startsWith(names(targs), "l")]), upper=unname(targs[startsWith(names(targs), "u")]))
         targs <- targs[sapply(targs, length) > 0L]
      }
   }

   if (any(sapply(targs, length) > 1L))
      stop("Elements of 'targs' arguments must be scalars.", call.=FALSE)

   if (is.null(targs$tau2))
      targs$tau2 <- 0

   return(targs)

}

############################################################################

transf.rtoz <- function(xi) {                 ### resulting value between -Inf (for -1) and +Inf (for +1)
   xi[xi >  1] <-  1
   xi[xi < -1] <- -1
   atanh(xi)
}

transf.ztor <- function(xi)
   tanh(xi)

transf.ztor.int <- function(xi, targs=NULL) {

   if (is.na(xi))
      return(NA_real_)

   targs <- .chktargsint(targs)

   if (is.null(targs$lower))
      targs$lower <- xi-10*sqrt(targs$tau2)
   if (is.null(targs$upper))
      targs$upper <- xi+10*sqrt(targs$tau2)

   toint <- function(zval, xi, tau2)
      tanh(zval) * dnorm(zval, mean=xi, sd=sqrt(tau2))

   cfunc <- function(xi, tau2, lower, upper)
      integrate(toint, lower=lower, upper=upper, xi=xi, tau2=tau2)$value

   if (targs$tau2 == 0) {
      zi <- transf.ztor(xi)
   } else {
      zi <- mapply(xi, FUN=cfunc, tau2=targs$tau2, lower=targs$lower, upper=targs$upper)
   }

   return(c(zi))

}

############################################################################

transf.exp.int <- function(xi, targs=NULL) {

   if (is.na(xi))
      return(NA_real_)

   targs <- .chktargsint(targs)

   if (is.null(targs$lower))
      targs$lower <- xi-10*sqrt(targs$tau2)
   if (is.null(targs$upper))
      targs$upper <- xi+10*sqrt(targs$tau2)

   toint <- function(zval, xi, tau2)
      exp(zval) * dnorm(zval, mean=xi, sd=sqrt(tau2))

   cfunc <- function(xi, tau2, lower, upper)
      integrate(toint, lower=lower, upper=upper, xi=xi, tau2=tau2)$value

   if (targs$tau2 == 0) {
      zi <- exp(xi)
   } else {
      zi <- mapply(xi, FUN=cfunc, tau2=targs$tau2, lower=targs$lower, upper=targs$upper)
   }

   return(c(zi))

}

############################################################################

transf.logit <- function(xi)                       ### resulting value between -Inf (for 0) and +Inf (for +1)
   qlogis(xi)

transf.ilogit <- function(xi)
   plogis(xi)

transf.ilogit.int <- function(xi, targs=NULL) {

   if (is.na(xi))
      return(NA_real_)

   targs <- .chktargsint(targs)

   if (is.null(targs$lower))
      targs$lower <- xi-10*sqrt(targs$tau2)
   if (is.null(targs$upper))
      targs$upper <- xi+10*sqrt(targs$tau2)

   toint <- function(zval, xi, tau2)
      plogis(zval) * dnorm(zval, mean=xi, sd=sqrt(tau2))

   cfunc <- function(xi, tau2, lower, upper)
      integrate(toint, lower=lower, upper=upper, xi=xi, tau2=tau2)$value

   if (targs$tau2 == 0) {
      zi <- transf.ilogit(xi)
   } else {
      zi <- mapply(xi, FUN=cfunc, tau2=targs$tau2, lower=targs$lower, upper=targs$upper)
   }

   return(c(zi))

}

############################################################################

transf.arcsin <- function(xi)                      ### resulting value between 0 (for 0) and asin(1) = pi/2 (for 1)
   asin(sqrt(xi))

transf.iarcsin <- function(xi) {
   zi <- sin(xi)^2
   zi[xi < 0] <- 0                                 ### if xi value is below 0 (e.g., CI bound), return 0
   zi[xi > asin(1)] <- 1                           ### if xi value is above maximum possible value, return 1
   return(c(zi))
}

# transf.iarcsin.int <- function(xi, targs=NULL) {
#
#   if (is.na(xi))
#      return(NA_real_)
#
#   targs <- .chktargsint(targs)
#
#   if (is.null(targs$lower))
#      targs$lower <- 0
#   if (is.null(targs$upper))
#      targs$upper <- asin(1)
#
#   toint <- function(zval, xi, tau2)
#      transf.iarcsin(zval) * dnorm(zval, mean=xi, sd=sqrt(tau2))
#
#   cfunc <- function(xi, tau2, lower, upper)
#      integrate(toint, lower=lower, upper=upper, xi=xi, tau2=tau2)$value
#
#   if (targs$tau2 == 0) {
#      zi <- transf.iarcsin(xi)
#   } else {
#      zi <- mapply(xi, FUN=cfunc, tau2=targs$tau2, lower=targs$lower, upper=targs$upper)
#   }
#
#   return(c(zi))
#
# }

############################################################################

transf.pft <- function(xi, ni) {                   ### Freeman-Tukey transformation for proportions
   xi <- xi*ni
   zi <- 1/2*(asin(sqrt(xi/(ni+1))) + asin(sqrt((xi+1)/(ni+1))))
   return(c(zi))
}

transf.ipft <- function(xi, ni) {                  ### inverse of Freeman-Tukey transformation for individual proportions
   zi <- suppressWarnings(1/2 * (1 - sign(cos(2*xi)) * sqrt(1 - (sin(2*xi)+(sin(2*xi)-1/sin(2*xi))/ni)^2)))
   zi <- ifelse(is.nan(zi), NA_real_, zi)
   zi[xi > transf.pft(1,ni)] <- 1                  ### if xi is above upper limit, return 1
   zi[xi < transf.pft(0,ni)] <- 0                  ### if xi is below lower limit, return 0
   return(c(zi))
}

transf.ipft.hm <- function(xi, targs) {            ### inverse of Freeman-Tukey transformation for a collection of proportions
   if (is.null(targs) || (is.list(targs) && is.null(targs$ni)))
      stop("Must specify the sample sizes via the 'targs' argument.", call.=FALSE)
   if (is.list(targs)) {
      ni <- targs$ni
   } else {
      ni <- ni
   }
   nhm <- 1/(mean(1/ni, na.rm=TRUE))               ### calculate harmonic mean of the ni's
   zi <- suppressWarnings(1/2 * (1 - sign(cos(2*xi)) * sqrt(1 - (sin(2*xi)+(sin(2*xi)-1/sin(2*xi))/nhm)^2)))
   zi <- ifelse(is.nan(zi), NA_real_, zi)          ### it may not be possible to calculate zi
   zi[xi > transf.pft(1,nhm)] <- 1                 ### if xi is above upper limit, return 1
   zi[xi < transf.pft(0,nhm)] <- 0                 ### if xi is below lower limit, return 0
   return(c(zi))
}

############################################################################

transf.isqrt <- function(xi) {
   zi <- xi*xi
   zi[xi < 0] <- 0                                 ### if xi value is below 0 (e.g., CI bound), return 0
   return(c(zi))
}

############################################################################

transf.irft <- function(xi, ti) {                  ### Freeman-Tukey transformation for incidence rates
   zi <- 1/2*(sqrt(xi) + sqrt(xi + 1/ti))          ### xi is the incidence rate (not the number of events!)
   return(c(zi))
}

transf.iirft <- function(xi, ti) {                 ### inverse of Freeman-Tukey transformation for incidence rates (see Freeman-Tukey_incidence.r in code directory)
   #zi <- (1/ti - 2*xi^2 + ti*xi^4)/(4*xi^2*ti)    ### old version where transf.irft was not multiplied by 1/2
   zi <- (1/ti - 8*xi^2 + 16*ti*xi^4)/(16*xi^2*ti) ### xi is the incidence rate (not the number of events!)
   zi <- ifelse(is.nan(zi), NA_real_, zi)
   zi[xi < transf.irft(0,ti)] <- 0                 ### if xi is below lower limit, return 0
   zi[zi <= .Machine$double.eps] <- 0              ### avoid finite precision errors in back-transformed values (transf.iirft(transf.irft(0, 1:200), 1:200))
   return(c(zi))
}

############################################################################

transf.ahw <- function(xi) {                       ### resulting value between 0 (for alpha=0) and 1 (for alpha=1)
   #zi <- (1-xi)^(1/3)
   zi <- 1 - (1-xi)^(1/3)
   return(c(zi))
}

transf.iahw <- function(xi) {
   #zi <- 1-xi^3
   zi <- 1 - (1-xi)^3
   zi <- ifelse(is.nan(zi), NA_real_, zi)
   zi[xi > 1] <- 1                                 ### if xi is above upper limit, return 1
   zi[xi < 0] <- 0                                 ### if xi is below lower limit, return 0
   return(c(zi))
}

transf.abt <- function(xi) {                       ### Bonett (2002) transformation of alphas (without bias correction)
#transf.abt <- function(xi, ni) {                  ### resulting value between 0 (for alpha=0) to Inf (for alpha=1)
   #zi <- log(1-xi) - log(ni/(ni-1))
   #zi <- log(1-xi)
   zi <- -log(1-xi)
   return(c(zi))
}

transf.iabt <- function(xi) {                      ### inverse of Bonett (2002) transformation
#transf.iabt <- function(xi, ni) {
   #zi <- 1 - exp(xi) * ni / (ni-1)
   #zi <- 1 - exp(xi)
   zi <- 1 - exp(-xi)
   zi <- ifelse(is.nan(zi), NA_real_, zi)
   zi[xi < 0] <- 0                                 ### if xi is below lower limit, return 0
   return(c(zi))
}

############################################################################

transf.dtou1 <- function(xi) {
   u2i <- pnorm(abs(xi)/2)
   return((2*u2i - 1) / u2i)
}

transf.dtou2 <- function(xi)
   pnorm(xi/2)

transf.dtou3 <- function(xi)
   pnorm(xi)

transf.dtocles <- function(xi)
   pnorm(xi/sqrt(2))

transf.dtorpb <- function(xi, n1i, n2i) {
   if (missing(n1i) || missing(n2i)) {
      hi <- 4
   } else {
      if (length(n1i) != length(n2i))
         stop("Length of 'n1i' does not match length of 'n2i'.", call.=FALSE)
      if (length(n1i) != length(xi))
         stop("Length of 'n1i' and 'n2i' does not match length of 'xi'.", call.=FALSE)
      mi <- n1i + n2i - 2
      hi <- mi / n1i + mi / n2i
   }
   return(xi / sqrt(xi^2 + hi))
}

transf.dtobesd <- function(xi) {
   rpbi <- xi / sqrt(xi^2 + 4)
   return(0.50 + rpbi/2)
}

transf.dtomd <- function(xi, targs=NULL) {
   if (is.null(targs) || (is.list(targs) && is.null(targs$sd)))
      stop("Must specify a standard deviation value via the 'targs' argument.", call.=FALSE)
   if (is.list(targs)) {
      sd <- targs$sd
   } else {
      sd <- targs
   }
   if (length(sd) != 1L)
      stop("Specify a single standard deviation value via the 'targs' argument.", call.=FALSE)
   return(xi * sd)
}

transf.rtorpb <- function(xi, pi) {
   if (missing(pi)) {
      pi <- 0.5
   } else {
      if (length(pi) == 1L)
         pi <- rep(pi, length(xi))
      if (length(xi) != length(pi))
         stop("Length of 'xi' does not match length of 'pi'.", call.=FALSE)
   }
   if (any(pi < 0 | pi > 1, na.rm=TRUE))
      stop("One or more 'pi' values are < 0 or > 1.", call.=FALSE)
   return(xi * dnorm(qnorm(pi)) / sqrt(pi*(1-pi)))
}

transf.rtod <- function(xi, n1i, n2i) {
   if (missing(n1i) || missing(n2i)) {
      hi <- 4
      pi <- 0.5
      n1i <- 1
      n2i <- 1
   } else {
      if (length(n1i) != length(n2i))
         stop("Length of 'n1i' does not match length of 'n2i'.", call.=FALSE)
      if (length(n1i) != length(xi))
         stop("Length of 'n1i' and 'n2i' does not match length of 'xi'.", call.=FALSE)
      mi <- n1i + n2i - 2
      hi <- mi / n1i + mi / n2i
      pi <- n1i / (n1i + n2i)
   }
   if (any(c(n1i < 0, n2i < 0), na.rm=TRUE))
      stop("One or more values specified via the 'n1i' or 'n2i' arguments are negative.")
   rpbi <- xi * dnorm(qnorm(pi)) / sqrt(pi*(1-pi))
   return(sqrt(hi) * rpbi / sqrt(1 - rpbi^2))
}

transf.lnortord <- function(xi, pc) {
   if (length(pc) == 1L)
      pc <- rep(pc, length(xi))
   if (length(xi) != length(pc))
      stop("Length of 'xi' does not match length of 'pc'.", call.=FALSE)
   if (any(pc < 0) || any(pc > 1))
      stop("The control group risk 'pc' must be between 0 and 1.", call.=FALSE)
   return(exp(xi)*pc / (1 - pc + pc * exp(xi)) - pc)
}

transf.lnortorr <- function(xi, pc) {
   if (length(pc) == 1L)
      pc <- rep(pc, length(xi))
   if (length(xi) != length(pc))
      stop("Length of 'xi' does not match length of 'pc'.", call.=FALSE)
   if (any(pc < 0) || any(pc > 1))
      stop("The control group risk 'pc' must be between 0 and 1.", call.=FALSE)
   return(exp(xi) / (pc * (exp(xi) - 1) + 1))
}

############################################################################

transf.lnortod.norm <- function(xi)
   xi / 1.65

transf.lnortod.logis <- function(xi)
   sqrt(3) / base::pi * xi

transf.dtolnor.norm <- function(xi)
   xi * 1.65

transf.dtolnor.logis <- function(xi)
   xi / sqrt(3) * base::pi

transf.lnortortet.pearson <- function(xi)
   cos(base::pi / (1 + sqrt(exp(xi))))

transf.lnortortet.digby <- function(xi)
   (exp(xi)^(3/4) - 1) / (exp(xi)^(3/4) + 1)

############################################################################