dansari <- function(x, m ,n)
{
    .C("dansari", as.integer(length(x)), d = as.double(x), as.integer(m),
       as.integer(n), PACKAGE="ctest")$d
}
# $Id: dpqperm.R,v 1.12 2001/12/08 14:24:34 hothorn Exp $

toltest <- function(x, scores, m) 
{
    a <- sort(scores*x - round(scores*x))
    upper <- a[1:m]       
    upper <- sum(abs(upper[upper < 0]))
    lower <- a[(length(a) - m):length(a)]
    lower <- sum(abs(lower[lower > 0]))
    max(upper, lower)/x   
}

dperm <- function(x, scores, m, paired = NULL, tol = 0.01, fact = NULL)
{
    if (is.null(x)) stop("Non-numeric argument to mathematical function")
    eq <- equiscores(scores, m, tol, fact)
    cp <- cperm(eq, m, paired)
    RVAL <- rep(0, length(x))
    RVAL[x %in% cp$T ] <- cp$Prob[cp$T %in% x]
    return(RVAL) 
}


pperm <- function(q, scores, m, paired = NULL, tol = 0.01, fact = NULL,
                     alternative=c("less", "greater", "two.sided"),
                     pprob=FALSE)
{
    if(is.null(q)) stop("Non-numeric argument to mathematical function")
    alternative <- match.arg(alternative)
    if (is.null(paired))  paired <- (length(scores) == m)
    eq <- equiscores(scores, m, tol, fact)
    cp <- cperm(eq, m, paired)
    PVALUE <- c()
    PPROB <- c()
    for (i in q) {
        if(alternative=="less")
            prob <- cp$Prob[cp$T <= i]
        if(alternative=="greater")
            prob <- cp$Prob[cp$T >= i]
        if(alternative=="two.sided" & paired) {
            if(2*sum(cp$Prob[cp$T <= i]) < 1)
                prob <- 2*cp$Prob[cp$T <= i]
            else
                prob <- 2*cp$Prob[cp$T >= i]
        }
        if(alternative=="two.sided" & !paired) {
            expect <- m/length(scores)*sum(scores)
            TH <- cp$T - expect
            ih <- i - expect
            prob <- c(cp$Prob[TH <= ifelse(ih > 0, -ih, ih)],
                      cp$Prob[TH >= ifelse(ih >= 0, ih, -ih)])
        }
    PVALUE <- c(PVALUE, sum(prob))
    PPROB <- c(PPROB, ifelse(!is.null(cp$Prob[cp$T == i]),
                                      cp$Prob[cp$T ==i], 0)) 
    }
    if(pprob)
        return(list(PVALUE = PVALUE, PPROB = PPROB))
    else
        return(PVALUE)
} 


qperm <- function(p, scores, m, paired = NULL, tol = 0.01, fact = NULL)
{
    if (is.null(p)) stop("Non-numeric argument to mathematical function")
    if (any(p < 0) || any(p > 1)) {
        warning("p is not a probability")
        return(NaN)
    }
    eq <- equiscores(scores, m, tol, fact)                          
    cp <- cperm(eq, m, paired)
    cs <- cumsum(cp$Prob)
    RVAL <- c()
    for (i in p) {
        quant <- which(cs < i)
        if (length(quant) == 0) quant <- 0
        RVAL <- c(RVAL, cp$T[max(quant) + 1])
    }
    return(RVAL)
}

rperm <- function(n, scores, m)
    sapply(1:n, dummy <- function(x) sum(sample(scores,m)))

