R/f_design_sample_size_calculator.R

Defines functions .addNumberOfSubjectsToPowerResult .addStudyDurationToDesignPlan .hideFutilityStopsIfNotApplicable getPowerSurvival .getNumberOfSubjects .getNumberOfSubjectsInner getPowerRates getPowerMeans .createDesignPlanRates .createDesignPlanMeans .getSampleSizeSequentialSurvival .getSampleSizeFixedSurvival .getEventsFixed .getEventProbabilities .getEventProbabilitiesOverall .getEventProbabilitiesGroupwise .getLambda getNumberOfSubjects getEventProbabilities .getEventProbabilityFunctionVec .getEventProbabilityFunction .getPiecewiseExpStartTimesWithoutLeadingZero .getSampleSizeSequentialRates .getSampleSizeFixedRates .getFarringtonManningValues .getFarringtonManningValuesRatio .getFarringtonManningValuesDiff .getColumnCumSum .getSampleSizeSequentialMeans .getSampleSizeFixedMeans .getSampleSize .warnInCaseOfDefinedPiValue .initDesignPlanSurvival .initDesignPlanSurvivalByPiecewiseSurvivalTimeObject .isUserDefinedMaxNumberOfSubjects .createDesignPlanSurvival .getSampleSizeSurvival getSampleSizeSurvival getSampleSizeRates .warnInCaseOfTwoSidedPowerIsDisabled .warnInCaseOfTwoSidedPowerArgument getSampleSizeMeans .getEffectScaleBoundaryDataSurvival .getEffectScaleBoundaryDataRates .getEffectScaleBoundaryDataMeans .addEffectScaleBoundaryDataToDesignPlan

Documented in getEventProbabilities getNumberOfSubjects getPowerMeans getPowerRates getPowerSurvival getSampleSizeMeans getSampleSizeRates getSampleSizeSurvival

## |
## |  *Sample size calculator*
## |
## |  This file is part of the R package rpact:
## |  Confirmatory Adaptive Clinical Trial Design and Analysis
## |
## |  Author: Gernot Wassmer, PhD, and Friedrich Pahlke, PhD
## |  Licensed under "GNU Lesser General Public License" version 3
## |  License text can be found here: https://www.r-project.org/Licenses/LGPL-3
## |
## |  RPACT company website: https://www.rpact.com
## |  rpact package website: https://www.rpact.org
## |
## |  Contact us for information about our services: info@rpact.com
## |
## |  File version: $Revision: 7126 $
## |  Last changed: $Date: 2023-06-23 14:26:39 +0200 (Fr, 23 Jun 2023) $
## |  Last changed by: $Author: pahlke $
## |

#' @include f_core_utilities.R
NULL

.addEffectScaleBoundaryDataToDesignPlan <- function(designPlan) {
    .assertIsTrialDesignPlan(designPlan)

    design <- designPlan$.design
    if (.isTrialDesignPlanMeans(designPlan)) {
        if (design$kMax == 1 && designPlan$.isSampleSizeObject()) {
            designPlan$maxNumberOfSubjects <- designPlan$nFixed
        }

        boundaries <- .getEffectScaleBoundaryDataMeans(designPlan)
    } else if (.isTrialDesignPlanRates(designPlan)) {
        if (designPlan$.isSampleSizeObject()) { # comes from getSampleSize
            if (designPlan$groups == 1) {
                designPlan$directionUpper <- (designPlan$pi1 > designPlan$thetaH0)
            } else {
                if (designPlan$riskRatio) {
                    designPlan$directionUpper <- (designPlan$pi1 / designPlan$pi2 > designPlan$thetaH0)
                } else {
                    designPlan$directionUpper <- (designPlan$pi1 - designPlan$pi2 > designPlan$thetaH0)
                }
            }
            designPlan$.setParameterType("directionUpper", C_PARAM_GENERATED)
        }
        if (design$kMax == 1 && designPlan$.isSampleSizeObject()) {
            designPlan$maxNumberOfSubjects <- designPlan$nFixed
        }
        boundaries <- .getEffectScaleBoundaryDataRates(designPlan)
    } else if (.isTrialDesignPlanSurvival(designPlan)) {
        if (designPlan$.isSampleSizeObject()) { # comes from getSampleSize
            designPlan$directionUpper <- (designPlan$hazardRatio > designPlan$thetaH0)
            designPlan$.setParameterType("directionUpper", C_PARAM_GENERATED)
        }

        if (design$kMax == 1 && designPlan$.isSampleSizeObject()) {
            designPlan$eventsPerStage <- matrix(designPlan$eventsFixed, nrow = 1)
        }
        boundaries <- .getEffectScaleBoundaryDataSurvival(designPlan)
    }

    if (designPlan$.design$sided == 1) {
        designPlan$criticalValuesEffectScale <- boundaries$criticalValuesEffectScaleUpper
        designPlan$.setParameterType("criticalValuesEffectScale", C_PARAM_GENERATED)
    } else {
        if (all(boundaries$criticalValuesEffectScaleLower < boundaries$criticalValuesEffectScaleUpper, na.rm = TRUE)) {
            designPlan$criticalValuesEffectScaleLower <- boundaries$criticalValuesEffectScaleLower
            designPlan$criticalValuesEffectScaleUpper <- boundaries$criticalValuesEffectScaleUpper
        } else {
            designPlan$criticalValuesEffectScaleLower <- boundaries$criticalValuesEffectScaleUpper
            designPlan$criticalValuesEffectScaleUpper <- boundaries$criticalValuesEffectScaleLower
        }
        designPlan$.setParameterType("criticalValuesEffectScaleUpper", C_PARAM_GENERATED)
        designPlan$.setParameterType("criticalValuesEffectScaleLower", C_PARAM_GENERATED)
    }

    if (!.isTrialDesignFisher(design) && any(design$futilityBounds > C_FUTILITY_BOUNDS_DEFAULT)) {
        if (design$sided == 1) {
            designPlan$futilityBoundsEffectScale <- round(boundaries$futilityBoundsEffectScaleUpper, 8)
            designPlan$.setParameterType("futilityBoundsEffectScale", C_PARAM_GENERATED)
        } else {
            if (all(designPlan$futilityBoundsEffectScaleLower < designPlan$futilityBoundsEffectScaleUpper, na.rm = TRUE)) {
                designPlan$futilityBoundsEffectScaleLower <- round(boundaries$futilityBoundsEffectScaleLower, 8)
                designPlan$futilityBoundsEffectScaleUpper <- round(boundaries$futilityBoundsEffectScaleUpper, 8)
            } else {
                designPlan$futilityBoundsEffectScaleLower <- round(boundaries$futilityBoundsEffectScaleUpper, 8)
                designPlan$futilityBoundsEffectScaleUpper <- round(boundaries$futilityBoundsEffectScaleLower, 8)
            }
            designPlan$.setParameterType("futilityBoundsEffectScaleLower", C_PARAM_GENERATED)
            designPlan$.setParameterType("futilityBoundsEffectScaleUpper", C_PARAM_GENERATED)
        }
    }
}

.getEffectScaleBoundaryDataMeans <- function(designPlan) {
    design <- designPlan$.design
    thetaH0 <- designPlan$thetaH0
    stDev <- designPlan$stDev
    maxNumberOfSubjects <- designPlan$maxNumberOfSubjects
    allocationRatioPlanned <- designPlan$allocationRatioPlanned
    directionUpper <- designPlan$directionUpper

    # initialize effect scale matrix
    futilityBoundsEffectScaleUpper <- rep(NA_real_, design$kMax - 1)
    futilityBoundsEffectScaleLower <- rep(NA_real_, design$kMax - 1)

    if (designPlan$normalApproximation) {
        criticalValues <- design$criticalValues
        futilityBounds <- design$futilityBounds
    } else {
        criticalValues <- stats::qt(
            1 - design$stageLevels,
            design$informationRates %*% t(maxNumberOfSubjects) - designPlan$groups
        )

        # outside validated range
        numberOfNAs <- sum(as.vector(criticalValues) > 50, na.rm = TRUE)
        criticalValues[criticalValues > 50] <- NA_real_
        if (any(is.na(criticalValues) & (design$criticalValues < 8))) {
            warning("The computation of ", .integerToWrittenNumber(numberOfNAs),
                " efficacy boundar", ifelse(numberOfNAs == 1, "y", "ies"), " on ",
                "treatment effect scale not performed presumably due to too small df",
                call. = FALSE
            )
        }

        futilityBounds <- stats::qt(
            stats::pnorm(design$futilityBounds),
            design$informationRates[1:(design$kMax - 1)] %*% t(maxNumberOfSubjects) - designPlan$groups
        )

        # outside validated range
        futilityBounds[futilityBounds < -50] <- NA_real_
    }
    futilityBounds[!is.na(futilityBounds) & futilityBounds <= C_FUTILITY_BOUNDS_DEFAULT] <- NA_real_

    if (designPlan$groups == 1) {
        criticalValuesEffectScaleUpper <- thetaH0 + criticalValues * stDev /
            sqrt(design$informationRates %*% t(maxNumberOfSubjects))
        criticalValuesEffectScaleLower <- thetaH0 - criticalValues * stDev /
            sqrt(design$informationRates %*% t(maxNumberOfSubjects))
        if (!.isTrialDesignFisher(design) && !all(is.na(futilityBounds))) {
            futilityBoundsEffectScaleUpper <- thetaH0 + futilityBounds * stDev /
                sqrt(design$informationRates[1:(design$kMax - 1)] %*% t(maxNumberOfSubjects))
        }
        if (!.isTrialDesignFisher(design) && design$sided == 2 && design$kMax > 1 &&
                (design$typeOfDesign == C_TYPE_OF_DESIGN_PT || !is.null(design$typeBetaSpending) && design$typeBetaSpending != "none")) {
            futilityBoundsEffectScaleLower <- thetaH0 - futilityBounds * stDev /
                sqrt(design$informationRates[1:(design$kMax - 1)] %*% t(maxNumberOfSubjects))
        }
    } else if (!designPlan$meanRatio) {
        criticalValuesEffectScaleUpper <- thetaH0 + criticalValues * stDev *
            (1 + allocationRatioPlanned) / (sqrt(allocationRatioPlanned *
                design$informationRates %*% t(maxNumberOfSubjects)))
        criticalValuesEffectScaleLower <- thetaH0 - criticalValues * stDev *
            (1 + allocationRatioPlanned) / (sqrt(allocationRatioPlanned *
                design$informationRates %*% t(maxNumberOfSubjects)))
        if (!.isTrialDesignFisher(design) && !all(is.na(futilityBounds))) {
            futilityBoundsEffectScaleUpper <- thetaH0 + futilityBounds * stDev *
                (1 + allocationRatioPlanned) / (sqrt(allocationRatioPlanned *
                    design$informationRates[1:(design$kMax - 1)] %*% t(maxNumberOfSubjects)))
        }
        if (!.isTrialDesignFisher(design) && design$sided == 2 && design$kMax > 1 &&
                (design$typeOfDesign == C_TYPE_OF_DESIGN_PT || !is.null(design$typeBetaSpending) && design$typeBetaSpending != "none")) {
            futilityBoundsEffectScaleLower <- thetaH0 - futilityBounds * stDev *
                (1 + allocationRatioPlanned) / (sqrt(allocationRatioPlanned *
                    design$informationRates[1:(design$kMax - 1)] %*% t(maxNumberOfSubjects)))
        }
    } else {
        criticalValuesEffectScaleUpper <- thetaH0 + criticalValues * stDev *
            sqrt(1 + 1 / allocationRatioPlanned + thetaH0^2 * (1 + allocationRatioPlanned)) /
            (sqrt(design$informationRates %*% t(maxNumberOfSubjects)))
        criticalValuesEffectScaleLower <- thetaH0 - criticalValues * stDev *
            sqrt(1 + 1 / allocationRatioPlanned + thetaH0^2 * (1 + allocationRatioPlanned)) /
            (sqrt(design$informationRates %*% t(maxNumberOfSubjects)))


        if (!.isTrialDesignFisher(design) && !all(is.na(futilityBounds))) {
            futilityBoundsEffectScaleUpper <- thetaH0 + futilityBounds * stDev *
                sqrt(1 + 1 / allocationRatioPlanned + thetaH0^2 * (1 + allocationRatioPlanned)) /
                (sqrt(design$informationRates[1:(design$kMax - 1)] %*% t(maxNumberOfSubjects)))
        }
        if (!.isTrialDesignFisher(design) && design$sided == 2 && design$kMax > 1 &&
                (design$typeOfDesign == C_TYPE_OF_DESIGN_PT || !is.null(design$typeBetaSpending) && design$typeBetaSpending != "none")) {
            futilityBoundsEffectScaleLower <- thetaH0 - futilityBounds * stDev *
                sqrt(1 + 1 / allocationRatioPlanned + thetaH0^2 * (1 + allocationRatioPlanned)) /
                (sqrt(design$informationRates[1:(design$kMax - 1)] %*% t(maxNumberOfSubjects)))
        }
    }

    directionUpper[is.na(directionUpper)] <- TRUE
    if (length(directionUpper) > 0 && all(!directionUpper)) {
        criticalValuesEffectScaleUpper <- -criticalValuesEffectScaleUpper + 2 * thetaH0
        criticalValuesEffectScaleLower <- -criticalValuesEffectScaleLower + 2 * thetaH0
        if (!all(is.na(futilityBoundsEffectScaleUpper))) {
            futilityBoundsEffectScaleUpper <- -futilityBoundsEffectScaleUpper + 2 * thetaH0
            futilityBoundsEffectScaleLower <- -futilityBoundsEffectScaleLower + 2 * thetaH0
        }
    }
    if (designPlan$meanRatio) {
        criticalValuesEffectScaleUpper[!is.na(criticalValuesEffectScaleUpper) & criticalValuesEffectScaleUpper <= 0] <- NA_real_
        criticalValuesEffectScaleLower[!is.na(criticalValuesEffectScaleLower) & criticalValuesEffectScaleLower <= 0] <- NA_real_
        futilityBoundsEffectScaleUpper[!is.na(futilityBoundsEffectScaleUpper) & futilityBoundsEffectScaleUpper <= 0] <- NA_real_
        futilityBoundsEffectScaleLower[!is.na(futilityBoundsEffectScaleLower) & futilityBoundsEffectScaleLower <= 0] <- NA_real_
    }

    return(list(
        criticalValuesEffectScaleUpper = matrix(criticalValuesEffectScaleUpper, nrow = design$kMax),
        criticalValuesEffectScaleLower = matrix(criticalValuesEffectScaleLower, nrow = design$kMax),
        futilityBoundsEffectScaleUpper = matrix(futilityBoundsEffectScaleUpper, nrow = design$kMax - 1),
        futilityBoundsEffectScaleLower = matrix(futilityBoundsEffectScaleLower, nrow = design$kMax - 1)
    ))
}

.getEffectScaleBoundaryDataRates <- function(designPlan) {
    design <- designPlan$.design
    thetaH0 <- designPlan$thetaH0
    pi2 <- designPlan$pi2
    maxNumberOfSubjects <- designPlan$maxNumberOfSubjects
    allocationRatioPlanned <- designPlan$allocationRatioPlanned
    directionUpper <- designPlan$directionUpper

    nParameters <- length(maxNumberOfSubjects)

    directionUpper[is.na(directionUpper)] <- TRUE

    criticalValuesEffectScaleUpper <- matrix(, nrow = design$kMax, ncol = nParameters)
    criticalValuesEffectScaleLower <- matrix(, nrow = design$kMax, ncol = nParameters)
    futilityBoundsEffectScaleUpper <- matrix(, nrow = design$kMax - 1, ncol = nParameters)
    futilityBoundsEffectScaleLower <- matrix(, nrow = design$kMax - 1, ncol = nParameters)
    if (length(allocationRatioPlanned) == 1) {
        allocationRatioPlanned <- rep(allocationRatioPlanned, nParameters)
    }
    futilityBounds <- design$futilityBounds
    futilityBounds[!is.na(futilityBounds) & futilityBounds <= C_FUTILITY_BOUNDS_DEFAULT] <- NA_real_

    if (designPlan$groups == 1) {
        n1 <- design$informationRates %*% t(maxNumberOfSubjects)
        for (j in (1:nParameters)) {
            criticalValuesEffectScaleUpper[, j] <- thetaH0 + (2 * directionUpper[j] - 1) *
                design$criticalValues * sqrt(thetaH0 * (1 - thetaH0)) / sqrt(n1[, j])
            if (design$sided == 2) {
                criticalValuesEffectScaleLower[, j] <- thetaH0 - (2 * directionUpper[j] - 1) *
                    design$criticalValues * sqrt(thetaH0 * (1 - thetaH0)) / sqrt(n1[, j])
            }
            if (!.isTrialDesignFisher(design) && !all(is.na(futilityBounds))) {
                futilityBoundsEffectScaleUpper[, j] <- thetaH0 + (2 * directionUpper[j] - 1) *
                    futilityBounds * sqrt(thetaH0 * (1 - thetaH0)) /
                    sqrt(n1[1:(design$kMax - 1), j])
            }
            if (!.isTrialDesignFisher(design) && design$sided == 2 && design$kMax > 1 &&
                    (design$typeOfDesign == C_TYPE_OF_DESIGN_PT || !is.null(design$typeBetaSpending) && design$typeBetaSpending != "none")) {
                futilityBoundsEffectScaleLower[, j] <- thetaH0 - (2 * directionUpper[j] - 1) *
                    futilityBounds * sqrt(thetaH0 * (1 - thetaH0)) /
                    sqrt(n1[1:(design$kMax - 1), j])
            }
        }
    } else if (!designPlan$riskRatio) {
        boundaries <- design$criticalValues

        # calculate pi1 that solves (pi1 - pi2 - thetaH0) / SE(pi1 - pi2 - thetaH0)
        # = crit by using Farrington & Manning approach
        for (j in (1:nParameters)) {
            n1 <- allocationRatioPlanned[j] * design$informationRates *
                maxNumberOfSubjects[j] / (1 + allocationRatioPlanned[j])
            n2 <- n1 / allocationRatioPlanned[j]

            for (i in (1:length(boundaries))) {
                tryCatch(
                    {
                        pi1Bound <- uniroot(
                            function(x) {
                                fm <- .getFarringtonManningValues(
                                    rate1 = x, rate2 = pi2, theta = thetaH0,
                                    allocation = allocationRatioPlanned[j], method = "diff"
                                )
                                (x - pi2 - thetaH0) / sqrt(fm$ml1 * (1 - fm$ml1) /
                                    n1[i] + fm$ml2 * (1 - fm$ml2) / n2[i]) -
                                    (2 * directionUpper[j] - 1) * boundaries[i]
                            },
                            lower = 0, upper = 1, tol = .Machine$double.eps^0.5
                        )$root
                    },
                    error = function(e) {
                        pi1Bound <<- NA_real_
                    }
                )

                # difference to pi2
                criticalValuesEffectScaleUpper[i, j] <- pi1Bound - pi2
            }
            if (design$sided == 2) {
                for (i in (1:length(boundaries))) {
                    tryCatch(
                        {
                            pi1Bound <- uniroot(
                                function(x) {
                                    fm <- .getFarringtonManningValues(
                                        rate1 = x, rate2 = pi2, theta = thetaH0,
                                        allocation = allocationRatioPlanned[j], method = "diff"
                                    )
                                    (x - pi2 - thetaH0) / sqrt(fm$ml1 * (1 - fm$ml1) /
                                        n1[i] + fm$ml2 * (1 - fm$ml2) / n2[i]) +
                                        (2 * directionUpper[j] - 1) * boundaries[i]
                                },
                                lower = 0, upper = 1, tol = .Machine$double.eps^0.5
                            )$root
                        },
                        error = function(e) {
                            pi1Bound <<- NA_real_
                        }
                    )

                    # difference to pi2
                    criticalValuesEffectScaleLower[i, j] <- pi1Bound - pi2
                }
            }
        }
        if (!.isTrialDesignFisher(design) && !all(is.na(futilityBounds))) {
            boundaries <- futilityBounds
            for (j in (1:nParameters)) {
                n1 <- allocationRatioPlanned[j] * design$informationRates *
                    maxNumberOfSubjects[j] / (1 + allocationRatioPlanned[j])
                n2 <- n1 / allocationRatioPlanned[j]
                for (i in (1:length(boundaries))) {
                    tryCatch(
                        {
                            pi1Bound <- uniroot(
                                function(x) {
                                    fm <- .getFarringtonManningValues(
                                        rate1 = x, rate2 = pi2, theta = thetaH0,
                                        allocation = allocationRatioPlanned[j], method = "diff"
                                    )
                                    (x - pi2 - thetaH0) / sqrt(fm$ml1 * (1 - fm$ml1) / n1[i] +
                                        fm$ml2 * (1 - fm$ml2) / n2[i]) -
                                        (2 * directionUpper[j] - 1) * boundaries[i]
                                },
                                lower = 0, upper = 1, tol = .Machine$double.eps^0.5
                            )$root
                        },
                        error = function(e) {
                            pi1Bound <<- NA_real_
                        }
                    )

                    # difference to pi2
                    futilityBoundsEffectScaleUpper[i, j] <- pi1Bound - pi2
                }
            }
        }

        if (!.isTrialDesignFisher(design) && design$sided == 2 && design$kMax > 1 &&
                (design$typeOfDesign == C_TYPE_OF_DESIGN_PT || !is.null(design$typeBetaSpending) && design$typeBetaSpending != "none")) {
            boundaries <- -futilityBounds
            for (j in (1:nParameters)) {
                n1 <- allocationRatioPlanned[j] * design$informationRates *
                    maxNumberOfSubjects[j] / (1 + allocationRatioPlanned[j])
                n2 <- n1 / allocationRatioPlanned[j]
                for (i in (1:length(boundaries))) {
                    tryCatch(
                        {
                            pi1Bound <- uniroot(
                                function(x) {
                                    fm <- .getFarringtonManningValues(
                                        rate1 = x, rate2 = pi2, theta = thetaH0,
                                        allocation = allocationRatioPlanned[j], method = "diff"
                                    )
                                    (x - pi2 - thetaH0) / sqrt(fm$ml1 * (1 - fm$ml1) / n1[i] +
                                        fm$ml2 * (1 - fm$ml2) / n2[i]) -
                                        (2 * directionUpper[j] - 1) * boundaries[i]
                                },
                                lower = 0, upper = 1, tol = .Machine$double.eps^0.5
                            )$root
                        },
                        error = function(e) {
                            pi1Bound <<- NA_real_
                        }
                    )
                    futilityBoundsEffectScaleLower[i, j] <- pi1Bound - pi2 # difference to pi2
                }
            }
        }
    } else {
        boundaries <- design$criticalValues
        # calculate pi1 that solves (pi1 - thetaH0 * pi2) / SE(pi1 - thetaH0 * pi2)
        # = crit by using Farrington & Manning approach
        for (j in (1:nParameters)) {
            n1 <- allocationRatioPlanned[j] * design$informationRates * maxNumberOfSubjects[j] /
                (1 + allocationRatioPlanned[j])
            n2 <- n1 / allocationRatioPlanned[j]
            for (i in (1:length(boundaries))) {
                tryCatch(
                    {
                        pi1Bound <- uniroot(
                            function(x) {
                                fm <- .getFarringtonManningValues(
                                    rate1 = x, rate2 = pi2, theta = thetaH0,
                                    allocation = allocationRatioPlanned[j], method = "ratio"
                                )
                                (x - thetaH0 * pi2) / sqrt(fm$ml1 * (1 - fm$ml1) / n1[i] +
                                    thetaH0^2 * fm$ml2 * (1 - fm$ml2) / n2[i]) -
                                    (2 * directionUpper[j] - 1) * boundaries[i]
                            },
                            lower = 0, upper = 1, tol = .Machine$double.eps^0.5
                        )$root
                    },
                    error = function(e) {
                        pi1Bound <<- NA_real_
                    }
                )

                # ratio to pi2
                criticalValuesEffectScaleUpper[i, j] <- pi1Bound / pi2
            }
            if (design$sided == 2) {
                for (i in (1:length(boundaries))) {
                    tryCatch(
                        {
                            pi1Bound <- uniroot(
                                function(x) {
                                    fm <- .getFarringtonManningValues(
                                        rate1 = x, rate2 = pi2, theta = thetaH0,
                                        allocation = allocationRatioPlanned[j], method = "ratio"
                                    )
                                    (x - thetaH0 * pi2) / sqrt(fm$ml1 * (1 - fm$ml1) / n1[i] +
                                        thetaH0^2 * fm$ml2 * (1 - fm$ml2) / n2[i]) +
                                        (2 * directionUpper[j] - 1) * boundaries[i]
                                },
                                lower = 0, upper = 1, tol = .Machine$double.eps^0.5
                            )$root
                        },
                        error = function(e) {
                            pi1Bound <<- NA_real_
                        }
                    )

                    # ratio to pi2
                    criticalValuesEffectScaleLower[i, j] <- pi1Bound / pi2
                }
            }
        }
        if (!.isTrialDesignFisher(design) && !all(is.na(futilityBounds))) {
            boundaries <- futilityBounds
            for (j in (1:nParameters)) {
                n1 <- allocationRatioPlanned[j] * design$informationRates * maxNumberOfSubjects[j] /
                    (1 + allocationRatioPlanned[j])
                n2 <- n1 / allocationRatioPlanned[j]
                for (i in (1:length(boundaries))) {
                    tryCatch(
                        {
                            pi1Bound <- uniroot(
                                function(x) {
                                    fm <- .getFarringtonManningValues(
                                        rate1 = x, rate2 = pi2, theta = thetaH0,
                                        allocation = allocationRatioPlanned[j], method = "ratio"
                                    )
                                    (x - thetaH0 * pi2) / sqrt(fm$ml1 * (1 - fm$ml1) / n1[i] +
                                        thetaH0^2 * fm$ml2 * (1 - fm$ml2) / n2[i]) -
                                        (2 * directionUpper[j] - 1) * boundaries[i]
                                },
                                lower = 0, upper = 1, tol = .Machine$double.eps^0.5
                            )$root
                        },
                        error = function(e) {
                            pi1Bound <<- NA_real_
                        }
                    )

                    # ratio to pi2
                    futilityBoundsEffectScaleUpper[i, j] <- pi1Bound / pi2
                }
            }
        }
        if (!.isTrialDesignFisher(design) && design$sided == 2 && design$kMax > 1 &&
                (design$typeOfDesign == C_TYPE_OF_DESIGN_PT || !is.null(design$typeBetaSpending) && design$typeBetaSpending != "none")) {
            boundaries <- -futilityBounds
            for (j in (1:nParameters)) {
                n1 <- allocationRatioPlanned[j] * design$informationRates * maxNumberOfSubjects[j] /
                    (1 + allocationRatioPlanned[j])
                n2 <- n1 / allocationRatioPlanned[j]
                for (i in (1:length(boundaries))) {
                    tryCatch(
                        {
                            pi1Bound <- uniroot(
                                function(x) {
                                    fm <- .getFarringtonManningValues(
                                        rate1 = x, rate2 = pi2, theta = thetaH0,
                                        allocation = allocationRatioPlanned[j], method = "ratio"
                                    )
                                    (x - thetaH0 * pi2) / sqrt(fm$ml1 * (1 - fm$ml1) / n1[i] +
                                        thetaH0^2 * fm$ml2 * (1 - fm$ml2) / n2[i]) -
                                        (2 * directionUpper[j] - 1) * boundaries[i]
                                },
                                lower = 0, upper = 1, tol = .Machine$double.eps^0.5
                            )$root
                        },
                        error = function(e) {
                            pi1Bound <<- NA_real_
                        }
                    )

                    # ratio to pi2
                    futilityBoundsEffectScaleLower[i, j] <- pi1Bound / pi2
                }
            }
        }
    }
    return(list(
        criticalValuesEffectScaleUpper = matrix(criticalValuesEffectScaleUpper, nrow = design$kMax),
        criticalValuesEffectScaleLower = matrix(criticalValuesEffectScaleLower, nrow = design$kMax),
        futilityBoundsEffectScaleUpper = matrix(futilityBoundsEffectScaleUpper, nrow = design$kMax - 1),
        futilityBoundsEffectScaleLower = matrix(futilityBoundsEffectScaleLower, nrow = design$kMax - 1)
    ))
}

