#univariate, greater than alternative
PermAllOS<- function(data, ind, mvec=0, nrep = 1000)
{
  require(mvtnorm)
  require(boot)
  require(emulator)
  require(parallel)
  require(CompQuadForm)
  require(multcomp)
  DiMStuOS <- function(data, ind, mvec=0, perm = F)
  {
    require(mvtnorm)
    require(boot)
    require(emulator)
    require(parallel)
    require(CompQuadForm)
    y = data[ind,2]
    strat = data[ind,1]
    if(perm == T)
    {
      strat = data[,1]
    }
    y1 = y[strat==1]
    y2 = y[strat==2]
    n1 = length(y1)
    n2 = length(y2)
    n = n1+n2
    stu = (mean(y1) - mean(y2) - mvec)/sqrt(var(y1)/length(y1) + var(y2)/length(y2))
    sighat = sqrt(var(y1)/n1 + var(y2)/n2)*sqrt(n)
    mu3hat = ((n/n1)^3*sum((y1 - mean(y1))^3) - (n/n2)^3*sum((y2 - mean(y2))^3))/n
    edge = pnorm(stu) + ((mu3hat)/(6*sqrt(n)*sighat^3))*dnorm(stu)*(2*stu^2+1)
    c((mean(y1) - mean(y2) - mvec), (mean(y1) - mean(y2) - mvec)/sqrt(var(y1)/length(y1) + var(y2)/length(y2)), edge)
  }

  y = data[ind,2]
  strat = data[,1]
  y1 = y[strat==1]
  y2 = y[strat==2]
  n1 = length(y1)
  n2 = length(y2)
  n = n1+n2
  dim = (mean(y1) - mean(y2) - mvec)
  stu = (mean(y1) - mean(y2) - mvec)/sqrt(var(y1)/n1 + var(y2)/n2)
  sighat = sqrt(var(y1)/n1 + var(y2)/n2)*sqrt(n)
  mu3hat = ((n/n1)^3*sum((y1 - mean(y1))^3) - (n/n2)^3*sum((y2 - mean(y2))^3))/n
  edge = pnorm(stu) + ((mu3hat)/(6*sqrt(n)*sighat^3))*dnorm(stu)*(2*stu^2+1)
  datnew = cbind(strat, y)
  predim = 0
  preedge = 0
  prestu = 0
  if(nrep > 0)
  {
    aall = boot(data = datnew, statistic = DiMStuOS, R = nrep, strata = datnew[,1],  mvec = mean(y1)-mean(y2))$t
    adim = aall[,1]
    astu = aall[,2]
    aedge = aall[,3]
    predim = mean(adim <= dim)
    prestu = mean(astu <= stu)
    preedge = mean(aedge <= edge)
  }
  c(dim = dim, stu = stu, edge = edge, predim = predim, prestu = prestu, preedge = preedge)
}



