library(ggplot2)
library(perm)

testWith <- function(n, mu, tst)
{
        pvals <- replicate(400, tst(rnorm(n, 0, 1), rnorm(n, mu, 1))$p.value)
        data.frame(mean(pvals), mean(pvals < 0.05))
}

tbl <- data.frame()
for(mu in seq(-2, 2, 0.02)) for(n in seq(10, 100, 5)) { tbl <- rbind(tbl, data.frame(mu=mu, n=n)) }
wres <- cbind(tbl, do.call(rbind, apply(tbl, c(1), function(row) testWith(row[2], row[1], wilcox.test))))
tres <- cbind(tbl, do.call(rbind, apply(tbl, c(1), function(row) testWith(row[2], row[1], permTS))))

diff <- tbl
diff$pvals <- tres$mean.pvals. - wres$mean.pvals.
diff$power <- tres$mean.pvals...0.05. - wres$mean.pvals...0.05.

svg("tres_pvals.svg")
ggplot(tres, aes(x=n, y=mu)) + geom_tile(aes(fill=mean.pvals.)) + scale_fill_gradientn(colours = rainbow(7))
dev.off()
svg("tres_power.svg")
ggplot(tres, aes(x=n, y=mu)) + geom_tile(aes(fill=mean.pvals...0.05.)) + scale_fill_gradientn(colours = rainbow(7))
dev.off()

svg("wres_pvals.svg")
ggplot(wres, aes(x=n, y=mu)) + geom_tile(aes(fill=mean.pvals.)) + scale_fill_gradientn(colours = rainbow(7))
dev.off()
svg("wres_power.svg")
ggplot(wres, aes(x=n, y=mu)) + geom_tile(aes(fill=mean.pvals...0.05.)) + scale_fill_gradientn(colours = rainbow(7))
dev.off()

svg("diff_pvals.svg")
ggplot(diff, aes(x=n, y=mu)) + geom_tile(aes(fill=pvals)) + scale_fill_gradientn(colours = rainbow(7))
dev.off()
svg("diff_power.svg")
ggplot(diff, aes(x=n, y=mu)) + geom_tile(aes(fill=power)) + scale_fill_gradientn(colours = rainbow(7))
dev.off()

print(sprintf("diff pvals: max %f, avg %f; diff power: max %f, avg %f", max(diff$pvals), mean(diff$pvals), max(diff$power), mean(diff$power)))