.packageName <- "exactRankTests"
ansari.exact <- function(x, ...) UseMethod("ansari.exact")

ansari.exact.default <-
function(x, y, alternative = c("two.sided", "less", "greater"),
         exact = NULL, conf.int = FALSE, conf.level = 0.95, ...) 
{
    alternative <- match.arg(alternative)
    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")
    }
    DNAME <- paste(deparse(substitute(x)), "and", deparse(substitute(y)))

    if(!is.numeric(x)) stop("`x' must be numeric")
    if(!is.numeric(y)) stop("`y' must be numeric")

    x <- x[complete.cases(x)]
    y <- y[complete.cases(y)]
    m <- length(x)
    if(m < 1)
        stop("not enough x observations")
    n <- length(y)
    if(n < 1)
        stop("not enough y observations")
    N <- m + n

    r <- cscores(c(x, y), type="Ansari")
    STATISTIC <- sum(r[seq(along = x)])

    if(is.null(exact))
        exact <- ((m < 50) && (n < 50))

    if(exact) {
        PVAL <-
            switch(alternative,
                   two.sided = pperm(STATISTIC, r, m, alternative="two.sided"),
                   less = pperm(STATISTIC, r, m, alternative="less"),
                   greater = pperm(STATISTIC, r, m, alternative="greater"))
        if (conf.int) {
            alpha <- 1 - conf.level
            x <- sort(x)
            y <- sort(y)
            ab <- function(sig) {
                rab <- rank(c(x/sig, y))
                sum(pmin(rab, N - rab + 1)[seq(along = x)])
            }
            ratio <- outer(x,y,"/")
            aratio <- ratio[ratio >= 0]
            sigma <- sort(aratio)

            cci <- function(alpha) {
              u <- absigma - qperm(alpha/2,  r, m) 
              l <- absigma - qperm(1 - alpha/2, r, m) 
              ## Check if the statistic exceeds both quantiles first.
              uci <- NULL
              lci <- NULL                    
              if(length(u[u >= 0]) == 0 || length(l[l > 0]) == 0) {
                  warning(paste("Samples differ in location: Cannot",
                                "compute confidence set, returning NA"))
                  return(c(NA, NA))
              } 
              if (is.null(uci)) {
                  u[u < 0] <- NA
                  uci <- min(sigma[which(u == min(u, na.rm=TRUE))])
              }
              if (is.null(lci)) {
                  l[l <= 0] <- NA
                  lci <- max(sigma[which(l == min(l, na.rm=TRUE))])
              }
              ## The process of the statistics does not need to be
              ## monotone in sigma: check this and interchange quantiles.
              if (uci > lci) {
                  l <- absigma - qperm(alpha/2,  r, m)
                  u <- absigma - qperm(1 - alpha/2, r, m)
                  u[u < 0] <- NA
                  uci <- min(sigma[which(u == min(u, na.rm=TRUE))])
                  l[l <= 0] <- NA
                  lci <- max(sigma[which(l == min(l, na.rm=TRUE))])
               }
               c(uci, lci)
            }

            cint <- if(length(sigma) < 1) {
                warning("Cannot compute confidence set, returning NA")
                c(NA, NA)
            }
            else {
                ## Compute statistics directly: computing the steps is
                ## not faster.
                absigma <-
                    sapply(sigma + c(diff(sigma)/2,
                                     sigma[length(sigma)]*1.01), ab)
                switch(alternative, two.sided = {
                    cci(alpha)
                }, greater= {
                    c(cci(alpha*2)[1], Inf)
                }, less= {
                    c(0, cci(alpha*2)[2])
                })
            }
            attr(cint, "conf.level") <- conf.level
            u <- absigma - qperm(0.5, r, m)
            sgr <- sigma[u <= 0]
            if (length(sgr) == 0) sgr <- NA
            else sgr <- max(sgr)
            sle <- sigma[u > 0]
            if (length(sle) == 0) sle <- NA
            else sle <- min(sle)
            ESTIMATE <- mean(c(sle, sgr))
        }
    }
    else {
        normalize <- function(s, r) {
          ss <- sum(r)
          N <- m + n 
          E <- m/N*ss
          V <- m*(N-m)/(N^2*(N-1))*(N*sum(r^2) - ss^2)
          (s - E)/sqrt(V)
        }
        p <- pnorm(normalize(STATISTIC, r))
        PVAL <- switch(alternative,
                       two.sided = 2 * min(p, 1 - p),
                       less = 1 - p,
                       greater = p)
    
        if(conf.int && !exact) {
            alpha <- 1 - conf.level
            ab <- function(sig, zq) {
                r <- rank(c(x / sig, y))
                s <- sum(pmin(r, N -r + 1)[seq(along = x)])
                normalize(s, r) - zq
            }
            ## Use uniroot here.
            ## Compute the range of sigma first.
            srangepos <- NULL
            srangeneg <- NULL
            if (any(x[x > 0]) && any(y[y > 0]))
                srangepos <-
                    c(min(x[x>0], na.rm=TRUE)/max(y[y>0], na.rm=TRUE), 
                      max(x[x>0], na.rm=TRUE)/min(y[y>0], na.rm=TRUE))
            if (any(x[x <= 0]) && any(y[y < 0]))
                srangeneg <-
                    c(min(x[x<=0], na.rm=TRUE)/max(y[y<0], na.rm=TRUE), 
                      max(x[x<=0], na.rm=TRUE)/min(y[y<0], na.rm=TRUE))
            if (any(is.infinite(c(srangepos, srangeneg)))) {
                warning(paste("Cannot compute asymptotic confidence",
                              "set or estimator"))
                conf.int <- FALSE
            } else {
                ccia <- function(alpha) {
                    ## Check if the statistic exceeds both quantiles
                    ## first.
                    statu <- ab(srange[1], zq=qnorm(alpha/2))
                    statl <- ab(srange[2], zq=qnorm(alpha/2, lower=FALSE))
                    if (statu > 0 || statl < 0) {
                        warning(paste("Samples differ in location:",
                                      "Cannot compute confidence set,",
                                      "returning NA"))
                        return(c(NA, NA))
                    }
                    u <- uniroot(ab, srange, tol=1e-4,
                                 zq=qnorm(alpha/2))$root
                    l <- uniroot(ab, srange, tol=1e-4,
                                 zq=qnorm(alpha/2, lower=FALSE))$root
                    ## The process of the statistics does not need to be
                    ## monotone: sort is ok here.
                    sort(c(u, l))
                }
                srange <- range(c(srangepos, srangeneg), na.rm=FALSE)
                cint <- switch(alternative, two.sided = {
                    ccia(alpha)
                }, greater= {
                    c(ccia(alpha*2)[1], Inf)
                }, less= {
                    c(0, ccia(alpha*2)[2])
                })
                attr(cint, "conf.level") <- conf.level
                ## Check if the statistic exceeds both quantiles first.
                statu <- ab(srange[1], zq=0)
                statl <- ab(srange[2], zq=0)
                if (statu > 0 || statl < 0) {
                    ESTIMATE <- NA
                    warning("Cannot compute estimate, returning NA")
                } else
                    ESTIMATE <- uniroot(ab, srange, tol=1e-4, zq=0)$root
            }
        }
    }
    
    names(STATISTIC) <- "AB"
    RVAL <- list(statistic = STATISTIC,
                 p.value = PVAL,
                 null.value = c("ratio of scales" = 1),
                 alternative = alternative,
                 method = "Ansari-Bradley test",
                 data.name = DNAME)
    if(conf.int)
        RVAL <- c(RVAL,
                  list(conf.int = cint,
                       estimate = c("ratio of scales" = ESTIMATE)))
    class(RVAL) <- "htest"
    return(RVAL)
}

