.packageName <- "ump"
arpv.binom <- function(x, n, p, plot = TRUE, ...) {

    if (! is.numeric(x)) stop("x not numeric")
    if (! is.numeric(n)) stop("n not numeric")
    if (! is.numeric(p)) stop("p not numeric")

    if (length(x) != 1) stop("x not scalar")
    if (length(n) != 1) stop("n not scalar")
    if (length(p) != 1) stop("p not scalar")

    if (as.integer(x) != x) stop("x not integer")
    if (as.integer(n) != n) stop("n not integer")

    if (! (n > 0)) stop("n not positive")
    if (! (0 <= x & x <= n)) stop("x not in 0, ..., n")
    if (! (0 < p & p < 1)) stop("p not in (0, 1)")

    mu <- n * p

    foo <- sign(x - mu)

    if (foo == 0) {

        p <- dbinom(x, n, p)
        alpha <- c(1 - p, 1)
        phi <- c(0, 1)

    } else {

    if (foo > 0) {
        c1 <- seq(0, floor(mu))
        c2 <- x
    } else {
        c1 <- x
        c2 <- seq(n, ceiling(mu))
    }

    p1 <- dbinom(c1, n, p)
    p2 <- dbinom(c2, n, p)
    P1 <- pbinom(c1 - 1, n, p)
    P2 <- pbinom(c2, n, p, lower.tail = FALSE)
    M1 <- mu * pbinom(c1 - 2, n - 1, p)
    M2 <- mu * pbinom(c2 - 1, n - 1, p, lower.tail = FALSE)

    alpha.max.1 <- (p1 * (c2 - c1) - M1 - M2 + c2 * P1 + c2 * P2) / (c2 - mu)
    alpha.max.2 <- (p2 * (c1 - c2) - M1 - M2 + c1 * P2 + c1 * P1) / (c1 - mu)
    alpha.min.1 <- (- M1 - M2 + c2 * P1 + c2 * P2) / (c2 - mu)
    alpha.min.2 <- (- M1 - M2 + c1 * P1 + c1 * P2) / (c1 - mu)
    alpha.max <- pmin(alpha.max.1, alpha.max.2)
    alpha.min <- pmax(alpha.min.1, alpha.min.2)
    alpha.max[c2 - c1 == 1] <- 1
    alpha.min[c2 - c1 == n] <- 0

    gamma.max.1 <- (alpha.max * (c2 - mu) + (M1 - c2 * P1) + (M2 - c2 * P2)) /
        (p1 * (c2 - c1))
    gamma.max.2 <- (alpha.max * (c1 - mu) + (M2 - c1 * P2) + (M1 - c1 * P1)) /
        (p2 * (c1 - c2))
    gamma.min.1 <- (alpha.min * (c2 - mu) + (M1 - c2 * P1) + (M2 - c2 * P2)) /
        (p1 * (c2 - c1))
    gamma.min.2 <- (alpha.min * (c1 - mu) + (M2 - c1 * P2) + (M1 - c1 * P1)) /
        (p2 * (c1 - c2))

    inies <- alpha.max > alpha.min

    alpha.min <- alpha.min[inies]
    alpha.max <- alpha.max[inies]

    alpha.min[alpha.min < 0] <- 0
    alpha.max[alpha.max > 1] <- 1

    if (foo > 0) {
        ##### c2 == x #####
        phi.min <- gamma.min.2[inies]
        phi.max <- gamma.max.2[inies]
    } else {
        ##### c1 == x #####
        phi.min <- gamma.min.1[inies]
        phi.max <- gamma.max.1[inies]
    }

    phi.min[phi.min < 0] <- 0
    phi.max[phi.max > 1] <- 1

    alpha <- c(alpha.min[1], alpha.max)
    phi <- c(phi.min[1], phi.max)

    }

    if (plot)
        arpv.plot(alpha, phi, ...)

    return(invisible(list(alpha = alpha, phi = phi)))

}
arpv.plot <- function(alpha, phi, df = TRUE, verticals = TRUE) {

    if (! is.numeric(alpha)) stop("alpha not numeric")
    if (! is.numeric(phi)) stop("phi not numeric")
    if (! is.logical(df)) stop("df not logical")

    if (length(alpha) != length(phi)) stop("alpha and phi not same length")

    if (! all(0 <= alpha & alpha <= 1)) stop("alpha not in [0, 1]")
    if (! all(0 <= phi & phi <= 1)) stop("phi not in [0, 1]")

    if (df) {
        plot(alpha, phi, xlab = "significance level",
            ylab = "fuzzy P-value", type = "l")
        u <- par("usr")
        lines(c(u[1], alpha[1]), c(0, 0))
        lines(c(alpha[length(alpha)], u[2]), c(1, 1))
    } else {
        dens <- diff(phi) / diff(alpha)
        plot(range(alpha), range(0, dens), type = "n",
            xlab = "significance level",
            ylab = "density of randomized P-value")
        nalpha <- length(alpha)
        ndens <- length(dens)
        segments(alpha[-nalpha], dens, alpha[-1], dens)
        if (verticals) {
            # outer verticals
            jalpha <- c(1, nalpha)
            jdens <- c(1, ndens)
            segments(alpha[jalpha], rep(0, 2), alpha[jalpha], dens[jdens],
                lty = 2)
            # inner verticals
            if (nalpha > 2)
                segments(alpha[-jalpha], pmin(dens[-1], dens[-ndens]),
                    alpha[-jalpha], pmax(dens[-1], dens[-ndens]), lty = 2)
        }
    }
}
fci.binom <- function(x, n, alpha = 0.05, p = seq(0, 1, length = 10001),
    flat = 1 / 4) {

    if (! is.numeric(x)) stop("x not numeric")
    if (! is.numeric(n)) stop("n not numeric")
    if (! is.numeric(alpha)) stop("alpha not numeric")
    if (! is.numeric(p)) stop("p not numeric")

    if (length(x) != 1) stop("x not scalar")
    if (length(n) != 1) stop("n not scalar")
    if (length(alpha) != 1) stop("alpha not scalar")

    if (as.integer(x) != x) stop("x not integer")
    if (as.integer(n) != n) stop("n not integer")

    if (! (n > 0)) stop("n not positive")
    if (! (0 <= x & x <= n)) stop("x not in 0, ..., n")
    if (! (0 < alpha & alpha < 1)) stop("alpha not in (0, 1)")
    if (! all(0 <= p & p <= 1)) stop("p not in [0, 1]")

    phi <- umpu.binom(x, n, p, alpha)
    support <- range(p[phi < 1])

    cat(100 * (1 - alpha), "percent fuzzy confidence interval\n")

    if (x == 0 || x == n) {

        cat("core is empty\n")
        if (x == 0) {
            cat("support is [", support[1], ", ", support[2], ")\n", sep = "")
            xlim <- c(0, support[2] * (1 + flat))
        } else {
            cat("support is (", support[1], ", ", support[2], "]\n", sep = "")
            xlim <- c(support[1] - (1 - support[1]) * flat, 1)
        }

        plot(p, 1 - phi, xlim = xlim, ylim = c(0, 1),
            ylab = expression(1 - phi(x, alpha, p)), type = "l")

    } else {

        core <- range(p[phi == 0])
        cat("core is [", core[1], ", ", core[2], "]\n", sep = "")
        cat("support is (", support[1], ", ", support[2], ")\n", sep = "")

        lim.low <- c(support[1], core[1])
        lim.hig <- c(core[2], support[2])

        extra <- flat * max(diff(lim.low), diff(lim.hig))
        lim.low <- range(lim.low - extra, lim.low + extra)
        lim.hig <- range(lim.hig - extra, lim.hig + extra)
        lim.low <- pmax(0, lim.low)
        lim.hig <- pmax(0, lim.hig)
        lim.low <- pmin(1, lim.low)
        lim.hig <- pmin(1, lim.hig)

        if (lim.low[2] < lim.hig[1]) {
            oldpar <- par(mfrow = c(1, 2))
            plot(p, 1 - phi, xlim = lim.low, ylim = c(0, 1),
                ylab = expression(1 - phi(x, alpha, p)), type = "l")
            plot(p, 1 - phi, xlim = lim.hig, ylim = c(0, 1),
                ylab = expression(1 - phi(x, alpha, p)), type = "l")
            par(mfrow = oldpar)
        } else {
            plot(p, 1 - phi, xlim = c(lim.low[1], lim.hig[2]), ylim = c(0, 1),
                ylab = expression(1 - phi(x, alpha, p)), type = "l")
        }
    }
}