#multivariate
PermAllMult <- function(data, ind, mvec=0, nrep = 1000)
{
  require(mvtnorm)
  require(boot)
  require(emulator)
  require(parallel)
  require(CompQuadForm)
  require(multcomp)
  MaxPre <- function(data, ind, mvec=0, perm = F)
  {
    require(mvtnorm)
    require(boot)
    require(emulator)
    require(parallel)
    require(CompQuadForm)
    require(multcomp)
    y = data[ind,-1]
    strat = data[ind,1]
    if(perm == T)
    {
      strat = data[,1]
    }
    y1 = y[strat==1,]
    y2 = y[strat==2,]
    n1 = nrow(y1)
    n2 = nrow(y2)
    cov1 = (cov(y1)/n1 + cov(y2)/n2)
    covpool = ((cov(y1)*(n1-1)+ cov(y2)*(n2-1))/(n1+n2-2))*(1/n1+1/n2)
    
    ses = sqrt(diag(cov1))
    covstu = cov2cor(cov1)
    dimm = (colMeans(y1) - colMeans(y2) - mvec)
    stu = abs(dimm)/ses
    mm = max(stu)
    hot <- quad.form.inv(cov1, dimm)
    pool <- quad.form.inv(covpool, dimm)
    datnew = cbind(strat, y)
    #this would use pmvnorm to prepivot, but it is computationally intensive
    #we replace it with 1000 monte-carlo simulations from the relevant Gaussian
    #premm  = as.numeric(pmvnorm(lower = -rep(mm, ncol(y1)), upper = rep(mm, ncol(y1)), sigma = covstu))
    BB = rmvnorm(1000, sigma = covstu)
    premm = mean(apply(abs(BB), 1, max) <= mm)
    vchalf =(t(chol(cov1)))
    A1 = quad.form.inv(covpool, vchalf)
    eig = eigen(A1)
    lambda = eig$values
    prepool= as.numeric((1 - davies(pool, lambda, h = rep(1, length(lambda)), delta = rep(0, nrow(cov1)), lim = 100000, acc = .000005)$Qq))
    c(hot, pool, mm, prepool, premm)
  }
  y = data[ind,-1]
  strat = data[,1]
  y1 = y[strat==1,]
  y2 = y[strat==2,]
  n1 = nrow(y1)
  n2 = nrow(y2)
  cov1 = (cov(y1)/n1 + cov(y2)/n2)
  covpool = ((cov(y1)*(n1-1)+ cov(y2)*(n2-1))/(n1+n2-2))*(1/n1+1/n2)
  
  ses = sqrt(diag(cov1))
  covstu = cov2cor(cov1)
  dimm = (colMeans(y1) - colMeans(y2) - mvec)
  stu = abs(dimm)/ses
  mm = max(stu)
  hot <- quad.form.inv(cov1, dimm)
  pool <- quad.form.inv(covpool, dimm)
  #this would use pmvnorm to prepivot, but it is computationally intensive
  #we replace it with 1000 monte-carlo simulations from the relevant Gaussian
  #premm = as.numeric(pmvnorm(lower = rep(-mm, ncol(y1)), upper = rep(mm, ncol(y1)), sigma = covstu))
  BB = rmvnorm(1000, sigma = covstu)
  premm = mean(apply(abs(BB), 1, max) <= mm)
  vchalf =(t(chol(cov1)))
  A1 = quad.form.inv(covpool, vchalf)
  eig = eigen(A1)
  lambda = eig$values
  prepool= as.numeric((1 - davies(pool, lambda, h = rep(1, length(lambda)), delta = rep(0, nrow(cov1)), lim = 100000, acc = .000005)$Qq))
  
  datnew = cbind(strat, y)
  boothot = 0
  bootpool = 0
  bootmm = 0
  bootprepool = 0
  bootpremm = 0
  if(nrep > 0)
  {
    aall = boot(data = datnew, statistic = MaxPre, R = nrep, strata = datnew[,1], mvec = colMeans(y1)-colMeans(y2))$t
    ahot= aall[,1]
    apool = aall[,2]
    amm = aall[,3]
    aprepool = aall[,4]
    apremm = aall[,5]
    boothot = mean(ahot <= hot)
    bootpool = mean(apool <= pool)
    bootmm = mean(amm <= mm)
    bootprepool = mean(aprepool <= prepool)
    bootpremm = mean(apremm <= premm)
  }
  
  
  c(hot = hot, pool = pool, mm = mm, prepool = prepool, premm = premm, boothot = boothot, bootpool = bootpool, bootmm = bootmm, bootprepool = bootprepool, bootpremm = bootpremm)
}