ansari.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("ansari.exact", c(DATA, list(...)))
    y$data.name <- DNAME
    y
}
# $Id: cscores.R,v 1.14.4.1 2003/11/03 15:16:50 hothorn Exp $

cscores <- function(y, ...) UseMethod("cscores")

cscores.default <- function(y, type=c("Data", "Wilcoxon", 
                            "NormalQuantile", "AnsariBradley", "Median",
                            "Savage", "ConSal"), int = FALSE, 
                            maxs=length(y), ... ) {
  type <- match.arg(type)
  if (!(all.equal(floor(maxs),maxs))  || maxs < 1) 
    stop("maxs is not an positiv integer")
  N <- length(y)
  RET <- switch(type, "Data" = y,
           "Wilcoxon" = rank(y),
           "NormalQuantile" = qnorm(rank(y)/(N+1)),
           "AnsariBradley" = { 
             r <- rank(y)
             pmin(r, N - r + 1) },
           "Median" = {
             r <- rank(y)
             r[r <= (N+1)/2] <- 0
             r[r > 0] <- 1
             r},
           "Savage" = {
             cscores.Surv(cbind(y, 1)) },
           "ConSal" = {
             (rank(y)/(N+1))^4
           }
  )
  attr(RET, "scores") <- type
  if (int) {
    fscore <- RET - floor(RET)
    if (all(fscore[fscore != 0] == 0.5)) 
      RET <- 2*RET
    else RET <- round(RET*maxs/max(RET))
  }
  RET
}