.getEffectScaleBoundaryDataSurvival <- function(designPlan) {
    design <- designPlan$.design
    thetaH0 <- designPlan$thetaH0
    eventsPerStage <- designPlan$eventsPerStage
    allocationRatioPlanned <- designPlan$allocationRatioPlanned
    directionUpper <- designPlan$directionUpper

    if (design$kMax == 1) {
        nParameters <- length(eventsPerStage)
    } else {
        nParameters <- ncol(eventsPerStage)
    }

    directionUpper[is.na(directionUpper)] <- TRUE

    if (length(allocationRatioPlanned) == 1) {
        allocationRatioPlanned <- rep(allocationRatioPlanned, nParameters)
    }

    futilityBounds <- design$futilityBounds
    futilityBounds[!is.na(futilityBounds) & futilityBounds <= C_FUTILITY_BOUNDS_DEFAULT] <- NA_real_

    criticalValues <- design$criticalValues

    criticalValuesEffectScaleUpper <- matrix(, nrow = design$kMax, ncol = nParameters)
    criticalValuesEffectScaleLower <- matrix(, nrow = design$kMax, ncol = nParameters)
    futilityBoundsEffectScaleUpper <- matrix(, nrow = design$kMax - 1, ncol = nParameters)
    futilityBoundsEffectScaleLower <- matrix(, nrow = design$kMax - 1, ncol = nParameters)

    for (j in (1:nParameters)) {
        if (design$sided == 1) {
            criticalValuesEffectScaleUpper[, j] <- thetaH0 * (exp((2 * directionUpper[j] - 1) * criticalValues *
                (1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
                    eventsPerStage[, j])))
        } else {
            criticalValuesEffectScaleUpper[, j] <- thetaH0 * (exp((2 * directionUpper[j] - 1) * criticalValues *
                (1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
                    eventsPerStage[, j])))
            criticalValuesEffectScaleLower[, j] <- thetaH0 * (exp(-(2 * directionUpper[j] - 1) * criticalValues *
                (1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
                    eventsPerStage[, j])))
        }

        if (!.isTrialDesignFisher(design) && !all(is.na(futilityBounds))) {
            futilityBoundsEffectScaleUpper[, j] <- thetaH0 * (exp((2 * directionUpper[j] - 1) * futilityBounds *
                (1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
                    eventsPerStage[1:(design$kMax - 1), j])))
        }
        if (!.isTrialDesignFisher(design) && design$sided == 2 && design$kMax > 1 &&
                (design$typeOfDesign == C_TYPE_OF_DESIGN_PT || !is.null(design$typeBetaSpending) && design$typeBetaSpending != "none")) {
            futilityBoundsEffectScaleLower[, j] <- thetaH0 * (exp(-(2 * directionUpper[j] - 1) * futilityBounds *
                (1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
                    eventsPerStage[1:(design$kMax - 1), j])))
        }
    }

    return(list(
        criticalValuesEffectScaleUpper = matrix(criticalValuesEffectScaleUpper, nrow = design$kMax),
        criticalValuesEffectScaleLower = matrix(criticalValuesEffectScaleLower, nrow = design$kMax),
        futilityBoundsEffectScaleUpper = matrix(futilityBoundsEffectScaleUpper, nrow = design$kMax - 1),
        futilityBoundsEffectScaleLower = matrix(futilityBoundsEffectScaleLower, nrow = design$kMax - 1)
    ))
}


#' @title
#' Get Sample Size Means
#'
#' @description
#' Returns the sample size for testing means in one or two samples.
#'
#' @inheritParams param_design_with_default
#' @inheritParams param_groups
#' @param normalApproximation The type of computation of the p-values. If \code{TRUE}, the variance is
#'        assumed to be known, default is \code{FALSE}, i.e., the calculations are performed
#'        with the t distribution.
#' @param meanRatio If \code{TRUE}, the sample size for
#'        one-sided testing of H0: \code{mu1 / mu2 = thetaH0} is calculated, default is \code{FALSE}.
#' @inheritParams param_thetaH0
#' @inheritParams param_alternative
#' @inheritParams param_stDev
#' @inheritParams param_allocationRatioPlanned_sampleSize
#' @inheritParams param_three_dots
#'
#' @details
#' At given design the function calculates the stage-wise (non-cumulated) and maximum
#' sample size for testing means.
#' In a two treatment groups design, additionally, an allocation ratio = n1/n2 can be specified.
#' A null hypothesis value thetaH0 != 0 for testing the difference of two means or
#' thetaH0 != 1 for testing the ratio of two means can be specified.
#' Critical bounds and stopping for futility bounds are provided at the effect scale
#' (mean, mean difference, or mean ratio, respectively) for each sample size calculation separately.
#'
#' @template return_object_trial_design_plan
#' @template how_to_get_help_for_generics
#'
#' @family sample size functions
#'
#' @template examples_get_sample_size_means
#'
#' @export
#'
getSampleSizeMeans <- function(design = NULL, ...,
        groups = 2,
        normalApproximation = FALSE,
        meanRatio = FALSE,
        thetaH0 = ifelse(meanRatio, 1, 0),
        alternative = seq(0.2, 1, 0.2), # C_ALTERNATIVE_DEFAULT
        stDev = 1, # C_STDEV_DEFAULT
        allocationRatioPlanned = NA_real_ # C_ALLOCATION_RATIO_DEFAULT
        ) {
    if (is.null(design)) {
        design <- .getDefaultDesign(..., type = "sampleSize")
        .warnInCaseOfUnknownArguments(
            functionName = "getSampleSizeMeans",
            ignore = .getDesignArgumentsToIgnoreAtUnknownArgumentCheck(design, powerCalculationEnabled = FALSE), ...
        )
    } else {
        .assertIsTrialDesign(design)
        .warnInCaseOfUnknownArguments(functionName = "getSampleSizeMeans", ...)
        .warnInCaseOfTwoSidedPowerArgument(...)
    }

    designPlan <- .createDesignPlanMeans(
        objectType = "sampleSize",
        design = design, normalApproximation = normalApproximation, meanRatio = meanRatio,
        thetaH0 = thetaH0, alternative = alternative, stDev = stDev, groups = groups,
        allocationRatioPlanned = allocationRatioPlanned, ...
    )

    return(.getSampleSize(designPlan))
}

.warnInCaseOfTwoSidedPowerArgument <- function(...) {
    args <- list(...)
    argNames <- names(args)
    if ("twoSidedPower" %in% argNames) {
        warning("'twoSidedPower' can only be defined in 'design'", call. = FALSE)
    }
}

.warnInCaseOfTwoSidedPowerIsDisabled <- function(design) {
    if (design$sided == 2 && !is.na(design$twoSidedPower) && !design$twoSidedPower &&
            design$.getParameterType("twoSidedPower") == C_PARAM_USER_DEFINED) {
        warning("design$twoSidedPower = FALSE will be ignored because design$sided = 2", call. = FALSE)
    }
}

#' @title
#' Get Sample Size Rates
#'
#' @description
#' Returns the sample size for testing rates in one or two samples.
#'
#' @inheritParams param_design_with_default
#' @inheritParams param_groups
#' @param normalApproximation If \code{FALSE}, the sample size
#'        for the case of one treatment group is calculated exactly using the binomial distribution,
#'        default is \code{TRUE}.
#' @param riskRatio If \code{TRUE}, the sample size for one-sided
#'        testing of H0: \code{pi1 / pi2 = thetaH0} is calculated, default is \code{FALSE}.
#' @inheritParams param_thetaH0
#' @inheritParams param_pi1_rates
#' @inheritParams param_pi2_rates
#' @inheritParams param_allocationRatioPlanned_sampleSize
#' @inheritParams param_three_dots
#'
#' @details
#' At given design the function calculates the stage-wise (non-cumulated) and maximum sample size for testing rates.
#' In a two treatment groups design, additionally, an allocation ratio = n1/n2 can be specified.
#' If a null hypothesis value thetaH0 != 0 for testing the difference of two rates
#' thetaH0 != 1 for testing the risk ratio is specified, the sample size
#' formula according to Farrington & Manning (Statistics in Medicine, 1990) is used.
#' Critical bounds and stopping for futility bounds are provided at the effect scale
#' (rate, rate difference, or rate ratio, respectively) for each sample size calculation separately.
#' For the two-sample case, the calculation here is performed at fixed pi2 as given as argument
#' in the function.
#'
#' @template return_object_trial_design_plan
#' @template how_to_get_help_for_generics
#'
#' @family sample size functions
#'
#' @template examples_get_sample_size_rates
#'
#' @export
#'
getSampleSizeRates <- function(design = NULL, ...,
        groups = 2,
        normalApproximation = TRUE,
        riskRatio = FALSE,
        thetaH0 = ifelse(riskRatio, 1, 0),
        pi1 = c(0.4, 0.5, 0.6), # C_PI_1_SAMPLE_SIZE_DEFAULT
        pi2 = 0.2, # C_PI_2_DEFAULT
        allocationRatioPlanned = NA_real_ # C_ALLOCATION_RATIO_DEFAULT
        ) {
    if (is.null(design)) {
        design <- .getDefaultDesign(..., type = "sampleSize")
        .warnInCaseOfUnknownArguments(
            functionName = "getSampleSizeRates",
            ignore = .getDesignArgumentsToIgnoreAtUnknownArgumentCheck(design, powerCalculationEnabled = FALSE), ...
        )
    } else {
        .assertIsTrialDesign(design)
        .warnInCaseOfUnknownArguments(functionName = "getSampleSizeRates", ...)
        .warnInCaseOfTwoSidedPowerArgument(...)
    }

    designPlan <- .createDesignPlanRates(
        objectType = "sampleSize",
        design = design, normalApproximation = normalApproximation, riskRatio = riskRatio,
        thetaH0 = thetaH0, pi1 = pi1, pi2 = pi2, groups = groups,
        allocationRatioPlanned = allocationRatioPlanned, ...
    )

    return(.getSampleSize(designPlan))
}

# Hidden parameter:
# @param accountForObservationTimes If \code{accountForObservationTimes = TRUE}, the number of
#        subjects is calculated assuming specific accrual and follow-up time, default is \code{TRUE}
#        (see details).
# If \code{accountForObservationTimes = FALSE}, only the event rates are used for the calculation
# of the maximum number of subjects.
# \code{accountForObservationTimes} can be selected as \code{FALSE}. In this case,
# the number of subjects is calculated from the event probabilities only.
# This kind of computation does not account for the specific accrual pattern and survival distribution.

#' @title
#' Get Sample Size Survival
#'
#' @description
#' Returns the sample size for testing the hazard ratio in a two treatment groups survival design.
#'
#' @inheritParams param_design_with_default
#' @inheritParams param_typeOfComputation
#' @inheritParams param_allocationRatioPlanned_sampleSize
#' @inheritParams param_thetaH0
#' @inheritParams param_lambda1
#' @inheritParams param_lambda2
#' @inheritParams param_pi1_survival
#' @inheritParams param_pi2_survival
#' @inheritParams param_median1
#' @inheritParams param_median2
#' @inheritParams param_piecewiseSurvivalTime
#' @inheritParams param_accrualTime
#' @inheritParams param_accrualIntensity
#' @inheritParams param_accrualIntensityType
#' @inheritParams param_eventTime
#' @inheritParams param_hazardRatio
#' @inheritParams param_kappa
#' @inheritParams param_dropoutRate1
#' @inheritParams param_dropoutRate2
#' @inheritParams param_dropoutTime
#' @param followUpTime The assumed (additional) follow-up time for the study, default is \code{6}.
#'        The total study duration is \code{accrualTime + followUpTime}.
#' @param maxNumberOfSubjects If \code{maxNumberOfSubjects > 0} is specified,
#'        the follow-up time for the required number of events is determined.
#' @inheritParams param_three_dots
#'
#' @details
#' At given design the function calculates the number of events and an estimate for the
#' necessary number of subjects for testing the hazard ratio in a survival design.
#' It also calculates the time when the required events are expected under the given
#' assumptions (exponentially, piecewise exponentially, or Weibull distributed survival times
#' and constant or non-constant piecewise accrual).
#' Additionally, an allocation ratio = \code{n1 / n2} can be specified where \code{n1} and \code{n2} are the number
#' of subjects in the two treatment groups.
#'
#' Optional argument \code{accountForObservationTimes}: if \code{accountForObservationTimes = TRUE}, the number of
#' subjects is calculated assuming specific accrual and follow-up time, default is \code{TRUE}.
#'
#' The formula of Kim & Tsiatis (Biometrics, 1990)
#' is used to calculate the expected number of events under the alternative
#' (see also Lakatos & Lan, Statistics in Medicine, 1992). These formulas are generalized
#' to piecewise survival times and non-constant piecewise accrual over time.\cr
#'
#' Optional argument \code{accountForObservationTimes}: if \code{accountForObservationTimes = FALSE},
#' only the event rates are used for the calculation of the maximum number of subjects.
#'
#' @template details_piecewise_survival
#'
#' @template details_piecewise_accrual
#'
#' @template return_object_trial_design_plan
#' @template how_to_get_help_for_generics
#'
#' @family sample size functions
#'
#' @template examples_get_sample_size_survival
#'
#' @export
#'
getSampleSizeSurvival <- function(design = NULL, ...,
        typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
        thetaH0 = 1, # C_THETA_H0_SURVIVAL_DEFAULT
        pi1 = NA_real_,
        pi2 = NA_real_,
        lambda1 = NA_real_,
        lambda2 = NA_real_,
        median1 = NA_real_,
        median2 = NA_real_,
        kappa = 1,
        hazardRatio = NA_real_,
        piecewiseSurvivalTime = NA_real_,
        allocationRatioPlanned = NA_real_, # C_ALLOCATION_RATIO_DEFAULT
        eventTime = 12, # C_EVENT_TIME_DEFAULT
        accrualTime = c(0, 12), # C_ACCRUAL_TIME_DEFAULT
        accrualIntensity = 0.1, # C_ACCRUAL_INTENSITY_DEFAULT
        accrualIntensityType = c("auto", "absolute", "relative"),
        followUpTime = NA_real_,
        maxNumberOfSubjects = NA_real_,
        dropoutRate1 = 0, # C_DROP_OUT_RATE_1_DEFAULT
        dropoutRate2 = 0, # C_DROP_OUT_RATE_2_DEFAULT
        dropoutTime = 12 # C_DROP_OUT_TIME_DEFAULT
        ) {
    if (is.null(design)) {
        design <- .getDefaultDesign(..., type = "sampleSize", ignore = c("accountForObservationTimes"))
        .warnInCaseOfUnknownArguments(
            functionName = "getSampleSizeSurvival",
            ignore = c(.getDesignArgumentsToIgnoreAtUnknownArgumentCheck(
                design,
                powerCalculationEnabled = FALSE
            ), "accountForObservationTimes"), ...
        )
    } else {
        .assertIsTrialDesign(design)
        .warnInCaseOfUnknownArguments(
            functionName = "getSampleSizeSurvival", ...,
            ignore = c("accountForObservationTimes")
        )
        .warnInCaseOfTwoSidedPowerArgument(...)
    }

    if (!is.na(maxNumberOfSubjects) && maxNumberOfSubjects == 0) {
        maxNumberOfSubjects <- NA_real_
    }

    # identify accrual time case
    accrualSetup <- getAccrualTime(
        accrualTime = accrualTime,
        accrualIntensity = accrualIntensity,
        accrualIntensityType = accrualIntensityType,
        maxNumberOfSubjects = maxNumberOfSubjects, showWarnings = FALSE
    )
    accrualSetup$.validate()

    accountForObservationTimes <- .getOptionalArgument("accountForObservationTimes", ...)
    if (is.null(accountForObservationTimes)) {
        accountForObservationTimes <- TRUE
    }

    if (!accrualSetup$maxNumberOfSubjectsCanBeCalculatedDirectly &&
            accrualSetup$followUpTimeMustBeUserDefined) {
        if (is.na(followUpTime)) {
            if (accrualSetup$piecewiseAccrualEnabled && !accrualSetup$endOfAccrualIsUserDefined) {
                stop(
                    C_EXCEPTION_TYPE_MISSING_ARGUMENT,
                    "'followUpTime', 'maxNumberOfSubjects' or end of accrual must be defined"
                )
            }

            stop(
                C_EXCEPTION_TYPE_MISSING_ARGUMENT,
                "'followUpTime' or 'maxNumberOfSubjects' must be defined"
            )
        }

        if (followUpTime == Inf) {
            followUpTime <- 1e12
        }

        if (!any(is.na(hazardRatio)) && !is.na(thetaH0)) {
            .assertIsValidHazardRatio(hazardRatio, thetaH0)
        }

        pwst <- getPiecewiseSurvivalTime(
            piecewiseSurvivalTime = piecewiseSurvivalTime,
            lambda1 = lambda1, lambda2 = lambda2,
            pi1 = pi1, pi2 = pi2,
            median1 = median1, median2 = median2,
            hazardRatio = hazardRatio, eventTime = eventTime, kappa = kappa,
            .silent = TRUE
        )
        paramName <- NULL
        if (!pwst$piecewiseSurvivalEnabled) {
            if (pwst$.getParameterType("pi1") == C_PARAM_USER_DEFINED ||
                    pwst$.getParameterType("pi1") == C_PARAM_DEFAULT_VALUE ||
                    pwst$.getParameterType("pi2") == C_PARAM_USER_DEFINED) {
                paramName <- "pi1"
            } else if (pwst$.getParameterType("lambda1") == C_PARAM_USER_DEFINED ||
                    pwst$.getParameterType("lambda2") == C_PARAM_USER_DEFINED) {
                paramName <- "lambda1"
            } else if (pwst$.getParameterType("hazardRatio") == C_PARAM_USER_DEFINED) {
                paramName <- "hazardRatio"
            } else if (pwst$.getParameterType("median1") == C_PARAM_USER_DEFINED ||
                    pwst$.getParameterType("median2") == C_PARAM_USER_DEFINED) {
                paramName <- "median1"
            }
        } else if (pwst$.getParameterType("hazardRatio") == C_PARAM_USER_DEFINED) {
            paramName <- "hazardRatio"
        }
        if (!is.null(paramName)) {
            paramValue <- pwst[[paramName]]
            if (!is.null(paramValue) && length(paramValue) > 1) {
                stop(
                    C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                    "the calculation of 'maxNumberOfSubjects' for given 'followUpTime' ",
                    "is only available for a single '", paramName, "'; ",
                    paramName, " = ", .arrayToString(
                        paramValue,
                        vectorLookAndFeelEnabled = TRUE
                    )
                )
            }
        }

        hr <- hazardRatio
        if (all(is.na(hazardRatio))) {
            hr <- pwst$hazardRatio
        }
        if (all(is.na(hazardRatio))) {
            .assertIsValidHazardRatio(hr, thetaH0)
        }

        maxNumberOfSubjectsTarget <- NA_real_
        withCallingHandlers(
            {
                # search for accrual time that provides a result
                at <- accrualSetup$accrualTime
                additionalAccrual <- 1
                searchAccrualTimeEnabled <- TRUE
                maxSearchIterations <- 50
                maxNumberOfSubjectsLower <- NA_real_
                maxNumberOfSubjectsLowerBefore <- 0
                sampleSize <- NULL
                expectionMessage <- NA_character_

                while (searchAccrualTimeEnabled && maxSearchIterations >= 0 &&
                    (is.na(maxNumberOfSubjectsLower) ||
                        maxNumberOfSubjectsLower < maxNumberOfSubjectsLowerBefore ||
                        maxNumberOfSubjectsLower < 1e8)) {
                    tryCatch(
                        {
                            maxNumberOfSubjectsLowerBefore <- ifelse(is.na(maxNumberOfSubjectsLower),
                                0, maxNumberOfSubjectsLower
                            )
                            maxNumberOfSubjectsLower <- getAccrualTime(
                                accrualTime = c(at, at[length(at)] + additionalAccrual),
                                accrualIntensity = accrualSetup$accrualIntensity,
                                accrualIntensityType = accrualIntensityType
                            )$maxNumberOfSubjects
                            additionalAccrual <- 2 * additionalAccrual
                            sampleSize <- .getSampleSizeSurvival(
                                design = design,
                                typeOfComputation = typeOfComputation,
                                thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
                                allocationRatioPlanned = allocationRatioPlanned,
                                accountForObservationTimes = accountForObservationTimes,
                                eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
                                accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
                                piecewiseSurvivalTime = piecewiseSurvivalTime,
                                lambda2 = lambda2, lambda1 = lambda1, median1 = median1, median2 = median2,
                                followUpTime = NA_real_, maxNumberOfSubjects = maxNumberOfSubjectsLower,
                                dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2, dropoutTime = dropoutTime,
                                hazardRatio = hazardRatio
                            )
                            searchAccrualTimeEnabled <- FALSE
                        },
                        error = function(e) {
                            expectionMessage <<- e$message
                        }
                    )
                    maxSearchIterations <- maxSearchIterations - 1
                }

                if (is.null(sampleSize) || is.na(sampleSize$followUpTime)) {
                    if (!is.na(expectionMessage) && grepl("'allocationRatioPlanned' > 0", expectionMessage)) {
                        stop(expectionMessage, call. = FALSE)
                    }
                    stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE,
                        "'additionalAccrual' could not be found, change accrual time specification",
                        call. = FALSE
                    )
                }

                # define lower bound for maxNumberOfSubjects
                maxNumberOfSubjectsLower <- ceiling(max(na.omit(c(
                    sampleSize$eventsFixed,
                    as.vector(sampleSize$eventsPerStage)
                ))))
                if (is.na(maxNumberOfSubjectsLower)) {
                    stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                        "'maxNumberOfSubjectsLower' could not be found",
                        call. = FALSE
                    )
                }

                # check whether accrual time already fulfills requirement
                # (followUpTime < given value) or need to be increased,
                # then define upper bound for maxNumberOfSubjects
                maxSearchIterations <- 50
                maxNumberOfSubjectsUpper <- NA_real_
                fut <- sampleSize$followUpTime
                iterations <- 1
                while (fut <= followUpTime) {
                    fut <- 2 * abs(fut)
                    iterations <- iterations + 1
                    if (iterations > 50) {
                        stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE,
                            "search algorithm failed to end",
                            call. = FALSE
                        )
                    }
                }
                while (!is.na(fut) && fut > followUpTime && maxSearchIterations >= 0) {
                    maxNumberOfSubjectsUpper <- getAccrualTime(
                        accrualTime = c(at, at[length(at)] + additionalAccrual),
                        accrualIntensity = accrualSetup$accrualIntensity,
                        accrualIntensityType = accrualIntensityType
                    )$maxNumberOfSubjects
                    additionalAccrual <- 2 * additionalAccrual
                    fut <- .getSampleSizeSurvival(
                        design = design,
                        typeOfComputation = typeOfComputation,
                        thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
                        allocationRatioPlanned = allocationRatioPlanned,
                        accountForObservationTimes = accountForObservationTimes,
                        eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
                        accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
                        piecewiseSurvivalTime = piecewiseSurvivalTime,
                        lambda2 = lambda2, lambda1 = lambda1, median1 = median1, median2 = median2,
                        followUpTime = NA_real_, maxNumberOfSubjects = maxNumberOfSubjectsUpper,
                        dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2, dropoutTime = dropoutTime,
                        hazardRatio = hazardRatio
                    )$followUpTime
                    maxSearchIterations <- maxSearchIterations - 1
                }
                if (is.na(maxNumberOfSubjectsUpper)) {
                    stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE,
                        "'maxNumberOfSubjectsUpper' could not be found ",
                        "(fut = ", fut, ", followUpTime = ", followUpTime, ")",
                        call. = FALSE
                    )
                }

                # use maxNumberOfSubjectsLower and maxNumberOfSubjectsUpper to find end of accrual
                if (dropoutRate1 != 0 || dropoutRate2 != 0) {
                    #  Adjust lower bound for given dropouts assuming exponential distribution
                    if (is.na(allocationRatioPlanned)) {
                        allocationRatioPlanned <- C_ALLOCATION_RATIO_DEFAULT
                    }

                    maxNumberOfSubjectsLower <- maxNumberOfSubjectsLower /
                        ((allocationRatioPlanned * (1 - dropoutRate1)^(
                            accrualSetup$accrualTime[length(accrualSetup$accrualTime)] / dropoutTime) +
                            (1 - dropoutRate2)^(accrualSetup$accrualTime[length(accrualSetup$accrualTime)] / dropoutTime)) /
                            (allocationRatioPlanned + 1))

                    prec <- 1
                    maxSearchIterations <- 50
                    while (prec > 1e-04 && maxSearchIterations >= 0) {
                        maxNumberOfSubjectsTarget <- (maxNumberOfSubjectsLower + maxNumberOfSubjectsUpper) / 2
                        fut <- .getSampleSizeSurvival(
                            design = design,
                            typeOfComputation = typeOfComputation,
                            thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
                            allocationRatioPlanned = allocationRatioPlanned,
                            accountForObservationTimes = accountForObservationTimes,
                            eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
                            accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
                            piecewiseSurvivalTime = piecewiseSurvivalTime,
                            lambda2 = lambda2, lambda1 = lambda1,
                            median1 = median1, median2 = median2,
                            followUpTime = NA_real_, maxNumberOfSubjects = maxNumberOfSubjectsTarget,
                            dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2,
                            dropoutTime = dropoutTime,
                            hazardRatio = hazardRatio
                        )$followUpTime
                        ifelse(fut <= followUpTime,
                            maxNumberOfSubjectsUpper <- maxNumberOfSubjectsTarget,
                            maxNumberOfSubjectsLower <- maxNumberOfSubjectsTarget
                        )
                        prec <- maxNumberOfSubjectsUpper - maxNumberOfSubjectsLower
                        maxSearchIterations <- maxSearchIterations - 1
                    }
                } else {
                    maxNumberOfSubjectsTarget <- .getOneDimensionalRootBisectionMethod(
                        fun = function(x) {
                            fut <- .getSampleSizeSurvival(
                                design = design,
                                typeOfComputation = typeOfComputation,
                                thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
                                allocationRatioPlanned = allocationRatioPlanned,
                                accountForObservationTimes = accountForObservationTimes,
                                eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
                                accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
                                piecewiseSurvivalTime = piecewiseSurvivalTime,
                                lambda2 = lambda2, lambda1 = lambda1,
                                median1 = median1, median2 = median2,
                                followUpTime = NA_real_, maxNumberOfSubjects = x,
                                dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2,
                                dropoutTime = dropoutTime,
                                hazardRatio = hazardRatio
                            )$followUpTime
                            return(followUpTime - fut)
                        },
                        lower = maxNumberOfSubjectsLower,
                        upper = maxNumberOfSubjectsUpper,
                        tolerance = 1e-04,
                        acceptResultsOutOfTolerance = TRUE,
                        maxSearchIterations = 50,
                        direction = 0,
                        suppressWarnings = FALSE,
                        callingFunctionInformation = "getSampleSizeSurvival"
                    )
                }
            },
            warning = function(w) {
                invokeRestart("muffleWarning")
            }
        )

        if (is.na(maxNumberOfSubjectsTarget)) {
            stop(
                C_EXCEPTION_TYPE_RUNTIME_ISSUE,
                "failed to calculate 'maxNumberOfSubjects' by given 'followUpTime' ",
                "(lower = ", maxNumberOfSubjectsLower, ", upper = ", maxNumberOfSubjectsUpper, ")"
            )
        }

        sampleSizeSurvival <- .getSampleSizeSurvival(
            design = design,
            typeOfComputation = typeOfComputation,
            thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
            allocationRatioPlanned = allocationRatioPlanned,
            accountForObservationTimes = accountForObservationTimes,
            eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
            accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
            piecewiseSurvivalTime = piecewiseSurvivalTime, lambda2 = lambda2, lambda1 = lambda1,
            median1 = median1, median2 = median2,
            followUpTime = NA_real_, maxNumberOfSubjects = maxNumberOfSubjectsTarget,
            dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2,
            dropoutTime = dropoutTime,
            hazardRatio = hazardRatio
        )
        sampleSizeSurvival$.setParameterType("followUpTime", C_PARAM_USER_DEFINED)
        sampleSizeSurvival$.accrualTime <- accrualSetup

        if (!is.na(sampleSizeSurvival$followUpTime)) {
            if (followUpTime == 1e12) {
                followUpTime <- Inf
            }

            if (sampleSizeSurvival$followUpTime >= -1e-02 && sampleSizeSurvival$followUpTime <= 1e-02) {
                sampleSizeSurvival$followUpTime <- 0
            }

            if (sampleSizeSurvival$followUpTime < followUpTime - 1e-02 ||
                    sampleSizeSurvival$followUpTime > followUpTime + 1e-02) {
                sampleSizeSurvival$.setParameterType("followUpTime", C_PARAM_GENERATED)
                warning("User defined 'followUpTime' (", followUpTime, ") ignored because ",
                    "follow-up time is ", round(sampleSizeSurvival$followUpTime, 4),
                    call. = FALSE
                )
            }
        }

        sampleSizeSurvival$.setParameterType("maxNumberOfSubjects", C_PARAM_GENERATED)
        sampleSizeSurvival$.setParameterType("accrualTime", C_PARAM_GENERATED)

        return(sampleSizeSurvival)
    }

    return(.getSampleSizeSurvival(
        design = design,
        typeOfComputation = typeOfComputation,
        thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
        allocationRatioPlanned = allocationRatioPlanned,
        accountForObservationTimes = accountForObservationTimes,
        eventTime = eventTime, accrualTime = accrualTime,
        accrualIntensity = accrualIntensity, accrualIntensityType = accrualIntensityType,
        kappa = kappa, piecewiseSurvivalTime = piecewiseSurvivalTime,
        lambda2 = lambda2, lambda1 = lambda1,
        median1 = median1, median2 = median2,
        followUpTime = followUpTime, maxNumberOfSubjects = maxNumberOfSubjects,
        dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2,
        dropoutTime = dropoutTime,
        hazardRatio = hazardRatio
    ))
}