equiscores <- function(scores, m, tol = 0.01, fact=NULL)
{
    if (any(is.null(scores))) 
        stop("Non-numeric argument to mathematical function")
    if (is.null(m)) 
        stop("Non-numeric argument to mathematical function")
    if (m < 1) 
        stop("m less than 1")
    if (m > length(scores)) 
        stop("m greater length(scores)")

    fscore <- scores - floor(scores)
    
    if (all(fscore == 0))
    { 
        # integer valued scores
        fact <- 1
        add <- min(scores) - 1
	scores <- scores - add
    } else {
        if (all(fscore[fscore != 0] == 0.5))
        {
            # midranked scores
            fact <- 2
            scores <- scores*fact
            add <- min(scores) - 1
            scores <- scores - add
        } else {
            # rational or real scores
            ssc <- sort(scores)
            b <- min(ssc[2:length(ssc)] - ssc[1:(length(ssc)-1)])
            if (b > 0) b <- ceiling(1/b) else b <- 100
            if (is.null(fact) || fact < b ) {
                if (toltest(b, scores, m) <= tol)
                    fact <- b
                else {
                    # do not induce more than 20.000 columns
                    maxfact <- min(1000, round(20000/sum(scores -
                                                    min(scores))))   
                    if (maxfact < b)
                        fact <- maxfact
                    else {
                        test <- function(x, tol, sc)
                                ifelse(toltest(x, sc, m) - tol > 0, 1/x, x)
                        fact <- optim(10, test, tol=tol, sc=scores)$par
                        if (fact > maxfact) {
                            fact <- maxfact
                            warning(paste("cannot hold tol, tolerance:",
                                     round(toltest(maxfact, scores, m), 6)))
                        } else fact <- min(fact)
                    }
                }
            } 
            scores <- round(scores * fact)
            add <- min(scores)-1
            scores <- scores - add
        }
    }

    RVAL <- list(scores = scores, fact = fact, add = add)
    class(RVAL) <- "equis"
    return(RVAL)
}


cperm <- function(escores, m, paired = NULL)
{
    if (!(class(escores) == "equis"))
        stop("scores are not of class equis") 

    N <- length(escores$scores)

    prob <- rep(0, max(cumsum(escores$scores)))

    if (is.null(paired))
        paired <- (N == m)
    else 
        paired <- (N == m) && paired  

    if (paired) {
        # paired two sample situation
        prob <- c(0, prob)
        prob <- .C("cpermdist1", prob = as.double(prob),
                   as.integer(escores$scores), as.integer(N))$prob
        t <- which(prob != 0)
        prob <- prob[t]
        # 0 is possible
	t <- t - 1
    } else {
        # independent samples
        col <- sum(sort(escores$scores)[(N + 1 - m):N])
        scores <- rep(1, N)
        prob <- .C("cpermdist2", prob = as.double(prob), as.integer(m),
                as.integer(col), as.integer(scores), as.integer(escores$scores),
                as.integer(N), as.integer(1))$prob
        t <- which(prob != 0)
        prob <- prob[t]
    }
    t <- (t + escores$add*m)/escores$fact
    RVAL <- list(T = t, Prob = prob)
    class(RVAL) <- "cperm"
    return(RVAL)
}
    # $Id: perm.test.R,v 1.7 2001/12/08 14:43:37 hothorn Exp $

perm.test <- function(x, ...) UseMethod("perm.test")

