threshhyb<-function(ywd, lowlev = 2, verbose = F, seed = 0)
{
  ans <- ywd
  n <- length(ywd$D)
  nlev <- log(n + 1, base = 2) - 1
  i <- nlev
  iloc <- 1
  while(i > lowlev) {
#       Extract current level coefficients from all wavelet coefficients
    raw <- ywd$D[iloc:(iloc + 2^i - 1)]
    d <- length(raw)	
#       Test:  if the variance is small enough, just use threshold sqrt(2logd)
    if((sum(raw^2) - d)/d <= sqrt(i^3/2^i)) {
      if(verbose)
	cat(paste("At level ", i, " the threshhold is sqrt(2log(d)): ", 
		  sqrt(2 * log(d)), "\n", sep = ""))
      ans$D[iloc:(iloc + 2^i - 1)] <- shrinkit(ywd$D[iloc:(iloc + 2^i - 1)], 
					       sqrt(2 * log(d)))
    }
      else {
#       Generate random subset
	if(length(seed) != 1) .Random.seed <- seed
	Iset <- sort(sample(d, d/2))
	rawI <- raw[Iset]
	rawIp <- raw[ - Iset]
	ggI <- sort(abs(rawI))
	ggIp <- sort(abs(rawIp))	
#       Calculate SURE for all possible thresholds
	surevecI <- sapply(c(ggI[ggI < sqrt(2 * log(d))], 0,
			     sqrt(2 * log(d))), sure, ggI)
	surevecIp <- sapply(c(ggIp[ggI < sqrt(2 * log(d))], 0, 
			      sqrt(2 * log(d))), sure, ggIp)	
#       Threshold that minimizes risk
	llI <- length(surevecI)
	llIp <- length(surevecIp)	
#       The minimum occurs either at sqrt(2logd), 
	if(min(surevecI) == surevecI[llI])
	  threshI <- sqrt(2 * log(d))
	else if(min(surevecI) == surevecI[llI - 1])
	  threshI <- 0
	else threshI <- ggI[match(min(surevecI), surevecI)]	
	#       or at 0, 
	if(min(surevecIp) == surevecIp[llIp])
	  threshIp <- sqrt(2 * log(d))
	else if(min(surevecIp) == surevecI[llIp - 1])
	  threshIp <- 0
	else
	  threshIp <- ggIp[match(min(surevecIp), surevecIp)]
		#       or at 0, 
	if(verbose) {
	  cat(paste("At level ", i, ", threshold1 is ", threshI, "\n",
		    sep = ""))
	  cat(paste("At level ", i, ", threshold2 is ", threshIp,
		    "\n", sep = ""))
	}
#       Perform shrinking
	newI <- shrinkit(rawI, threshIp)
	newIp <- shrinkit(rawIp, threshI)
	new <- rep(0, d)
	new[Iset] <- newI
	new[ - Iset] <- newIp
	ans$D[iloc:(iloc + 2^i - 1)] <- new
      }
#       Otherwise, go through all this stuff
    iloc <- iloc + 2^i
    i <- i - 1
  }
  ans
}

shrinkit <- function(coeffs, thresh)
{
  sign(coeffs) * pmax(abs(coeffs) - thresh, 0)
}

sure <- function(t, x)
{
  ax <- sort(abs(x))	
					#  num gives the number of elements
					# in x that are <= t
  num <- match(F, ax <= t, nomatch = length(ax) + 1) - 1
  length(ax) - 2 * num + sum(pmin(ax, t)^2)
}