.getSampleSizeSurvival <- function(...,
        design = NULL,
        typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
        thetaH0 = 1,
        pi2 = NA_real_,
        pi1 = NA_real_,
        allocationRatioPlanned = NA_real_,
        accountForObservationTimes = TRUE,
        eventTime = C_EVENT_TIME_DEFAULT,
        accrualTime = C_ACCRUAL_TIME_DEFAULT,
        accrualIntensity = C_ACCRUAL_INTENSITY_DEFAULT,
        accrualIntensityType = c("auto", "absolute", "relative"),
        kappa = 1,
        piecewiseSurvivalTime = NA_real_,
        lambda2 = NA_real_,
        lambda1 = NA_real_,
        median1 = NA_real_,
        median2 = NA_real_,
        followUpTime = NA_real_,
        maxNumberOfSubjects = NA_real_,
        dropoutRate1 = 0,
        dropoutRate2 = dropoutRate1,
        dropoutTime = NA_real_,
        hazardRatio = NA_real_) {
    designPlan <- .createDesignPlanSurvival(
        objectType = "sampleSize",
        design = design,
        typeOfComputation = typeOfComputation,
        thetaH0 = thetaH0,
        pi2 = pi2,
        pi1 = pi1,
        allocationRatioPlanned = allocationRatioPlanned,
        accountForObservationTimes = accountForObservationTimes,
        eventTime = eventTime,
        accrualTime = accrualTime,
        accrualIntensity = accrualIntensity,
        accrualIntensityType = accrualIntensityType,
        kappa = kappa,
        piecewiseSurvivalTime = piecewiseSurvivalTime,
        lambda2 = lambda2,
        lambda1 = lambda1,
        median1 = median1,
        median2 = median2,
        followUpTime = followUpTime,
        maxNumberOfSubjects = maxNumberOfSubjects,
        dropoutRate1 = dropoutRate1,
        dropoutRate2 = dropoutRate2,
        dropoutTime = dropoutTime,
        hazardRatio = hazardRatio
    )
    return(.getSampleSize(designPlan))
}

.createDesignPlanSurvival <- function(..., objectType = c("sampleSize", "power"),
        design,
        typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
        thetaH0, pi2, pi1,
        allocationRatioPlanned,
        accountForObservationTimes,
        eventTime,
        accrualTime,
        accrualIntensity,
        accrualIntensityType,
        kappa,
        piecewiseSurvivalTime,
        lambda2,
        lambda1,
        median1,
        median2,
        followUpTime = NA_real_,
        directionUpper = NA,
        maxNumberOfEvents = NA_real_,
        maxNumberOfSubjects,
        dropoutRate1,
        dropoutRate2,
        dropoutTime,
        hazardRatio) {
    objectType <- match.arg(objectType)
    typeOfComputation <- .matchArgument(typeOfComputation, "Schoenfeld")

    .assertIsTrialDesignInverseNormalOrGroupSequential(design)
    .assertIsValidAlphaAndBeta(design$alpha, design$beta)
    .assertIsValidSidedParameter(design$sided)
    .assertIsSingleLogical(accountForObservationTimes, "accountForObservationTimes", naAllowed = TRUE)
    .assertIsSingleNumber(thetaH0, "thetaH0")
    .assertIsValidThetaH0(thetaH0, endpoint = "survival", groups = 2)
    .assertIsValidKappa(kappa)
    .assertIsValidDirectionUpper(directionUpper, design$sided, objectType, userFunctionCallEnabled = TRUE)

    if (objectType == "power") {
        .assertIsSingleNumber(maxNumberOfEvents, "maxNumberOfEvents")
        .assertIsInClosedInterval(maxNumberOfEvents, "maxNumberOfEvents",
            lower = 1, upper = maxNumberOfSubjects
        )
    }

    if (!any(is.na(pi1)) && (any(pi1 <= 0) || any(pi1 >= 1))) {
        stop(
            C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
            "event rate 'pi1' (", .arrayToString(pi1), ") is out of bounds (0; 1)"
        )
    }

    if (!any(is.na(pi2)) && (any(pi2 <= 0) || any(pi2 >= 1))) {
        stop(
            C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
            "event rate 'pi2' (", .arrayToString(pi2), ") is out of bounds (0; 1)"
        )
    }

    if (design$sided == 2 && thetaH0 != 1) {
        stop(
            C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
            "two-sided case is implemented for superiority testing only (i.e., thetaH0 = 1)"
        )
    }

    if (thetaH0 <= 0) {
        stop(
            C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
            "null hypothesis hazard ratio is not allowed be negative or zero"
        )
    }

    if (!(typeOfComputation %in% c("Schoenfeld", "Freedman", "HsiehFreedman"))) {
        stop(
            C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
            "computation type ('", typeOfComputation, "') must be one of the following: ",
            "'Schoenfeld', 'Freedman', or 'HsiehFreedman' "
        )
    }

    if (typeOfComputation != "Schoenfeld" && thetaH0 != 1) {
        stop(
            C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
            "Freedman test calculation is possible only for superiority testing (thetaH0 != 1)"
        )
    }

    if (is.numeric(accrualTime) && all(is.na(accrualTime))) {
        accrualTime <- C_ACCRUAL_TIME_DEFAULT
    }
    if (all(is.na(accrualIntensity))) {
        accrualIntensity <- C_ACCRUAL_INTENSITY_DEFAULT
    }

    accrualSetup <- getAccrualTime(
        accrualTime = accrualTime,
        accrualIntensity = accrualIntensity,
        accrualIntensityType = accrualIntensityType,
        maxNumberOfSubjects = maxNumberOfSubjects
    )
    accrualSetup$.validate()

    if (is.na(allocationRatioPlanned)) {
        allocationRatioPlanned <- C_ALLOCATION_RATIO_DEFAULT
    }

    .assertIsValidAllocationRatioPlannedSampleSize(allocationRatioPlanned, maxNumberOfSubjects)

    designPlan <- TrialDesignPlanSurvival(
        design = design,
        typeOfComputation = typeOfComputation,
        thetaH0 = thetaH0,
        allocationRatioPlanned = allocationRatioPlanned,
        accountForObservationTimes = accountForObservationTimes,
        eventTime = eventTime,
        accrualTime = accrualSetup$.getAccrualTimeWithoutLeadingZero(),
        accrualIntensity = accrualSetup$accrualIntensity,
        kappa = kappa,
        followUpTime = followUpTime,
        maxNumberOfSubjects = maxNumberOfSubjects,
        dropoutRate1 = dropoutRate1,
        dropoutRate2 = dropoutRate2,
        dropoutTime = dropoutTime,
        hazardRatio = hazardRatio
    )

    .setValueAndParameterType(
        designPlan, "allocationRatioPlanned",
        allocationRatioPlanned, C_ALLOCATION_RATIO_DEFAULT
    )
    .setValueAndParameterType(designPlan, "dropoutRate1", dropoutRate1, C_DROP_OUT_RATE_1_DEFAULT)
    .setValueAndParameterType(designPlan, "dropoutRate2", dropoutRate2, C_DROP_OUT_RATE_2_DEFAULT)
    .setValueAndParameterType(designPlan, "dropoutTime", dropoutTime, C_DROP_OUT_TIME_DEFAULT)
    .setValueAndParameterType(designPlan, "kappa", kappa, 1)

    designPlan$.setSampleSizeObject(objectType)

    designPlan$criticalValuesPValueScale <- matrix(design$stageLevels, ncol = 1)
    if (design$sided == 2) {
        designPlan$criticalValuesPValueScale <- designPlan$criticalValuesPValueScale * 2
        designPlan$.setParameterType("criticalValuesPValueScale", C_PARAM_GENERATED)
    }

    if (any(design$futilityBounds > C_FUTILITY_BOUNDS_DEFAULT)) {
        designPlan$futilityBoundsPValueScale <- matrix(1 - stats::pnorm(design$futilityBounds), ncol = 1)
        designPlan$.setParameterType("futilityBoundsPValueScale", C_PARAM_GENERATED)
    }

    designPlan$.accrualTime <- accrualSetup

    designPlan$totalAccrualTime <- accrualSetup$accrualTime[length(accrualSetup$accrualTime)]
    if (length(accrualSetup$accrualTime) > 2) {
        designPlan$.setParameterType("totalAccrualTime", C_PARAM_GENERATED)
    } else {
        designPlan$.setParameterType("totalAccrualTime", C_PARAM_NOT_APPLICABLE)
    }

    if (is.na(maxNumberOfSubjects)) {
        if (!is.na(accrualSetup$maxNumberOfSubjects)) {
            designPlan$maxNumberOfSubjects <- accrualSetup$maxNumberOfSubjects
            designPlan$.setParameterType(
                "maxNumberOfSubjects",
                accrualSetup$.getParameterType("maxNumberOfSubjects")
            )
        }
    } else if (maxNumberOfSubjects == 0) {
        designPlan$.setParameterType("maxNumberOfSubjects", C_PARAM_GENERATED)
    } else {
        designPlan$.setParameterType("maxNumberOfSubjects", C_PARAM_USER_DEFINED)
    }

    if (identical(as.integer(accrualSetup$accrualTime), C_ACCRUAL_TIME_DEFAULT) ||
            identical(
                as.integer(c(0L, accrualSetup$.getAccrualTimeWithoutLeadingZero())),
                C_ACCRUAL_TIME_DEFAULT
            )) {
        designPlan$.setParameterType("accrualTime", C_PARAM_DEFAULT_VALUE)
    } else {
        designPlan$.setParameterType("accrualTime", accrualSetup$.getParameterType("accrualTime"))
    }

    if (length(designPlan$accrualIntensity) == 1 &&
            designPlan$accrualIntensity == C_ACCRUAL_INTENSITY_DEFAULT) {
        designPlan$.setParameterType("accrualIntensity", C_PARAM_DEFAULT_VALUE)
    } else {
        designPlan$.setParameterType(
            "accrualIntensity",
            accrualSetup$.getParameterType("accrualIntensity")
        )
    }

    .assertIsSingleNumber(designPlan$eventTime, "eventTime")
    .assertIsSingleNumber(designPlan$allocationRatioPlanned, "allocationRatioPlanned")
    .assertIsSingleNumber(designPlan$kappa, "kappa")
    if (objectType == "power") {
        .assertIsValidMaxNumberOfSubjects(designPlan$maxNumberOfSubjects)
    }
    .assertIsSingleNumber(designPlan$dropoutRate1, "dropoutRate1")
    .assertIsSingleNumber(designPlan$dropoutRate2, "dropoutRate2")
    .assertIsSingleNumber(designPlan$dropoutTime, "dropoutTime")

    if (objectType == "power") {
        pi1Default <- C_PI_1_DEFAULT
    } else {
        pi1Default <- C_PI_1_SAMPLE_SIZE_DEFAULT
    }
    designPlan$.piecewiseSurvivalTime <- getPiecewiseSurvivalTime(
        piecewiseSurvivalTime = piecewiseSurvivalTime, lambda2 = lambda2, lambda1 = lambda1,
        median1 = median1, median2 = median2,
        hazardRatio = hazardRatio, pi1 = pi1, pi2 = pi2, eventTime = eventTime, kappa = kappa,
        .pi1Default = pi1Default
    )
    designPlan$.setParameterType("kappa", designPlan$.piecewiseSurvivalTime$.getParameterType("kappa"))

    if (designPlan$.piecewiseSurvivalTime$.getParameterType("pi1") == C_PARAM_DEFAULT_VALUE &&
            length(designPlan$.piecewiseSurvivalTime$pi1) > 1 &&
            length(accrualSetup$accrualIntensity) > 1 && all(accrualSetup$accrualIntensity < 1)) {
        designPlan$.piecewiseSurvivalTime$pi1 <- designPlan$.piecewiseSurvivalTime$pi1[1]
        warning("Only the first default 'pi1' (", designPlan$.piecewiseSurvivalTime$pi1, ") was used ",
            "because the accrual intensities (", .arrayToString(accrualSetup$accrualIntensity), ") ",
            "were defined relative (all accrual intensities are < 1)",
            call. = FALSE
        )
    }

    .initDesignPlanSurvival(designPlan)

    designPlan$.setParameterType("followUpTime", C_PARAM_NOT_APPLICABLE)
    if (designPlan$accountForObservationTimes) {
        .assertIsSingleNumber(dropoutRate1, "dropoutRate1")
        .assertIsSingleNumber(dropoutRate2, "dropoutRate2")
        .assertIsSingleNumber(dropoutTime, "dropoutTime")

        if (!is.na(dropoutTime) && dropoutTime <= 0) {
            stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'dropoutTime' (", dropoutTime, ") must be > 0")
        }

        if (dropoutRate1 < 0 || dropoutRate1 >= 1) {
            stop(
                C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
                "'dropoutRate1' (", dropoutRate1, ") is out of bounds [0; 1)"
            )
        }

        if (dropoutRate2 < 0 || dropoutRate2 >= 1) {
            stop(
                C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
                "'dropoutRate2' (", dropoutRate2, ") is out of bounds [0; 1)"
            )
        }

        if (!is.na(eventTime) && eventTime <= 0) {
            stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'eventTime' (", eventTime, ") must be > 0")
        }

        .assertIsValidAccrualTime(accrualSetup$.getAccrualTimeWithoutLeadingZero())

        .assertIsValidFollowUpTime(followUpTime)
        .setValueAndParameterType(designPlan, "followUpTime", followUpTime, C_FOLLOW_UP_TIME_DEFAULT)

        if (.isUserDefinedMaxNumberOfSubjects(designPlan) && !is.null(followUpTime) &&
                length(followUpTime) == 1 && !is.na(followUpTime)) {
            warning("Follow-up time will be calculated, value entered (",
                followUpTime, ") is not taken into account",
                call. = FALSE
            )
        } else if (is.na(followUpTime)) {
            designPlan$followUpTime <- C_FOLLOW_UP_TIME_DEFAULT
            designPlan$.setParameterType("followUpTime", C_PARAM_DEFAULT_VALUE)
        }

        if (objectType == "power") {
            designPlan$followUpTime <- NA_real_
            designPlan$.setParameterType("followUpTime", C_PARAM_NOT_APPLICABLE)
        }
    } else {
        for (p in c(
            "accrualTime", "accrualIntensity",
            "eventTime", "dropoutRate1", "dropoutRate2", "dropoutTime", "followUpTime",
            "analysisTime", "studyDuration"
        )) {
            designPlan$.setParameterType(p, C_PARAM_NOT_APPLICABLE)
        }
        if (designPlan$.getParameterType("accrualTime") == C_PARAM_USER_DEFINED ||
                !identical(accrualTime, C_ACCRUAL_TIME_DEFAULT)) {
            designPlan$.warnInCaseArgumentExists(accrualSetup$accrualTime, "accrualTime")
        }
        designPlan$.warnInCaseArgumentExists(dropoutRate1, "dropoutRate1")
        designPlan$.warnInCaseArgumentExists(dropoutRate2, "dropoutRate2")
        if (!identical(dropoutTime, C_DROP_OUT_TIME_DEFAULT)) {
            designPlan$.warnInCaseArgumentExists(dropoutTime, "dropoutTime")
        }
        designPlan$.warnInCaseArgumentExists(maxNumberOfSubjects, "maxNumberOfSubjects")
        if (!identical(followUpTime, C_FOLLOW_UP_TIME_DEFAULT)) {
            designPlan$.warnInCaseArgumentExists(followUpTime, "followUpTime")
        }
    }

    .setValueAndParameterType(designPlan, "directionUpper", directionUpper, TRUE)
    if (objectType == "power") {
        .setValueAndParameterType(designPlan, "maxNumberOfEvents", maxNumberOfEvents, NA_real_)
        designPlan$.setParameterType("accountForObservationTimes", C_PARAM_NOT_APPLICABLE)
    }

    return(designPlan)
}

.isUserDefinedMaxNumberOfSubjects <- function(designPlan) {
    if (!is.null(designPlan) && length(designPlan$.getParameterType("maxNumberOfSubjects")) > 0) {
        if (designPlan$.getParameterType("maxNumberOfSubjects") == C_PARAM_USER_DEFINED) {
            return(TRUE)
        }
    }

    return(!is.null(designPlan$maxNumberOfSubjects) && length(designPlan$maxNumberOfSubjects) == 1 &&
        !is.na(designPlan$maxNumberOfSubjects) && designPlan$maxNumberOfSubjects > 0)
}

.initDesignPlanSurvivalByPiecewiseSurvivalTimeObject <- function(designPlan, pwstSetup) {
    designPlan$pi1 <- pwstSetup$pi1
    designPlan$.setParameterType("pi1", pwstSetup$.getParameterType("pi1"))

    designPlan$pi2 <- pwstSetup$pi2
    designPlan$.setParameterType("pi2", pwstSetup$.getParameterType("pi2"))

    designPlan$hazardRatio <- pwstSetup$hazardRatio
    designPlan$.setParameterType("hazardRatio", pwstSetup$.getParameterType("hazardRatio"))

    designPlan$lambda1 <- pwstSetup$lambda1
    designPlan$.setParameterType("lambda1", pwstSetup$.getParameterType("lambda1"))

    designPlan$lambda2 <- pwstSetup$lambda2
    designPlan$.setParameterType("lambda2", pwstSetup$.getParameterType("lambda2"))

    designPlan$median1 <- pwstSetup$median1
    designPlan$.setParameterType("median1", pwstSetup$.getParameterType("median1"))

    designPlan$median2 <- pwstSetup$median2
    designPlan$.setParameterType("median2", pwstSetup$.getParameterType("median2"))

    designPlan$piecewiseSurvivalTime <- pwstSetup$piecewiseSurvivalTime
    designPlan$.setParameterType(
        "piecewiseSurvivalTime",
        pwstSetup$.getParameterType("piecewiseSurvivalTime")
    )

    designPlan$eventTime <- pwstSetup$eventTime
    designPlan$.setParameterType("eventTime", pwstSetup$.getParameterType("eventTime"))

    if (pwstSetup$.isLambdaBased()) {
        return(length(designPlan$hazardRatio))
    }

    return(length(designPlan$pi1))
}

.initDesignPlanSurvival <- function(designPlan) {
    numberOfResults <- .initDesignPlanSurvivalByPiecewiseSurvivalTimeObject(
        designPlan, designPlan$.piecewiseSurvivalTime
    )

    if (designPlan$.piecewiseSurvivalTime$.isLambdaBased()) {
        if (length(designPlan$accountForObservationTimes) == 0 ||
                is.na(designPlan$accountForObservationTimes) ||
                !designPlan$accountForObservationTimes) {
            designPlan$accountForObservationTimes <- TRUE
            designPlan$.setParameterType("accountForObservationTimes", C_PARAM_DEFAULT_VALUE)
        }

        if (!designPlan$accountForObservationTimes) {
            designPlan$accountForObservationTimes <- TRUE
            warning("'accountForObservationTimes' was set to TRUE ",
                "because piecewise exponential survival function is enabled",
                call. = FALSE
            )
        }
    } else {
        if (.isUserDefinedMaxNumberOfSubjects(designPlan)) {
            if (length(designPlan$accountForObservationTimes) != 0 &&
                    !is.na(designPlan$accountForObservationTimes) &&
                    !designPlan$accountForObservationTimes) {
                stop(
                    C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                    "'accountForObservationTimes' must be TRUE because 'maxNumberOfSubjects' is > 0"
                )
            }

            designPlan$.setParameterType("accountForObservationTimes", C_PARAM_GENERATED)
            designPlan$accountForObservationTimes <- TRUE
        } else {
            if (length(designPlan$accountForObservationTimes) == 0 ||
                    is.na(designPlan$accountForObservationTimes)) {
                designPlan$accountForObservationTimes <- FALSE
                designPlan$.setParameterType("accountForObservationTimes", C_PARAM_DEFAULT_VALUE)
            } else {
                designPlan$.setParameterType(
                    "accountForObservationTimes",
                    ifelse(designPlan$accountForObservationTimes,
                        C_PARAM_DEFAULT_VALUE, C_PARAM_USER_DEFINED
                    )
                )
            }
        }
    }

    designPlan$.setParameterType("chi", C_PARAM_NOT_APPLICABLE)

    if (designPlan$.isSampleSizeObject()) {
        designPlan$.setParameterType("directionUpper", C_PARAM_NOT_APPLICABLE)
        designPlan$.setParameterType("maxNumberOfEvents", C_PARAM_NOT_APPLICABLE)
    }

    return(numberOfResults = numberOfResults)
}

.warnInCaseOfDefinedPiValue <- function(designPlan, piValueName) {
    piValue <- designPlan[[piValueName]]
    if (!is.null(piValue) && !is.na(piValue) && length(piValue) > 0) {
        designPlan$.setParameterType(piValueName, C_PARAM_NOT_APPLICABLE)
        warning("'pi2' (", .arrayToString(piValue), ") will be ignored ",
            "because piecewise exponential survival function is enabled",
            call. = FALSE
        )
        designPlan[[piValueName]] <- NA_real_
    }
}

