"simulation1full" <- 
function(nreps = 100)
{
	kvals <- c(5, 50, 500)
	muvals <- c(3, 4, 5, 7)
	results <- array(NA, dim = c(18, 4, 3, nreps), dimnames = list(
		c("exponential", "cauchy", "postmean", "exphard", "a=1",
		"a=0.5", "a=0.2", "a=0.1", "SURE", "adapt", 
		"FDR q=0.01", "FDR q=0.1", "FDR q=0.4", "blockthresh",
		"neighblock", "neighcoeff", "universal soft", 
		"universal hard"), as.character(muvals), as.character(
		kvals), NULL))
	#  simulate matrix of errors
	#
	set.seed(153)
	errmat <- matrix(rnorm(nreps * 1000), nrow = 1000)
	jperm <- sample(1:1000)
	#
	#  now conduct individual simulations
	#
	# if you want to calculate errors in a different norm, use the perrnorm parameter
	#  in simcomp1
	#
	for(jk in (1:3)) {
		for(jmu in (1:4)) {
			zmean <- rep(0, 1000)
			zmean[1:kvals[jk]] <- muvals[jmu]
			z <- matrix(rep(zmean, nreps), nrow = 1000) +
				errmat
			results[1, jmu, jk,  ] <- simcomp1(
				"ebayesthresh", zmean, z, a = NA)
			results[2, jmu, jk,  ] <- simcomp1(
				"ebayesthresh", zmean, z, prior = 
				"cauchy")
			results[3, jmu, jk,  ] <- simcomp1(
				"ebayesthresh", zmean, z, threshrule = 
				"mean", a = NA)
			results[4, jmu, jk,  ] <- simcomp1(
				"ebayesthresh", zmean, z, threshrule = 
				"hard", a = NA)
			results[5, jmu, jk,  ] <- simcomp1(
				"ebayesthresh", zmean, z, a = 1)
			results[6, jmu, jk,  ] <- simcomp1(
				"ebayesthresh", zmean, z, a = 0.5)
			results[7, jmu, jk,  ] <- simcomp1(
				"ebayesthresh", zmean, z, a = 
				0.20000000000000001)
			results[8, jmu, jk,  ] <- simcomp1(
				"ebayesthresh", zmean, z, a = 
				0.10000000000000001)
			results[9, jmu, jk,  ] <- simcomp1("suresmooth",
				zmean, z, adapt = F)
			results[10, jmu, jk,  ] <- simcomp1(
				"suresmooth", zmean, z, adapt = T)
			results[11, jmu, jk,  ] <- simcomp1("fdrsmooth",
				zmean, z, q = 0.01)
			results[12, jmu, jk,  ] <- simcomp1("fdrsmooth",
				zmean, z, q = 0.10000000000000001)
			results[13, jmu, jk,  ] <- simcomp1("fdrsmooth",
				zmean, z, q = 0.40000000000000002)
			results[14, jmu, jk,  ] <- simcomp1(
				"blockthresh.fun", zmean[jperm], z[
				jperm,  ])
			results[15, jmu, jk,  ] <- simcomp1(
				"neighblock.fun", zmean[jperm], z[
				jperm,  ])
			results[16, jmu, jk,  ] <- simcomp1(
				"neighcoeff.fun", zmean[jperm], z[
				jperm,  ])
			results[17, jmu, jk,  ] <- simcomp1(
				"universalsmooth", zmean, z, soft = T)
			results[18, jmu, jk,  ] <- simcomp1(
				"universalsmooth", zmean, z, soft = F)
		}
	}
	return(results)
}

"simcomp1" <- 
function(method, zmean, z, fullcomp = T, perrnorm = 2, ...)
{
	nreps <- dim(z)[2]
	zs <- apply(z, 2, method, ...)
	mse <- apply(zs - zmean, 2, vecnorm, p = perrnorm)^perrnorm
	cat(method, "\n")
	return(mse)
}

"fdrsmooth" <- 
function(x, q = 0.050000000000000003)
{
	#  implement the FDR procedure
	#
	n <- length(x)
	y <- rev(sort(abs(x)))
	t <-  - qnorm(((1:n) * q)/(2 * n))
	tfdr <- min(t[y >= t], t[1])
	xfdr <- x
	xfdr[abs(x) <= tfdr] <- 0
	return(xfdr)
}

"blockthresh.fun" <- 
function(y, noise.level = 1)
{
	if(noise.level == 0)
		noise.level <- mad(diff(y))/sqrt(2)
	m <- length(y)
	L <- floor(logb(m))
	nb <- floor(m/L)
	thresh <- 4.5052409999999998 * L * noise.level^2
	SS <- rep(0, nb)
	for(b in (1:nb)) {
		SS[b] <- sum(y[((b - 1) * L + 1):(b * L)]^2)
	}
	factor <- c(SS > thresh)
	factor <- rep(factor, rep(L, nb))
	y[1:(nb * L)] <- y[1:(nb * L)] * factor
	if(m > nb * L) {
		Lb <- m - nb * L
		SSL <- sum(y[(nb * L + 1):m]^2)
		y[(nb * L + 1):m] <- y[(nb * L + 1):m] * c(SSL > thresh
			)
	}
	return(y)
}