umpu.binom <- function(x, n, p, alpha, maxiter = 10, tol = 1e-9) {

    if (! is.numeric(x)) stop("x not numeric")
    if (! is.numeric(n)) stop("n not numeric")
    if (! is.numeric(p)) stop("p not numeric")
    if (! is.numeric(alpha)) stop("alpha not numeric")
    if (! is.numeric(maxiter)) stop("maxiter not numeric")
    if (! is.numeric(tol)) stop("tol not numeric")

    if (length(n) != 1) stop("n not scalar")
    if (length(maxiter) != 1) stop("maxiter not scalar")
    if (length(tol) != 1) stop("tol not scalar")
    foo <- (length(x) > 1) + (length(p) > 1) + (length(alpha) > 1)
    if (foo > 1) stop("at most one of x, p, alpha can be non-scalar")

    if (as.integer(n) != n) stop("n not integer")
    if (as.integer(maxiter) != maxiter) stop("maxiter not integer")
    if (any(as.integer(x) != x)) stop("x not integer")

    if (! (n > 0)) stop("n not positive")
    if (! (maxiter > 0)) stop("maxiter not positive")
    if (! (tol > 0)) stop("tol not positive")
    if (! all(0 <= x & x <= n)) stop("x not in 0, ..., n")

    if (! all(0 <= p & p <= 1)) stop("p not in [0, 1]")
    if (! all(0 <= alpha & alpha <= 1)) stop("alpha not in [0, 1]")

    if (length(p) > 1) {
        out <- .C("umpubinomt",
            x = as.integer(x),
            n = as.integer(n),
            alpha = as.double(alpha),
            p = as.double(p),
            np = length(p),
            maxiter = as.integer(maxiter),
            result = double(length(p)),
            tol = as.double(tol),
            PACKAGE = "ump")
    } else if (length(alpha) > 1) {
        out <- .C("umpubinoma",
            x = as.integer(x),
            n = as.integer(n),
            alpha = as.double(alpha),
            nalpha = length(alpha),
            p = as.double(p),
            maxiter = as.integer(maxiter),
            result = double(length(alpha)),
            tol = as.double(tol),
            PACKAGE = "ump")
    } else {
        out <- .C("umpubinomx",
            x = as.integer(x),
            nx = length(x),
            n = as.integer(n),
            alpha = as.double(alpha),
            p = as.double(p),
            maxiter = as.integer(maxiter),
            result = double(length(x)),
            tol = as.double(tol),
            PACKAGE = "ump")
    }
    return(out$result)
}
.First.lib <- function(lib, pkg)
{
    library.dynam("ump", pkg, lib)
}