.getSampleSize <- function(designPlan) {
    if (.isTrialDesignPlanMeans(designPlan) || .isTrialDesignPlanRates(designPlan)) {
        if (identical(designPlan$allocationRatioPlanned, 0)) {
            designPlan$optimumAllocationRatio <- TRUE
            designPlan$.setParameterType("optimumAllocationRatio", C_PARAM_USER_DEFINED)
        }

        if (.isTrialDesignPlanMeans(designPlan)) {
            sampleSizeFixed <- .getSampleSizeFixedMeans(
                alpha = designPlan$getAlpha(),
                beta = designPlan$getBeta(),
                sided = designPlan$getSided(),
                twoSidedPower = designPlan$getTwoSidedPower(),
                normalApproximation = designPlan$normalApproximation,
                meanRatio = designPlan$meanRatio,
                thetaH0 = designPlan$thetaH0,
                alternative = designPlan$alternative,
                stDev = designPlan$stDev,
                groups = designPlan$groups,
                allocationRatioPlanned = designPlan$allocationRatioPlanned
            )
        } else {
            sampleSizeFixed <- .getSampleSizeFixedRates(
                alpha = designPlan$getAlpha(),
                beta = designPlan$getBeta(),
                sided = designPlan$getSided(),
                normalApproximation = designPlan$normalApproximation,
                riskRatio = designPlan$riskRatio,
                thetaH0 = designPlan$thetaH0,
                pi1 = designPlan$pi1,
                pi2 = designPlan$pi2,
                groups = designPlan$groups,
                allocationRatioPlanned = designPlan$allocationRatioPlanned
            )
        }

        # Fixed
        designPlan$nFixed <- sampleSizeFixed$nFixed
        designPlan$.setParameterType("nFixed", C_PARAM_GENERATED)
        if (designPlan$groups == 2) {
            designPlan$nFixed1 <- sampleSizeFixed$n1Fixed
            designPlan$nFixed2 <- sampleSizeFixed$n2Fixed
            designPlan$.setParameterType("nFixed1", C_PARAM_GENERATED)
            designPlan$.setParameterType("nFixed2", C_PARAM_GENERATED)
            designPlan$numberOfSubjects1 <- matrix(designPlan$nFixed1, nrow = 1)
            designPlan$numberOfSubjects2 <- matrix(designPlan$nFixed2, nrow = 1)
        }
        designPlan$numberOfSubjects <- matrix(designPlan$nFixed, nrow = 1)

        if (!is.null(sampleSizeFixed$allocationRatioPlanned) &&
                (length(designPlan$allocationRatioPlanned) !=
                    length(sampleSizeFixed$allocationRatioPlanned) ||
                    sum(designPlan$allocationRatioPlanned == sampleSizeFixed$allocationRatioPlanned) !=
                        length(designPlan$allocationRatioPlanned))) {
            designPlan$allocationRatioPlanned <- sampleSizeFixed$allocationRatioPlanned
            designPlan$.setParameterType("allocationRatioPlanned", C_PARAM_GENERATED)
        }

        # Sequential
        if (designPlan$.design$kMax > 1) {
            designCharacteristics <- getDesignCharacteristics(designPlan$.design)
            if (.isTrialDesignPlanMeans(designPlan)) {
                sampleSizeSequential <- .getSampleSizeSequentialMeans(
                    sampleSizeFixed, designCharacteristics
                )
            } else {
                sampleSizeSequential <- .getSampleSizeSequentialRates(
                    sampleSizeFixed, designCharacteristics
                )
            }

            designPlan$informationRates <- sampleSizeSequential$informationRates
            if (ncol(designPlan$informationRates) == 1 &&
                    identical(designPlan$informationRates[, 1], designPlan$.design$informationRates)) {
                designPlan$.setParameterType("informationRates", C_PARAM_NOT_APPLICABLE)
            } else {
                designPlan$.setParameterType("informationRates", C_PARAM_GENERATED)
            }

            designPlan$maxNumberOfSubjects <- sampleSizeSequential$maxNumberOfSubjects
            designPlan$.setParameterType("maxNumberOfSubjects", C_PARAM_GENERATED)
            if (designPlan$groups == 2) {
                designPlan$maxNumberOfSubjects1 <- .getNumberOfSubjects1(
                    designPlan$maxNumberOfSubjects, designPlan$allocationRatioPlanned
                )
                designPlan$maxNumberOfSubjects2 <- .getNumberOfSubjects2(
                    designPlan$maxNumberOfSubjects, designPlan$allocationRatioPlanned
                )
                designPlan$.setParameterType("maxNumberOfSubjects1", C_PARAM_GENERATED)
                designPlan$.setParameterType("maxNumberOfSubjects2", C_PARAM_GENERATED)
            }

            designPlan$numberOfSubjects <- sampleSizeSequential$numberOfSubjects
            designPlan$.setParameterType("numberOfSubjects", C_PARAM_GENERATED)

            if (designPlan$groups == 2) {
                designPlan$numberOfSubjects1 <- sampleSizeSequential$numberOfSubjects1
                designPlan$numberOfSubjects2 <- sampleSizeSequential$numberOfSubjects2
                designPlan$.setParameterType("numberOfSubjects1", C_PARAM_GENERATED)
                designPlan$.setParameterType("numberOfSubjects2", C_PARAM_GENERATED)
            }

            designPlan$expectedNumberOfSubjectsH0 <- sampleSizeSequential$expectedNumberOfSubjectsH0
            designPlan$expectedNumberOfSubjectsH01 <- sampleSizeSequential$expectedNumberOfSubjectsH01
            designPlan$expectedNumberOfSubjectsH1 <- sampleSizeSequential$expectedNumberOfSubjectsH1
            designPlan$.setParameterType("expectedNumberOfSubjectsH0", C_PARAM_GENERATED)
            designPlan$.setParameterType("expectedNumberOfSubjectsH01", C_PARAM_GENERATED)
            designPlan$.setParameterType("expectedNumberOfSubjectsH1", C_PARAM_GENERATED)

            designPlan$.setParameterType("eventsFixed", C_PARAM_NOT_APPLICABLE)
            designPlan$.setParameterType("nFixed1", C_PARAM_NOT_APPLICABLE)
            designPlan$.setParameterType("nFixed2", C_PARAM_NOT_APPLICABLE)
            designPlan$.setParameterType("nFixed", C_PARAM_NOT_APPLICABLE)

            if (designPlan$allocationRatioPlanned[1] == 1) {
                designPlan$.setParameterType("numberOfSubjects1", C_PARAM_NOT_APPLICABLE)
                designPlan$.setParameterType("numberOfSubjects2", C_PARAM_NOT_APPLICABLE)
            }

            if (!is.null(sampleSizeSequential$rejectPerStage)) {
                designPlan$rejectPerStage <- matrix(sampleSizeSequential$rejectPerStage,
                    nrow = designPlan$.design$kMax
                )
                designPlan$.setParameterType("rejectPerStage", C_PARAM_GENERATED)

                designPlan$earlyStop <- sum(designPlan$rejectPerStage[1:(designPlan$.design$kMax - 1), ])
                designPlan$.setParameterType("earlyStop", C_PARAM_GENERATED)
            }
            if (!is.null(sampleSizeSequential$futilityPerStage) &&
                    any(designPlan$.design$futilityBounds != C_FUTILITY_BOUNDS_DEFAULT)) {
                designPlan$futilityPerStage <- matrix(sampleSizeSequential$futilityPerStage,
                    nrow = designPlan$.design$kMax - 1
                )
                designPlan$.setParameterType("futilityPerStage", C_PARAM_GENERATED)

                designPlan$futilityStop <- sum(designPlan$futilityPerStage)
                designPlan$.setParameterType("futilityStop", C_PARAM_GENERATED)

                designPlan$earlyStop <- designPlan$earlyStop + sum(designPlan$futilityPerStage)
            }
        }

        .addEffectScaleBoundaryDataToDesignPlan(designPlan)

        return(designPlan)
    } else if (.isTrialDesignPlanSurvival(designPlan)) {
        # Fixed
        designPlan <- .getSampleSizeFixedSurvival(designPlan)

        # Sequential
        if (designPlan$.design$kMax > 1) {
            designCharacteristics <- getDesignCharacteristics(designPlan$.design)
            designPlan <- .getSampleSizeSequentialSurvival(designPlan, designCharacteristics)
        }

        if (designPlan$accountForObservationTimes && !any(is.na(designPlan$followUpTime)) &&
                all(designPlan$followUpTime == C_FOLLOW_UP_TIME_DEFAULT)) {
            designPlan$.setParameterType("followUpTime", C_PARAM_DEFAULT_VALUE)
        }

        .addEffectScaleBoundaryDataToDesignPlan(designPlan)

        if (designPlan$.getParameterType("maxNumberOfSubjects") == C_PARAM_GENERATED &&
                designPlan$.accrualTime$.getParameterType("maxNumberOfSubjects") != C_PARAM_GENERATED &&
                all(designPlan$accrualIntensity < 1)) {
            numberOfDefinedAccrualIntensities <- length(designPlan$accrualIntensity)

            accrualTime <- designPlan$accrualTime
            if (length(accrualTime) > 0 && accrualTime[1] != 0) {
                accrualTime <- c(0, accrualTime)
            }

            if (any(designPlan$accrualIntensity < 0)) {
                stop(
                    C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                    "'accrualIntensityRelative' (",
                    .arrayToString(designPlan$accrualIntensity), ") must be >= 0"
                )
            }

            designPlan$accrualIntensityRelative <- designPlan$accrualIntensity
            if (identical(designPlan$accrualIntensityRelative, C_ACCRUAL_INTENSITY_DEFAULT)) {
                designPlan$.setParameterType("accrualIntensityRelative", C_PARAM_NOT_APPLICABLE)
            } else {
                designPlan$.setParameterType(
                    "accrualIntensityRelative",
                    designPlan$.getParameterType("accrualIntensity")
                )
            }

            accrualIntensityAbsolute <- c()
            for (maxNumberOfSubjects in designPlan$maxNumberOfSubjects) {
                accrualSetup <- getAccrualTime(
                    accrualTime = accrualTime,
                    accrualIntensity = designPlan$accrualIntensityRelative,
                    accrualIntensityType = "relative",
                    maxNumberOfSubjects = maxNumberOfSubjects
                )
                accrualIntensityAbsolute <- c(accrualIntensityAbsolute, accrualSetup$accrualIntensity)
            }
            designPlan$accrualIntensity <- accrualIntensityAbsolute
            designPlan$.setParameterType("accrualIntensity", C_PARAM_GENERATED)

            if (numberOfDefinedAccrualIntensities > 1) {
                paramName <- NULL
                if (designPlan$.getParameterType("pi1") == C_PARAM_USER_DEFINED ||
                        designPlan$.getParameterType("pi1") == C_PARAM_DEFAULT_VALUE ||
                        designPlan$.getParameterType("pi2") == C_PARAM_USER_DEFINED) {
                    paramName <- "pi1"
                } else if (designPlan$.getParameterType("median1") == C_PARAM_USER_DEFINED ||
                        designPlan$.getParameterType("median2") == C_PARAM_USER_DEFINED) {
                    paramName <- "median1"
                }
                if (!is.null(paramName)) {
                    paramValue <- designPlan[[paramName]]
                    if (!is.null(paramValue) && length(paramValue) > 1) {
                        stop(
                            C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                            "the definition of relative accrual intensities ",
                            "(all 'accrualIntensity' values < 1) ",
                            "is only available for a single value ",
                            "(", paramName, " = ", .arrayToString(
                                paramValue,
                                vectorLookAndFeelEnabled = TRUE
                            ), ")"
                        )
                    }
                }
            }
        }

        designPlan$maxNumberOfEvents <- designPlan$eventsPerStage[designPlan$.design$kMax, ]
        designPlan$.setParameterType("maxNumberOfEvents", C_PARAM_GENERATED)

        if (!any(is.na(designPlan$followUpTime))) {
            if (any(designPlan$followUpTime < -1e-02)) {
                warning("Accrual duration longer than maximal study ",
                    "duration (time to maximal number of events); followUpTime = ",
                    .arrayToString(designPlan$followUpTime),
                    call. = FALSE
                )
            }
        } else {
            warning("Follow-up time could not be calculated for hazardRatio = ",
                .arrayToString(designPlan$hazardRatio[indices]),
                call. = FALSE
            )
        }

        if (designPlan$.getParameterType("accountForObservationTimes") != C_PARAM_USER_DEFINED) {
            designPlan$.setParameterType("accountForObservationTimes", C_PARAM_NOT_APPLICABLE)
        }
        designPlan$.setParameterType("chi", C_PARAM_NOT_APPLICABLE)

        .addStudyDurationToDesignPlan(designPlan)

        return(designPlan)
    }

    stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE, "unknown trial plan class '", .getClassName(designPlan), "'")
}

.getSampleSizeFixedMeans <- function(..., alpha = 0.025, beta = 0.2, sided = 1,
        twoSidedPower = C_TWO_SIDED_POWER_DEFAULT,
        normalApproximation = FALSE, meanRatio = FALSE,
        thetaH0 = 0, alternative = C_ALTERNATIVE_DEFAULT,
        stDev = C_STDEV_DEFAULT, groups = 2, allocationRatioPlanned = C_ALLOCATION_RATIO_DEFAULT) {
    nFixed <- rep(NA_real_, length(alternative))

    for (i in 1:length(alternative)) {
        theta <- alternative[i]

        if (groups == 1) {
            if (sided == 1 || !twoSidedPower) {
                if (normalApproximation == FALSE) {
                    up <- 2
                    while (stats::pt(
                        stats::qt(1 - alpha / sided, up - 1), max(0.001, up - 1),
                        sqrt(up) * abs(theta - thetaH0) / stDev
                    ) > beta) {
                        up <- 2 * up
                    }
                    nFixed[i] <- .getOneDimensionalRoot(
                        function(n) {
                            return(stats::pt(
                                stats::qt(1 - alpha / sided, max(0.001, n - 1)),
                                max(0.001, n - 1), sqrt(n) * abs(theta - thetaH0) / stDev
                            ) - beta)
                        },
                        lower = 0.001, upper = up, tolerance = 1e-04,
                        callingFunctionInformation = ".getSampleSizeFixedMeans"
                    )
                } else {
                    nFixed[i] <- (.getOneMinusQNorm(alpha / sided) +
                        .getOneMinusQNorm(beta))^2 / ((theta - thetaH0) / stDev)^2
                }
            } else {
                up <- 2
                while (stats::pt(
                    stats::qt(1 - alpha / 2, max(0.001, up - 1)), max(0.001, up - 1),
                    sqrt(up) * (theta - thetaH0) / stDev
                ) -
                    stats::pt(
                        -stats::qt(1 - alpha / 2, max(0.001, up - 1)),
                        max(0.001, up - 1), sqrt(up) * (theta - thetaH0) / stDev
                    ) > beta) {
                    up <- 2 * up
                }
                if (normalApproximation == FALSE) {
                    nFixed[i] <- .getOneDimensionalRoot(
                        function(n) {
                            return(stats::pt(
                                stats::qt(1 - alpha / 2, max(0.001, n - 1)), max(0.001, n - 1),
                                sqrt(n) * (theta - thetaH0) / stDev
                            ) -
                                stats::pt(
                                    -stats::qt(1 - alpha / 2, max(0.001, n - 1)),
                                    max(0.001, n - 1), sqrt(n) * (theta - thetaH0) / stDev
                                ) - beta)
                        },
                        lower = 0.001, upper = up, tolerance = 1e-04,
                        callingFunctionInformation = ".getSampleSizeFixedMeans"
                    )
                } else {
                    nFixed[i] <- .getOneDimensionalRoot(
                        function(n) {
                            return(stats::pnorm(.getOneMinusQNorm(alpha / 2) - sqrt(n) * (theta - thetaH0) / stDev) -
                                stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(n) * (theta - thetaH0) / stDev) - beta)
                        },
                        lower = 0.001, upper = up, tolerance = 1e-04,
                        callingFunctionInformation = ".getSampleSizeFixedMeans"
                    )
                }
            }
        } else if (groups == 2) {
            if (sided == 1 || !twoSidedPower) {
                if (!meanRatio) {
                    # allocationRatioPlanned = 0 provides optimum sample size
                    if (allocationRatioPlanned == 0) {
                        allocationRatioPlanned <- 1
                    }
                    if (normalApproximation == FALSE) {
                        up <- 2
                        while (stats::pt(
                            stats::qt(1 - alpha / sided, up *
                                (1 + allocationRatioPlanned) - 2),
                            up * (1 + allocationRatioPlanned) - 2,
                            sqrt(up) * sqrt(allocationRatioPlanned / (1 + allocationRatioPlanned)) *
                                abs(theta - thetaH0) / stDev
                        ) > beta) {
                            up <- 2 * up
                        }
                        n2Fixed <- .getOneDimensionalRoot(
                            function(x) {
                                return(stats::pt(
                                    stats::qt(1 - alpha / sided, max(
                                        0.001,
                                        x * (1 + allocationRatioPlanned) - 2
                                    )),
                                    max(0.001, x * (1 + allocationRatioPlanned) - 2),
                                    sqrt(x) * sqrt(allocationRatioPlanned /
                                        (1 + allocationRatioPlanned)) *
                                        abs(theta - thetaH0) / stDev
                                ) - beta)
                            },
                            lower = 0.001, upper = up, tolerance = 1e-04,
                            callingFunctionInformation = ".getSampleSizeFixedMeans"
                        )
                        nFixed[i] <- n2Fixed * (1 + allocationRatioPlanned)
                    } else {
                        nFixed[i] <- (1 + allocationRatioPlanned)^2 / allocationRatioPlanned *
                            (.getOneMinusQNorm(alpha / sided) + .getOneMinusQNorm(beta))^2 /
                            ((theta - thetaH0) / stDev)^2
                    }
                } else {
                    # allocationRatioPlanned = 0 provides optimum sample size
                    if (allocationRatioPlanned == 0) {
                        allocationRatioPlanned <- 1 / thetaH0
                    }
                    if (!normalApproximation) {
                        up <- 2
                        while (stats::pt(
                            stats::qt(
                                1 - alpha / sided,
                                up * (1 + allocationRatioPlanned) - 2
                            ),
                            up * (1 + allocationRatioPlanned) - 2,
                            sqrt(up * allocationRatioPlanned /
                                (1 + allocationRatioPlanned * thetaH0^2)) *
                                abs(theta - thetaH0) / stDev
                        ) > beta) {
                            up <- 2 * up
                        }
                        n2Fixed <- .getOneDimensionalRoot(
                            function(n2) {
                                return(stats::pt(
                                    stats::qt(1 - alpha / sided, max(
                                        0.001,
                                        n2 * (1 + allocationRatioPlanned) - 2
                                    )),
                                    max(0.001, n2 * (1 + allocationRatioPlanned) - 2),
                                    sqrt(n2 * allocationRatioPlanned /
                                        (1 + allocationRatioPlanned * thetaH0^2)) *
                                        abs(theta - thetaH0) / stDev
                                ) - beta)
                            },
                            lower = 0.001, upper = up, tolerance = 1e-04,
                            callingFunctionInformation = ".getSampleSizeFixedMeans"
                        )
                        nFixed[i] <- n2Fixed * (1 + allocationRatioPlanned)
                    } else {
                        nFixed[i] <- (1 + 1 / allocationRatioPlanned + thetaH0^2 *
                            (1 + allocationRatioPlanned)) *
                            (.getOneMinusQNorm(alpha / sided) + .getOneMinusQNorm(beta))^2 /
                            ((theta - thetaH0) / stDev)^2
                    }
                }
            } else {
                if (!normalApproximation) {
                    if (allocationRatioPlanned == 0) {
                        allocationRatioPlanned <- 1
                    }
                    up <- 2
                    while (stats::pt(
                        stats::qt(1 - alpha / 2, max(0.001, up * (1 + allocationRatioPlanned) - 2)),
                        max(0.001, up * (1 + allocationRatioPlanned) - 2),
                        sqrt(up) * sqrt(allocationRatioPlanned / (1 + allocationRatioPlanned)) *
                            (theta - thetaH0) / stDev
                    ) - stats::pt(
                        -stats::qt(
                            1 - alpha / 2,
                            up * (1 + allocationRatioPlanned) - 2
                        ),
                        up * (1 + allocationRatioPlanned) - 2,
                        sqrt(up) * sqrt(allocationRatioPlanned / (1 + allocationRatioPlanned)) *
                            (theta - thetaH0) / stDev
                    ) > beta) {
                        up <- 2 * up
                    }
                    n2Fixed <- .getOneDimensionalRoot(
                        function(n2) {
                            return(stats::pt(
                                stats::qt(1 - alpha / 2, max(0.001, n2 * (1 + allocationRatioPlanned) - 2)),
                                max(0.001, n2 * (1 + allocationRatioPlanned) - 2),
                                sqrt(n2) * sqrt(allocationRatioPlanned / (1 + allocationRatioPlanned)) *
                                    (theta - thetaH0) / stDev
                            ) - stats::pt(
                                -stats::qt(
                                    1 - alpha / 2,
                                    max(0.001, n2 * (1 + allocationRatioPlanned) - 2)
                                ),
                                max(0.001, n2 * (1 + allocationRatioPlanned) - 2),
                                sqrt(n2) * sqrt(allocationRatioPlanned / (1 + allocationRatioPlanned)) *
                                    (theta - thetaH0) / stDev
                            ) - beta)
                        },
                        lower = 0.001, upper = up, tolerance = 1e-04,
                        callingFunctionInformation = ".getSampleSizeFixedMeans"
                    )
                    nFixed[i] <- n2Fixed * (1 + allocationRatioPlanned)
                } else {
                    up <- 2
                    while (stats::pnorm(.getOneMinusQNorm(alpha / 2) - sqrt(up / 4) * (theta - thetaH0) / stDev) -
                        stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(up / 4) *
                            (theta - thetaH0) / stDev) > beta) {
                        up <- 2 * up
                    }

                    nFixed[i] <- (1 + allocationRatioPlanned)^2 / (4 * allocationRatioPlanned) *
                        .getOneDimensionalRoot(
                            function(n) {
                                return(stats::pnorm(.getOneMinusQNorm(alpha / 2) -
                                    sqrt(n / 4) * (theta - thetaH0) / stDev) -
                                    stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(n / 4) *
                                        (theta - thetaH0) / stDev) - beta)
                            },
                            lower = 0.001, upper = up, tolerance = 1e-04,
                            callingFunctionInformation = ".getSampleSizeFixedMeans"
                        )
                }
            }
        }
    }

    if (groups == 1) {
        return(list(
            alpha = alpha,
            beta = beta,
            sided = sided,
            groups = groups,
            thetaH0 = thetaH0,
            alternative = alternative,
            stDev = stDev,
            normalApproximation = normalApproximation,
            nFixed = nFixed
        ))
    } else if (groups == 2) {
        n1Fixed <- nFixed * allocationRatioPlanned / (1 + allocationRatioPlanned)
        n2Fixed <- n1Fixed / allocationRatioPlanned
        return(list(
            alpha = alpha,
            beta = beta,
            sided = sided,
            groups = groups,
            allocationRatioPlanned = allocationRatioPlanned,
            thetaH0 = thetaH0,
            meanRatio = meanRatio,
            alternative = alternative,
            stDev = stDev,
            normalApproximation = normalApproximation,
            n1Fixed = n1Fixed,
            n2Fixed = n2Fixed,
            nFixed = nFixed
        ))
    }
}

.getSampleSizeSequentialMeans <- function(fixedSampleSize, designCharacteristics) {
    kMax <- designCharacteristics$.design$kMax
    numberOfSubjects <- matrix(NA_real_, kMax, length(fixedSampleSize$alternative))
    numberOfSubjects1 <- matrix(NA_real_, kMax, length(fixedSampleSize$alternative))
    numberOfSubjects2 <- matrix(NA_real_, kMax, length(fixedSampleSize$alternative))
    maxNumberOfSubjects <- rep(NA_real_, length(fixedSampleSize$alternative))
    expectedNumberOfSubjectsH0 <- rep(NA_real_, length(fixedSampleSize$alternative))
    expectedNumberOfSubjectsH01 <- rep(NA_real_, length(fixedSampleSize$alternative))
    expectedNumberOfSubjectsH1 <- rep(NA_real_, length(fixedSampleSize$alternative))

    informationRates <- designCharacteristics$information / designCharacteristics$shift

    for (i in (1:length(fixedSampleSize$alternative))) {
        maxNumberOfSubjects[i] <- fixedSampleSize$nFixed[i] * designCharacteristics$inflationFactor

        numberOfSubjects[, i] <- maxNumberOfSubjects[i] *
            c(informationRates[1], (informationRates[2:kMax] - informationRates[1:(kMax - 1)]))

        expectedNumberOfSubjectsH0[i] <- designCharacteristics$averageSampleNumber0 *
            fixedSampleSize$nFixed[i]
        expectedNumberOfSubjectsH01[i] <- designCharacteristics$averageSampleNumber01 *
            fixedSampleSize$nFixed[i]
        expectedNumberOfSubjectsH1[i] <- designCharacteristics$averageSampleNumber1 *
            fixedSampleSize$nFixed[i]

        if (fixedSampleSize$groups == 2) {
            if (length(fixedSampleSize$allocationRatioPlanned) > 1) {
                allocationRatioPlanned <- fixedSampleSize$allocationRatioPlanned[i]
            } else {
                allocationRatioPlanned <- fixedSampleSize$allocationRatioPlanned
            }
            numberOfSubjects1[, i] <- numberOfSubjects[, i] * allocationRatioPlanned /
                (1 + allocationRatioPlanned)
            numberOfSubjects2[, i] <- numberOfSubjects[, i] / (1 + allocationRatioPlanned)
        }
    }

    if (fixedSampleSize$groups == 1) {
        return(list(
            alpha = fixedSampleSize$alpha,
            beta = fixedSampleSize$beta,
            sided = fixedSampleSize$sided,
            groups = fixedSampleSize$groups,
            thetaH0 = fixedSampleSize$thetaH0,
            alternative = fixedSampleSize$alternative,
            stDev = fixedSampleSize$stDev,
            normalApproximation = fixedSampleSize$normalApproximation,
            informationRates = matrix(informationRates, ncol = 1),
            maxNumberOfSubjects = maxNumberOfSubjects,
            numberOfSubjects = .getColumnCumSum(numberOfSubjects),
            expectedNumberOfSubjectsH0 = expectedNumberOfSubjectsH0,
            expectedNumberOfSubjectsH01 = expectedNumberOfSubjectsH01,
            expectedNumberOfSubjectsH1 = expectedNumberOfSubjectsH1,
            rejectPerStage = designCharacteristics$rejectionProbabilities,
            futilityPerStage = designCharacteristics$futilityProbabilities
        ))
    } else {
        return(list(
            alpha = fixedSampleSize$alpha,
            beta = fixedSampleSize$beta,
            sided = fixedSampleSize$sided,
            groups = fixedSampleSize$groups,
            allocationRatioPlanned = fixedSampleSize$allocationRatioPlanned,
            thetaH0 = fixedSampleSize$thetaH0,
            alternative = fixedSampleSize$alternative,
            stDev = fixedSampleSize$stDev,
            normalApproximation = fixedSampleSize$normalApproximation,
            meanRatio = fixedSampleSize$meanRatio,
            informationRates = matrix(informationRates, ncol = 1),
            maxNumberOfSubjects = maxNumberOfSubjects,
            numberOfSubjects = .getColumnCumSum(numberOfSubjects),
            numberOfSubjects1 = .getColumnCumSum(numberOfSubjects1),
            numberOfSubjects2 = .getColumnCumSum(numberOfSubjects2),
            expectedNumberOfSubjectsH0 = expectedNumberOfSubjectsH0,
            expectedNumberOfSubjectsH01 = expectedNumberOfSubjectsH01,
            expectedNumberOfSubjectsH1 = expectedNumberOfSubjectsH1,
            rejectPerStage = designCharacteristics$rejectionProbabilities,
            futilityPerStage = designCharacteristics$futilityProbabilities
        ))
    }
}

.getColumnCumSum <- function(x) {
    if (is.matrix(x)) {
        result <- x
        for (i in 1:ncol(x)) {
            result[, i] <- cumsum(x[, i])
        }
        return(result)
    }

    return(cumsum(x))
}

.getFarringtonManningValuesDiff <- function(..., rate1, rate2, theta, allocation) {
    if (theta == 0) {
        ml1 <- (allocation * rate1 + rate2) / (1 + allocation)
        ml2 <- ml1
        return(c(ml1, ml2))
    }

    a <- 1 + 1 / allocation
    b <- -(1 + 1 / allocation + rate1 + rate2 / allocation + theta * (1 / allocation + 2))
    c <- theta^2 + theta * (2 * rate1 + 1 / allocation + 1) + rate1 + rate2 / allocation
    d <- -theta * (1 + theta) * rate1

    v <- b^3 / (3 * a)^3 - b * c / (6 * a^2) + d / (2 * a)
    if (!is.na(v) && (v == 0)) {
        u <- sqrt(b^2 / (3 * a)^2 - c / (3 * a))
        w <- acos(-1) / 2
    } else {
        u <- sign(v) * sqrt(b^2 / (3 * a)^2 - c / (3 * a))
        w <- 1 / 3 * (acos(-1) + acos(v / u^3))
    }
    ml1 <- min(max(0, 2 * u * cos(w) - b / (3 * a)), 1)
    ml2 <- min(max(0, ml1 - theta), 1)

    return(c(ml1, ml2))
}

.getFarringtonManningValuesRatio <- function(..., rate1, rate2, theta, allocation) {
    if (theta == 1) {
        ml1 <- (allocation * rate1 + rate2) / (1 + allocation)
        ml2 <- ml1
        return(c(ml1, ml2))
    }

    a <- 1 + 1 / allocation
    b <- -((1 + rate2 / allocation) * theta + 1 / allocation + rate1)
    c <- (rate1 + rate2 / allocation) * theta
    ml1 <- (-b - sqrt(b^2 - 4 * a * c)) / (2 * a)
    ml2 <- ml1 / theta

    return(c(ml1, ml2))
}

#
# @title
# Get Farrington Manning Values
#
# @description
# Calculates and returns the maximum likelihood estimates under H0.
#
# @details
# Calculation of maximum likelihood estimates under
# H0: pi1 - pi2 = theta or H0: pi1 / pi2 = theta
#
# @references
# Farrington & Manning (1990)
# Wassmer (2003)
#
# @keywords internal
#
.getFarringtonManningValues <- function(rate1, rate2, theta, allocation, method = c("diff", "ratio")) {
    method <- match.arg(method)
    if (method == "diff") {
        ml <- .getFarringtonManningValuesDiff(rate1 = rate1, rate2 = rate2, theta = theta, allocation = allocation)
    } else {
        ml <- .getFarringtonManningValuesRatio(rate1 = rate1, rate2 = rate2, theta = theta, allocation = allocation)
    }
    return(list(theta = theta, method = method, ml1 = ml[1], ml2 = ml[2]))
}