#ANOVA
PermAnova <- function(data, ind, mvec=rep(0, length(unique(data[,1]))), C, nrep = 1000)
{
  require(mvtnorm)
  require(boot)
  require(emulator)
  require(parallel)
  require(CompQuadForm)
  require(multcomp)
  AnovaPre <- function(data, ind, mvec=rep(0, length(unique(data[,1]))), C, perm = F)
  {
    require(mvtnorm)
    require(boot)
    require(emulator)
    require(parallel)
    require(CompQuadForm)
    require(multcomp)
    y = data[ind,-1]
    strat = data[ind,1]
    if(perm == T)
    {
      strat = data[,1]
    }
    y = y - mvec[strat]
    k = length(unique(strat))
    n = table(strat)
    N = sum(n)
    SSE = sum((n-1)*tapply(y, strat, var))
    MSE = SSE/(N-k)
    muhat = tapply(y, strat, mean)
    Gammahat = diag(MSE/n)
    SSR = sum(n*(tapply(y, strat, mean) - mean(y))^2)
    Vhat = quad.form(Gammahat, C)
    mhat = muhat%*%C
    tukey = max(abs(mhat)/diag(Vhat)^(1/2))
    fstat = SSR/(k-1)/MSE
    sigma2 = tapply(y, strat, var)
    CR = sum((n/sigma2)*(muhat - sum((n/sigma2)*muhat)/sum((n/sigma2)))^2)
    c(fstat, tukey, CR)
  }
  y = data[ind,-1]
  strat = data[,1]
  k = length(unique(strat))
  n = table(strat)
  N = sum(n)
  SSE = sum((n-1)*tapply(y, strat, var))
  MSE = SSE/(N-k)
  muhat = tapply(y, strat, mean)
  Gammahat = diag(MSE/n)
  SSR = sum(n*(tapply(y, strat, mean) - mean(y))^2)
  Vhat = quad.form(Gammahat, C)
  mhat = muhat%*%C
  tukey = max(abs(mhat)/diag(Vhat)^(1/2))
  fstat = SSR/(k-1)/MSE
  sigma2 = tapply(y, strat, var)
  CR = sum((n/sigma2)*(muhat - sum((n/sigma2)*muhat)/sum((n/sigma2)))^2)
  mvec = muhat
  bootfstat = 0
  boottukey = 0
  bootCR = 0
  datnew = cbind(strat, y)
  if(nrep > 0)
  {
    aall = boot(data = datnew, statistic = AnovaPre, R = nrep, strata = datnew[,1], mvec = mvec, C = C)$t
    afstat= aall[,1]
    atukey = aall[,2]
    aCR = aall[,3]
    bootfstat = mean(afstat<= fstat)
    boottukey  = mean(atukey <= tukey)
    bootCR = mean(aCR <= CR)
  }
  c(fstat=fstat, tukey = tukey,CR = CR, bootfstat = bootfstat, boottukey = boottukey,  bootCR = bootCR)
}



#MANOVA
PermManova <- function(data, ind, mvec=matrix(0, length(unique(data[,1]), ncol(data[,-1]))), nrep = 1000)
{
  require(mvtnorm)
  require(boot)
  require(emulator)
  require(parallel)
  require(CompQuadForm)
  require(multcomp)
  ManovaPre <- function(data, ind, mvec=matrix(0, length(unique(data[,1]), ncol(data[,-1]))), perm = F)
  {
    require(mvtnorm)
    require(boot)
    require(emulator)
    require(parallel)
    require(CompQuadForm)
    require(multcomp)
    y = data[ind,-1]
    strat = data[ind,1]
    if(perm == T)
    {
      strat = data[,1]
    }
    y = y - mvec[strat,]
    fstrat = factor(strat)
    linmod = lm(y~fstrat)
    mm = manova(y~fstrat)
    eigs = summary(mm)$Eigenvalues
    tP = sum(eigs/(1+eigs))
    tL = sum(eigs)
    tW = -prod(1/(1+eigs))
    tR = max(eigs)
    c(tP = tP, tL = tL, tW = tW, tR = tR)
  }
  y = data[ind,-1]
  strat = data[,1]
  fstrat = factor(strat)
  linmod = lm(y~fstrat)
  mm = manova(y~fstrat)
  eigs = summary(mm)$Eigenvalues
  tP = sum(eigs/(1+eigs))
  tL = sum(eigs)
  tW = -prod(1/(1+eigs))
  tR = max(eigs)
  
  dd = data.frame(fstrat, y)
  mvec = as.matrix(aggregate(y~fstrat, data = dd, mean)[,-1])
  boottP = 0
  boottW = 0
  boottL= 0
  boottR = 0
  datnew = cbind(strat, y)
  if(nrep > 0)
  {
    aall = boot(data = datnew, statistic = ManovaPre, R = nrep, strata = datnew[,1],  mvec = mvec)$t
    atP= aall[,1]
    atL = aall[,2]
    atW= aall[,3]
    atR = aall[,4]
    boottP = mean(atP <= tP)
    boottL = mean(atL <= tL)
    boottW = mean(atW <= tW)
    boottR = mean(atR <= tR)
  }
  c(tP = tP, tL = tL, tW = tW, tR = tR, boottP = boottP, boottL = boottL, boottW = boottW, boottR = boottR)
}