"neighblock.fun" <- 
function(y, noise.level = 1)
{
	m <- length(y)
	L0 <- floor(log(m)/2)
	nb <- floor(m/L0)
	L1 <- max(1, floor(L0/2))
	L <- L0 + 2 * L1
	thresh <- 4.5049999999999999 * L * noise.level^2
	yy <- c(y[(m - L1 + 1):m], y, y[1:L1])
	SS <- rep(0, nb)
	for(b in 1:nb) {
		SS[b] <- sum(yy[((b - 1) * L0 + 1):((b - 1) * L0 + L)]^
			2)
	}
	factor <- (SS > thresh)
	factor <- rep(factor, rep(L0, nb))
	y[1:(nb * L0)] <- y[1:(nb * L0)] * factor
	if(m > nb * L0) {
		Lb <- floor((L - (m - nb * L0))/2)
		SSL <- sum(c(y[(m - L + Lb + 1):m], y[1:Lb])^2)
		y[(nb * L0 + 1):m] <- y[(nb * L0 + 1):m] * (SSL > 
			thresh)
	}
	return(y)
}

"neighcoeff.fun" <- 
function(y, noise.level = 1)
{
	m <- length(y)
	if(noise.level == 0)
		noise.level <- mad(diff(y))/sqrt(2)
	thresh <- 2 * logb(m) * noise.level^2
	yy <- c(y[m], y, y[1])
	SS <- rep(0, m)
	for(b in 1:m)
		SS[b] <- sum((yy[b:(b + 2)])^2)
	y[SS <= thresh] <- 0
	return(y)
}

"suresmooth" <- 
function(x, adapt = T, returnthresh = F)
{
	#
	#  find the SURE threshold and apply it to the sequence x
	#  also find carry out the ADAPT smoothing if adapt=T
	#
	n <- length(x)
	#
	#  first do the adapt test 
	#
	if(adapt & (sum(x^2)/n <= 1 + (log(n)/log(2))^1.5/sqrt(n))) 
			thresh <- sqrt(2 * log(n)) else {
		#
		#  the sure calculation
		#
		sx <- sort(x^2)
		sureseq <- n - 2 * (1:n) + seq(from = n - 1, to = 0,
			by = -1) * sx + cumsum(sx)
		thresh <- sqrt(min(sx[sureseq == min(sureseq)], 2 *
			log(n)))
	}
	xhat <- sign(x) * pmax(0, abs(x) - thresh)
	if(returnthresh)
		return(thresh)
	else return(xhat)
}

"universalsmooth" <- 
function(x, soft = T)
{
	n <- length(x)
	tt <- sqrt(2 * log(n))
	sx <- sign(x)
	ax <- abs(x)
	if(soft)
		xt <- sx * pmax(0, ax - tt)
	else {
		xt <- x
		xt[ax <= tt] <- 0
	}
	return(xt)
}

"table1ebayess" <- 
function(zs = ebayess.res, tex = F, avcal = T, relcal = F)
{
	#   takes results of simulation2() or simulation3() and
	#   produces Table 2.   Can either produce
	#   an array or else latex output
	#   If avcal, calculates means, if avcal=F calculates std errors
	#   If relcal=T does things relative to results for first method considered
	#
	if(relcal) {
		for(jj in (2:dim(zs)[1]))
			zs[jj,  ,  ,  ] <- zs[jj,  ,  ,  ] - zs[1,
				,  ,  ]
		zs[1,  ,  ,  ] <- 0
	}
	if(avcal) {
		zr <- apply(zs, c(1, 2, 3), mean)
	}
	else {
		zr <- sqrt(apply(zs, c(1, 2, 3), var)/dim(zs)[4])
	}
	if(!tex) {
		return(zr)
	}
	zt <- round(cbind(zr[,  , 1], zr[,  , 2]))
	zt <- apply(zt, 1, paste, collapse = " & ")
	zt <- paste(names(zt), zt, sep = " & ", collapse = " \\  ")
	return(zt)
}

"table2ebayess" <- 
function(zs = ebayess.res, tex = F)
{
	#   takes results and
	#   produces inefficiency plot.   Can either produce
	#   an array or else latex output
	#
	nmethods <- dim(zs)[1]
	methnames <- dimnames(zs)[[1]]
	zav <- apply(zs, c(1, 2, 3), mean)
	zav <- t(matrix(zav, nrow = nmethods, dimnames = list(methnames,
		NULL)))
	zbest <- apply(zav, 1, min)
	zineff <- 100 * (zav/zbest - 1)
	zineffmat <- matrix(NA, nrow = nmethods, ncol = 4, dimnames = 
		list(methnames, c("median", "mean", "10th", "max")))
	zineffmat[, 1] <- apply(zineff, 2, median)
	zineffmat[, 2] <- apply(zineff, 2, mean)
	zineffmat[, 3] <- apply(zineff, 2, quantile, probs = 9/11)
	zineffmat[, 4] <- apply(zineff, 2, max)
	if(!tex) {
		return(zineffmat)
	}
	zt <- round(zineffmat)
	zt <- apply(zt, 1, paste, collapse = " & ")
	zt <- paste(names(zt), zt, sep = " & ", collapse = " \\  ")
	return(zt)
}