.getSampleSizeFixedRates <- function(..., alpha = 0.025, beta = 0.2, sided = 1,
        normalApproximation = TRUE, riskRatio = FALSE,
        thetaH0 = 0, pi1 = seq(0.4, 0.6, 0.1), pi2 = 0.2,
        groups = 2, allocationRatioPlanned = 1) {
    if (groups == 1) {
        nFixed <- rep(NA_real_, length(pi1))

        for (i in 1:length(pi1)) {
            if (normalApproximation) {
                nFixed[i] <- (.getOneMinusQNorm(alpha / sided) * sqrt(thetaH0 * (1 - thetaH0)) +
                    .getOneMinusQNorm(beta) * sqrt(pi1[i] * (1 - pi1[i])))^2 /
                    (pi1[i] - thetaH0)^2
            } else {
                ifelse(pi1[i] > thetaH0, lower.tail <- FALSE, lower.tail <- TRUE)
                iterations <- 1
                if (lower.tail) {
                    nup <- 2
                    while ((stats::pbinom(stats::qbinom(alpha, nup, thetaH0, lower.tail = lower.tail) - 1,
                        nup, pi1[i],
                        lower.tail = lower.tail
                    ) < 1 - beta) && (iterations <= 50)) {
                        nup <- 2 * nup
                        iterations <- iterations + 1
                    }
                    if (iterations > 50) {
                        nFixed[i] <- Inf
                    } else {
                        prec <- 2
                        nlow <- 2
                        while (prec > 1) {
                            nFixed[i] <- round((nlow + nup) / 2)
                            ifelse(stats::pbinom(stats::qbinom(alpha, nFixed[i], thetaH0, lower.tail = lower.tail) - 1,
                                nFixed[i], pi1[i],
                                lower.tail = lower.tail
                            ) < 1 - beta,
                            nlow <- nFixed[i], nup <- nFixed[i]
                            )
                            prec <- nup - nlow
                        }
                        if (stats::pbinom(stats::qbinom(alpha, nFixed[i], thetaH0, lower.tail = lower.tail) - 1,
                                nFixed[i], pi1[i],
                                lower.tail = lower.tail
                            ) < 1 - beta) {
                            nFixed[i] <- nFixed[i] + 1
                        }
                    }
                } else {
                    nup <- 2
                    while ((stats::pbinom(stats::qbinom(alpha, nup, thetaH0, lower.tail = lower.tail),
                        nup, pi1[i],
                        lower.tail = lower.tail
                    ) < 1 - beta) && (iterations <= 50)) {
                        nup <- 2 * nup
                        iterations <- iterations + 1
                    }
                    if (iterations > 50) {
                        nFixed[i] <- Inf
                    } else {
                        prec <- 2
                        nlow <- 2
                        while (prec > 1) {
                            nFixed[i] <- round((nlow + nup) / 2)
                            ifelse(stats::pbinom(stats::qbinom(alpha, nFixed[i], thetaH0, lower.tail = lower.tail),
                                nFixed[i], pi1[i],
                                lower.tail = lower.tail
                            ) < 1 - beta,
                            nlow <- nFixed[i], nup <- nFixed[i]
                            )
                            prec <- nup - nlow
                        }
                        if (stats::pbinom(stats::qbinom(alpha, nFixed[i], thetaH0, lower.tail = lower.tail),
                                nFixed[i], pi1[i],
                                lower.tail = lower.tail
                            ) < 1 - beta) {
                            nFixed[i] <- nFixed[i] + 1
                        }
                    }
                }
            }
        }

        return(list(
            alpha = alpha,
            beta = beta,
            sided = sided,
            groups = groups,
            thetaH0 = thetaH0,
            pi1 = pi1,
            normalApproximation = normalApproximation,
            nFixed = nFixed
        ))
    }

    if (groups == 2) {
        n1Fixed <- rep(NA_real_, length(pi1))
        n2Fixed <- rep(NA_real_, length(pi1))
        nFixed <- rep(NA_real_, length(pi1))
        if (allocationRatioPlanned == 0) {
            allocationRatioPlannedVec <- rep(NA_real_, length(pi1))
        }

        for (i in 1:length(pi1)) {
            if (!riskRatio) {
                # allocationRatioPlanned = 0 provides optimum sample size
                if (allocationRatioPlanned == 0) {
                    allocationRatioPlannedVec[i] <- stats::optimize(function(x) {
                        fm <- .getFarringtonManningValues(
                            rate1 = pi1[i], rate2 = pi2,
                            theta = thetaH0, allocation = x, method = "diff"
                        )
                        n1 <- (.getOneMinusQNorm(alpha / sided) * sqrt(fm$ml1 * (1 - fm$ml1) + fm$ml2 * (1 - fm$ml2) * x) +
                            .getOneMinusQNorm(beta) * sqrt(pi1[i] * (1 - pi1[i]) + pi2 * (1 - pi2) * x))^2 /
                            (pi1[i] - pi2 - thetaH0)^2
                        return((1 + x) / x * n1)
                    }, interval = c(0, 5), tol = 0.0001)$minimum
                    fm <- .getFarringtonManningValues(
                        rate1 = pi1[i], rate2 = pi2, theta = thetaH0,
                        allocation = allocationRatioPlannedVec[i], method = "diff"
                    )
                    n1Fixed[i] <- (.getOneMinusQNorm(alpha / sided) * sqrt(fm$ml1 * (1 - fm$ml1) +
                        fm$ml2 * (1 - fm$ml2) * allocationRatioPlannedVec[i]) +
                        .getOneMinusQNorm(beta) * sqrt(pi1[i] * (1 - pi1[i]) + pi2 * (1 - pi2) *
                            allocationRatioPlannedVec[i]))^2 / (pi1[i] - pi2 - thetaH0)^2
                } else {
                    fm <- .getFarringtonManningValues(
                        rate1 = pi1[i], rate2 = pi2,
                        theta = thetaH0, allocation = allocationRatioPlanned, method = "diff"
                    )
                    n1Fixed[i] <- (.getOneMinusQNorm(alpha / sided) * sqrt(fm$ml1 * (1 - fm$ml1) +
                        fm$ml2 * (1 - fm$ml2) * allocationRatioPlanned) +
                        .getOneMinusQNorm(beta) * sqrt(pi1[i] * (1 - pi1[i]) + pi2 * (1 - pi2) *
                            allocationRatioPlanned))^2 / (pi1[i] - pi2 - thetaH0)^2
                }
            } else {
                if (allocationRatioPlanned == 0) {
                    # allocationRatioPlanned = 0 provides optimum sample size
                    allocationRatioPlannedVec[i] <- stats::optimize(function(x) {
                        fm <- .getFarringtonManningValues(
                            rate1 = pi1[i], rate2 = pi2,
                            theta = thetaH0, allocation = x, method = "ratio"
                        )
                        n1 <- (.getOneMinusQNorm(alpha / sided) * sqrt(fm$ml1 * (1 - fm$ml1) +
                            fm$ml2 * (1 - fm$ml2) * x * thetaH0^2) +
                            .getOneMinusQNorm(beta) * sqrt(pi1[i] * (1 - pi1[i]) + pi2 *
                                (1 - pi2) * x * thetaH0^2))^2 / (pi1[i] - thetaH0 * pi2)^2
                        return((1 + x) / x * n1)
                    }, interval = c(0, 5), tol = 0.0001)$minimum
                    fm <- .getFarringtonManningValues(
                        rate1 = pi1[i], rate2 = pi2, theta = thetaH0,
                        allocation = allocationRatioPlannedVec[i], method = "ratio"
                    )
                    n1Fixed[i] <- (.getOneMinusQNorm(alpha / sided) * sqrt(fm$ml1 * (1 - fm$ml1) +
                        fm$ml2 * (1 - fm$ml2) * allocationRatioPlannedVec[i] * thetaH0^2) +
                        .getOneMinusQNorm(beta) * sqrt(pi1[i] * (1 - pi1[i]) + pi2 * (1 - pi2) *
                            allocationRatioPlannedVec[i] * thetaH0^2))^2 / (pi1[i] - thetaH0 * pi2)^2
                } else {
                    fm <- .getFarringtonManningValues(
                        rate1 = pi1[i], rate2 = pi2,
                        theta = thetaH0, allocation = allocationRatioPlanned, method = "ratio"
                    )
                    n1Fixed[i] <- (.getOneMinusQNorm(alpha / sided) * sqrt(fm$ml1 * (1 - fm$ml1) +
                        fm$ml2 * (1 - fm$ml2) * allocationRatioPlanned * thetaH0^2) +
                        .getOneMinusQNorm(beta) * sqrt(pi1[i] * (1 - pi1[i]) + pi2 * (1 - pi2) *
                            allocationRatioPlanned * thetaH0^2))^2 / (pi1[i] - thetaH0 * pi2)^2
                }
            }
        }
        if (allocationRatioPlanned == 0) {
            allocationRatioPlanned <- allocationRatioPlannedVec
        }

        n2Fixed <- n1Fixed / allocationRatioPlanned
        nFixed <- n1Fixed + n2Fixed

        return(list(
            alpha = alpha,
            beta = beta,
            sided = sided,
            groups = groups,
            allocationRatioPlanned = allocationRatioPlanned,
            thetaH0 = thetaH0,
            pi1 = pi1,
            pi2 = pi2,
            normalApproximation = normalApproximation,
            riskRatio = riskRatio,
            n1Fixed = n1Fixed,
            n2Fixed = n2Fixed,
            nFixed = nFixed
        ))
    }
}

.getSampleSizeSequentialRates <- function(fixedSampleSize, designCharacteristics) {
    kMax <- designCharacteristics$.design$kMax
    numberOfSubjects <- matrix(NA_real_, kMax, length(fixedSampleSize$pi1))
    numberOfSubjects1 <- matrix(NA_real_, kMax, length(fixedSampleSize$pi1))
    numberOfSubjects2 <- matrix(NA_real_, kMax, length(fixedSampleSize$pi1))
    maxNumberOfSubjects <- rep(NA_real_, length(fixedSampleSize$pi1))
    expectedNumberOfSubjectsH0 <- rep(NA_real_, length(fixedSampleSize$pi1))
    expectedNumberOfSubjectsH01 <- rep(NA_real_, length(fixedSampleSize$pi1))
    expectedNumberOfSubjectsH1 <- rep(NA_real_, length(fixedSampleSize$pi1))

    informationRates <- designCharacteristics$information / designCharacteristics$shift

    for (i in 1:length(fixedSampleSize$pi1)) {
        maxNumberOfSubjects[i] <- fixedSampleSize$nFixed[i] * designCharacteristics$inflationFactor

        numberOfSubjects[, i] <- maxNumberOfSubjects[i] * c(
            informationRates[1],
            (informationRates[2:kMax] - informationRates[1:(kMax - 1)])
        )

        expectedNumberOfSubjectsH0[i] <- designCharacteristics$averageSampleNumber0 * fixedSampleSize$nFixed[i]
        expectedNumberOfSubjectsH01[i] <- designCharacteristics$averageSampleNumber01 * fixedSampleSize$nFixed[i]
        expectedNumberOfSubjectsH1[i] <- designCharacteristics$averageSampleNumber1 * fixedSampleSize$nFixed[i]

        if (fixedSampleSize$groups == 2) {
            if (length(fixedSampleSize$allocationRatioPlanned) > 1) {
                allocationRatioPlanned <- fixedSampleSize$allocationRatioPlanned[i]
            } else {
                allocationRatioPlanned <- fixedSampleSize$allocationRatioPlanned
            }
            numberOfSubjects1[, i] <- numberOfSubjects[, i] * allocationRatioPlanned /
                (1 + allocationRatioPlanned)
            numberOfSubjects2[, i] <- numberOfSubjects[, i] / (1 + allocationRatioPlanned)
        }
    }

    if (fixedSampleSize$groups == 1) {
        return(list(
            alpha = fixedSampleSize$alpha,
            beta = fixedSampleSize$beta,
            sided = fixedSampleSize$sided,
            groups = fixedSampleSize$groups,
            thetaH0 = fixedSampleSize$thetaH0,
            pi1 = fixedSampleSize$pi1,
            normalApproximation = fixedSampleSize$normalApproximation,
            informationRates = matrix(informationRates, ncol = 1),
            maxNumberOfSubjects = maxNumberOfSubjects,
            numberOfSubjects = .getColumnCumSum(numberOfSubjects),
            expectedNumberOfSubjectsH0 = expectedNumberOfSubjectsH0,
            expectedNumberOfSubjectsH01 = expectedNumberOfSubjectsH01,
            expectedNumberOfSubjectsH1 = expectedNumberOfSubjectsH1,
            rejectPerStage = designCharacteristics$rejectionProbabilities,
            futilityPerStage = designCharacteristics$futilityProbabilities
        ))
    } else {
        return(list(
            alpha = fixedSampleSize$alpha,
            beta = fixedSampleSize$beta,
            sided = fixedSampleSize$sided,
            groups = fixedSampleSize$groups,
            allocationRatioPlanned = fixedSampleSize$allocationRatioPlanned,
            thetaH0 = fixedSampleSize$thetaH0,
            pi1 = fixedSampleSize$pi1,
            pi2 = fixedSampleSize$pi2,
            normalApproximation = fixedSampleSize$normalApproximation,
            riskRatio = fixedSampleSize$riskRatio,
            informationRates = matrix(informationRates, ncol = 1),
            maxNumberOfSubjects = maxNumberOfSubjects,
            numberOfSubjects = .getColumnCumSum(numberOfSubjects),
            numberOfSubjects1 = .getColumnCumSum(numberOfSubjects1),
            numberOfSubjects2 = .getColumnCumSum(numberOfSubjects2),
            expectedNumberOfSubjectsH0 = expectedNumberOfSubjectsH0,
            expectedNumberOfSubjectsH01 = expectedNumberOfSubjectsH01,
            expectedNumberOfSubjectsH1 = expectedNumberOfSubjectsH1,
            rejectPerStage = designCharacteristics$rejectionProbabilities,
            futilityPerStage = designCharacteristics$futilityProbabilities
        ))
    }
}

.getPiecewiseExpStartTimesWithoutLeadingZero <- function(piecewiseSurvivalTime) {
    if (is.null(piecewiseSurvivalTime) || length(piecewiseSurvivalTime) == 0 ||
            all(is.na(piecewiseSurvivalTime))) {
        return(NA_real_)
    }

    if (piecewiseSurvivalTime[1] != 0) {
        stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
            "the first value of 'piecewiseSurvivalTime' (",
            .arrayToString(piecewiseSurvivalTime), ") must be 0",
            call. = FALSE
        )
    }

    if (length(piecewiseSurvivalTime) == 1) {
        return(numeric(0))
    }

    if (length(piecewiseSurvivalTime) < 2) {
        stop(
            C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "length of 'piecewiseSurvivalTime' (",
            length(piecewiseSurvivalTime), ") must be > 1"
        )
    }

    return(piecewiseSurvivalTime[2:length(piecewiseSurvivalTime)])
}

.getEventProbabilityFunction <- function(..., time, piecewiseLambda, piecewiseSurvivalTime, phi, kappa) {
    if (length(piecewiseLambda) == 1) {
        if (kappa != 1 && phi > 0) {
            stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "Weibull distribution cannot ",
                "be used together with specified dropout rate (use simulation instead)",
                call. = FALSE
            )
        }

        return(piecewiseLambda / (piecewiseLambda + phi) *
            pweibull(time, shape = kappa, scale = 1 / (piecewiseLambda + phi), lower.tail = TRUE, log.p = FALSE))
    }

    if (length(piecewiseSurvivalTime) != length(piecewiseLambda)) {
        stop(
            C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
            "length of 'piecewiseSurvivalTime' (", .arrayToString(piecewiseSurvivalTime),
            ") must be equal to length of 'piecewiseLambda' (", .arrayToString(piecewiseLambda), ")"
        )
    }

    piecewiseSurvivalTime <- .getPiecewiseExpStartTimesWithoutLeadingZero(piecewiseSurvivalTime)

    if (kappa != 1) {
        stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
            "Weibull distribution cannot be used for piecewise survival definition",
            call. = FALSE
        )
    }
    len <- length(piecewiseSurvivalTime)
    for (i in 1:len) {
        if (i == 1) {
            if (time <= piecewiseSurvivalTime[1]) {
                return(piecewiseLambda[1] / (piecewiseLambda[1] + phi) *
                    (1 - exp(-((piecewiseLambda[1] + phi) * time))))
            }
        } else if (i == 2) {
            cdfPart <- piecewiseLambda[1] / (piecewiseLambda[1] + phi) *
                (1 - exp(-((piecewiseLambda[1] + phi) * piecewiseSurvivalTime[1])))
            if (time <= piecewiseSurvivalTime[2]) {
                cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1]
                cdf <- cdfPart + piecewiseLambda[2] / (piecewiseLambda[2] + phi) * exp(-cdfFactor) * (
                    exp(-phi * piecewiseSurvivalTime[1]) - exp(-piecewiseLambda[2] *
                        (time - piecewiseSurvivalTime[1]) - phi * time))
                return(cdf)
            }
        } else if (i == 3) {
            cdfPart <- cdfPart + piecewiseLambda[2] / (piecewiseLambda[2] + phi) *
                exp(-piecewiseLambda[1] * piecewiseSurvivalTime[1]) * (
                    exp(-phi * piecewiseSurvivalTime[1]) - exp(-piecewiseLambda[2] *
                        (piecewiseSurvivalTime[2] - piecewiseSurvivalTime[1]) - phi * piecewiseSurvivalTime[2]))
            if (time <= piecewiseSurvivalTime[3]) {
                cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1] +
                    piecewiseLambda[2] * (piecewiseSurvivalTime[2] - piecewiseSurvivalTime[1])
                cdf <- cdfPart + piecewiseLambda[3] / (piecewiseLambda[3] + phi) * exp(-cdfFactor) * (
                    exp(-phi * piecewiseSurvivalTime[2]) - exp(-piecewiseLambda[3] *
                        (time - piecewiseSurvivalTime[2]) - phi * time))
                return(cdf)
            }
        } else if (i > 3) {
            cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1] +
                sum(piecewiseLambda[2:(i - 2)] * (piecewiseSurvivalTime[2:(i - 2)] -
                    piecewiseSurvivalTime[1:(i - 3)]))
            cdfPart <- cdfPart + piecewiseLambda[i - 1] / (piecewiseLambda[i - 1] + phi) * exp(-cdfFactor) * (
                exp(-phi * piecewiseSurvivalTime[i - 2]) - exp(-piecewiseLambda[i - 1] *
                    (piecewiseSurvivalTime[i - 1] - piecewiseSurvivalTime[i - 2]) - phi * piecewiseSurvivalTime[i - 1]))
            if (time <= piecewiseSurvivalTime[i]) {
                cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1] +
                    sum(piecewiseLambda[2:(i - 1)] * (piecewiseSurvivalTime[2:(i - 1)] - piecewiseSurvivalTime[1:(i - 2)]))
                cdf <- cdfPart + piecewiseLambda[i] / (piecewiseLambda[i] + phi) * exp(-cdfFactor) * (
                    exp(-phi * piecewiseSurvivalTime[i - 1]) - exp(-piecewiseLambda[i] *
                        (time - piecewiseSurvivalTime[i - 1]) - phi * time))
                return(cdf)
            }
        }
    }

    if (len == 1) {
        cdfPart <- piecewiseLambda[1] / (piecewiseLambda[1] + phi) *
            (1 - exp(-((piecewiseLambda[1] + phi) * piecewiseSurvivalTime[1])))
    } else if (len == 2) {
        cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1]
        cdfPart <- cdfPart + piecewiseLambda[len] / (piecewiseLambda[len] + phi) * exp(-cdfFactor) * (
            exp(-phi * piecewiseSurvivalTime[len - 1]) - exp(-piecewiseLambda[len] *
                (piecewiseSurvivalTime[len] - piecewiseSurvivalTime[len - 1]) - phi * piecewiseSurvivalTime[len]))
    } else {
        cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1] +
            sum(piecewiseLambda[2:(len - 1)] * (piecewiseSurvivalTime[2:(len - 1)] - piecewiseSurvivalTime[1:(len - 2)]))
        cdfPart <- cdfPart + piecewiseLambda[len] / (piecewiseLambda[len] + phi) * exp(-cdfFactor) * (
            exp(-phi * piecewiseSurvivalTime[len - 1]) - exp(-piecewiseLambda[len] *
                (piecewiseSurvivalTime[len] - piecewiseSurvivalTime[len - 1]) - phi * piecewiseSurvivalTime[len]))
    }

    if (len == 1) {
        cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1]
    } else {
        cdfFactor <- cdfFactor + piecewiseLambda[len] * (piecewiseSurvivalTime[len] - piecewiseSurvivalTime[len - 1])
    }

    cdf <- cdfPart + piecewiseLambda[len + 1] / (piecewiseLambda[len + 1] + phi) * exp(-cdfFactor) * (
        exp(-phi * piecewiseSurvivalTime[len]) - exp(-piecewiseLambda[len + 1] *
            (time - piecewiseSurvivalTime[len]) - phi * time))

    return(cdf)
}

.getEventProbabilityFunctionVec <- function(..., timeVector, piecewiseLambda, piecewiseSurvivalTime, phi, kappa) {
    result <- c()
    for (time in timeVector) {
        result <- c(result, .getEventProbabilityFunction(
            time = time, piecewiseLambda = piecewiseLambda,
            piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa
        ))
    }
    return(result)
}

#' @title
#' Get Event Probabilities
#'
#' @description
#' Returns the event probabilities for specified parameters at given time vector.
#'
#' @param time A numeric vector with time values.
#' @inheritParams param_lambda1
#' @inheritParams param_lambda2
#' @inheritParams param_piecewiseSurvivalTime
#' @inheritParams param_hazardRatio
#' @inheritParams param_kappa
#' @inheritParams param_allocationRatioPlanned_sampleSize
#' @inheritParams param_accrualTime
#' @inheritParams param_accrualIntensity
#' @inheritParams param_accrualIntensityType
#' @inheritParams param_dropoutRate1
#' @inheritParams param_dropoutRate2
#' @inheritParams param_dropoutTime
#' @param maxNumberOfSubjects If \code{maxNumberOfSubjects > 0} is specified,
#'        the end of accrual at specified \code{accrualIntensity} for the specified
#'        number of subjects is determined or \code{accrualIntensity} is calculated
#'        at fixed end of accrual.
#' @inheritParams param_three_dots
#'
#' @details
#' The function computes the overall event probabilities in a two treatment groups design.
#' For details of the parameters see \code{\link[=getSampleSizeSurvival]{getSampleSizeSurvival()}}.
#'
#' @return Returns a \code{\link{EventProbabilities}} object.
#' The following generics (R generic functions) are available for this result object:
#' \itemize{
#'   \item \code{\link[=names.FieldSet]{names()}} to obtain the field names,
#'   \item \code{\link[=print.FieldSet]{print()}} to print the object,
#'   \item \code{\link[=summary.ParameterSet]{summary()}} to display a summary of the object,
#'   \item \code{\link[=plot.EventProbabilities]{plot()}} to plot the object,
#'   \item \code{\link[=as.data.frame.ParameterSet]{as.data.frame()}} to coerce the object to a \code{\link[base]{data.frame}},
#'   \item \code{\link[=as.matrix.FieldSet]{as.matrix()}} to coerce the object to a \code{\link[base]{matrix}}.
#' }
#' @template how_to_get_help_for_generics
#'
#' @template examples_get_event_probabilities
#'
#' @export
#'
getEventProbabilities <- function(time, ...,
        accrualTime = c(0, 12), # C_ACCRUAL_TIME_DEFAULT
        accrualIntensity = 0.1, # C_ACCRUAL_INTENSITY_DEFAULT
        accrualIntensityType = c("auto", "absolute", "relative"),
        kappa = 1,
        piecewiseSurvivalTime = NA_real_,
        lambda2 = NA_real_,
        lambda1 = NA_real_,
        allocationRatioPlanned = 1,
        hazardRatio = NA_real_,
        dropoutRate1 = 0, # C_DROP_OUT_RATE_1_DEFAULT
        dropoutRate2 = 0, # C_DROP_OUT_RATE_2_DEFAULT
        dropoutTime = 12, # C_DROP_OUT_TIME_DEFAULT
        maxNumberOfSubjects = NA_real_) {
    .warnInCaseOfUnknownArguments(functionName = "getEventProbabilities", ...)

    .assertIsNumericVector(time, "time")
    .assertIsValidMaxNumberOfSubjects(maxNumberOfSubjects, naAllowed = TRUE)
    .assertIsValidAllocationRatioPlannedSampleSize(allocationRatioPlanned, maxNumberOfSubjects)
    .assertIsValidKappa(kappa)
    .assertIsSingleNumber(hazardRatio, "hazardRatio", naAllowed = TRUE)

    if (!is.na(dropoutTime) && dropoutTime <= 0) {
        stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'dropoutTime' (", dropoutTime, ") must be > 0", call. = FALSE)
    }

    if (dropoutRate1 < 0 || dropoutRate1 >= 1) {
        stop(
            C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
            "'dropoutRate1' (", dropoutRate1, ") is out of bounds [0; 1)"
        )
    }

    if (dropoutRate2 < 0 || dropoutRate2 >= 1) {
        stop(
            C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
            "'dropoutRate2' (", dropoutRate2, ") is out of bounds [0; 1)"
        )
    }

    accrualSetup <- getAccrualTime(
        accrualTime = accrualTime,
        accrualIntensity = accrualIntensity,
        accrualIntensityType = accrualIntensityType,
        maxNumberOfSubjects = maxNumberOfSubjects
    )
    accrualTime <- accrualSetup$.getAccrualTimeWithoutLeadingZero()
    accrualIntensity <- accrualSetup$accrualIntensity
    maxNumberOfSubjects <- accrualSetup$maxNumberOfSubjects

    setting <- getPiecewiseSurvivalTime(
        piecewiseSurvivalTime = piecewiseSurvivalTime,
        lambda2 = lambda2, lambda1 = lambda1,
        hazardRatio = hazardRatio, kappa = kappa,
        delayedResponseAllowed = TRUE,
        .lambdaBased = TRUE
    )

    if (!setting$delayedResponseEnabled && length(setting$lambda1) > 1 &&
            setting$.getParameterType("lambda1") == C_PARAM_USER_DEFINED) {
        warning("Only the first 'lambda1' (", lambda1[1], ") was used to calculate event probabilities", call. = FALSE)
        setting <- getPiecewiseSurvivalTime(
            piecewiseSurvivalTime = piecewiseSurvivalTime,
            lambda2 = lambda2, lambda1 = lambda1[1],
            hazardRatio = hazardRatio, kappa = kappa,
            delayedResponseAllowed = TRUE,
            .lambdaBased = TRUE
        )
    }

    piecewiseSurvivalTime <- setting$piecewiseSurvivalTime
    lambda2 <- setting$lambda2
    lambda1 <- setting$lambda1
    hazardRatio <- setting$hazardRatio

    phi <- -log(1 - c(dropoutRate1, dropoutRate2)) / dropoutTime

    if (length(accrualTime) != length(accrualIntensity)) {
        stop(
            C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
            "length of 'accrualTime' (", (length(accrualTime) + 1),
            ") must be equal to length of 'accrualIntensity' (", length(accrualIntensity), ")"
        )
    }

    if (any(accrualIntensity <= 0)) {
        stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all values of 'accrualIntensity' must be > 0")
    }

    if (any(accrualTime <= 0)) {
        stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all values of 'accrualTime' must be > 0")
    }

    if (kappa != 1 && any(phi > 0)) {
        stop(
            C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
            "for Weibull distribution (kappa != 1) drop-out rates (phi) cannot be specified"
        )
    }

    if (any(phi < 0)) {
        stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all drop-out rates (phi) must be >= 0")
    }

    .assertIsNumericVector(lambda2, "lambda2")
    if (any(lambda2 <= 0)) {
        stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all rates (lambda2) must be > 0")
    }

    eventProbabilities <- EventProbabilities(
        .piecewiseSurvivalTime = setting,
        .accrualTime           = accrualSetup,
        time                   = time,
        accrualTime            = accrualTime,
        accrualIntensity       = accrualIntensity,
        kappa                  = kappa,
        piecewiseSurvivalTime  = piecewiseSurvivalTime,
        lambda1                = lambda1,
        lambda2                = lambda2,
        allocationRatioPlanned = allocationRatioPlanned,
        hazardRatio            = hazardRatio,
        dropoutRate1           = dropoutRate1,
        dropoutRate2           = dropoutRate2,
        dropoutTime            = dropoutTime,
        maxNumberOfSubjects    = maxNumberOfSubjects
    )

    eventProbabilities$.setParameterType("time", C_PARAM_USER_DEFINED)
    eventProbabilities$.setParameterType(
        "accrualTime",
        accrualSetup$.getParameterType("accrualTime")
    )
    eventProbabilities$.setParameterType(
        "accrualIntensity",
        accrualSetup$.getParameterType("accrualIntensity")
    )
    eventProbabilities$.setParameterType("kappa", setting$.getParameterType("kappa"))
    eventProbabilities$.setParameterType(
        "piecewiseSurvivalTime",
        setting$.getParameterType("piecewiseSurvivalTime")
    )
    eventProbabilities$.setParameterType("lambda1", setting$.getParameterType("lambda1"))
    eventProbabilities$.setParameterType("lambda2", setting$.getParameterType("lambda2"))
    .setValueAndParameterType(eventProbabilities, "allocationRatioPlanned", allocationRatioPlanned, 1)
    eventProbabilities$.setParameterType("hazardRatio", setting$.getParameterType("hazardRatio"))
    .setValueAndParameterType(eventProbabilities, "dropoutRate1", dropoutRate1, C_DROP_OUT_RATE_1_DEFAULT)
    .setValueAndParameterType(eventProbabilities, "dropoutRate2", dropoutRate2, C_DROP_OUT_RATE_2_DEFAULT)
    .setValueAndParameterType(eventProbabilities, "dropoutTime", dropoutTime, C_DROP_OUT_TIME_DEFAULT)
    eventProbabilities$.setParameterType(
        "maxNumberOfSubjects",
        accrualSetup$.getParameterType("maxNumberOfSubjects")
    )

    eventProbabilities$cumulativeEventProbabilities <- numeric(0)
    eventProbabilities$overallEventProbabilities <- numeric(0) # deprecated
    eventProbabilities$eventProbabilities1 <- numeric(0)
    eventProbabilities$eventProbabilities2 <- numeric(0)

    for (timeValue in time) {
        eventProbs <- .getEventProbabilitiesGroupwise(
            time = timeValue,
            accrualTimeVector = accrualSetup$.getAccrualTimeWithoutLeadingZero(),
            accrualIntensity = accrualSetup$accrualIntensity, lambda2 = lambda2,
            lambda1 = lambda1, piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi,
            kappa = kappa, allocationRatioPlanned = allocationRatioPlanned, hazardRatio = hazardRatio
        )

        eventProbabilities$cumulativeEventProbabilities <- c(
            eventProbabilities$cumulativeEventProbabilities,
            .getEventProbabilitiesOverall(eventProbs, allocationRatioPlanned)
        )

        eventProbabilities$overallEventProbabilities <-
            eventProbabilities$cumulativeEventProbabilities # deprecated

        eventProbabilities$eventProbabilities1 <- c(
            eventProbabilities$eventProbabilities1,
            eventProbs[1]
        )
        eventProbabilities$eventProbabilities2 <- c(
            eventProbabilities$eventProbabilities2,
            eventProbs[2]
        )
    }

    eventProbabilities$.setParameterType("cumulativeEventProbabilities", C_PARAM_GENERATED)
    eventProbabilities$.setParameterType("eventProbabilities1", C_PARAM_GENERATED)
    eventProbabilities$.setParameterType("eventProbabilities2", C_PARAM_GENERATED)

    return(eventProbabilities)
}