perm.test.default <-
function(x,y=NULL, paired = FALSE, alternative = c("two.sided", "less", "greater"),
         mu = 0, exact=NULL, conf.int = FALSE, conf.level = 0.95, tol=NULL, ...) {

    if (is.null(x)) stop("x is missing")
    if (is.null(y)) paired <- TRUE

    alternative <- match.arg(alternative)
    if(!missing(mu) && ((length(mu) > 1) || !is.finite(mu)))
        stop("mu must be a single number")
    if(conf.int) {
        if(!((length(conf.level) == 1)
           && is.finite(conf.level)
           && (conf.level < 1)))
           stop("conf.level must be a single number between 0 and 1")
    }

    MIDP <- NULL
    if (!is.null(y))
        DNAME <- paste(deparse(substitute(x)), "and",
                     deparse(substitute(y)))
    else DNAME <- deparse(substitute(x))

    if (paired) {
        if (!is.null(y)) x <- x - y
        m <- length(x)
        METHOD <- "1-sample Permutation Test"
        if (is.null(exact)) exact <- (m <= 50)
        if (any(x != floor(x)) && exact) {
            if (is.null(tol)) {
                x <- x - min(x)
                x <- round(x*m/max(x))
                METHOD <- paste(METHOD, "(scores mapped into 1:m)")
            }
            METHOD <- paste(METHOD, "using rounded scores")
        }

        STATISTIC <- sum(x[x > 0])
        if (exact) {
            PVAL <- switch(alternative,
               "two.sided" = {
                    pperm(STATISTIC, abs(x), m, alternative="two.sided",
                          pprob=TRUE, tol=tol)
                }, "greater" = {
                    pperm(STATISTIC, abs(x), m, alternative="greater",
                          pprob=TRUE, tol=tol)
                }, "less" = {
                    pperm(STATISTIC, abs(x), m, alternative="less", pprob=TRUE,
                          tol=tol)
                })
            MIDP <- PVAL$PPROB
            PVAL <- PVAL$PVALUE
        } else {
            METHOD <- paste("Asymptotic", METHOD)
            wmean <- sum(abs(x))/2
            wvar <- sum(abs(x)^2)/4
            PVAL <- pnorm((STATISTIC - wmean) / sqrt(wvar))
            PVAL <- min(PVAL, 1-PVAL)
            if(alternative == "two.sided")
                PVAL <- 2 * PVAL
            if (conf.int) { 
                conf.int <- FALSE
                warning("cannot compute asymptotic confidence intervals")
            }
        }
        if (conf.int) {
        warning("Cannot compute confidence interval for paired data!")
        # I do not find any related reference, can R"ohmel 96 be adapted?
        if (FALSE)
            alpha <- 1-conf.level
            xscores <- equiscores(x,m)

            Hx <- rep(0, sum(xscores$scores)*m)

            Hx <- .C("cpermdist2", H = as.double(Hx), as.integer(m),
                as.integer(sum(xscores$scores)), as.integer(rep(1,m)),
                as.integer(xscores$scores), as.integer(m),
                as.integer(length(Hx)+1))$H

            Hx <- matrix(Hx, nrow=m, byrow=TRUE)
            Hx <- rbind(Hx, 0)
            Hx[nrow(Hx), ncol(Hx)] <- 1

            D <- c()
  
            for (i in 1:m) {
                xmean <- which(Hx[i+1,] > 0) - 1
                xm <- length(xmean)
                xmeans <- ((xmean + xscores$add*i)/xscores$fact)/i
                diff <- outer(xmeans, xmeans, "+")
                diff <- sort(diff[!lower.tri(diff)]) / 2
                count <- as.vector(outer(Hx[i+1, xmean+1], Hx[i+1, xmean+1], "*")) 
                od <- order(diff)
                count <- count[od]
                diff <- diff[od]
                du <- duplicated(diff)
                for (k in length(count):1)
                    if (du[k]) count[k-1] <- count[k-1] + count[k]

                diff <- unique(diff)
                count <- count[!du]

                old <- diff %in% D[,1]
                Dold <- D[,1] %in% diff
 
                if (any(Dold))
                    D[Dold,2] <- D[Dold,2] + count[old]
  
                if (any(!old))
                    D <- rbind(D, cbind(diff[!old], count[!old]))

                D <- D[order(D[,1]),]
            } 

            od <- order(D[,1])
            Dcount <- cumsum(D[od,2])
            Dsort <- D[od,1]

            L <- sum(D[,2])

            cint <- switch(alternative,
                "two.sided" = {
                    qu <- floor((1-alpha/2)*L)             
                    qu <- which(abs(Dcount - qu) == min(abs(Dcount - qu)))
                    ql <- ceiling(alpha/2*L)
                    ql <- which(abs(Dcount - ql) == min(abs(Dcount - ql)))
                    c(Dsort[ql], Dsort[qu])
                }, "greater" = {
                    qu <- floor((1-alpha)*L)
                    qu <- which(abs(Dcount - qu) == min(abs(Dcount - qu)))
                    c(-Inf, Dsort[qu])
                }, "less" = {
                    ql <- ceiling(alpha*L)
                    ql <- which(abs(Dcount - ql) == min(abs(Dcount - ql)))
                    c(Dsort[ql], Inf)
                })
            attr(cint, "conf.level") <- conf.level
        }
#    }
     } else {
        m <- length(x)
        n <- length(y)
        x <- x - mu
        if (is.null(exact)) exact <- (m <= 50 && n <= 50)
        cxy <- c(x,y)

        METHOD <- "2-sample Permutation Test"
        if (any(cxy != floor(cxy)) && exact) {
            if (is.null(tol)) {
                cxy <- cxyx - min(cxyx)
                cxy <- round(cxy*(n+m)/max(cxy))
                METHOD <- paste(METHOD, "(scores mapped into 1:(m+n))")
	        x <- cxy[seq(along=x)]
	        y <- cxy[-seq(along=x)]
            }
            METHOD <- paste(METHOD, "using rounded scores")
        }

        STATISTIC <- sum(x)
        if (exact) {
            PVAL <- switch(alternative,
                "two.sided" = {
                    pperm(STATISTIC, c(x,y), m, alternative="two.sided",
                          pprob=TRUE, tol=tol)
                }, "greater" = {
                    pperm(STATISTIC, c(x,y), m, alternative="greater",
                          pprob=TRUE, tol=tol)
                }, "less" = {
                    pperm(STATISTIC, c(x,y), m, alternative="less", pprob=TRUE,
                          tol=tol)
                })
            MIDP <- PVAL$PPROB
            PVAL <- PVAL$PVALUE
        } else {
            METHOD <- paste("Asymptotic", METHOD)
            N <- m + n
            wmean <- m/N*sum(x)
            wvar <- m*n/(N*(N-1))*sum((x - wmean/m)^2)
            PVAL <- pnorm((STATISTIC - wmean)/sqrt(wvar))
            PVAL <- min(PVAL, 1-PVAL)
            if(alternative == "two.sided")
                PVAL <- 2 * min(PVAL, 1 - PVAL)
            if (conf.int) {
                conf.int <- FALSE
                    warning("cannot compute asymptotic confidence intervals")
            }
        }
        if (conf.int) {
            alpha <- 1-conf.level
            xscores <- equiscores(x,m)
            yscores <- equiscores(y,n)

            Hx <- rep(0, sum(xscores$scores)*m)
            Hy <- rep(0, sum(yscores$scores)*n)

            Hx <- .C("cpermdist2", H = as.double(Hx), as.integer(m),
                as.integer(sum(xscores$scores)), as.integer(rep(1,m)),
                as.integer(xscores$scores), as.integer(m),
                as.integer(length(Hx)+1))$H

            Hx <- matrix(Hx, nrow=m, byrow=TRUE)
            Hx <- rbind(Hx, 0)
            Hx[nrow(Hx), ncol(Hx)] <- 1

            Hy <- .C("cpermdist2", H = as.double(Hy), as.integer(n),
                as.integer(sum(yscores$scores)), as.integer(rep(1,n)),
                as.integer(yscores$scores), as.integer(n),
                as.integer(length(Hy)+1))$H

            Hy <- matrix(Hy, nrow=n, byrow=TRUE)
            Hy <- rbind(Hy, 0)
            Hy[nrow(Hy), ncol(Hy)] <- 1

            L <- sum(choose(m, 1:m)*choose(n, 1:m))
            D <- c()

            for (i in 1:min(m,n)) {
                xmean <- which(Hx[i+1,] > 0) - 1
                ymean <- which(Hy[i+1,] > 0) - 1
                xm <- length(xmean)
                ym <- length(ymean)

                diff <- as.vector(outer(((xmean + xscores$add*i)/xscores$fact)/i,
                       ((ymean + yscores$add*i)/yscores$fact)/i, "-"))
                count <- as.vector(outer(Hx[i+1, xmean+1], Hy[i+1, ymean+1], "*")) 
                od <- order(diff)
                count <- count[od]
                diff <- diff[od]
                du <- duplicated(diff)
                for (k in length(count):1)
                    if (du[k]) count[k-1] <- count[k-1] + count[k]

                diff <- unique(diff)
                count <- count[!du]

                old <- diff %in% D[,1]
                Dold <- D[,1] %in% diff

                if (any(Dold))
                     D[Dold,2] <- D[Dold,2] + count[old]

                if (any(!old))
                    D <- rbind(D, cbind(diff[!old], count[!old]))

                D <- D[order(D[,1]),]
            }
            od <- order(D[,1])
            Dcount <- cumsum(D[od,2])
            Dsort <- D[od,1]
 
            cint <- switch(alternative,
                "two.sided" = {
                    qu <- floor((1-alpha/2)*L)             
                    qu <- which(abs(Dcount - qu) == min(abs(Dcount - qu)))
                    ql <- ceiling(alpha/2*L)
                    ql <- which(abs(Dcount - ql) == min(abs(Dcount - ql)))
                    c(Dsort[ql], Dsort[qu])
                }, "greater" = {
                    qu <- floor((1-alpha)*L)
                    qu <- which(abs(Dcount - qu) == min(abs(Dcount - qu)))
                    c(-Inf, Dsort[qu])
                }, "less" = {
                    ql <- ceiling(alpha*L)
                    ql <- which(abs(Dcount - ql) == min(abs(Dcount - ql)))
                    c(Dsort[ql], Inf)
                })
            attr(cint, "conf.level") <- conf.level
        }
    }
    if (exact) {
        names(MIDP) <- "point prob"
        RVAL <- list(statistic = STATISTIC,
                     p.value = PVAL,
                     pointprob = MIDP,
                     null.value = c(mu = mu),
                     alternative = alternative,
                     method = METHOD,
                     data.name = DNAME)
        if(conf.int)
             RVAL$conf.int <- cint
    } else {
        RVAL <- list(statistic = STATISTIC,
                     p.value = PVAL,
                     null.value = c(mu = mu),
                     alternative = alternative,
                     method = METHOD,
                     data.name = DNAME)
    }
    class(RVAL) <- "htest"
    return(RVAL)
}

