testGaussien <- function(x,y,a){
  m <- length(x); n <- length(y);
  sqrt(m*n/(m+n))*(mean(x)-mean(y)) > qnorm(1-a)
}

testMW <- function(x,y,a){
  te <- wilcox.test(x,y, paired=FALSE, alternative="greater");
  te$p.value<a
}

estimPuissance <- function(m,n,test,a){
  N <- 1000;
  mu1 <- 1; mu2 <- 0; sig2 <- 1; alpha <- 0.025;
  res <- array(T, N)
  for (j in 1:N){
    x <- mu1 + sqrt(sig2)*rnorm(m);
    y <- mu2 + sqrt(sig2)*rnorm(n);
    res[j] <- test(x,y,a);
  }
  mean(res);
}

puissanceOfSize <- function(mMax, nMax, test, a){
  res <- matrix(0, ncol=mMax, nrow=nMax);
  for (m in 1:mMax){
    print(m/mMax)
    for (n in 1:nMax){
      res[n, m] <- estimPuissance(m, n, testMW, 0.025)
    }
  }
  res
}

plotPuissance <- function(P){
  n <- dim(P)[1]; m <- dim(P)[2];
  lx <- c(); ly <- c(); lp <- c();
  for (j in 1:n){
    for (k in 1:m){
      lx <- c(lx, j);
      ly <- c(ly, k);
      lp <- c(lp,  P[j,k]);
    }
  }
  puis <- 0.8
  plot(lx[lp>puis], ly[lp>puis], col="red", xlim=c(0,n+1), ylim = c(0, m+1));
  points(lx[lp<=puis], ly[lp<=puis], col="blue");
}

m <- 50; n <- 50; a <- 0.025
# PG <- puissanceOfSize(m, n, testGaussien, a);
# write.table(as.table(PG), file="PuissanceTestGaussien.html")
PG <- as.matrix(read.table("PuissanceTestGaussien.html"))

# PMW <- puissanceOfSize(m, n, testMW, a);
# write.table(as.table(PMW), file="PuissanceTestMW.html")
PMW <- as.matrix(read.table("PuissanceTestMW.html"))

par(mfrow=c(1,2));
plotPuissance(PG)
plotPuissance(PMW)
Im <- seq(1,m,0.1)
# lines(Im, 4*Im/(Im - 4), col="green", lwd=3)