#' @title
#' Get Number Of Subjects
#'
#' @description
#' Returns the number of recruited subjects at given time vector.
#'
#' @param time A numeric vector with time values.
#' @inheritParams param_accrualTime
#' @inheritParams param_accrualIntensity
#' @inheritParams param_accrualIntensityType
#' @param maxNumberOfSubjects If \code{maxNumberOfSubjects > 0} is specified,
#'        the end of accrual at specified \code{accrualIntensity} for the specified number of
#'        subjects is determined or \code{accrualIntensity} is calculated at fixed end of accrual.
#' @inheritParams param_three_dots
#'
#' @details
#' Calculate number of subjects over time range at given accrual time vector
#' and accrual intensity. Intensity can either be defined in absolute or
#' relative terms (for the latter, \code{maxNumberOfSubjects} needs to be defined)\cr
#' The function is used by \code{\link[=getSampleSizeSurvival]{getSampleSizeSurvival()}}.
#'
#' @return Returns a \code{\link{NumberOfSubjects}} object.
#' The following generics (R generic functions) are available for this result object:
#' \itemize{
#'   \item \code{\link[=names.FieldSet]{names()}} to obtain the field names,
#'   \item \code{\link[=print.FieldSet]{print()}} to print the object,
#'   \item \code{\link[=summary.ParameterSet]{summary()}} to display a summary of the object,
#'   \item \code{\link[=plot.NumberOfSubjects]{plot()}} to plot the object,
#'   \item \code{\link[=as.data.frame.ParameterSet]{as.data.frame()}} to coerce the object to a \code{\link[base]{data.frame}},
#'   \item \code{\link[=as.matrix.FieldSet]{as.matrix()}} to coerce the object to a \code{\link[base]{matrix}}.
#' }
#' @template how_to_get_help_for_generics
#'
#' @seealso \code{\link{AccrualTime}} for defining the accrual time.
#'
#' @template examples_get_number_of_subjects
#'
#' @export
#'
getNumberOfSubjects <- function(time, ...,
        accrualTime = c(0, 12), # C_ACCRUAL_TIME_DEFAULT
        accrualIntensity = 0.1, # C_ACCRUAL_INTENSITY_DEFAULT
        accrualIntensityType = c("auto", "absolute", "relative"),
        maxNumberOfSubjects = NA_real_) {
    .warnInCaseOfUnknownArguments(functionName = "getNumberOfSubjects", ...)

    .assertIsNumericVector(time, "time")

    accrualSetup <- getAccrualTime(
        accrualTime = accrualTime,
        accrualIntensity = accrualIntensity,
        accrualIntensityType = accrualIntensityType,
        maxNumberOfSubjects = maxNumberOfSubjects
    )
    accrualTime <- accrualSetup$.getAccrualTimeWithoutLeadingZero()
    accrualIntensity <- accrualSetup$accrualIntensity
    maxNumberOfSubjects <- accrualSetup$maxNumberOfSubjects

    if (length(accrualTime) != length(accrualIntensity)) {
        stop(
            C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
            "length of 'accrualTime' (", length(accrualTime),
            ") must be equal to length of 'accrualIntensity' (", length(accrualIntensity), ")"
        )
    }

    if (any(accrualIntensity < 0)) {
        stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all values of 'accrualIntensity' must be >= 0")
    }

    if (all(accrualIntensity < 1)) {
        stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "at least one value of 'accrualIntensity' must be >= 1")
    }

    if (any(accrualTime <= 0)) {
        stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all values of 'accrualTime' must be > 0")
    }

    numberOfSubjects <- .getNumberOfSubjects(
        time = time, accrualTime = accrualTime,
        accrualIntensity = accrualIntensity, maxNumberOfSubjects = maxNumberOfSubjects
    )

    result <- NumberOfSubjects(
        .accrualTime = accrualSetup,
        time = time,
        accrualTime = accrualTime,
        accrualIntensity = accrualIntensity,
        maxNumberOfSubjects = maxNumberOfSubjects,
        numberOfSubjects = numberOfSubjects
    )

    result$.setParameterType("time", C_PARAM_USER_DEFINED)
    result$.setParameterType("accrualTime", accrualSetup$.getParameterType("accrualTime"))
    result$.setParameterType("accrualIntensity", accrualSetup$.getParameterType("accrualIntensity"))
    result$.setParameterType("maxNumberOfSubjects", accrualSetup$.getParameterType("maxNumberOfSubjects"))
    result$.setParameterType("numberOfSubjects", C_PARAM_GENERATED)

    return(result)
}


.getLambda <- function(..., groupNumber, lambda2, lambda1, hazardRatio, kappa) {
    if (groupNumber == 1) {
        if (!any(is.na(lambda1))) {
            return(lambda1)
        }

        lambda2 <- lambda2 * hazardRatio^(1 / kappa)
    }
    return(lambda2)
}

.getEventProbabilitiesGroupwise <- function(..., time, accrualTimeVector, accrualIntensity, lambda2,
        lambda1, piecewiseSurvivalTime, phi, kappa, allocationRatioPlanned, hazardRatio) {
    .assertIsSingleNumber(time, "time")

    if (length(accrualTimeVector) > 1 && accrualTimeVector[1] == 0) {
        accrualTimeVector <- accrualTimeVector[2:length(accrualTimeVector)]
    }

    accrualTimeVectorLength <- length(accrualTimeVector)
    densityIntervals <- accrualTimeVector
    if (accrualTimeVectorLength > 1) {
        densityIntervals[2:accrualTimeVectorLength] <-
            accrualTimeVector[2:accrualTimeVectorLength] -
            accrualTimeVector[1:(accrualTimeVectorLength - 1)]
    }

    if (length(densityIntervals) > 1 && length(accrualIntensity) > 1 &&
            length(densityIntervals) != length(accrualIntensity)) {
        stop(
            C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'densityIntervals' (", .arrayToString(densityIntervals),
            ") and 'accrualIntensity' (", .arrayToString(accrualIntensity), ") must have same length"
        )
    }

    densityVector <- accrualIntensity / sum(densityIntervals * accrualIntensity)

    eventProbs <- rep(NA_real_, 2)

    for (k in 1:accrualTimeVectorLength) {
        if (time <= accrualTimeVector[k]) {
            for (groupNumber in c(1, 2)) { # two groups: 1 = treatment, 2 = control

                lambdaTemp <- .getLambda(
                    groupNumber           = groupNumber,
                    lambda2               = lambda2,
                    lambda1               = lambda1,
                    hazardRatio           = hazardRatio,
                    kappa                 = kappa
                )

                inner <- function(x) {
                    .getEventProbabilityFunctionVec(
                        timeVector = x, piecewiseLambda = lambdaTemp,
                        piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi[groupNumber], kappa = kappa
                    )
                }
                timeValue1 <- 0
                if (k > 1) {
                    timeValue1 <- time - accrualTimeVector[1]
                }

                eventProbs[groupNumber] <- densityVector[1] * integrate(inner, timeValue1, time)$value

                if (k > 2) {
                    for (j in 2:(k - 1)) {
                        eventProbs[groupNumber] <- eventProbs[groupNumber] +
                            densityVector[j] * integrate(
                                inner, time - accrualTimeVector[j],
                                time - accrualTimeVector[j - 1]
                            )$value
                    }
                }
                if (k > 1) {
                    eventProbs[groupNumber] <- eventProbs[groupNumber] +
                        densityVector[k] * integrate(inner, 0, time - accrualTimeVector[k - 1])$value
                }
            }

            return(eventProbs)
        }
    }

    for (groupNumber in c(1, 2)) {
        lambdaTemp <- .getLambda(
            groupNumber           = groupNumber,
            lambda2               = lambda2,
            lambda1               = lambda1,
            hazardRatio           = hazardRatio,
            kappa                 = kappa
        )

        inner <- function(x) {
            .getEventProbabilityFunctionVec(
                timeVector = x, piecewiseLambda = lambdaTemp,
                piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi[groupNumber], kappa = kappa
            )
        }

        eventProbs[groupNumber] <- densityVector[1] *
            integrate(inner, time - accrualTimeVector[1], time)$value
        if (accrualTimeVectorLength > 1) {
            for (j in (2:accrualTimeVectorLength)) {
                eventProbs[groupNumber] <- eventProbs[groupNumber] +
                    densityVector[j] * integrate(
                        inner, time - accrualTimeVector[j],
                        time - accrualTimeVector[j - 1]
                    )$value
            }
        }
    }

    return(eventProbs)
}

.getEventProbabilitiesOverall <- function(eventProbs, allocationRatioPlanned) {
    return((allocationRatioPlanned * eventProbs[1] + eventProbs[2]) / (1 + allocationRatioPlanned))
}

.getEventProbabilities <- function(..., time, accrualTimeVector, accrualIntensity, lambda2,
        lambda1, piecewiseSurvivalTime, phi, kappa, allocationRatioPlanned, hazardRatio) {
    eventProbs <- .getEventProbabilitiesGroupwise(
        time = time, accrualTimeVector = accrualTimeVector,
        accrualIntensity = accrualIntensity, lambda2 = lambda2,
        lambda1 = lambda1, piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi,
        kappa = kappa, allocationRatioPlanned = allocationRatioPlanned, hazardRatio = hazardRatio
    )

    return(.getEventProbabilitiesOverall(eventProbs, allocationRatioPlanned))
}

.getEventsFixed <- function(..., typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
        twoSidedPower, alpha, beta, sided, hazardRatio, thetaH0, allocationRatioPlanned) {
    typeOfComputation <- match.arg(typeOfComputation)

    if (typeOfComputation == "Schoenfeld") {
        eventsFixed <- (.getOneMinusQNorm(alpha / sided) + .getOneMinusQNorm(beta))^2 /
            (log(hazardRatio) - log(thetaH0))^2 *
            (1 + allocationRatioPlanned)^2 / allocationRatioPlanned
        if (twoSidedPower && (sided == 2)) {
            up <- 2 * eventsFixed
            eventsFixed <- .getOneDimensionalRoot(
                function(n) {
                    return(stats::pnorm(.getOneMinusQNorm(alpha / 2) - sqrt(n) *
                        (log(hazardRatio) - log(thetaH0)) * sqrt(allocationRatioPlanned) /
                        (1 + allocationRatioPlanned)) -
                        stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(n) *
                            (log(hazardRatio) - log(thetaH0)) * sqrt(allocationRatioPlanned) /
                            (1 + allocationRatioPlanned)) - beta)
                },
                lower = 0.001, upper = up, tolerance = 1e-04,
                callingFunctionInformation = ".getEventsFixed"
            )
        }
        return(eventsFixed)
    }

    if (typeOfComputation == "Freedman") {
        eventsFixed <- (.getOneMinusQNorm(alpha / sided) + .getOneMinusQNorm(beta))^2 *
            (1 + hazardRatio * allocationRatioPlanned)^2 / (1 - hazardRatio)^2 /
            allocationRatioPlanned
        if (twoSidedPower && (sided == 2)) {
            up <- 2 * eventsFixed
            eventsFixed <- .getOneDimensionalRoot(
                function(n) {
                    return(stats::pnorm(.getOneMinusQNorm(alpha / 2) - sqrt(n) *
                        sqrt(allocationRatioPlanned) * (1 - hazardRatio) /
                        (1 + allocationRatioPlanned * hazardRatio)) -
                        stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(n) *
                            sqrt(allocationRatioPlanned) * (1 - hazardRatio) /
                            (1 + allocationRatioPlanned * hazardRatio)) - beta)
                },
                lower = 0.001, upper = up, tolerance = 1e-04,
                callingFunctionInformation = ".getEventsFixed"
            )
        }
        return(eventsFixed)
    }

    if (typeOfComputation == "HsiehFreedman") {
        eventsFixed <- (.getOneMinusQNorm(alpha / sided) + .getOneMinusQNorm(beta))^2 *
            (1 + hazardRatio)^2 / (1 - hazardRatio)^2 *
            (1 + allocationRatioPlanned)^2 / (4 * allocationRatioPlanned)
        if (twoSidedPower && sided == 2) {
            up <- 2 * eventsFixed
            eventsFixed <- .getOneDimensionalRoot(
                function(n) {
                    return(stats::pnorm(.getOneMinusQNorm(alpha / 2) - sqrt(n) *
                        2 * sqrt(allocationRatioPlanned) / (1 + allocationRatioPlanned) *
                        (1 - hazardRatio) / (1 + hazardRatio)) -
                        stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(n) *
                            2 * sqrt(allocationRatioPlanned) / (1 + allocationRatioPlanned) *
                            (1 - hazardRatio) / (1 + hazardRatio)) - beta)
                },
                lower = 0.001, upper = up, tolerance = 1e-04,
                callingFunctionInformation = ".getEventsFixed"
            )
        }
        return(eventsFixed)
    }
}

.getSampleSizeFixedSurvival <- function(designPlan) {
    alpha <- designPlan$getAlpha()
    beta <- designPlan$getBeta()
    sided <- designPlan$getSided()
    twoSidedPower <- designPlan$getTwoSidedPower()
    typeOfComputation <- designPlan$typeOfComputation
    thetaH0 <- designPlan$thetaH0
    pi1 <- designPlan$pi1
    pi2 <- designPlan$pi2
    allocationRatioPlanned <- designPlan$allocationRatioPlanned
    accountForObservationTimes <- designPlan$accountForObservationTimes
    accrualTime <- designPlan$accrualTime
    kappa <- designPlan$kappa
    piecewiseSurvivalTime <- designPlan$piecewiseSurvivalTime
    maxNumberOfSubjects <- designPlan$maxNumberOfSubjects
    hazardRatio <- designPlan$hazardRatio

    .assertIsValidHazardRatio(hazardRatio, thetaH0)

    if (designPlan$.piecewiseSurvivalTime$.isLambdaBased()) {
        numberOfResults <- length(hazardRatio)
    } else {
        numberOfResults <- length(pi1)
    }

    designPlan$eventsFixed <- rep(NA_real_, numberOfResults) # number of events
    designPlan$nFixed <- rep(NA_real_, numberOfResults) # number of subjects
    designPlan$chi <- rep(NA_real_, numberOfResults) # probability of an event

    calculateAllocationRatioPlanned <- FALSE
    if (allocationRatioPlanned == 0) {
        allocationRatioPlannedVec <- rep(NA_real_, numberOfResults)
        calculateAllocationRatioPlanned <- TRUE
        designPlan$optimumAllocationRatio <- TRUE
        designPlan$.setParameterType("optimumAllocationRatio", C_PARAM_USER_DEFINED)
    }

    userDefinedMaxNumberOfSubjects <- .isUserDefinedMaxNumberOfSubjects(designPlan)
    if (userDefinedMaxNumberOfSubjects && allocationRatioPlanned == 0) {
        stop(
            C_EXCEPTION_TYPE_CONFLICTING_ARGUMENTS,
            "determination of optimum allocation ('allocationRatioPlanned' = 0) not possible ",
            "for given 'maxNumberOfSubjects' (", designPlan$maxNumberOfSubjects, ")"
        )
    }

    if (userDefinedMaxNumberOfSubjects) {
        timeVector <- rep(NA_real_, numberOfResults)
    }

    designPlan$.calculateFollowUpTime <- FALSE

    lambda1 <- designPlan$lambda1
    if (designPlan$.piecewiseSurvivalTime$piecewiseSurvivalEnabled) {
        lambda1 <- rep(NA_real_, numberOfResults)
    }

    for (i in 1:numberOfResults) {
        phi <- -c(
            log(1 - designPlan$dropoutRate1),
            log(1 - designPlan$dropoutRate2)
        ) / designPlan$dropoutTime

        if (!userDefinedMaxNumberOfSubjects) {
            if (calculateAllocationRatioPlanned) {
                # allocationRatioPlanned = 0 provides optimum sample size
                allocationRatioPlanned <- stats::optimize(function(x) {
                    numberEvents <- .getEventsFixed(
                        typeOfComputation = typeOfComputation, twoSidedPower = twoSidedPower,
                        alpha = alpha, beta = beta, sided = sided, hazardRatio = hazardRatio[i],
                        thetaH0 = thetaH0, allocationRatioPlanned = x
                    )

                    if (!accountForObservationTimes) {
                        probEvent <- (x * pi1[i] + pi2) / (1 + x)
                    } else {
                        probEvent <- .getEventProbabilities(
                            time = accrualTime[length(accrualTime)] + designPlan$followUpTime,
                            accrualTimeVector = accrualTime,
                            accrualIntensity = designPlan$accrualIntensity,
                            lambda2 = designPlan$lambda2, lambda1 = lambda1[i],
                            piecewiseSurvivalTime = piecewiseSurvivalTime,
                            phi = phi, kappa = kappa, allocationRatioPlanned = x,
                            hazardRatio = hazardRatio[i]
                        )
                    }
                    return(numberEvents / probEvent)
                }, interval = c(0, 5), tol = 0.0001)$minimum
                allocationRatioPlannedVec[i] <- allocationRatioPlanned
            }

            designPlan$eventsFixed[i] <- .getEventsFixed(
                typeOfComputation = typeOfComputation, twoSidedPower = twoSidedPower,
                alpha = alpha, beta = beta, sided = sided, hazardRatio = hazardRatio[i],
                thetaH0 = thetaH0, allocationRatioPlanned = allocationRatioPlanned
            )

            if (!accountForObservationTimes) {
                designPlan$chi[i] <- (allocationRatioPlanned * pi1[i] + pi2) /
                    (1 + allocationRatioPlanned)
            } else {
                designPlan$chi[i] <- .getEventProbabilities(
                    time = accrualTime[length(accrualTime)] + designPlan$followUpTime,
                    accrualTimeVector = accrualTime,
                    accrualIntensity = designPlan$accrualIntensity,
                    lambda2 = designPlan$lambda2,
                    lambda1 = lambda1[i],
                    piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa,
                    allocationRatioPlanned = allocationRatioPlanned, hazardRatio = hazardRatio[i]
                )
            }
            designPlan$.setParameterType("chi", C_PARAM_GENERATED)

            designPlan$nFixed[i] <- designPlan$eventsFixed[i] / designPlan$chi[i]
        } else {
            if (length(maxNumberOfSubjects) > 1) {
                stop(
                    C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                    "length of user defined 'maxNumberOfSubjects' (",
                    .arrayToString(maxNumberOfSubjects), ") must be 1"
                )
            }

            designPlan$.calculateFollowUpTime <- TRUE

            designPlan$eventsFixed[i] <- .getEventsFixed(
                typeOfComputation = typeOfComputation, twoSidedPower = twoSidedPower,
                alpha = alpha, beta = beta, sided = sided, hazardRatio = hazardRatio[i],
                thetaH0 = thetaH0, allocationRatioPlanned = allocationRatioPlanned
            )

            designPlan$nFixed[i] <- maxNumberOfSubjects
            if (designPlan$eventsFixed[i] > maxNumberOfSubjects) {
                if (length(hazardRatio) > 1) {
                    stop(
                        C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                        sprintf(
                            paste0(
                                "'maxNumberOfSubjects' (%s) is smaller than the number ",
                                "of events (%.3f) at index %s (hazard ratio = %.3f)"
                            ),
                            maxNumberOfSubjects, designPlan$eventsFixed[i], i, hazardRatio[i]
                        )
                    )
                } else {
                    stop(
                        C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                        sprintf(
                            paste0(
                                "'maxNumberOfSubjects' (%s) is smaller than the number ",
                                "of events (%.3f)"
                            ),
                            maxNumberOfSubjects, designPlan$eventsFixed[i]
                        )
                    )
                }
            }

            up <- 2
            iterate <- 1
            while (designPlan$eventsFixed[i] / .getEventProbabilities(
                time = up, accrualTimeVector = accrualTime, accrualIntensity = designPlan$accrualIntensity,
                lambda2 = designPlan$lambda2, lambda1 = lambda1[i],
                piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa,
                allocationRatioPlanned = allocationRatioPlanned,
                hazardRatio = hazardRatio[i]
            ) > maxNumberOfSubjects) {
                up <- 2 * up
                iterate <- iterate + 1
                if (iterate > 50) {
                    stop(
                        C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                        "the number of subjects is too small to reach maximum number of events ",
                        "(presumably due to drop-out rates), search algorithm failed"
                    )
                }
            }

            timeVector[i] <- .getOneDimensionalRoot(
                function(x) {
                    designPlan$eventsFixed[i] / .getEventProbabilities(
                        time = x, accrualTimeVector = accrualTime,
                        accrualIntensity = designPlan$accrualIntensity, lambda2 = designPlan$lambda2,
                        lambda1 = lambda1[i],
                        piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa,
                        allocationRatioPlanned = allocationRatioPlanned,
                        hazardRatio = hazardRatio[i]
                    ) - maxNumberOfSubjects
                },
                lower = 0, upper = up, tolerance = 1e-06, acceptResultsOutOfTolerance = TRUE,
                callingFunctionInformation = ".getSampleSizeSequentialSurvival"
            )

            if (!is.na(timeVector[i])) {
                designPlan$chi[i] <- .getEventProbabilities(
                    time = timeVector[i],
                    accrualTimeVector = accrualTime,
                    accrualIntensity = designPlan$accrualIntensity, lambda2 = designPlan$lambda2,
                    lambda1 = lambda1[i],
                    piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa,
                    allocationRatioPlanned = allocationRatioPlanned, hazardRatio = hazardRatio[i]
                )
                designPlan$.setParameterType("chi", C_PARAM_GENERATED)
            }
        }
    }

    if (calculateAllocationRatioPlanned) {
        allocationRatioPlanned <- allocationRatioPlannedVec
        designPlan$allocationRatioPlanned <- allocationRatioPlanned
        designPlan$.setParameterType("allocationRatioPlanned", C_PARAM_GENERATED)
    }

    if (userDefinedMaxNumberOfSubjects) {
        designPlan$followUpTime <- timeVector - accrualTime[length(accrualTime)]
        designPlan$.setParameterType("followUpTime", C_PARAM_GENERATED)
    }

    designPlan$nFixed2 <- designPlan$nFixed / (1 + allocationRatioPlanned)
    designPlan$nFixed1 <- designPlan$nFixed2 * allocationRatioPlanned

    if (designPlan$.design$kMax == 1 &&
            designPlan$.accrualTime$.isRelativeAccrualIntensity(designPlan$accrualIntensity)) {
        designPlan$accrualIntensity <- designPlan$nFixed / designPlan$accrualTime
        designPlan$.setParameterType("accrualIntensity", C_PARAM_GENERATED)
    }

    designPlan$numberOfSubjects1 <- matrix(designPlan$nFixed1, nrow = 1)
    designPlan$numberOfSubjects2 <- matrix(designPlan$nFixed2, nrow = 1)

    if (!designPlan$.piecewiseSurvivalTime$piecewiseSurvivalEnabled) {
        eventRatio <- allocationRatioPlanned * pi1 / pi2
    } else {
        eventRatio <- NA_real_
    }

    # Fixed
    designPlan$hazardRatio <- hazardRatio
    designPlan$expectedEventsH1 <- designPlan$eventsFixed
    designPlan$maxNumberOfSubjects <- designPlan$nFixed
    designPlan$numberOfSubjects <- matrix(designPlan$nFixed, nrow = 1)

    designPlan$.setParameterType("eventsFixed", C_PARAM_GENERATED)
    designPlan$.setParameterType("nFixed1", C_PARAM_GENERATED)
    designPlan$.setParameterType("nFixed2", C_PARAM_GENERATED)
    designPlan$.setParameterType("nFixed", C_PARAM_GENERATED)

    if (designPlan$accountForObservationTimes) {
        designPlan$analysisTime <- matrix(accrualTime[length(accrualTime)] +
            designPlan$followUpTime, nrow = 1)
        designPlan$.setParameterType("analysisTime", C_PARAM_GENERATED)
    }
    return(designPlan)
}