cscores.factor <- function(y, ...) {
  if (nlevels(y) > 2) stop("cannot compute scores for more than 2 levels")
  RET <- as.integer(y)
  attr(RET, "scores") <- "Median"
  RET
}

irank <- function(x, ox=NULL) {
  if (is.null(ox))
    .Call("irank", as.double(x), as.integer(order(x)-1), 
          PACKAGE="exactRankTests")
  else 
    .Call("irank", as.double(x), as.integer(ox-1), 
          PACKAGE="exactRankTests")
}

cscores.Surv <- function(y, type="LogRank", int=FALSE, maxs=nrow(y), ...) {
  type <- match.arg(type)
  if (!(all.equal(floor(maxs),maxs))  || maxs < 1)
    stop("maxs is not an positiv integer")
  time <- y[,1]
  event <- y[,2]
  N <- length(time)
  ot <- order(time)
  rt <- irank(time, ot)
  fact <- event/(N - rt + 1)
  scores <- event - cumsum(fact[ot])[rt]
  if (int)
    RET <- round(scores*maxs/max(scores))
  else
    RET <- scores
  attr(RET, "scores") <- "LogRank"
  RET
}

# $Id: dpqperm.R,v 1.24.4.1 2004/09/27 06:07:11 hothorn Exp $

findfact <- function(scores, m, tol) {
    sc <- rev(sort(scores))
    foo <- function(fact) {
      abs(sum(abs(fact*sc[1:m] - round(fact*sc[1:m],0)))/fact - tol)
    }
    grid <- seq(from = 1, to = 10000/(sum(sc) -min(sc)), by = 1)
    fact <- which.min(sapply(grid, foo))
    ## fact <- optimize(foo, c(1, 100000/(sum(sc) -min(sc))))$minimum
    thistol <- sum(abs(fact*sc[1:m] - round(fact*sc[1:m],0)))/fact
    if (thistol > tol)
      warning(paste("cannot hold tol, tolerance:", round(thistol, 6)))
    fact
}