perm.test.formula <-
function(formula, data, subset, na.action, ...)
{
    if(missing(formula)
       || (length(formula) != 3)  
       || (length(attr(terms(formula[-2]), "term.labels")) != 1)
       || (length(attr(terms(formula[-3]), "term.labels")) != 1))
        stop("formula missing or incorrect")
    if(missing(na.action))
        na.action <- getOption("na.action")
    m <- match.call(expand.dots = FALSE)
    if(is.matrix(eval(m$data, parent.frame())))
        m$data <- as.data.frame(data) 
    m[[1]] <- as.name("model.frame")
    m$... <- NULL
    mf <- eval(m, parent.frame())
    DNAME <- paste(names(mf), collapse = " by ")
    names(mf) <- NULL
    response <- attr(attr(mf, "terms"), "response")
    g <- factor(mf[[-response]])
    if(nlevels(g) != 2)
        stop("grouping factor must have exactly 2 levels")
    DATA <- split(mf[[response]], g)
    names(DATA) <- c("x", "y")
    y <- do.call("perm.test", c(DATA, list(...)))
    y$data.name <- DNAME
    y
}
# $Id: wilcox.exact.R,v 1.15 2001/12/08 13:45:14 hothorn Exp $

wilcox.exact <- function(x, ...) UseMethod("wilcox.exact")