# note that fixed sample size must be calculated before on 'designPlan'
.getSampleSizeSequentialSurvival <- function(designPlan, designCharacteristics) {
    if (designPlan$.piecewiseSurvivalTime$.isLambdaBased()) {
        numberOfResults <- length(designPlan$hazardRatio)
    } else {
        numberOfResults <- length(designPlan$pi1)
    }

    kMax <- designCharacteristics$.design$kMax
    designPlan$eventsPerStage <- matrix(NA_real_, kMax, numberOfResults)
    analysisTime <- matrix(NA_real_, kMax, numberOfResults)
    numberOfSubjects <- matrix(NA_real_, kMax, numberOfResults)
    designPlan$expectedEventsH0 <- rep(NA_real_, numberOfResults)
    designPlan$expectedEventsH01 <- rep(NA_real_, numberOfResults)
    designPlan$expectedEventsH1 <- rep(NA_real_, numberOfResults)
    expectedNumberOfSubjectsH1 <- rep(NA_real_, numberOfResults)
    studyDuration <- rep(NA_real_, numberOfResults)
    designPlan$chi <- rep(NA_real_, numberOfResults)

    informationRates <- designCharacteristics$information / designCharacteristics$shift

    lambda1 <- designPlan$lambda1
    if (designPlan$.piecewiseSurvivalTime$piecewiseSurvivalEnabled) {
        lambda1 <- rep(NA_real_, numberOfResults)
    }

    if (designPlan$accountForObservationTimes && designPlan$.calculateFollowUpTime) {
        designPlan$followUpTime <- rep(NA_real_, numberOfResults)
    }

    for (i in 1:numberOfResults) {
        designPlan$eventsPerStage[, i] <- designPlan$eventsFixed[i] * informationRates *
            designCharacteristics$inflationFactor

        if (!designPlan$accountForObservationTimes) {
            if (length(designPlan$allocationRatioPlanned) > 1) {
                allocationRatioPlanned <- designPlan$allocationRatioPlanned[i]
            } else {
                allocationRatioPlanned <- designPlan$allocationRatioPlanned
            }
            designPlan$chi[i] <- (allocationRatioPlanned * designPlan$pi1[i] + designPlan$pi2) /
                (1 + allocationRatioPlanned)
            designPlan$.setParameterType("chi", C_PARAM_GENERATED)
            numberOfSubjects[kMax, i] <- designPlan$eventsPerStage[kMax, i] / designPlan$chi[i]
        } else {
            phi <- -c(log(1 - designPlan$dropoutRate1), log(1 - designPlan$dropoutRate2)) /
                designPlan$dropoutTime

            if (designPlan$.calculateFollowUpTime) {
                if (designPlan$eventsPerStage[kMax, i] > designPlan$maxNumberOfSubjects[i]) {
                    stop(
                        C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                        sprintf(
                            paste0(
                                "the number of subjects (%s) is smaller than the number ",
                                "of events (%s) at stage %s"
                            ),
                            designPlan$maxNumberOfSubjects[i],
                            designPlan$eventsPerStage[kMax, i], i
                        )
                    )
                }

                up <- 2
                iterate <- 1
                while (designPlan$eventsPerStage[kMax, i] / .getEventProbabilities(
                    time = up,
                    accrualTimeVector = designPlan$accrualTime,
                    accrualIntensity = designPlan$accrualIntensity,
                    lambda2 = designPlan$lambda2,
                    lambda1 = lambda1[i],
                    piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
                    phi = phi, kappa = designPlan$kappa,
                    allocationRatioPlanned = designPlan$allocationRatioPlanned,
                    hazardRatio = designPlan$hazardRatio[i]
                ) > designPlan$maxNumberOfSubjects[i]) {
                    up <- 2 * up
                    iterate <- iterate + 1
                    if (iterate > 50) {
                        stop(
                            C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                            "the number of subjects is too small to reach maximum number of events ",
                            "(presumably due to drop-out rates)"
                        )
                    }
                }

                totalTime <- .getOneDimensionalRoot(
                    function(x) {
                        designPlan$eventsPerStage[kMax, i] / designPlan$maxNumberOfSubjects[i] -
                            .getEventProbabilities(
                                time = x, accrualTimeVector = designPlan$accrualTime,
                                accrualIntensity = designPlan$accrualIntensity,
                                lambda2 = designPlan$lambda2,
                                lambda1 = lambda1[i],
                                piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
                                phi = phi, kappa = designPlan$kappa,
                                allocationRatioPlanned = designPlan$allocationRatioPlanned,
                                hazardRatio = designPlan$hazardRatio[i]
                            )
                    },
                    lower = 0, upper = up, tolerance = 1e-06,
                    callingFunctionInformation = ".getSampleSizeSequentialSurvival"
                )

                # analysis times
                for (j in 1:kMax) {
                    analysisTime[j, i] <- .getOneDimensionalRoot(
                        function(x) {
                            designPlan$eventsPerStage[j, i] / designPlan$maxNumberOfSubjects[i] -
                                .getEventProbabilities(
                                    time = x, accrualTimeVector = designPlan$accrualTime,
                                    accrualIntensity = designPlan$accrualIntensity,
                                    lambda2 = designPlan$lambda2,
                                    lambda1 = lambda1[i],
                                    piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
                                    phi = phi, kappa = designPlan$kappa,
                                    allocationRatioPlanned = designPlan$allocationRatioPlanned,
                                    hazardRatio = designPlan$hazardRatio[i]
                                )
                        },
                        lower = 0, upper = totalTime, tolerance = 1e-06, acceptResultsOutOfTolerance = TRUE,
                        callingFunctionInformation = ".getSampleSizeSequentialSurvival"
                    )
                }
                analysisTime[kMax, i] <- totalTime

                designPlan$followUpTime[i] <- totalTime -
                    designPlan$accrualTime[length(designPlan$accrualTime)]

                numberOfSubjects[, i] <- .getNumberOfSubjects(
                    time = analysisTime[, i],
                    accrualTime = designPlan$accrualTime,
                    accrualIntensity = designPlan$accrualIntensity,
                    maxNumberOfSubjects = designPlan$maxNumberOfSubjects[i]
                )
            } else {
                if (length(designPlan$allocationRatioPlanned) > 1) {
                    allocationRatioPlanned <- designPlan$allocationRatioPlanned[i]
                } else {
                    allocationRatioPlanned <- designPlan$allocationRatioPlanned
                }

                if (is.na(designPlan$followUpTime)) {
                    stop(
                        C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                        "'followUpTime' must be defined because 'designPlan$.calculateFollowUpTime' = FALSE"
                    )
                }

                designPlan$chi[i] <- .getEventProbabilities(
                    time = designPlan$accrualTime[length(designPlan$accrualTime)] + designPlan$followUpTime,
                    accrualTimeVector = designPlan$accrualTime,
                    accrualIntensity = designPlan$accrualIntensity,
                    lambda2 = designPlan$lambda2,
                    lambda1 = lambda1[i],
                    piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
                    phi = phi, kappa = designPlan$kappa,
                    allocationRatioPlanned = allocationRatioPlanned,
                    hazardRatio = designPlan$hazardRatio[i]
                )
                designPlan$.setParameterType("chi", C_PARAM_GENERATED)
                numberOfSubjects[kMax, i] <- designPlan$eventsPerStage[kMax, i] / designPlan$chi[i]

                #   Analysis times
                for (j in 1:(kMax - 1)) {
                    analysisTime[j, i] <- .getOneDimensionalRoot(
                        function(x) {
                            designPlan$eventsPerStage[j, i] / numberOfSubjects[kMax, i] -
                                .getEventProbabilities(
                                    time = x, accrualTimeVector = designPlan$accrualTime,
                                    accrualIntensity = designPlan$accrualIntensity,
                                    lambda2 = designPlan$lambda2,
                                    lambda1 = lambda1[i],
                                    piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
                                    phi = phi, kappa = designPlan$kappa,
                                    allocationRatioPlanned = allocationRatioPlanned,
                                    hazardRatio = designPlan$hazardRatio[i]
                                )
                        },
                        lower = 0, upper = designPlan$accrualTime[length(designPlan$accrualTime)] +
                            designPlan$followUpTime, tolerance = 1e-06,
                        callingFunctionInformation = ".getSampleSizeSequentialSurvival"
                    )
                }
                analysisTime[kMax, i] <- designPlan$accrualTime[length(designPlan$accrualTime)] +
                    designPlan$followUpTime

                numberOfSubjects[, i] <- .getNumberOfSubjects(
                    time = analysisTime[, i],
                    accrualTime = designPlan$accrualTime, accrualIntensity = designPlan$accrualIntensity,
                    maxNumberOfSubjects = numberOfSubjects[kMax, i]
                )
            }

            stoppingProbs <- designCharacteristics$rejectionProbabilities +
                c(designCharacteristics$futilityProbabilities, 0)

            if (all(is.na(designCharacteristics$futilityProbabilities))) {
                warning("Expected number of subjects H1 and study duration H1 ",
                    "cannot be calculated because the futility probabilities ",
                    "are not applicable for the specified design",
                    call. = FALSE
                )
            }

            stoppingProbs[kMax] <- 1 - sum(stoppingProbs[1:(kMax - 1)])

            studyDuration[i] <- analysisTime[, i] %*% stoppingProbs

            expectedNumberOfSubjectsH1[i] <- numberOfSubjects[, i] %*% stoppingProbs
        }

        designPlan$expectedEventsH0[i] <- designCharacteristics$averageSampleNumber0 *
            designPlan$eventsFixed[i]
        designPlan$expectedEventsH01[i] <- designCharacteristics$averageSampleNumber01 *
            designPlan$eventsFixed[i]
        designPlan$expectedEventsH1[i] <- designCharacteristics$averageSampleNumber1 *
            designPlan$eventsFixed[i]
        designPlan$.setParameterType("expectedEventsH0", C_PARAM_GENERATED)
        designPlan$.setParameterType("expectedEventsH01", C_PARAM_GENERATED)
        designPlan$.setParameterType("expectedEventsH1", C_PARAM_GENERATED)

        designPlan$numberOfSubjects2 <- numberOfSubjects / (1 + designPlan$allocationRatioPlanned)
        designPlan$numberOfSubjects1 <-
            designPlan$numberOfSubjects2 * designPlan$allocationRatioPlanned
        designPlan$.setParameterType("numberOfSubjects1", C_PARAM_GENERATED)
        designPlan$.setParameterType("numberOfSubjects2", C_PARAM_GENERATED)
    }

    if (!is.null(designCharacteristics$rejectionProbabilities)) {
        designPlan$rejectPerStage <- matrix(designCharacteristics$rejectionProbabilities,
            nrow = designPlan$.design$kMax
        )
        designPlan$.setParameterType("rejectPerStage", C_PARAM_GENERATED)

        designPlan$earlyStop <- sum(designPlan$rejectPerStage[1:(designPlan$.design$kMax - 1), ])
        designPlan$.setParameterType("earlyStop", C_PARAM_GENERATED)
    }

    if (!is.null(designCharacteristics$futilityProbabilities) &&
            any(designPlan$.design$futilityBounds != C_FUTILITY_BOUNDS_DEFAULT)) {
        designPlan$futilityPerStage <- matrix(designCharacteristics$futilityProbabilities,
            nrow = designPlan$.design$kMax - 1
        )
        designPlan$.setParameterType("futilityPerStage", C_PARAM_GENERATED)

        designPlan$futilityStop <- sum(designPlan$futilityPerStage)
        designPlan$.setParameterType("futilityStop", C_PARAM_GENERATED)

        designPlan$earlyStop <- designPlan$earlyStop + sum(designPlan$futilityPerStage)
    }

    designPlan$informationRates <- matrix(informationRates, ncol = 1)
    if (!is.matrix(numberOfSubjects)) {
        designPlan$numberOfSubjects <- matrix(numberOfSubjects[kMax, ], nrow = 1)
    } else {
        designPlan$numberOfSubjects <- numberOfSubjects
    }

    designPlan$maxNumberOfSubjects <- designPlan$numberOfSubjects[kMax, ]
    if (designPlan$.getParameterType("maxNumberOfSubjects") == C_PARAM_NOT_APPLICABLE ||
            length(designPlan$maxNumberOfSubjects) > 1) {
        designPlan$.setParameterType("maxNumberOfSubjects", C_PARAM_GENERATED)
    }

    designPlan$maxNumberOfSubjects1 <- .getNumberOfSubjects1(
        designPlan$maxNumberOfSubjects, designPlan$allocationRatioPlanned
    )
    designPlan$maxNumberOfSubjects2 <- .getNumberOfSubjects2(
        designPlan$maxNumberOfSubjects, designPlan$allocationRatioPlanned
    )
    designPlan$.setParameterType("maxNumberOfSubjects1", C_PARAM_GENERATED)
    designPlan$.setParameterType("maxNumberOfSubjects2", C_PARAM_GENERATED)

    if (ncol(designPlan$informationRates) == 1 &&
            identical(designPlan$informationRates[, 1], designPlan$.design$informationRates)) {
        designPlan$.setParameterType("informationRates", C_PARAM_NOT_APPLICABLE)
    } else {
        designPlan$.setParameterType("informationRates", C_PARAM_GENERATED)
    }
    designPlan$.setParameterType("numberOfSubjects", C_PARAM_GENERATED)
    designPlan$.setParameterType("eventsPerStage", C_PARAM_GENERATED)

    if (designPlan$accountForObservationTimes) {
        designPlan$analysisTime <- analysisTime
        designPlan$expectedNumberOfSubjectsH1 <- expectedNumberOfSubjectsH1
        designPlan$studyDuration <- studyDuration
        designPlan$studyDurationH1 <- studyDuration # deprecated

        designPlan$.setParameterType("analysisTime", C_PARAM_GENERATED)
        designPlan$.setParameterType("expectedNumberOfSubjectsH1", C_PARAM_GENERATED)
        designPlan$.setParameterType("studyDuration", C_PARAM_GENERATED)
    }

    designPlan$.setParameterType("eventsFixed", C_PARAM_NOT_APPLICABLE)
    designPlan$.setParameterType("nFixed1", C_PARAM_NOT_APPLICABLE)
    designPlan$.setParameterType("nFixed2", C_PARAM_NOT_APPLICABLE)
    designPlan$.setParameterType("nFixed", C_PARAM_NOT_APPLICABLE)

    if (designPlan$allocationRatioPlanned[1] == 1) {
        designPlan$.setParameterType("numberOfSubjects1", C_PARAM_NOT_APPLICABLE)
        designPlan$.setParameterType("numberOfSubjects2", C_PARAM_NOT_APPLICABLE)
    }

    designPlan$.calculateFollowUpTime <- NA

    return(designPlan)
}

# Note that 'directionUpper' and 'maxNumberOfSubjects' are only applicable
# for 'objectType' = "power"
.createDesignPlanMeans <- function(..., objectType = c("sampleSize", "power"),
        design, normalApproximation = FALSE, meanRatio = FALSE,
        thetaH0 = ifelse(meanRatio, 1, 0), alternative = NA_real_,
        stDev = C_STDEV_DEFAULT, directionUpper = NA,
        maxNumberOfSubjects = NA_real_, groups = 2, allocationRatioPlanned = NA_real_) {
    objectType <- match.arg(objectType)

    .assertIsTrialDesignInverseNormalOrGroupSequential(design)
    .assertIsValidAlphaAndBeta(design$alpha, design$beta)
    .assertIsValidSidedParameter(design$sided)
    .assertIsValidStandardDeviation(stDev)
    .assertIsValidGroupsParameter(groups)
    .assertIsSingleNumber(thetaH0, "thetaH0")
    .assertIsSingleLogical(meanRatio, "meanRatio")
    .assertIsValidThetaH0(thetaH0, endpoint = "means", groups = groups, ratioEnabled = meanRatio)
    .assertIsSingleLogical(normalApproximation, "normalApproximation")

    if (meanRatio) {
        if (identical(alternative, C_ALTERNATIVE_POWER_SIMULATION_DEFAULT)) {
            alternative <- C_ALTERNATIVE_POWER_SIMULATION_MEAN_RATIO_DEFAULT
        }
        .assertIsInOpenInterval(alternative, "alternative", 0, NULL, naAllowed = TRUE)
    }

    .assertIsValidDirectionUpper(directionUpper, design$sided, objectType, userFunctionCallEnabled = TRUE)

    if (objectType == "sampleSize" && !any(is.na(alternative))) {
        if (design$sided == 1 && any(alternative - thetaH0 <= 0)) {
            stop(
                C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                "any 'alternative' (", .arrayToString(alternative),
                ") must be > 'thetaH0' (", thetaH0, ")"
            )
        }

        if (any(alternative - thetaH0 == 0)) {
            stop(
                C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                "any 'alternative' (", .arrayToString(alternative),
                ") must be != 'thetaH0' (", thetaH0, ")"
            )
        }
    }

    designPlan <- TrialDesignPlanMeans(design = design, meanRatio = meanRatio)
    designPlan$.setSampleSizeObject(objectType)

    designPlan$criticalValuesPValueScale <- matrix(design$stageLevels, ncol = 1)
    if (design$sided == 2) {
        designPlan$criticalValuesPValueScale <- designPlan$criticalValuesPValueScale * 2
        designPlan$.setParameterType("criticalValuesPValueScale", C_PARAM_GENERATED)
    }

    if (any(design$futilityBounds > C_FUTILITY_BOUNDS_DEFAULT)) {
        designPlan$futilityBoundsPValueScale <- matrix(1 - stats::pnorm(design$futilityBounds), ncol = 1)
        designPlan$.setParameterType("futilityBoundsPValueScale", C_PARAM_GENERATED)
    }

    if (groups == 2) {
        if (design$sided == 2 && ((thetaH0 != 0 && !meanRatio) || (thetaH0 != 1 && meanRatio))) {
            stop(
                C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                "two-sided case is implemented only for superiority testing (i.e., thetaH0 = ", ifelse(meanRatio, 1, 0), ")"
            )
        }

        if (is.na(allocationRatioPlanned)) {
            allocationRatioPlanned <- C_ALLOCATION_RATIO_DEFAULT
        }

        if (allocationRatioPlanned < 0) {
            stop(
                C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                "'allocationRatioPlanned' (", allocationRatioPlanned, ") must be >= 0"
            )
        }

        .setValueAndParameterType(designPlan, "allocationRatioPlanned", allocationRatioPlanned, 1)

        if (meanRatio && thetaH0 <= 0) {
            stop(
                C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                "null hypothesis mean ratio is not allowed be negative or zero, ",
                "i.e., 'thetaH0' must be > 0 if 'meanRatio' = TRUE"
            )
        }
    }

    .setValueAndParameterType(designPlan, "normalApproximation", normalApproximation, FALSE)
    .setValueAndParameterType(designPlan, "meanRatio", meanRatio, FALSE)
    .setValueAndParameterType(designPlan, "thetaH0", thetaH0, 0)
    if (objectType == "power") {
        .setValueAndParameterType(
            designPlan, "alternative", alternative,
            C_ALTERNATIVE_POWER_SIMULATION_DEFAULT
        )
    } else {
        .setValueAndParameterType(designPlan, "alternative", alternative, C_ALTERNATIVE_DEFAULT)
    }
    .setValueAndParameterType(designPlan, "stDev", stDev, C_STDEV_DEFAULT)
    if (objectType == "power") {
        .assertIsValidMaxNumberOfSubjects(maxNumberOfSubjects)
        .setValueAndParameterType(designPlan, "maxNumberOfSubjects", maxNumberOfSubjects, NA_real_)
        .setValueAndParameterType(designPlan, "directionUpper", directionUpper, TRUE)

        designPlan$.setParameterType("effect", C_PARAM_GENERATED)
    }
    .setValueAndParameterType(designPlan, "groups", groups, 2)
    if (groups == 1) {
        if (isTRUE(meanRatio)) {
            warning("'meanRatio' (", meanRatio, ") will be ignored ",
                "because it is not applicable for 'groups' = 1",
                call. = FALSE
            )
        }
        designPlan$.setParameterType("meanRatio", C_PARAM_NOT_APPLICABLE)

        if (length(allocationRatioPlanned) == 1 && !is.na(allocationRatioPlanned)) {
            warning("'allocationRatioPlanned' (", allocationRatioPlanned,
                ") will be ignored because it is not applicable for 'groups' = 1",
                call. = FALSE
            )
        }
        designPlan$.setParameterType("allocationRatioPlanned", C_PARAM_NOT_APPLICABLE)
    }

    return(designPlan)
}

#
# note that 'directionUpper' and 'maxNumberOfSubjects' are
# only applicable for 'objectType' = "power"
#
.createDesignPlanRates <- function(..., objectType = c("sampleSize", "power"),
        design, normalApproximation = TRUE, riskRatio = FALSE,
        thetaH0 = ifelse(riskRatio, 1, 0), pi1 = C_PI_1_SAMPLE_SIZE_DEFAULT,
        pi2 = C_PI_2_DEFAULT, directionUpper = NA,
        maxNumberOfSubjects = NA_real_, groups = 2, allocationRatioPlanned = NA_real_) {
    objectType <- match.arg(objectType)

    .assertIsTrialDesignInverseNormalOrGroupSequential(design)
    .assertIsValidAlphaAndBeta(design$alpha, design$beta)
    .assertIsValidSidedParameter(design$sided)
    .assertIsValidGroupsParameter(groups)
    .assertIsSingleLogical(normalApproximation, "normalApproximation")
    .assertIsSingleLogical(riskRatio, "riskRatio")
    .assertIsValidDirectionUpper(directionUpper, design$sided, objectType, userFunctionCallEnabled = TRUE)

    if (groups == 1) {
        if (!any(is.na(pi1)) && any(pi1 == thetaH0) && (objectType == "sampleSize")) {
            stop(
                C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                "any 'pi1' (", .arrayToString(pi1), ") must be != 'thetaH0' (", thetaH0, ")"
            )
        }

        if (any(is.na(pi1)) || any(pi1 <= 0) || any(pi1 >= 1)) {
            stop(
                C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
                "probability 'pi1' (", .arrayToString(pi1), ") is out of bounds (0; 1)"
            )
        }

        if (thetaH0 >= 1 || thetaH0 <= 0) {
            stop(C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS, "'thetaH0' (", thetaH0, ") is out of bounds (0; 1)")
        }

        if (!normalApproximation && design$sided == 2 && (objectType == "sampleSize")) {
            stop(
                C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                "exact sample size calculation not available for two-sided testing"
            )
        }
    } else if (groups == 2) {
        if (!any(is.na(c(pi1, pi2))) && any(abs(pi1 - pi2 - thetaH0) < 1E-12) &&
                (objectType == "sampleSize") && !riskRatio) {
            stop(
                C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                "any 'pi1 - pi2' (", .arrayToString(pi1 - pi2), ") must be != 'thetaH0' (", thetaH0, ")"
            )
        }

        if (!any(is.na(c(pi1, pi2))) && any(abs(pi1 / pi2 - thetaH0) < 1E-12) &&
                (objectType == "sampleSize") && riskRatio) {
            stop(
                C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                "any 'pi1 / pi2' (", .arrayToString(pi1 / pi2), ") must be != 'thetaH0' (", thetaH0, ")"
            )
        }

        if (any(is.na(pi1)) || any(pi1 <= 0) || any(pi1 >= 1)) {
            stop(
                C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
                "probability 'pi1' (", .arrayToString(pi1), ") is out of bounds (0; 1)"
            )
        }

        if (any(is.na(pi2)) || any(pi2 <= 0) || any(pi2 >= 1)) {
            stop(
                C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
                "probability 'pi2' (", .arrayToString(pi2), ") is out of bounds (0; 1)"
            )
        }

        if (design$sided == 2 && ((thetaH0 != 0 && !riskRatio) || (thetaH0 != 1 && riskRatio))) {
            stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "two-sided case is implemented only for superiority testing")
        }

        if (!normalApproximation) {
            stop(
                C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                "only normal approximation case is implemented for two groups"
            )
        }

        if (is.na(allocationRatioPlanned)) {
            allocationRatioPlanned <- C_ALLOCATION_RATIO_DEFAULT
        }

        if (allocationRatioPlanned < 0) {
            stop(
                C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                "'allocationRatioPlanned' (", allocationRatioPlanned, ") must be >= 0"
            )
        }

        if (riskRatio && thetaH0 <= 0) {
            stop(
                C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                "null hypothesis risk ratio is not allowed be negative or zero, ",
                "i.e., 'thetaH0' must be > 0 if 'riskRatio' = TRUE"
            )
        }
    }

    designPlan <- TrialDesignPlanRates(design = design)
    designPlan$.setSampleSizeObject(objectType)

    designPlan$criticalValuesPValueScale <- matrix(design$stageLevels, ncol = 1)
    if (design$sided == 2) {
        designPlan$criticalValuesPValueScale <- designPlan$criticalValuesPValueScale * 2
        designPlan$.setParameterType("criticalValuesPValueScale", C_PARAM_GENERATED)
    }

    if (any(design$futilityBounds > C_FUTILITY_BOUNDS_DEFAULT)) {
        designPlan$futilityBoundsPValueScale <- matrix(1 - stats::pnorm(design$futilityBounds), ncol = 1)
        designPlan$.setParameterType("futilityBoundsPValueScale", C_PARAM_GENERATED)
    }

    if (objectType == "power") {
        .assertIsValidMaxNumberOfSubjects(maxNumberOfSubjects)
        .setValueAndParameterType(designPlan, "maxNumberOfSubjects", maxNumberOfSubjects, NA_real_)
        .setValueAndParameterType(designPlan, "directionUpper", directionUpper, TRUE)

        designPlan$.setParameterType("effect", C_PARAM_GENERATED)
    }

    .setValueAndParameterType(designPlan, "normalApproximation", normalApproximation, TRUE)
    .setValueAndParameterType(designPlan, "thetaH0", thetaH0, ifelse(riskRatio, 1, 0))
    .assertIsValidThetaH0(thetaH0, endpoint = "rates", groups = groups, ratioEnabled = riskRatio)
    if (objectType == "power") {
        .setValueAndParameterType(designPlan, "pi1", pi1, C_PI_1_DEFAULT)
    } else {
        .setValueAndParameterType(designPlan, "pi1", pi1, C_PI_1_SAMPLE_SIZE_DEFAULT)
    }
    .setValueAndParameterType(designPlan, "pi2", pi2, 0.2)
    if (groups == 1) {
        if (designPlan$.getParameterType("pi2") == C_PARAM_USER_DEFINED) {
            warning("'pi2' (", pi2, ") will be ignored ",
                "because it is not applicable for 'groups' = 1",
                call. = FALSE
            )
        }
        designPlan$.setParameterType("pi2", C_PARAM_NOT_APPLICABLE)

        if (isTRUE(riskRatio)) {
            warning("'riskRatio' (", riskRatio, ") will be ignored ",
                "because it is not applicable for 'groups' = 1",
                call. = FALSE
            )
        }
        designPlan$.setParameterType("riskRatio", C_PARAM_NOT_APPLICABLE)

        if (length(allocationRatioPlanned) == 1 && !is.na(allocationRatioPlanned)) {
            warning("'allocationRatioPlanned' (", allocationRatioPlanned,
                ") will be ignored because it is not applicable for 'groups' = 1",
                call. = FALSE
            )
        }
        designPlan$.setParameterType("allocationRatioPlanned", C_PARAM_NOT_APPLICABLE)
    } else {
        .setValueAndParameterType(designPlan, "riskRatio", riskRatio, FALSE)
        .setValueAndParameterType(
            designPlan, "allocationRatioPlanned",
            allocationRatioPlanned, C_ALLOCATION_RATIO_DEFAULT
        )
    }
    .setValueAndParameterType(designPlan, "groups", groups, 2)

    return(designPlan)
}

#' @title
#' Get Power Means
#'
#' @description
#' Returns the power, stopping probabilities, and expected sample size for
#' testing means in one or two samples at given sample size.
#'
#' @inheritParams param_design_with_default
#' @inheritParams param_groups
#' @param normalApproximation The type of computation of the p-values. If \code{TRUE}, the variance is
#'        assumed to be known, default is \code{FALSE}, i.e., the calculations are performed
#'        with the t distribution.
#' @param meanRatio If \code{TRUE}, the sample size for
#'        one-sided testing of H0: \code{mu1 / mu2 = thetaH0} is calculated, default is \code{FALSE}.
#' @inheritParams param_thetaH0
#' @inheritParams param_alternative
#' @inheritParams param_stDev
#' @inheritParams param_allocationRatioPlanned
#' @inheritParams param_directionUpper
#' @inheritParams param_maxNumberOfSubjects
#' @inheritParams param_three_dots
#'
#' @details
#' At given design the function calculates the power, stopping probabilities,
#' and expected sample size, for testing means at given sample size.
#' In a two treatment groups design, additionally, an
#' allocation ratio = \code{n1 / n2} can be specified.
#' A null hypothesis value thetaH0 != 0 for testing the difference of two means
#' or \code{thetaH0 != 1} for testing the ratio of two means can be specified.
#' For the specified sample size, critical bounds and stopping for futility
#' bounds are provided at the effect scale (mean, mean difference, or
#' mean ratio, respectively)
#'
#' @template return_object_trial_design_plan
#' @template how_to_get_help_for_generics
#'
#' @family power functions
#'
#' @template examples_get_power_means
#'
#' @export
#'
getPowerMeans <- function(design = NULL, ...,
        groups = 2L,
        normalApproximation = FALSE,
        meanRatio = FALSE,
        thetaH0 = ifelse(meanRatio, 1, 0),
        alternative = seq(0, 1, 0.2), # C_ALTERNATIVE_POWER_SIMULATION_DEFAULT
        stDev = 1, # C_STDEV_DEFAULT
        directionUpper = NA,
        maxNumberOfSubjects = NA_real_,
        allocationRatioPlanned = NA_real_ # C_ALLOCATION_RATIO_DEFAULT
        ) {
    .assertIsValidMaxNumberOfSubjects(maxNumberOfSubjects)

    if (is.null(design)) {
        design <- .getDefaultDesign(..., type = "power")
        .warnInCaseOfUnknownArguments(
            functionName = "getPowerMeans",
            ignore = .getDesignArgumentsToIgnoreAtUnknownArgumentCheck(design, powerCalculationEnabled = TRUE), ...
        )
    } else {
        .warnInCaseOfUnknownArguments(functionName = "getPowerMeans", ...)
        .assertIsTrialDesign(design)
        .warnInCaseOfTwoSidedPowerArgument(...)
        .warnInCaseOfTwoSidedPowerIsDisabled(design)
    }

    designPlan <- .createDesignPlanMeans(
        objectType = "power",
        design = design, normalApproximation = normalApproximation, meanRatio = meanRatio,
        thetaH0 = thetaH0, alternative = alternative, stDev = stDev, directionUpper = directionUpper,
        maxNumberOfSubjects = maxNumberOfSubjects, groups = groups,
        allocationRatioPlanned = allocationRatioPlanned, ...
    )

    if (designPlan$groups == 1) {
        theta <- (designPlan$alternative - designPlan$thetaH0) / designPlan$stDev
        if (!is.na(designPlan$directionUpper) && !designPlan$directionUpper) {
            theta <- -theta
        }
        if (designPlan$normalApproximation) {
            powerAndAverageSampleNumber <- getPowerAndAverageSampleNumber(
                design, theta, maxNumberOfSubjects
            )
        } else {
            thetaAdj <- (sign(theta) * .getOneMinusQNorm(design$alpha / design$sided) -
                .getQNorm(stats::pt(
                    sign(theta) * stats::qt(1 - design$alpha / design$sided, maxNumberOfSubjects - 1),
                    maxNumberOfSubjects - 1,
                    theta * sqrt(maxNumberOfSubjects)
                ))) / sqrt(maxNumberOfSubjects)
            powerAndAverageSampleNumber <- getPowerAndAverageSampleNumber(
                design, thetaAdj, maxNumberOfSubjects
            )
        }
    } else {
        if (!designPlan$meanRatio) {
            theta <- sqrt(designPlan$allocationRatioPlanned) / (1 + designPlan$allocationRatioPlanned) *
                (designPlan$alternative - designPlan$thetaH0) / designPlan$stDev
        } else {
            theta <- sqrt(designPlan$allocationRatioPlanned) /
                sqrt((1 + designPlan$allocationRatioPlanned * designPlan$thetaH0^2) *
                    (1 + designPlan$allocationRatioPlanned)) *
                (designPlan$alternative - designPlan$thetaH0) / designPlan$stDev
        }
        if (!is.na(designPlan$directionUpper) && !designPlan$directionUpper) {
            theta <- -theta
        }
        if (designPlan$normalApproximation) {
            powerAndAverageSampleNumber <- getPowerAndAverageSampleNumber(
                design, theta, maxNumberOfSubjects
            )
        } else {
            thetaAdj <- (sign(theta) * .getOneMinusQNorm(design$alpha / design$sided) -
                .getQNorm(stats::pt(
                    sign(theta) * stats::qt(1 - design$alpha / design$sided, maxNumberOfSubjects - 2),
                    maxNumberOfSubjects - 2,
                    theta * sqrt(maxNumberOfSubjects)
                ))) / sqrt(maxNumberOfSubjects)
            powerAndAverageSampleNumber <- getPowerAndAverageSampleNumber(
                design, thetaAdj, maxNumberOfSubjects
            )
        }
    }

    designPlan$effect <- designPlan$alternative - designPlan$thetaH0

    designPlan$expectedNumberOfSubjects <- powerAndAverageSampleNumber$averageSampleNumber
    designPlan$overallReject <- powerAndAverageSampleNumber$overallReject
    designPlan$rejectPerStage <- powerAndAverageSampleNumber$rejectPerStage
    designPlan$futilityStop <- powerAndAverageSampleNumber$overallFutility
    designPlan$futilityPerStage <- powerAndAverageSampleNumber$futilityPerStage
    designPlan$earlyStop <- powerAndAverageSampleNumber$overallEarlyStop

    parameterNames <- c("overallReject")
    if (design$kMax > 1) {
        parameterNames <- c(
            parameterNames,
            "expectedNumberOfSubjects",
            "rejectPerStage",
            "futilityStop",
            "futilityPerStage",
            "earlyStop"
        )
    }
    for (parameterName in parameterNames) {
        designPlan$.setParameterType(parameterName, C_PARAM_GENERATED)
    }

    .addNumberOfSubjectsToPowerResult(designPlan)
    .addEffectScaleBoundaryDataToDesignPlan(designPlan)
    .hideFutilityStopsIfNotApplicable(designPlan)

    return(designPlan)
}