dperm <- function(x, scores, m, paired = NULL, tol = 0.01, fact = NULL,
                  density=FALSE, simulate=FALSE, B=10000)
{
    if (is.null(x)) stop("Non-numeric argument to mathematical function")
    # use raw scores for Monte-Carle simulations
    eq <- equiscores(scores, m, tol, fact, simulate)
    if (simulate) 
      cp <- cperm(eq, m, paired, B=B)
    else 
      cp <- cperm(eq, m, paired)
    # return everything up to x
    if (density & length(x) == 1) {
      indx <- which(cp$T <= x)
      RVAL <- data.frame(cp$T[indx], cp$Prob[indx])
    } else {
      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, simulate=FALSE, B=10000)
{
    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, simulate)
    if (simulate) 
      cp <- cperm(eq, m, paired, B=B)
    else 
      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
            if (ih != 0) 
              prob <- c(cp$Prob[TH <= ifelse(ih > 0, -ih, ih)],
                        cp$Prob[TH >= ifelse(ih >= 0, ih, -ih)])
            else
              prob <- 1
        }
    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,
                  simulate = FALSE, B=10000)
{
    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, simulate)
    if (simulate) 
      cp <- cperm(eq, m, paired, B=B)
    else 
      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=NULL, tol = 0.01, fact=NULL, simulate=FALSE)
{
  if (any(is.null(scores))) 
      stop("Non-numeric argument to mathematical function")

  if (!is.null(m)) { 
    if (m < 1) 
        stop("m less than 1")
    if (m > length(scores)) 
        stop("m greater length(scores)")
    paired <- (length(scores) == m)
  } else {
    paired <- FALSE
  }

  # use raw scores for simulation based computation of the permutation 
  # distribution
  if (simulate) {
    if(paired & any(scores < 0)) stop("scores must be positive")
    RVAL <- list(scores = scores, fact = 1, add = 0)
    class(RVAL) <- "equis"
    return(RVAL)
  }  

  if (is.null(fact)) {
    # first, handle integer and midranked scores, 
    # m is not needed here
    fact <- 1
    fscore <- scores - floor(scores)
    if (any(fscore != 0)) {
      if (all(fscore[fscore != 0] == 0.5)) # midranks
        fact <- 2
    }
    scores <- scores*fact

    fscore <- scores - floor(scores)

    # ok, lets face the problems    
    if (!all(fscore == 0)) { # rational or real scores
      fact <- ceiling(findfact(scores - min(scores) + 1, m, tol))
      scores <- round(scores * fact)
    }
  } else {
    # the user knows what to do
    scores <- round(scores * fact)
  } 
   
  if(!paired) 
    add <- min(scores) - 1
  else 
    add <- 0
  scores <- scores - add

  if(any(scores < 0)) stop("scores must be positive")

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


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

    N <- length(escores$scores)

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

    if (!is.null(B)) {
        if (B < 1) stop("B must be positive integer")
        RVAL <- .Call("sim2is", scores=as.double(escores$scores),
                                mfirst=as.integer(m), Nsim=as.integer(B), 
                                PACKAGE="exactRankTests")
        names(RVAL) <- c("T", "Prob")
        if (escores$add != 0 & paired) 
          warning("escores$add not zero for paired samples!")
        RVAL$T <- (RVAL$T + escores$add*m)/escores$fact
        return(RVAL)
    }

    if (paired) {
        # paired two sample situation
        prob <- .Call("cpermdist1", scores = as.integer(escores$scores),
                   PACKAGE="exactRankTests")
        t <- which(prob != 0)
        prob <- prob[t]
        # 0 is possible
	t <- t - 1
        if (escores$add != 0) 
          warning("escores$add not zero for paired samples!")
    } else {
        # independent samples
        col <- sum(sort(escores$scores)[(N + 1 - m):N])
        scores <- rep(1, N)
        prob <- .Call("cpermdist2", ma = as.integer(m),
                mb = as.integer(col), scorea = as.integer(scores), 
                scoreb = as.integer(escores$scores),
                retProb = as.logical(TRUE), PACKAGE="exactRankTests")
        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.17.4.1 2003/09/18 11:24:15 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")
    }

    if(!is.numeric(x)) stop("`x' must be numeric")

    MIDP <- NULL
    if (!is.null(y)) {
        if(!is.numeric(y)) stop("`y' must be numeric")
        DNAME <- paste(deparse(substitute(x)), "and",
                     deparse(substitute(y)))
    } else DNAME <- deparse(substitute(x))

    if (paired) {
        if (!is.null(y)) {
            OK <- complete.cases(x, y)
            x <- x[OK] - y[OK]
        } else {
            x <- x[complete.cases(x)]
        }
        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))
            if (alternative == "greater")
                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) {
        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)

            dummy <- rep(1,m)
            Hx <- .Call("cpermdist2", as.integer(m),
                as.integer(sum(xscores$scores)), as.integer(dummy),
                as.integer(xscores$scores), 
                as.logical(FALSE), PACKAGE="exactRankTests")

            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 {
        x <- x[complete.cases(x)]
        y <- y[complete.cases(y)]
        m <- length(x)
        n <- length(y)
        if(n < m & conf.int) {
            warning("x has more observations than y, returning perm.test(y, x, ...)")
            return(perm.test(y, x, paired, alternative,
                   mu, exact, conf.int, conf.level, tol=NULL,...))
        }
        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 <- cxy - min(cxy)
                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(c(x,y))
            wvar <- m*n/(N*(N-1))*sum((c(x,y) - wmean/m)^2)
            PVAL <- pnorm((STATISTIC - wmean)/sqrt(wvar))
            if (alternative == "greater")
                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)

            tsx <- sum(xscores$scores)
            tsy <- sum(yscores$scores)

            dummy <- rep(1,m)
            Hx <- .Call("cpermdist2", as.integer(m),
                as.integer(tsx), as.integer(dummy),
                as.integer(xscores$scores), 
                as.logical(FALSE), PACKAGE="exactRankTests")

            Hx <- matrix(Hx, nrow=m+1, byrow=TRUE)
           
            dummy <- rep(1,n)
            Hy <- .Call("cpermdist2", as.integer(n),
                as.integer(tsy), as.integer(dummy),
                as.integer(yscores$scores), 
                as.logical(FALSE), PACKAGE="exactRankTests")

            Hy <-  matrix(Hy, nrow=n+1, byrow=TRUE)

            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)
                if (any(is.na(du))) print("NA in du")
                if (length(du) < length(count)) print("du < count")

                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
        }
    }
    names(STATISTIC) <- "T"
    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.17 2003/01/16 10:44:25 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.numeric(x)) stop("`x' must be numeric")

    if(!is.null(y)) {
        if(!is.numeric(y)) stop("`y' must be numeric")
        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))
                if(alternative == "greater")
                    PVAL <- 1 - PVAL
                if(alternative == "two.sided") 
                    PVAL <- 2 * min(PVAL, 1-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))
            if (alternative == "greater")
               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
}

.onLoad <- function(lib, pkg) {
    vg190 = compareVersion(paste(version$major, 
                                 version$minor, sep="."), "1.9.0")
    
    if (vg190 < 0) {
        if(!require(ctest))
            warning("Could not load package ctest")
    } else {
        if(!require(stats))                 
            warning("Could not load package stats")      
    }
}