wilcox.exact.default <-
function(x, y = NULL, alternative = c("two.sided", "less", "greater"), 
         mu = 0, paired = FALSE, exact = NULL, 
         conf.int = FALSE, conf.level = 0.95, ...) 
{
    alternative <- match.arg(alternative)
    if(!missing(mu) && ((length(mu) > 1) || !is.finite(mu))) 
        stop("mu must be a single number")
    if(conf.int) {
        if(!((length(conf.level) == 1)
             && is.finite(conf.level)
             && (conf.level > 0)
             && (conf.level < 1)))
            stop("conf.level must be a single number between 0 and 1")
    }
    MIDP <- NULL

    if(!is.null(y)) {
        DNAME <- paste(deparse(substitute(x)), "and",
                       deparse(substitute(y)))
        if(paired) {
            if(length(x) != length(y)) 
                stop("x and y must have the same length")
            OK <- complete.cases(x, y)
            x <- x[OK] - y[OK]
            y <- NULL
        }
        else {
            x <- x[is.finite(x)]
            y <- y[is.finite(y)]
        }
    } else {
        DNAME <- deparse(substitute(x))
        if(paired) 
            stop("y missing for paired test")
        x <- x[is.finite(x)]
    }

    if(length(x) < 1) 
        stop("not enough (finite) x observations")
    CORRECTION <- 0
    if(is.null(y)) {
        METHOD <- "Exact Wilcoxon signed rank test"
        x <- x - mu
        ZEROES <- any(x == 0)
        if(ZEROES) 
            x <- x[x != 0]
        n <- length(x)
        if(is.null(exact)) 
            exact <- (n < 50)
        r <- rank(abs(x))
        STATISTIC <- sum(r[x > 0])
        names(STATISTIC) <- "V"
        TIES <- (length(r) != length(unique(r)))

        if (exact) {
            PVAL <-
               switch(alternative,
                      "two.sided" = pperm(STATISTIC, r, n,
                                          alternative="two.sided", pprob=TRUE),
                      "greater" = pperm(STATISTIC, r, n,
                                        alternative="greater", pprob=TRUE),
                      "less" = pperm(STATISTIC, r, n,
                                     alternative="less", pprob=TRUE))
            MIDP <- PVAL$PPROB
            PVAL <- PVAL$PVALUE
            if(conf.int && !is.na(x)) {
                ## Exact confidence intervale for the median in the
                ## one-sample case.  When used with paired values this
                ## gives a confidence interval for mean(x) - mean(y).
                x <- x + mu             # we want a conf.int for the median
                alpha <- 1 - conf.level
                diffs <- outer(x, x, "+")
                diffs <- sort(diffs[!lower.tri(diffs)]) / 2
                if (TIES) {
                  fs <- function(d)
                    xx <- x - d; sum(rank(abs(xx))[xx > 0])
                  w <- sapply(diffs, fs)
                } else {
                  w <- sum(rank(abs(x))):1
                }
                cint <-
                    switch(alternative,
                           "two.sided" = {
                               qu <- qperm(alpha/2, r, n) 
                               ql <- qperm(1-alpha/2, r, n)
                               if (qu <= min(w)) lci <- max(diffs)
                               else lci <- min(diffs[w <= qu])
                               if (ql >= max(w)) uci <- min(diffs)
                               else uci <- max(diffs[w > ql])
                               c(uci, lci)
                           },
                           "greater"= {
                               ql <- qperm(1-alpha, r, n)
                               if (ql >= max(w)) uci <- min(diffs)
                               else uci <- max(diffs[w > ql])
                               c(uci, Inf)
                           },
                           "less"= {
                               qu <- qperm(alpha, r, n)
                               if (qu <= min(w)) lci <- max(diffs)
                               else lci <- min(diffs[w <= qu])
                               c(-Inf, lci)
                })
                attr(cint, "conf.level") <- conf.level    
                wmean <- sum(r)/2
                ESTIMATE <- mean(c(min(diffs[w <= ceiling(wmean)]),
                                   max(diffs[w > wmean])))
                names(ESTIMATE) <- "(pseudo)median"
                }
            } else {
                METHOD <- "Asymptotic Wilcoxon signed rank test"
                wmean <- sum(r)/2
                wvar <- sum(r^2)/4
                PVAL <- pnorm((STATISTIC - wmean) / sqrt(wvar))
                PVAL <- min(PVAL, 1-PVAL)
                if(alternative == "two.sided") 
                    PVAL <- 2 * PVAL 

                if(conf.int && !is.na(x)) {
                    ## Asymptotic confidence intervale for the median in the
                    ## one-sample case.  When used with paired values this
                    ## gives a confidence interval for mean(x) - mean(y).
                    ## Algorithm not published, thus better documented here.
                    x <- x + mu
                    alpha <- 1 - conf.level
                    ## These are sample based limits for the median
                    mumin <- min(x)
                    mumax <- max(x)
                    ## wdiff(d, zq) returns the absolute difference between  
                    ## the asymptotic Wilcoxon statistic of x - mu - d and
                    ## the quantile zq.
                    CORRECTION.CI <- 0
                    wdiff <- function(d, zq) {
                        xd <- x - d
                        xd <- xd[xd != 0]
                        nx <- length(xd)
                        dr <- rank(abs(xd))
                        zd <- sum(dr[xd > 0])
                        zd <- (zd - wmean) / sqrt(wvar)
                        zd - zq
                    }
                    ## Here we search the root of the function wdiff on the set
                    ## c(mumin, mumax).
                    ##
                    ## This returns a value from c(mumin, mumax) for which
                    ## the asymptotic Wilcoxon statistic is equal to the
                    ## quantile zq.  This means that the statistic is not
                    ## within the critical region, and that implies that d
                    ## is a confidence limit for the median.
                    ##
                    ## As in the exact case, interchange quantiles.
                    cint <- switch(alternative, "two.sided" = {
                        l <- uniroot(wdiff, c(mumin, mumax), tol=1e-4,
                                     zq=qnorm(alpha/2, lower=FALSE))$root
                        u <- uniroot(wdiff, c(mumin, mumax), tol=1e-4,
                                     zq=qnorm(alpha/2))$root
                        c(l, u)
                    }, "greater"= {
                        l <- uniroot(wdiff, c(mumin, mumax), tol=1e-4,
                                     zq=qnorm(alpha, lower=FALSE))$root
                        c(l, +Inf)
                    }, "less"= {
                        u <- uniroot(wdiff, c(mumin, mumax), tol=1e-4,
                                     zq=qnorm(alpha))$root
                        c(-Inf, u)
                    })
                    attr(cint, "conf.level") <- conf.level    
                    ESTIMATE <- uniroot(wdiff, c(mumin, mumax), tol=1e-4,
                                    zq=0)$root
                    names(ESTIMATE) <- "(pseudo)median"
                 }

	}
    }
    else {
        if(length(y) < 1) 
            stop("not enough y observations")
        METHOD <- "Exact Wilcoxon rank sum test"
        r <- rank(c(x - mu, y))
        n.x <- length(x)
        n.y <- length(y)
        if(is.null(exact)) 
            exact <- (n.x < 50) && (n.y < 50)
        STATISTIC <- sum(r[seq(along = x)]) - n.x * (n.x + 1) / 2
        names(STATISTIC) <- "W"
        TIES <- (length(r) != length(unique(r)))
        if(exact) {
            PVAL <-
                switch(alternative,
                       "two.sided" = pperm(STATISTIC + n.x*(n.x +1)/2, r,
                                           n.x, alternative="two.sided", pprob=TRUE),
                       "greater" = pperm(STATISTIC + n.x*(n.x+1)/2, r,
                                       n.x, alternative="greater", pprob=TRUE), 
                           "less" = pperm(STATISTIC+ n.x*(n.x +1)/2, r, n.x,
                                          alternative="less", pprob=TRUE))
             MIDP <- PVAL$PPROB    
             PVAL <- PVAL$PVALUE
             if(conf.int) {
                 ## Exact confidence interval for the location parameter 
                 ## mean(y) - mean(x) in the two-sample case (cf. the
                 ## one-sample case).
                 if (mu != 0) r <- rank(c(x,y))
                 alpha <- 1 - conf.level
                 diffs <- sort(outer(x, y, "-"))
                 if (TIES) {
                   fs <- function(d)
                     sum(rank(c(x-d,y))[seq(along = x)]) - n.x * (n.x + 1) / 2
                   w <- sapply(diffs, fs)
                 } else {
                   w <- (n.x*n.y):1
                 }
                 cint <-
                     switch(alternative,
                            "two.sided" = {
                                qu <- qperm(alpha/2, r, n.x) - n.x*(n.x+1)/2
                                ql <- qperm(1-alpha/2, r, n.x) - n.x*(n.x+1)/2
                                if (qu <= min(w)) lci <- max(diffs)
                                else lci <- min(diffs[w <= qu])
                                if (ql >= max(w)) uci <- min(diffs)
                                else uci <- max(diffs[w > ql])
                                c(uci, lci)
                            },
                            "greater"= {
                                ql <- qperm(1-alpha, r, n.x) - n.x*(n.x+1)/2
                                if (ql >= max(w)) uci <- min(diffs)
                                else uci <- max(diffs[w > ql])
                                c(uci, +Inf)
                             },
                             "less"= {
                                 qu <- qperm(alpha, r, n.x) - n.x*(n.x+1)/2
                                 if (qu <= min(w)) lci <- max(diffs)
                                 else lci <- min(diffs[w <= qu])
                                 c(-Inf, lci)
                              })
                  attr(cint, "conf.level") <- conf.level    
                  wmean <- n.x/(n.x+n.y)*sum(r) - n.x*(n.x+1)/2
                  ESTIMATE <- mean(c(min(diffs[w <= ceiling(wmean)]),
                                     max(diffs[w > wmean])))
                  names(ESTIMATE) <- "difference in location"
             }
        } else {
            METHOD <- "Asymptotic Wilcoxon rank sum test"
            N <- n.x + n.y
            wmean <- n.x/N*sum(r)
            wvar <- n.x*n.y/(N*(N-1))*sum((r - wmean/n.x)^2)
            PVAL <- pnorm((STATISTIC + n.x*(n.x+1)/2 - wmean)/sqrt(wvar))
            PVAL <- min(PVAL, 1-PVAL)
            if(alternative == "two.sided") 
                PVAL <- 2 * min(PVAL, 1 - PVAL)
            if(conf.int) {
                ## Asymptotic confidence interval for the location
                ## parameter mean(x) - mean(y) in the two-sample case
                ## (cf. one-sample case).
                ##
                ## Algorithm not published, for a documentation see the
                ## one-sample case.
                alpha <- 1 - conf.level
                mumin <- min(x) - max(y)
                mumax <- max(x) - min(y)
                CORRECTION.CI <- 0
                wdiff <- function(d, zq) {
                    dr <- rank(c(x - d, y))
                    dz <- sum(dr[seq(along = x)])
                    dz <- (dz - wmean) / sqrt(wvar)
                    dz - zq
                }
                cint <- switch(alternative, "two.sided" = {
                    l <- uniroot(wdiff, c(mumin, mumax), tol=1e-4,
                                  zq=qnorm(alpha/2, lower=FALSE))$root
                    u <- uniroot(wdiff, c(mumin, mumax), tol=1e-4,
                                  zq=qnorm(alpha/2))$root
                    c(l, u)
                }, "greater"= {
                    l <- uniroot(wdiff, c(mumin, mumax), tol=1e-4,
                                  zq=qnorm(alpha, lower=FALSE))$root
                    c(l, +Inf)
                }, "less"= {
                    u <- uniroot(wdiff, c(mumin, mumax), tol=1e-4,
                                  zq=qnorm(alpha))$root
                    c(-Inf, u)
                })
                attr(cint, "conf.level") <- conf.level    
                ESTIMATE <- uniroot(wdiff, c(mumin, mumax), tol=1e-4,
                                    zq=0)$root
                names(ESTIMATE) <- "difference in location"
            }
        }
    }

    if (!is.null(MIDP)) {
        names(MIDP) <- "point prob"
        RVAL <- list(statistic = STATISTIC,
                     pointprob = MIDP,
                     p.value = PVAL, 
                     null.value = c(mu = mu),
                     alternative = alternative,
                     method = METHOD, 
                     data.name = DNAME)
    } else {
        RVAL <- list(statistic = STATISTIC,
                     p.value = PVAL, 
                     null.value = c(mu = mu),  
                     alternative = alternative,
                     method = METHOD, 
                     data.name = DNAME)
    }
    if(conf.int) {
        RVAL$conf.int <- cint
        RVAL$estimate <- ESTIMATE
    }
    class(RVAL) <- "htest"
    return(RVAL)
}

wilcox.exact.formula <-
function(formula, data, subset, na.action, ...)
{
    if(missing(formula)
       || (length(formula) != 3)  
       || (length(attr(terms(formula[-2]), "term.labels")) != 1)
       || (length(attr(terms(formula[-3]), "term.labels")) != 1))
        stop("formula missing or incorrect")
    if(missing(na.action))
        na.action <- getOption("na.action")
    m <- match.call(expand.dots = FALSE)
    if(is.matrix(eval(m$data, parent.frame())))
        m$data <- as.data.frame(data) 
    m[[1]] <- as.name("model.frame")
    m$... <- NULL
    mf <- eval(m, parent.frame())
    DNAME <- paste(names(mf), collapse = " by ")
    names(mf) <- NULL
    response <- attr(attr(mf, "terms"), "response")
    g <- factor(mf[[-response]])
    if(nlevels(g) != 2)
        stop("grouping factor must have exactly 2 levels")
    DATA <- split(mf[[response]], g)
    names(DATA) <- c("x", "y")
    y <- do.call("wilcox.exact", c(DATA, list(...)))
    y$data.name <- DNAME
    y
}
.First.lib <- function(lib, pkg) {
  require(ctest)
  library.dynam("exactRankTests", pkg, lib)
}