#' @title
#' Get Power Rates
#'
#' @description
#' Returns the power, stopping probabilities, and expected sample size for testing rates
#' in one or two samples at given sample sizes.
#'
#' @inheritParams param_design_with_default
#' @inheritParams param_groups
#' @param riskRatio If \code{TRUE}, the power for one-sided
#'        testing of H0: \code{pi1 / pi2 = thetaH0} is calculated, default is \code{FALSE}.
#' @inheritParams param_thetaH0
#' @inheritParams param_pi1_rates
#' @inheritParams param_pi2_rates
#' @inheritParams param_directionUpper
#' @inheritParams param_maxNumberOfSubjects
#' @inheritParams param_allocationRatioPlanned
#' @inheritParams param_three_dots
#'
#' @details
#' At given design the function calculates the power, stopping probabilities, and expected sample size,
#' for testing rates for given maximum sample size.
#' The sample sizes over the stages are calculated according to the specified information rate in the design.
#' In a two treatment groups design, additionally, an allocation ratio = n1/n2 can be specified.
#' If a null hypothesis value thetaH0 != 0 for testing the difference of two rates
#' or \code{thetaH0 != 1} for testing the risk ratio is specified, the
#' formulas according to Farrington & Manning (Statistics in Medicine, 1990) are used (only one-sided testing).
#' Critical bounds and stopping for futility bounds are provided at the effect scale
#' (rate, rate difference, or rate ratio, respectively).
#' For the two-sample case, the calculation here is performed at fixed pi2 as given as argument in the function.
#' Note that the power calculation for rates is always based on the normal approximation.
#'
#' @template return_object_trial_design_plan
#' @template how_to_get_help_for_generics
#'
#' @family power functions
#'
#' @template examples_get_power_rates
#'
#' @export
#'
getPowerRates <- function(design = NULL, ...,
        groups = 2L,
        riskRatio = FALSE,
        thetaH0 = ifelse(riskRatio, 1, 0),
        pi1 = seq(0.2, 0.5, 0.1), # C_PI_1_DEFAULT
        pi2 = 0.2, # C_PI_2_DEFAULT
        directionUpper = NA,
        maxNumberOfSubjects = NA_real_,
        allocationRatioPlanned = NA_real_ # C_ALLOCATION_RATIO_DEFAULT
        ) {
    if (is.null(design)) {
        design <- .getDefaultDesign(..., type = "power")
        .warnInCaseOfUnknownArguments(
            functionName = "getPowerRates",
            ignore = .getDesignArgumentsToIgnoreAtUnknownArgumentCheck(design, powerCalculationEnabled = TRUE), ...
        )
    } else {
        .warnInCaseOfUnknownArguments(functionName = "getPowerRates", ...)
        .assertIsTrialDesign(design)
        .warnInCaseOfTwoSidedPowerArgument(...)
        .warnInCaseOfTwoSidedPowerIsDisabled(design)
    }

    designPlan <- .createDesignPlanRates(
        objectType = "power",
        design = design, riskRatio = riskRatio,
        thetaH0 = thetaH0, pi1 = pi1, pi2 = pi2, directionUpper = directionUpper,
        maxNumberOfSubjects = maxNumberOfSubjects, groups = groups,
        allocationRatioPlanned = allocationRatioPlanned, ...
    )

    if (!is.na(allocationRatioPlanned) && allocationRatioPlanned <= 0) {
        stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "allocation ratio must be > 0")
    }

    allocationRatioPlanned <- designPlan$allocationRatioPlanned

    theta <- rep(NA_real_, length(pi1))
    if (groups == 1) {
        designPlan$effect <- pi1 - thetaH0
        theta <- (pi1 - thetaH0) / sqrt(pi1 * (1 - pi1)) + sign(pi1 - thetaH0) *
            .getOneMinusQNorm(design$alpha / design$sided) *
            (1 - sqrt(thetaH0 * (1 - thetaH0) / (pi1 * (1 - pi1)))) / sqrt(maxNumberOfSubjects)
    } else {
        if (!riskRatio) {
            designPlan$effect <- pi1 - pi2 - thetaH0
            for (i in (1:length(pi1))) {
                fm <- .getFarringtonManningValues(
                    rate1 = pi1[i], rate2 = pi2,
                    theta = thetaH0, allocation = allocationRatioPlanned, method = "diff"
                )
                theta[i] <- sqrt(allocationRatioPlanned) / (1 + allocationRatioPlanned) *
                    (pi1[i] - pi2 - thetaH0) * sqrt(1 + allocationRatioPlanned) /
                    sqrt(pi1[i] * (1 - pi1[i]) + allocationRatioPlanned * pi2 * (1 - pi2)) +
                    sign(pi1[i] - pi2 - thetaH0) * .getOneMinusQNorm(design$alpha / design$sided) *
                        (1 - sqrt(fm$ml1 * (1 - fm$ml1) + allocationRatioPlanned * fm$ml2 * (1 - fm$ml2)) /
                            sqrt(pi1[i] * (1 - pi1[i]) + allocationRatioPlanned * pi2 * (1 - pi2))) /
                        sqrt(maxNumberOfSubjects)
            }
        } else {
            designPlan$effect <- pi1 / pi2 - thetaH0
            for (i in (1:length(pi1))) {
                fm <- .getFarringtonManningValues(
                    rate1 = pi1[i], rate2 = pi2,
                    theta = thetaH0, allocation = allocationRatioPlanned, method = "ratio"
                )
                theta[i] <- sqrt(allocationRatioPlanned) / (1 + allocationRatioPlanned) *
                    (pi1[i] - thetaH0 * pi2) * sqrt(1 + allocationRatioPlanned) /
                    sqrt(pi1[i] * (1 - pi1[i]) + allocationRatioPlanned * thetaH0^2 * pi2 * (1 - pi2)) +
                    sign(pi1[i] - thetaH0 * pi2) * .getOneMinusQNorm(design$alpha / design$sided) *
                        (1 - sqrt(fm$ml1 * (1 - fm$ml1) + allocationRatioPlanned * thetaH0^2 *
                            fm$ml2 * (1 - fm$ml2)) / sqrt(pi1[i] * (1 - pi1[i]) + allocationRatioPlanned *
                            thetaH0^2 * pi2 * (1 - pi2))) /
                        sqrt(maxNumberOfSubjects)
            }
        }
    }

    if (!is.na(designPlan$directionUpper) && !designPlan$directionUpper) {
        theta <- -theta
    }

    powerAndAverageSampleNumber <- getPowerAndAverageSampleNumber(design, theta, maxNumberOfSubjects)

    designPlan$expectedNumberOfSubjects <- powerAndAverageSampleNumber$averageSampleNumber
    designPlan$overallReject <- powerAndAverageSampleNumber$overallReject
    designPlan$rejectPerStage <- powerAndAverageSampleNumber$rejectPerStage
    designPlan$futilityStop <- powerAndAverageSampleNumber$overallFutility
    designPlan$futilityPerStage <- powerAndAverageSampleNumber$futilityPerStage
    designPlan$earlyStop <- powerAndAverageSampleNumber$overallEarlyStop

    parameterNames <- c("overallReject")
    if (design$kMax > 1) {
        parameterNames <- c(
            parameterNames,
            "expectedNumberOfSubjects",
            "rejectPerStage",
            "futilityStop",
            "futilityPerStage",
            "earlyStop"
        )
    }
    for (parameterName in parameterNames) {
        designPlan$.setParameterType(parameterName, C_PARAM_GENERATED)
    }

    .addNumberOfSubjectsToPowerResult(designPlan)
    .addEffectScaleBoundaryDataToDesignPlan(designPlan)
    .hideFutilityStopsIfNotApplicable(designPlan)

    return(designPlan)
}

.getNumberOfSubjectsInner <- function(..., timeValue, accrualTime, accrualIntensity, maxNumberOfSubjects) {
    .assertIsSingleNumber(timeValue, "timeValue")
    if (length(accrualTime) != length(accrualIntensity)) {
        stop(
            C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "length of 'accrualTime' (", length(accrualIntensity), ") ",
            "must be equel to length of 'accrualIntensity' (", length(accrualIntensity), ")"
        )
    }

    densityIntervals <- accrualTime
    if (length(accrualTime) > 1) {
        densityIntervals[2:length(accrualTime)] <- accrualTime[2:length(accrualTime)] -
            accrualTime[1:(length(accrualTime) - 1)]
    }
    densityVector <- accrualIntensity / sum(densityIntervals * accrualIntensity)
    for (l in 1:length(densityVector)) {
        if (timeValue <= accrualTime[l]) {
            if (l == 1) {
                return(timeValue * densityVector[l] * maxNumberOfSubjects)
            } else {
                return((sum(densityVector[1:(l - 1)] * densityIntervals[1:(l - 1)]) +
                    (timeValue - accrualTime[l - 1]) * densityVector[l]) * maxNumberOfSubjects)
            }
        }
    }
    return(maxNumberOfSubjects)
}

.getNumberOfSubjects <- function(..., time, accrualTime, accrualIntensity, maxNumberOfSubjects) {
    subjectNumbers <- c()
    for (timeValue in time) {
        if (is.na(timeValue)) {
            return(NA_real_)
        }

        subjectNumbers <- c(
            subjectNumbers,
            .getNumberOfSubjectsInner(
                timeValue = timeValue, accrualTime = accrualTime,
                accrualIntensity = accrualIntensity, maxNumberOfSubjects = maxNumberOfSubjects
            )
        )
    }
    return(subjectNumbers)
}

#' @title
#' Get Power Survival
#'
#' @description
#' Returns the power, stopping probabilities, and expected sample size for testing
#' the hazard ratio in a two treatment groups survival design.
#'
#' @inheritParams param_design_with_default
#' @inheritParams param_typeOfComputation
#' @inheritParams param_allocationRatioPlanned
#' @inheritParams param_thetaH0
#' @inheritParams param_lambda1
#' @inheritParams param_lambda2
#' @inheritParams param_pi1_survival
#' @inheritParams param_pi2_survival
#' @inheritParams param_median1
#' @inheritParams param_median2
#' @inheritParams param_piecewiseSurvivalTime
#' @inheritParams param_directionUpper
#' @inheritParams param_accrualTime
#' @inheritParams param_accrualIntensity
#' @inheritParams param_accrualIntensityType
#' @inheritParams param_eventTime
#' @inheritParams param_hazardRatio
#' @inheritParams param_kappa
#' @inheritParams param_dropoutRate1
#' @inheritParams param_dropoutRate2
#' @inheritParams param_dropoutTime
#' @param maxNumberOfEvents \code{maxNumberOfEvents > 0} is the maximum number of events, it determines
#'        the power of the test and needs to be specified.
#' @inheritParams param_maxNumberOfSubjects_survival
#' @inheritParams param_three_dots
#'
#' @details
#' At given design the function calculates the power, stopping probabilities, and expected
#' sample size at given number of events and number of subjects.
#' It also calculates the time when the required events are expected under the given
#' assumptions (exponentially, piecewise exponentially, or Weibull distributed survival times
#' and constant or non-constant piecewise accrual).
#' Additionally, an allocation ratio = n1/n2 can be specified where n1 and n2 are the number
#' of subjects in the two treatment groups.
#'
#' The formula of Kim & Tsiatis (Biometrics, 1990)
#' is used to calculate the expected number of events under the alternative
#' (see also Lakatos & Lan, Statistics in Medicine, 1992). These formulas are generalized to piecewise survival times and
#' non-constant piecewise accrual over time.\cr
#'
#' @template details_piecewise_survival
#'
#' @template details_piecewise_accrual
#'
#' @template return_object_trial_design_plan
#' @template how_to_get_help_for_generics
#'
#' @family power functions
#'
#' @template examples_get_power_survival
#'
#' @export
#'
getPowerSurvival <- function(design = NULL, ...,
        typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
        thetaH0 = 1, # C_THETA_H0_SURVIVAL_DEFAULT
        directionUpper = NA,
        pi1 = NA_real_,
        pi2 = NA_real_,
        lambda1 = NA_real_,
        lambda2 = NA_real_,
        median1 = NA_real_,
        median2 = NA_real_,
        kappa = 1,
        hazardRatio = NA_real_,
        piecewiseSurvivalTime = NA_real_,
        allocationRatioPlanned = 1, # C_ALLOCATION_RATIO_DEFAULT
        eventTime = 12, # C_EVENT_TIME_DEFAULT
        accrualTime = c(0, 12), # C_ACCRUAL_TIME_DEFAULT
        accrualIntensity = 0.1, # C_ACCRUAL_INTENSITY_DEFAULT
        accrualIntensityType = c("auto", "absolute", "relative"),
        maxNumberOfSubjects = NA_real_,
        maxNumberOfEvents = NA_real_,
        dropoutRate1 = 0, # C_DROP_OUT_RATE_1_DEFAULT
        dropoutRate2 = 0, # C_DROP_OUT_RATE_2_DEFAULT
        dropoutTime = 12 # C_DROP_OUT_TIME_DEFAULT
        ) {
    if (is.null(design)) {
        design <- .getDefaultDesign(..., type = "power")
        .warnInCaseOfUnknownArguments(
            functionName = "getPowerSurvival",
            ignore = .getDesignArgumentsToIgnoreAtUnknownArgumentCheck(design, powerCalculationEnabled = TRUE), ...
        )
    } else {
        .assertIsTrialDesign(design)
        .warnInCaseOfUnknownArguments(functionName = "getPowerSurvival", ...)
        .warnInCaseOfTwoSidedPowerArgument(...)
        .warnInCaseOfTwoSidedPowerIsDisabled(design)
    }

    designPlan <- .createDesignPlanSurvival(
        objectType = "power",
        design = design,
        typeOfComputation = typeOfComputation,
        thetaH0 = thetaH0,
        pi2 = pi2,
        pi1 = pi1,
        allocationRatioPlanned = allocationRatioPlanned,
        accountForObservationTimes = TRUE,
        eventTime = eventTime,
        accrualTime = accrualTime,
        accrualIntensity = accrualIntensity,
        accrualIntensityType = accrualIntensityType,
        kappa = kappa,
        piecewiseSurvivalTime = piecewiseSurvivalTime,
        lambda2 = lambda2,
        lambda1 = lambda1,
        median1 = median1,
        median2 = median2,
        directionUpper = directionUpper,
        maxNumberOfEvents = maxNumberOfEvents,
        maxNumberOfSubjects = maxNumberOfSubjects,
        dropoutRate1 = dropoutRate1,
        dropoutRate2 = dropoutRate2,
        dropoutTime = dropoutTime,
        hazardRatio = hazardRatio
    )

    .assertIsSingleNumber(allocationRatioPlanned, "allocationRatioPlanned")
    .assertIsInOpenInterval(allocationRatioPlanned, "allocationRatioPlanned", 0, C_ALLOCATION_RATIO_MAXIMUM)

    if (designPlan$typeOfComputation == "Schoenfeld") {
        theta <- sqrt(allocationRatioPlanned) / (1 + allocationRatioPlanned) *
            (log(designPlan$hazardRatio / thetaH0))
    } else if (designPlan$typeOfComputation == "Freedman") {
        theta <- sqrt(allocationRatioPlanned) * (designPlan$hazardRatio - 1) /
            (allocationRatioPlanned * designPlan$hazardRatio + 1)
    } else if (designPlan$typeOfComputation == "HsiehFreedman") {
        theta <- sqrt(4 * allocationRatioPlanned) / (1 + allocationRatioPlanned) *
            (designPlan$hazardRatio - 1) / (designPlan$hazardRatio + 1)
    }

    if (!is.na(designPlan$directionUpper) && !designPlan$directionUpper) {
        theta <- -theta
    }

    powerAndAverageSampleNumber <- getPowerAndAverageSampleNumber(
        design = design, theta = theta, nMax = maxNumberOfEvents
    )

    kMax <- design$kMax
    sided <- design$sided

    if (designPlan$.piecewiseSurvivalTime$.isLambdaBased()) {
        numberOfResults <- length(designPlan$hazardRatio)
    } else {
        numberOfResults <- length(designPlan$pi1)
    }

    stoppingProbs <- matrix(NA_real_, kMax, numberOfResults)
    designPlan$analysisTime <- matrix(NA_real_, kMax, numberOfResults)
    designPlan$numberOfSubjects <- matrix(NA_real_, kMax, numberOfResults)
    designPlan$studyDuration <- rep(NA_real_, numberOfResults)
    designPlan$expectedNumberOfSubjects <- rep(NA_real_, numberOfResults)
    eventsPerStage <- maxNumberOfEvents * design$informationRates
    parameterNames <- c(
        "analysisTime",
        "numberOfSubjects",
        "studyDuration",
        "expectedNumberOfSubjects",
        "eventsPerStage"
    )

    phi <- -c(log(1 - dropoutRate1), log(1 - dropoutRate2)) / dropoutTime

    lambda1 <- designPlan$lambda1
    if (designPlan$.piecewiseSurvivalTime$piecewiseSurvivalEnabled) {
        lambda1 <- rep(NA_real_, numberOfResults)
    }

    for (i in 1:numberOfResults) {
        # Analysis times
        up <- 2
        iterate <- 1
        while (eventsPerStage[kMax] / designPlan$maxNumberOfSubjects > .getEventProbabilities(
            time = up, accrualTimeVector = designPlan$accrualTime,
            accrualIntensity = designPlan$accrualIntensity,
            lambda2 = designPlan$lambda2,
            lambda1 = lambda1[i], phi = phi,
            piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
            kappa = kappa, allocationRatioPlanned = allocationRatioPlanned,
            hazardRatio = designPlan$hazardRatio[i]
        )) {
            up <- 2 * up
            iterate <- iterate + 1
            if (iterate > 50) {
                stop(
                    C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
                    "'maxNumberOfSubjects' (", designPlan$maxNumberOfSubjects, ") ",
                    "is too small to reach maximum number of events ",
                    "(presumably due to drop-out rates)"
                )
            }
        }
        for (j in 1:kMax) {
            designPlan$analysisTime[j, i] <- .getOneDimensionalRoot(
                function(x) {
                    eventsPerStage[j] / designPlan$maxNumberOfSubjects -
                        .getEventProbabilities(
                            time = x,
                            accrualTimeVector = designPlan$accrualTime,
                            accrualIntensity = designPlan$accrualIntensity,
                            lambda2 = designPlan$lambda2,
                            lambda1 = lambda1[i], phi = phi,
                            piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
                            kappa = kappa, allocationRatioPlanned = allocationRatioPlanned,
                            hazardRatio = designPlan$hazardRatio[i]
                        )
                },
                lower = 0, upper = up, tolerance = 1e-06,
                callingFunctionInformation = "getPowerSurvival"
            )
            if (is.na(designPlan$analysisTime[j, i])) {
                warning("Cannot calculate analysis time at stage ", j, ": ",
                    "'maxNumberOfSubjects' (", designPlan$maxNumberOfSubjects, ") is too ",
                    "small to reach maximum number of events",
                    call. = FALSE
                )
            }
        }
        if (kMax > 1) {
            designPlan$numberOfSubjects[, i] <- .getNumberOfSubjects(
                time = designPlan$analysisTime[, i],
                accrualTime = designPlan$accrualTime,
                accrualIntensity = designPlan$accrualIntensity,
                maxNumberOfSubjects = designPlan$maxNumberOfSubjects
            )

            powerAndAverageSampleNumber$futilityPerStage[is.na(
                powerAndAverageSampleNumber$futilityPerStage[, i]
            ), i] <- 0

            stoppingProbs[, i] <- powerAndAverageSampleNumber$rejectPerStage[, i] +
                c(powerAndAverageSampleNumber$futilityPerStage[, i], 0)

            stoppingProbs[kMax, i] <- 1 - sum(stoppingProbs[1:(kMax - 1), i])
            designPlan$studyDuration[i] <- designPlan$analysisTime[, i] %*% stoppingProbs[, i]
            designPlan$.setParameterType("studyDuration", C_PARAM_GENERATED)
            designPlan$expectedNumberOfSubjects[i] <-
                designPlan$numberOfSubjects[, i] %*% stoppingProbs[, i]
        }
    }

    if (kMax == 1) {
        designPlan$expectedNumberOfSubjects <- .getNumberOfSubjects(
            time = designPlan$analysisTime[1, ],
            accrualTime = designPlan$accrualTime,
            accrualIntensity = designPlan$accrualIntensity,
            maxNumberOfSubjects = designPlan$maxNumberOfSubjects
        )

        designPlan$numberOfSubjects <- matrix(designPlan$expectedNumberOfSubjects, nrow = 1)
    }

    designPlan$eventsPerStage <- matrix(eventsPerStage, ncol = 1)
    designPlan$.setParameterType("eventsPerStage", C_PARAM_GENERATED)

    designPlan$expectedNumberOfEvents <- powerAndAverageSampleNumber$averageSampleNumber
    designPlan$overallReject <- powerAndAverageSampleNumber$overallReject
    designPlan$rejectPerStage <- powerAndAverageSampleNumber$rejectPerStage
    designPlan$futilityStop <- powerAndAverageSampleNumber$overallFutility
    designPlan$futilityPerStage <- powerAndAverageSampleNumber$futilityPerStage
    designPlan$earlyStop <- powerAndAverageSampleNumber$overallEarlyStop

    parameterNames <- c(
        parameterNames,
        "expectedNumberOfEvents",
        "overallReject",
        "rejectPerStage",
        "futilityStop",
        "futilityPerStage",
        "earlyStop"
    )

    for (parameterName in parameterNames) {
        designPlan$.setParameterType(parameterName, C_PARAM_GENERATED)
    }

    if (kMax == 1L) {
        designPlan$.setParameterType("numberOfSubjects", C_PARAM_NOT_APPLICABLE)
        designPlan$.setParameterType("eventsPerStage", C_PARAM_NOT_APPLICABLE)

        designPlan$.setParameterType("futilityStop", C_PARAM_NOT_APPLICABLE)
        designPlan$.setParameterType("earlyStop", C_PARAM_NOT_APPLICABLE)
        designPlan$.setParameterType("rejectPerStage", C_PARAM_NOT_APPLICABLE)
    }

    if (!any(is.na(designPlan$analysisTime)) && !any(is.na(designPlan$accrualTime))) {
        designPlan$followUpTime <- designPlan$analysisTime[kMax, ] -
            designPlan$accrualTime[length(designPlan$accrualTime)]
        designPlan$.setParameterType("followUpTime", C_PARAM_GENERATED)
    }

    .addEffectScaleBoundaryDataToDesignPlan(designPlan)
    .addStudyDurationToDesignPlan(designPlan)
    .hideFutilityStopsIfNotApplicable(designPlan)

    return(designPlan)
}

.hideFutilityStopsIfNotApplicable <- function(designPlan) {
    if (all(designPlan$.design$futilityBounds == C_FUTILITY_BOUNDS_DEFAULT)) {
        designPlan$.setParameterType("futilityStop", C_PARAM_NOT_APPLICABLE)
        designPlan$.setParameterType("futilityPerStage", C_PARAM_NOT_APPLICABLE)
    }
}

.addStudyDurationToDesignPlan <- function(designPlan) {
    if (!designPlan$accountForObservationTimes) {
        return(invisible())
    }

    kMax <- designPlan$.design$kMax
    if (kMax == 1) {
        designPlan$studyDuration <- designPlan$analysisTime[1, ]
        designPlan$.setParameterType("studyDuration", C_PARAM_GENERATED)
        designPlan$maxStudyDuration <- designPlan$studyDuration
    } else {
        designPlan$maxStudyDuration <- designPlan$analysisTime[kMax, ]
        designPlan$.setParameterType("maxStudyDuration", C_PARAM_GENERATED)
    }
}

.addNumberOfSubjectsToPowerResult <- function(designPlan) {
    design <- designPlan$.design

    designPlan$numberOfSubjects <- matrix(rep(NA_real_, design$kMax), ncol = 1)
    designPlan$numberOfSubjects[1, 1] <- design$informationRates[1] * designPlan$maxNumberOfSubjects
    if (design$kMax > 1) {
        designPlan$numberOfSubjects[2:design$kMax, 1] <- (design$informationRates[2:design$kMax] -
            design$informationRates[1:(design$kMax - 1)]) * designPlan$maxNumberOfSubjects
    }

    designPlan$numberOfSubjects <- .getColumnCumSum(designPlan$numberOfSubjects)

    designPlan$numberOfSubjects1 <- .getNumberOfSubjects1(
        designPlan$numberOfSubjects, designPlan$allocationRatioPlanned
    )
    designPlan$numberOfSubjects2 <- .getNumberOfSubjects2(
        designPlan$numberOfSubjects, designPlan$allocationRatioPlanned
    )

    if (designPlan$.design$kMax == 1) {
        designPlan$nFixed <- as.numeric(designPlan$numberOfSubjects)
        designPlan$.setParameterType("nFixed", C_PARAM_GENERATED)
        if (designPlan$groups == 2) {
            designPlan$nFixed1 <- as.numeric(designPlan$numberOfSubjects1)
            designPlan$nFixed2 <- as.numeric(designPlan$numberOfSubjects2)
            designPlan$.setParameterType("nFixed1", C_PARAM_GENERATED)
            designPlan$.setParameterType("nFixed2", C_PARAM_GENERATED)
        }
    } else {
        designPlan$.setParameterType("numberOfSubjects", C_PARAM_GENERATED)

        if ((designPlan$groups == 1) ||
                designPlan$allocationRatioPlanned == 1) {
            designPlan$.setParameterType("numberOfSubjects1", C_PARAM_NOT_APPLICABLE)
            designPlan$.setParameterType("numberOfSubjects2", C_PARAM_NOT_APPLICABLE)
        } else {
            designPlan$.setParameterType("numberOfSubjects1", C_PARAM_GENERATED)
            designPlan$.setParameterType("numberOfSubjects2", C_PARAM_GENERATED)
        }
    }
}

Try the rpact package in your browser

Any scripts or data that you put into this service are public.

rpact documentation built on July 9, 2023, 6:30 p.m.