#Median

PermMedian <- function(data, ind, mvec=rep(0, length(unique(data[,1]))), nrep = 1000)
{
  require(mvtnorm)
  require(boot)
  require(emulator)
  require(parallel)
  require(CompQuadForm)
  require(multcomp)
  MedPre <- function(data, ind, mvec = rep(0, length(unique(data[,1]))), perm = F)
  {
    require(mvtnorm)
    require(boot)
    require(emulator)
    require(parallel)
    require(CompQuadForm)
    require(multcomp)
    y = data[ind,-1]
    strat = data[ind,1]
    if(perm == T)
    {
      strat = data[,1]
    }
    y = y - mvec[strat]
    meds = tapply(y, strat, median)
    dim = meds[1] - meds[2]
    
    y1s = sort(y[strat==1])
    y2s = sort(y[strat==2])
    
    n1 = length(y1s)
    n2 = length(y2s)
    l1 = 1:n1
    l2 = 1:n2
    
    p1 = pbinom((n1-1)/2, n1, (l1-1)/n1) - pbinom((n1-1)/2, n1, l1/n1)
    p2 = pbinom((n2-1)/2, n2, (l2-1)/n2) - pbinom((n2-1)/2, n2, l2/n2)
    
    v1 = n1*sum((y1s - meds[1])^2*p1)
    v2 = n2*sum((y2s-meds[2])^2*p2)
    stu = dim/sqrt(v1/n1 +v2/n2)
    c(dim = dim, stu = stu)
  }
  y = data[ind,-1]
  strat = data[,1]
  meds = tapply(y, strat, median)
  dim = meds[1] - meds[2]
  y1s = sort(y[strat==1])
  y2s = sort(y[strat==2])
  
  n1 = length(y1s)
  n2 = length(y2s)
  l1 = 1:n1
  l2 = 1:n2
  
  p1 = pbinom((n1-1)/2, n1, (l1-1)/n1) - pbinom((n1-1)/2, n1, l1/n1)
  p2 = pbinom((n2-1)/2, n2, (l2-1)/n2) - pbinom((n2-1)/2, n2, l2/n2)
  
  v1 = n1*sum((y1s - meds[1])^2*p1)
  v2 = n2*sum((y2s-meds[2])^2*p2)
  stu = dim/sqrt(v1/n1 +v2/n2)
  mvec = meds
  datnew = cbind(strat, y)
  bootdim = 0
  bootstu = 0
  if(nrep > 0)
  {
    aall = boot(data = datnew, statistic = MedPre, R = nrep, strata = datnew[,1],  mvec = mvec)$t
    adim= aall[,1]
    astu = aall[,2]
    bootdim = mean(adim<= dim)
    bootstu  = mean(astu <= stu)
  }
  c(dim = dim, stu = stu, bootdim = bootdim, bootstu = bootstu)
}