R/waverify.R

Defines functions detailer energizer mowaverify2d.default mowaverify2d.SpatialVx mowaverify2d plot.waverify2d print.waverify2d waverify2d.default waverify2d.SpatialVx waverify2d makeWaveNames makeEmptyWave

Documented in detailer energizer makeEmptyWave makeWaveNames mowaverify2d mowaverify2d.default mowaverify2d.SpatialVx plot.waverify2d print.waverify2d waverify2d waverify2d.default waverify2d.SpatialVx

makeEmptyWave <- function(x) {
   n <- length(x)
   for(i in 1:n) x[[i]] <- x[[i]]*0
   return(x)
} # end of 'makeEmptyWave' function.

makeWaveNames <- function(J) {
   out <- list()
   out$HH <- paste( "HH", 1:J, sep="")
   out$HL <- paste( "HL", 1:J, sep="")
   out$LH <- paste( "LH", 1:J, sep="")
   out$LL <- paste( "LL", J, sep="")
   return(out)
} # end of 'makeWaveNames' function.

waverify2d <- function(X, ..., Clim=NULL, wavelet.type="haar", J=NULL) {
    UseMethod("waverify2d", X)
} # end of 'waverify2d' function.

waverify2d.SpatialVx <- function(X, ..., Clim=NULL, wavelet.type="haar", J=NULL,
    useLL=FALSE, compute.shannon=FALSE, which.space="field", time.point=1, obs = 1, model=1, verbose=FALSE) {

    out <- list()
    a <- attributes(X)
    attributes(out) <- a

    ## Begin: Get the data sets
    dat <- datagrabber(X, time.point=time.point, obs = obs, model=model)

    x <- dat$X
    y <- dat$Xhat
    ## End: Get the data sets

    out2 <- waverify2d.default(X=x, Y=y, Clim=Clim, wavelet.type=wavelet.type, J=J,
			        useLL=useLL, compute.shannon=compute.shannon,
				which.space=which.space, verbose=verbose)

    out$A <- out2$A
    out$J <- out2$J
    out$X.wave <- out2$X.wave
    out$Y.wave <- out2$Y.wave
    out$Clim.wave <- out2$Clim.wave
    out$Shannon.entropy <- out2$Shannon.entropy
    out$energy <- out2$energy
    out$mse <- out2$mse
    out$rmse <- out2$rmse
    out$acc <- out2$acc
    out$wave.method <- out2$wave.method

    attr(out, "time.point") <- time.point
    attr( out, "obs" ) <- obs
    attr(out, "model") <- model
    class(out) <- "waverify2d"
    return(out)
} # end of 'waverify2d.SpatialVx' function.

waverify2d.default <- function(X, ..., Y, Clim=NULL,
			wavelet.type="haar", J=NULL, useLL=FALSE,
			compute.shannon=FALSE, which.space="field",
			verbose=FALSE) {
   ##
   ## Function to apply wavelet decomposition to each field of a verification
   ## set, V=(Y, X), containing the forecast and observation/analysis fields,
   ## respectively.  Various verification metrics are applied to each wave number.
   ##
   ## 'X', 'Y', 'Clim'  an observed and forecast field, and optional climate field (to calculate the ACC).
   ## 'wavelet.type' character name giving the wavelet decomposition function to be used
   ##	(see help file for 'dwt.2d' from 'waveslim' package).
   ## 'J' numeric giving the number of levels of the wavelet to use.  If NULL, uses log2(dim( Y)[1]).
   ## 'useLL' logical, if TRUE the average component "LL" is used for the inverse wavelet transform, which ensures positive energy.
   ## 'compute.shannon' logical whether or not to compute the Shannon entropy.
   ## 'verbose' logical, should progress information be printed to the screen?
   ##
   ## Details: A diadic wavelet (in future, may also allow nondiadic fields) decomposition is applied to both forecast and observed
   ##	fields.  The rmse is taken for each scale of the decomposition.  If a climatology field is also supplied,
   ##	then the ACC is also take for each scale.
   ##
   ## Value: a list object containing the components:
   ##
   ## 'Y.wave', 'X.wave' objects of class "dwt.2d" or "modwt.2d" (if fields are not dyadic) giving the
   ##	wavelet decomposition coefficients for each field.
   ## 'J' the number of levels used for the wavelet decomposition.
   ## 'Shannon.entropy' 2X1 matrix (if compute.shannon is TRUE) giving the Shannon entropy score for the chosen wavelet function applied to each image.
   ## 
   ## Depends: pacakges 'waveslim' 
   ##
   if( verbose)  begin.time <- Sys.time()
   out <- list()
    attributes(out) <- atmp <- list(...)
    if(is.null(atmp$loc.byrow)) attr(out, "loc.byrow") <- FALSE

   climate <- !is.null( Clim)
   
   N <- dim( Y)[1]
   
   ## Do the 2-d wavelet transforms.
   if( is.null( J)) J <- log2(N)
   out$J <- J
   
   wv.X <- dwt.2d( X, wf=wavelet.type, J=J)
   wv.Y <- dwt.2d( Y, wf=wavelet.type, J=J)
   if( climate) wv.Clim <- dwt.2d( Clim, wf=wavelet.type, J=J)
   out$Y.wave <- wv.Y
   out$X.wave <- wv.X
   if( climate) out$Clim.wave <- wv.Clim
   
   if(compute.shannon) {
      Shannon.entropy <- matrix( NA, nrow=1, ncol=2)
      colnames( Shannon.entropy) <- c("Verification", "Forecast")
      tmp <- abs(unlist(wv.X))
      tmp <- tmp/sum(tmp,na.rm=TRUE)
      id0 <- tmp==0
      tmp[!id0] <- tmp[!id0]*log(tmp[!id0])
      Shannon.entropy[,1] <- -sum(tmp,na.rm=TRUE)
      tmp <- abs(unlist(wv.Y))
      tmp <- tmp/sum(tmp,na.rm=TRUE)
      id0 <- tmp==0
      tmp[!id0] <- tmp[!id0]*log(tmp[!id0])
      Shannon.entropy[,2] <- -sum(tmp, na.rm=TRUE)
      out$Shannon.entropy <- Shannon.entropy
   }
   
   # List names for coefficients of each wavelet object.
   tmp <- makeWaveNames(J=J)
   HH <- tmp$HH
   HL <- tmp$HL
   LH <- tmp$LH
   LL <- tmp$LL

   mse.i <- rmse.i <- numeric(J) + NA
   if( climate) acc.i <- rmse.i
   
   if( !climate) {
      energy.i <- matrix( NA, nrow=J, ncol=2)
      colnames( energy.i) <- c("Verification", "Forecast")
   } else {
      energy.i <- matrix( NA, nrow=J, ncol=3)
      colnames( energy.i) <- c("Verification", "Forecast", "Climatology")
   }
  
   if(verbose) cat("\n", "calculating statistics at each level.\n")
   for( i in 1:J) {
     if( verbose) cat(i, " ")
     energy.i[i,1] <- sum(wv.X[[HH[i]]]^2, na.rm=TRUE) + sum(wv.X[[HL[i]]]^2, na.rm=TRUE) + sum(wv.X[[LH[i]]]^2, na.rm=TRUE)
     energy.i[i,2] <- sum(wv.Y[[HH[i]]]^2, na.rm=TRUE) + sum(wv.Y[[HL[i]]]^2, na.rm=TRUE) + sum(wv.Y[[LH[i]]]^2, na.rm=TRUE)
     if( climate) energy.i[i,3] <- sum(wv.Clim[[HH[i]]]^2, na.rm=TRUE) + sum(wv.Clim[[HL[i]]]^2, na.rm=TRUE) + sum(wv.Clim[[LH[i]]]^2, na.rm=TRUE)
   
      if(verbose) {
	   if(which.space=="field") cat("Finding inverse transforms for this scale (i.e., detail reconstruction).\n")
	   else cat("Finding coefficient detail for this scale.\n")
      }
      inv.hold0.i <- detailer(wv.X, level=i, which.space=which.space, useLL=useLL, lnames=tmp)
      inv.hold1.i <- detailer(wv.Y, level=i, which.space=which.space, useLL=useLL, lnames=tmp)
      if( climate) inv.hold2.i <- detailer(wv.Clim, level=i, which.space=which.space, useLL=useLL, lnames=tmp)
      if(verbose) cat("Inverse transforms for present scale found.  Calculating vx stats.\n")
      mse.i[i] <- mean( (c( inv.hold0.i) - c( inv.hold1.i))^2, na.rm=TRUE)
      rmse.i[i] <- sqrt( mse.i[i])
   
      if( climate) {
   	if(verbose) cat("Calculating ACC.\n")
   	denom <- sqrt( sum( (c( inv.hold0.i) - c( inv.hold2.i))^2, na.rm=TRUE))*sqrt( sum( (c( inv.hold1.i) - c( inv.hold2.i))^2, na.rm=TRUE))
   	numer <- sum( diag( t( inv.hold0.i - inv.hold2.i) %*% (inv.hold1.i - inv.hold2.i)), na.rm=TRUE)
   	acc.i[i] <- numer/denom
   	if(verbose) cat("ACC calculated.\n")
      }
   } # end of for 'i' loop.
   
   out$energy <- energy.i
   out$mse <- mse.i
   out$rmse <- rmse.i
   if( climate) out$acc <- acc.i
   out$type <- wavelet.type
   out$wave.method <- "DWT"
   if(verbose) print( Sys.time() - begin.time)
   class(out) <- "waverify2d"
   return(out)
} # end of 'waverify2d.default' function.

print.waverify2d <- function(x, ...) {
    a <- attributes(x)
    if(!is.null(a$data.name)) {
	cat(a$data.name, "\n")
	cat(a$msg)
    }

    cat("\n", "Wavelet method = ", x$wave.method)
    cat("\n", " number of levels (J) = ", x$J, "\n")
    cat("Energy:\n")
    print(x$energy)
    cat("\n", "MSE:\n")
    print(x$mse)
    cat("\n", "RMSE:\n")
    print(x$rmse)

    invisible()
} # end of 'print.waverify2d' function.

plot.waverify2d <- function(x, ..., main1="X", main2="Y", main3="Climate",
    which.plots=c("all", "dwt2d", "details", "energy", "mse", "rmse", "acc"),
    separate=FALSE, col, horizontal=TRUE) {

    if(missing(col)) col <- c("gray", tim.colors(64))

    which.plots <- tolower(which.plots)
    which.plots <- match.arg(which.plots)

    a <- attributes(x)
    loc.byrow <- a$loc.byrow

    if(!is.null(a$projection) && any(is.element(c("all","details"), which.plots))) {
        proj <- a$projection
	domap <- a$map
	loc <- a$loc
	xd <- a$xdim
	if(proj) loc2 <- list(x=matrix(loc[,1], xd[1], xd[2], byrow=loc.byrow), y=matrix(loc[,2], xd[1], xd[2], byrow=loc.byrow))
    } else proj <- domap <- FALSE

    if(domap) {
	locr <- apply(loc, 2, range, finite=TRUE)
	ax <- list(x=pretty(round(loc[,1], digits=2)), y=pretty(round(loc[,2], digits=2)))
    }

    if(!is.null(a$data.name)) {
	if(length(a$data.name) == a$nforecast + 2) dn <- a$data.name[-1]
	else dn <- a$data.name
	if(is.numeric(a$model)) model.num <- a$model
	else model.num <- (1:a$nforecast)[a$model == dn[-1]]
	dn <- c(dn[1], dn[model.num + 1])
    }

    if(missing(main1) && !is.null(a$data.name)) main1 <- dn[1]
    if(missing(main2) && !is.null(a$data.name)) main2 <- dn[2]


    if(!is.null(a$msg)) par(oma=c(0,0,2,0))

    X <- x$X.wave
    Y <- x$Y.wave
    if(is.null(X) | is.null(Y)) do.detail <- FALSE
    else do.detail <- TRUE
    if(climate <- !is.null(x$Clim.wave)) Clim <- x$Clim.wave

    if(any(is.element(c("all","dwt2d2"),which.plots)) & (x$wave.method == "DWT")) {
	if(!do.detail) stop("plot.waverify2d: sorry, no detail information provided for one or both of the fields.")
	if(!separate) {
	   if(climate) nc <- 3
	   else nc <- 2
	   par(mfrow=c(1,nc), mar=c(4,4,4,4))
	}
	plot(X, main=main1, ...)
	plot(Y, main=main2, ...)
	if(climate) plot(Clim, main=main3, ...)
    }
    J <- x$J
    if(any(is.element(c("all", "details"), which.plots))) {
      if(!do.detail) stop("plot.waverify2d: sorry, no detail information provided for one or both of the fields.")
      tmp <- makeWaveNames(J=J)
      HH <- tmp$HH
      HL <- tmp$HL
      LH <- tmp$LH
      LL <- tmp$LL
      if(!separate) {
         nr <- ceiling(J/2)
         if(climate) nc <- 6
         else nc <- 4
         par(mfrow=c(nr,nc), mar=c(1,2,4,1))
      }
      holdX <- makeEmptyWave(X)
      holdY <- holdX
      if(climate) holdClim <- holdY
      for(i in 1:J) {
	ni <- length(c(holdX[[HH[i]]]))*3
	holdX[[HH[i]]] <- X[[HH[i]]]
	holdX[[HL[i]]] <- X[[HL[i]]]
	holdX[[LH[i]]] <- X[[LH[i]]]
	holdY[[HH[i]]] <- Y[[HH[i]]]
        holdY[[HL[i]]] <- Y[[HL[i]]]
        holdY[[LH[i]]] <- Y[[LH[i]]]
	if(climate) {
	   holdClim[[HH[i]]] <- Clim[[HH[i]]]
           holdClim[[HL[i]]] <- Clim[[HL[i]]]
           holdClim[[LH[i]]] <- Clim[[LH[i]]]
	}

	lamX <- sqrt( sum( c( holdX[[HH[i]]]^2) + c(holdX[[HL[i]]]^2) + c(holdX[[LH[i]]]^2), na.rm=TRUE)/ni)
	holdX[[HH[i]]] <- sign(X[[HH[i]]])*pmax(c(abs(X[[HH[i]]]) - lamX), 0)
        holdX[[HL[i]]] <- sign(X[[HL[i]]])*pmax(c(abs(X[[HL[i]]]) - lamX), 0)
        holdX[[LH[i]]] <- sign(X[[LH[i]]])*pmax(c(abs(X[[LH[i]]]) - lamX), 0)

	lamY <- sqrt(sum(c(holdY[[HH[i]]]^2) + c(holdY[[HL[i]]]^2) + c(holdY[[LH[i]]]^2), na.rm=TRUE)/ni)
        holdY[[HH[i]]] <- sign(Y[[HH[i]]])*pmax(c(abs(Y[[HH[i]]]) - lamY), 0)
        holdY[[HL[i]]] <- sign(Y[[HL[i]]])*pmax(c(abs(Y[[HL[i]]]) - lamY), 0)
        holdY[[LH[i]]] <- sign(Y[[LH[i]]])*pmax(c(abs(Y[[LH[i]]]) - lamY), 0)

	if(climate) {
	   lamClim <- sqrt( sum( c( holdClim[[HH[i]]]^2) + c(holdClim[[HL[i]]]^2) + c(holdClim[[LH[i]]]^2), na.rm=TRUE)/ni)
           holdClim[[HH[i]]] <- sign(Clim[[HH[i]]])*pmax(c(abs(Clim[[HH[i]]]) - lamClim), 0)
           holdClim[[HL[i]]] <- sign(Clim[[HL[i]]])*pmax(c(abs(Clim[[HL[i]]]) - lamClim), 0)
           holdClim[[LH[i]]] <- sign(Clim[[LH[i]]])*pmax(c(abs(Clim[[LH[i]]]) - lamClim), 0)
	}
	if(x$wave.method=="DWT") {
	   invX <- idwt.2d(holdX)
	   invY <- idwt.2d(holdY)
	   if(climate) invClim <- idwt.2d(holdClim)
	} else {
	   invX <- imodwt.2d(holdX)
	   invY <- imodwt.2d(holdY)
	   if(climate) invClim <- imodwt.2d(holdClim)
	}

	if(!is.null(x$loc)) loc <- x$loc
	else {
	   xdim <- dim(invX)
	   loc <- cbind(rep(1:xdim[1],xdim[2]),rep(1:xdim[2],each=xdim[1]))
	}
	if(!climate) zl <- range(c(c(invX),c(invY)),na.rm=TRUE)
	else zl <- range(c(c(invX),c(invY),c(invClim)),na.rm=TRUE)

	if(domap) {
	    if(proj) {
		map(xlim=locr[,1], ylim=locr[,2], type="n")
		# axis(1, at=ax$x, labels=ax$x)
		# axis(2, at=ax$y, labels=ax$y)
		poly.image(loc2$x, loc2$y, invX, main=paste(main1, " (level ", i, ")", sep=""),
				zlim=zl, col=col, axes=FALSE, add=TRUE, ...)
		map(add=TRUE, lwd=2)
		map(add=TRUE, database="state")
		map.axes()

		map(xlim=locr[,1], ylim=locr[,2], type="n")
		# axis(1, at=ax$x, labels=ax$x)
		# axis(2, at=ax$y, labels=ax$y)
		poly.image(loc2$x, loc2$y, invY, main=paste(main2, " (level ", i, ")", sep=""), zlim=zl, col=col, axes=FALSE, add=TRUE, ...)
		map(add=TRUE, lwd=2)
                map(add=TRUE, database="state")
		map.axes()

		if(climate) {
                    map(xlim=locr[,1], ylim=locr[,2], type="n")
		    # axis(1, at=ax$x, labels=ax$x)
		    # axis(2, at=ax$y, labels=ax$y)
		    poly.image(loc2$x, loc2$y, invClim, main=paste(main3, " (level ", i, ")", sep=""), zlim=zl, axes=FALSE, add=TRUE, ...)
		    map(add=TRUE, lwd=2)
                    map(add=TRUE, database="state")
		    map.axes()
		}
		title("")
	    } else {
		map(xlim=locr[,1], ylim=locr[,2], type="n")
		# axis(1, at=ax$x, labels=ax$x)
		# axis(2, at=ax$y, labels=ax$y)
		image(as.image(invX, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main=paste(main1, " (level ", i, ")", sep=""),
                                zlim=zl, col=col, axes=FALSE, add=TRUE, ...)
		map(add=TRUE, lwd=2)
                map(add=TRUE, database="state")
		map.axes()

		map(xlim=locr[,1], ylim=locr[,2], type="n")
		# axis(1, at=ax$x, labels=ax$x)
		# axis(2, at=ax$y, labels=ax$y)
		image(as.image(invY, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main=paste(main2, " (level ", i, ")", sep=""),
				zlim=zl, col=col, axes=FALSE, add=TRUE, ...)
		map(add=TRUE, lwd=2)
                map(add=TRUE, database="state")
		map.axes()

		if(climate) {
		    map(xlim=locr[,1], ylim=locr[,2], type="n")
		    # axis(1, at=ax$x, labels=ax$x)
		    # axis(2, at=ax$y, labels=ax$y)
		    image(as.image(invClim, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE),
				main=paste(main3, " (level ", i, ")", sep=""), zlim=zl, axes=FALSE, add=TRUE, ...)
		    map(add=TRUE, lwd=2)
                    map(add=TRUE, database="state")
		    map.axes()
		}
		title("")
	    }

	} else {

	    if(proj) {
		poly.image(loc2$x, loc2$y, invX, main=paste(main1, " (level ", i, ")", sep=""), zlim=zl, col=col, axes=FALSE, ...)
		poly.image(loc2$x, loc2$y, invY, main=paste(main2, " (level ", i, ")", sep=""), zlim=zl, col=col, axes=FALSE, ...)
		if(climate) poly.image(loc2$x, loc2$y, invClim, main=paste(main3, " (level ", i, ")", sep=""), zlim=zl, col=col, axes=FALSE, ...)
	    } else {
	        image(invX, main=paste(main1, " (level ", i, ")", sep=""), zlim=zl, col=col, axes=FALSE, ...)
	        image(invY, main=paste(main2, " (level ", i, ")", sep=""), zlim=zl, col=col, axes=FALSE, ...)
	        if(climate) image(invX, main=paste(main3, " (level ", i, ")", sep=""), zlim=zl, col=col, axes=FALSE, ...) 
	    }
	    title("")
	} # end of if else 'domap' stmts.
	image.plot(invX, legend.only=TRUE, horizontal=horizontal, col=col, zlim=zl, ...)

	holdX[[HH[i]]] <- holdX[[HH[i]]]*0
	holdX[[HL[i]]] <- holdX[[HL[i]]]*0
	holdX[[LH[i]]] <- holdX[[LH[i]]]*0
	holdY[[HH[i]]] <- holdY[[HH[i]]]*0
        holdY[[HL[i]]] <- holdY[[HL[i]]]*0
        holdY[[LH[i]]] <- holdY[[LH[i]]]*0
	if(climate) {
	   holdClim[[HH[i]]] <- holdClim[[HH[i]]]*0
           holdClim[[HL[i]]] <- holdClim[[HL[i]]]*0
           holdClim[[LH[i]]] <- holdClim[[LH[i]]]*0
	}
      } # end of for 'i' loop.
    } # end of if do detail plots.
    par(mfrow=c(1,1), mar=c(5,4,4,2)+0.1)
    if(is.element("energy",which.plots)) {
	yl <- range(c(x$energy),na.rm=TRUE)
	plot(1:J, x$energy[,1], type="l", xlab="level", ylab="Energy", lwd=1.5, col="darkblue", ylim=yl, ...)
	lines(1:J, x$energy[,2], lty=2, col="darkblue", lwd=1.5, ylim=yl, ...)
	if(climate) {
	   lines(1:J, x$energy[,3], lty=3, col="green", lwd=1, ylim=yl, ...)
	   legend("topleft", legend=c(main1,main2,main3), col=c("darkblue", "darkblue", "green"), lty=1:3, lwd=c(1.5, 1.5, 1), bty="n")
	} else legend("topleft", legend=c(main1, main2), col="darkblue", lty=1:2, lwd=1.5, bty="n")
    } # end of plot energies stmts.
    if(any(is.element(c("all", "mse"),which.plots))) plot(1:J, x$mse, type="l", xlab="level", ylab="MSE", lwd=1.5, col="darkblue", ...)
    if(any(is.element(c("all", "rmse"),which.plots))) plot(1:J, x$rmse, type="l", xlab="level", ylab="RMSE", lwd=1.5, col="darkblue", ...)
    if(any(is.element(c("all", "acc"),which.plots)) & climate) plot(1:J, x$acc, type="l", xlab="level", ylab="ACC", lwd=1.5, col="darkblue", ...)
    if(!is.null(a$msg)) mtext(a$msg, line=0.05, outer=TRUE)
    invisible()
} # end of 'plot.waverify2d' function.

mowaverify2d <- function(X, ..., Clim=NULL, wavelet.type="haar", J=4) {
    UseMethod("mowaverify2d",X)
} # end of 'mowaverify2d' function.

mowaverify2d.SpatialVx <- function(X, ..., Clim=NULL, wavelet.type="haar", J=4,
				    useLL=FALSE, compute.shannon=FALSE, which.space="field",
				    time.point=1, obs = 1, model=1, verbose=FALSE) {

    out <- list()
    a <- attributes(X)
    attributes(out) <- a

    ## Begin: Get the data sets
    dat <- datagrabber(X, time.point=time.point, obs = obs, model=model)
    x <- dat$X
    y <- dat$Xhat
    ## End: Get the data sets

    out2 <- mowaverify2d.default(X=x, Clim=Clim, Y=y, wavelet.type=wavelet.type, J=J,
				    useLL=useLL, compute.shannon=compute.shannon, which.space=which.space, verbose=verbose)

    out$J <- out2$J
    out$X.wave <- out2$X.wave
    out$Y.wave <- out2$Y.wave
    out$energy <- out2$energy
    out$mse <- out2$mse
    out$rmse <- out2$rmse
    if(!is.null(out2$acc)) out$acc <- out2$acc
    out$wave.type <- out2$wave.type
    out$wave.method <- out2$wave.method

    attr(out, "time.point") <- time.point
    attr( out, "obs" ) <- obs
    attr(out, "model") <- model
    class(out) <- "waverify2d"
    return(out)
} # end of 'mowaverify2d.SpatialVx' function.

mowaverify2d.default <- function(X, ..., Clim=NULL, Y, wavelet.type="haar", J=4,
				    useLL=FALSE, compute.shannon=FALSE, which.space="field", verbose=FALSE) {
   out <- list()
   attributes(out) <- atmp <- list(...)
   if(is.null(atmp$loc.byrow)) attr(out, "loc.byrow") <- FALSE

   out$J <- J
   climate <- !is.null(Clim)
   xdim <- dim(X)
   if(verbose) cat("\n", "Finding MO wavelet decompositions.\n")
   wv.X <- modwt.2d(X, wf=wavelet.type, J=J)
   wv.Y <- modwt.2d(Y, wf=wavelet.type, J=J)
   if( climate) wv.Clim <- modwt.2d( Clim, wf=wavelet.type, J=J)
   if(verbose) cat("\n", "MO wavelet decompositions found.\n")
   out$X.wave <- wv.X
   out$Y.wave <- wv.Y
   if( climate) out$Clim.wave <- wv.Clim
   lnames <- makeWaveNames(J=J)
   HH <- lnames$HH
   HL <- lnames$HL
   LH <- lnames$LH
   LL <- lnames$LL

   mse.i <- rmse.i <- numeric(J) + NA
   if( climate) {
	acc.i <- rmse.i
	energy.i <- matrix( NA, nrow=J, ncol=3)
	colnames(energy.i) <- c("Observed", "Forecast", "Climatology")
   } else {
	energy.i <- matrix( NA, nrow=J, ncol=2)
	colnames( energy.i) <- c("Observed", "Forecast")
   }

   if( compute.shannon) {
      if(verbose) cat("\n", "Cacluating Shannon entropy for each field.\n")
      Shannon.entropy <- matrix( NA, nrow=1, ncol=2)
      colnames( Shannon.entropy) <- c("Verification", "Forecast")
      tmp <- abs(unlist(wv.X))
      tmp <- tmp/sum(tmp,na.rm=TRUE)
      id0 <- tmp==0
      tmp[id0] <- 0
      tmp[!id0] <- tmp[!id0]*log(tmp[!id0])
      Shannon.entropy[,1] <- -sum(tmp, na.rm=TRUE)
      tmp <- abs(unlist(wv.Y))
      tmp <- tmp/sum(tmp,na.rm=TRUE)
      id0 <- tmp==0
      tmp[id0] <- 0
      tmp[!id0] <- tmp[!id0]*log(tmp[!id0])
      Shannon.entropy[,2] <- -sum(tmp, na.rm=TRUE)
      out$Shannon.entropy <- Shannon.entropy
      if(verbose) cat("\n", "Shannon entropy found.\n")
   }

   if(verbose) cat("n", "Calculating statistics for each detail field.\n")
   for(i in 1:J) {
	energy.i[i,1] <- sum( wv.X[[HH[i]]]^2, na.rm=TRUE) + sum( wv.X[[HL[i]]]^2, na.rm=TRUE) + sum( wv.X[[LH[i]]]^2, na.rm=TRUE)
	energy.i[i,2] <- sum( wv.Y[[HH[i]]]^2, na.rm=TRUE) + sum( wv.Y[[HL[i]]]^2, na.rm=TRUE) + sum( wv.Y[[LH[i]]]^2, na.rm=TRUE)
	if( climate) energy.i[i,3] <- sum( wv.Clim[[HH[i]]]^2, na.rm=TRUE) + sum( wv.Clim[[HL[i]]]^2, na.rm=TRUE) + sum( wv.Clim[[LH[i]]]^2, na.rm=TRUE)

	invX <- detailer(wv.X, level=i, which.space=which.space, trans.type="MODWT", useLL=useLL, lnames=lnames, J=J)
	invY <- detailer(wv.Y, level=i, which.space=which.space, trans.type="MODWT", useLL=useLL, lnames=lnames, J=J)
	if( climate) invClim <- detailer(wv.Clim, level=i, which.space=which.space, trans.type="MODWT", useLL=useLL, lnames=lnames, J=J)

	mse.i[i] <- mean( (c(invX) - c(invY))^2, na.rm=TRUE)
   	rmse.i[i] <- sqrt( mse.i[i])

	if( climate) {
           denom <- sqrt( sum( colSums((c( invX) - c( invClim))^2,na.rm=TRUE), na.rm=TRUE))*sqrt( sum( colSums((c(invY) - c(invClim))^2,na.rm=TRUE), na.rm=TRUE))
           numer <- sum( diag( t(invX - invClim) %*% (invY - invClim)), na.rm=TRUE)
           acc.i[i] <- numer/denom
	}

   } # end of for 'i' loop.
   out$energy <- energy.i
   out$mse <- mse.i
   out$rmse <- rmse.i
   if(climate) out$acc <- acc.i
   out$wave.type <- wavelet.type
   out$wave.method <- "MODWT"
   class(out) <- "waverify2d"
   return(out)
} # end of 'mowaverify2d.default' function.

energizer <- function(x,lnames=NULL,J=NULL) {
   if(is.null(J)) J <- (length(x)-1)/3
   if(is.null(lnames)) lnames <- makeWaveNames(J=J)
   HH <- lnames$HH
   HL <- lnames$HL
   LH <- lnames$LH
   LL <- lnames$LL
   out <- numeric(J)+NA
   for(i in 1:J) {
	good <- !(is.na(x[[HH[i]]]) & is.na(x[[HL[i]]]) & is.na(x[[LH[i]]]))
	N <- length(x[[HH[i]]][good])
	out[i] <- sum(colSums(x[[HH[i]]]^2 + x[[HL[i]]]^2 + x[[LH[i]]]^2,na.rm=TRUE), na.rm=TRUE)/N
   } # end of for 'i' loop.
   return(out)
} # end of 'energizer' function.

detailer <- function(x, level, which.space="field", trans.type="DWT", useLL=FALSE, lnames=NULL, J=NULL) {
   if(is.null(lnames)) {
	if(is.null(J)) stop("detailer: must supply J if lnames is NULL.")
	lnames <- makeWaveNames(J=J)
   }
   HH <- lnames$HH
   HL <- lnames$HL
   LH <- lnames$LH
   LL <- lnames$LL

   if(which.space=="field") {
	hold <- makeEmptyWave(x)
	if(useLL) hold[[LL]] <- x[[LL]]
	hold[[HH[level]]] <- x[[HH[level]]]
	hold[[HL[level]]] <- x[[HL[level]]]
	hold[[LH[level]]] <- x[[LH[level]]]
	ni <- length(x[[HH[level]]])*3
	lam <- sqrt(sum(c(hold[[HH[level]]]^2)+c(hold[[HL[level]]]^2)+c(hold[[LH[level]]]^2),na.rm=TRUE)/ni)
	hold[[HH[level]]] <- sign(hold[[HH[level]]])*pmax(c(abs(hold[[HH[level]]])-lam),0)
	hold[[HL[level]]] <- sign(hold[[HL[level]]])*pmax(c(abs(hold[[HL[level]]])-lam),0)
	hold[[LH[level]]] <- sign(hold[[LH[level]]])*pmax(c(abs(hold[[LH[level]]])-lam),0)
	if(trans.type=="DWT") out <- idwt.2d(hold)
	else out <- imodwt.2d(hold)
   } else out <- (x[[HH[level]]]/(2^level))^2 + (x[[HL[level]]]/(2^level))^2 + (x[[LH[level]]]/(2^level))^2
   return(out)
} # end of 'detailer' function.

# wavePurifyVx <- function( x, y, object=NULL, climate=NULL,
#     which.stats=c("bias", "ts", "ets", "pod", "far", "f", "hk", "mse"),
#     thresholds=NULL, return.fields=FALSE, time.point=1, obs = 1, model=1, verbose=FALSE, ...) {
# 
#    if(verbose) begin.time <- Sys.time()
#    isnt.xy <- missing(x) | missing(y)
#    if(isnt.xy & !is.null(object)) {
# 	out <- list()
# 	out.attr <- attributes(object)
# 	attributes(out) <- out.attr
# 
# 	## Begin: Get the data sets
#          dat <- datagrabber(object, time.point=time.point, obs = obs, model=model)
# 
#         x <- dat$X
#         y <- dat$Xhat
#         ## End: Get the data sets
# 
# 	xdim <- out.attr$xdim
# 	u <- out.attr$thresholds
# 	thresholds <- cbind( u$X[, obs ], u$Xhat[, model ] )
# 
#    } else if(!isnt.xy) {
# 
# 	out <- list()
# 	attributes(out) <- atmp <- list(...)
# 	if(is.null(atmp$loc.byrow)) attr(out, "loc.byrow") <- FALSE
# 	xdim <- dim(x)
# 	if(is.null(thresholds) & any(is.element(c("bias", "ts", "ets", "pod", "far", "f", "hk"), which.stats))) {
# 	   thresholds <- cbind( quantile(c(y), probs=c(0, 0.1, 0.25, 0.33, 0.5, 0.66, 0.75, 0.9, 0.95)),
# 				quantile(c(x), probs=c(0, 0.1, 0.25, 0.33, 0.5, 0.66, 0.75, 0.9, 0.95)))
# 	   out$qs <- c("0", "0.1", "0.25", "0.33", "0.5", "0.66", "0.75", "0.9", "0.95")
# 	} else if(!is.matrix(thresholds)) thresholds <- cbind(thresholds,thresholds)
# 	colnames(thresholds) <- c("Forecast","Verification")
# 	out$thresholds <- thresholds
# 	out$Vx.name <- as.character(substitute(x))
# 	out$Fcst.name <- as.character(substitute(y))
# 
#    } else stop("wavePurifyVx: Must supply arguments x and y, or object.")
# 
#    args <- list(...)
#    if(isnt.J <- is.null(args$J)) {
# 	J <- floor(log2(min(xdim)))
# 	args$J <- J
#    } else {
# 	J <- args$J
#    }
#    out$args <- args
#    if(all(floor(log2(xdim))==ceiling(log2(xdim)))) dyadic <- TRUE
#    else dyadic <- FALSE
#    if(verbose) cat("\n", "Denoising the fields.\n")
#    if(dyadic) {
# 	if(isnt.J) {
# 	   X <- denoise.dwt.2d(x=x, J=J, ...)
# 	   Y <- denoise.dwt.2d(x=y, J=J, ...)
# 	   if(!is.null(climate)) Climate <- denoise.dwt.2d(x=climate, J=J, ...)
# 	} else {
# 	   X <- denoise.dwt.2d(x=x, ...)
# 	   Y <- denoise.dwt.2d(x=y, ...)
# 	   if(!is.null(climate)) Climate <- denoise.dwt.2d(x=climate, ...)
# 	} # end of if else 'isnt.J' stmts.
#    } else {
# 	if(isnt.J) {
#            X <- denoise.modwt.2d(x=x, J=J, ...)
#            Y <- denoise.modwt.2d(x=y, J=J, ...)
# 	   if(!is.null(climate)) Climate <- denoise.modwt.2d(x=climate, J=J, ...)
#         } else {
#            X <- denoise.modwt.2d(x=x, ...)
#            Y <- denoise.modwt.2d(x=y, ...)
# 	   if(!is.null(climate)) Climate <- denoise.modwt.2d(x=climate, ...)
#         } # end of if 'is.J' stmts.
#    } # end of if else 'dyadic' stmts.
#    if(return.fields) {
# 	out$X <- x
# 	out$Xhat <- y
# 	out$X2 <- X
# 	out$Y2 <- Y
# 	if(!is.null(climate)) {
# 	   out$Climate <- climate
# 	   out$Climate2 <- Climate
# 	}
#    }
#    if(!is.null(thresholds)) q <- dim(thresholds)[1]
#    else q <- 1
#    if(is.element("bias",which.stats)) out$bias <- numeric(q)+NA 
#    if(is.element("ts",which.stats)) out$ts <- numeric(q)+NA 
#    if(is.element("ets",which.stats)) out$ets <- numeric(q)+NA 
#    if(is.element("pod",which.stats)) out$pod <- numeric(q)+NA 
#    if(is.element("far",which.stats)) out$far <- numeric(q)+NA 
#    if(is.element("f",which.stats)) out$f <- numeric(q)+NA 
#    if(is.element("hk",which.stats)) out$hk <- numeric(q)+NA 
#    if(is.element("mse",which.stats)) out$mse <- numeric(q)+NA
#    if(!is.null(climate)) out$acc <- numeric(q)+NA
#    binmat <- matrix(0, xdim[1], xdim[2])
#    if(verbose) {
# 	cat("\n", "Looping through thresholds = \n") 
# 	print(thresholds)
# 	cat("\n")
#    }
#    for( threshold in 1:q) {
#       if(verbose) cat(threshold, " ")
#       if(is.element("mse",which.stats)) {
# 	if(!is.null(thresholds)) {
# 	   X2 <- X
# 	   Y2 <- Y
# 	   X2[X<thresholds[threshold,2]] <- 0
# 	   Y2[Y<thresholds[threshold,1]] <- 0
# 	   out$mse[threshold] <- vxstats(Y2,X2,which.stats="mse")$mse
# 	} else out$mse[threshold] <- vxstats(Y,X,which.stats="mse")$mse
#       } # end of if do mse stmts.
#       if(any(is.element(c("bias", "ts", "ets", "pod", "far", "f", "hk"), which.stats))) {
# 	Xbin <- Ybin <- binmat
# 	Xbin[X>=thresholds[threshold,2]] <- 1
# 	Ybin[Y>=thresholds[threshold,1]] <- 1
# 	if(threshold==1) {
# 	   dostats <- which.stats
# 	   if(is.element("mse",dostats)) dostats <- dostats[dostats != "mse"]
# 	}
# 	tmp <- vxstats(Ybin,Xbin,which.stats=dostats)
# 	if(is.element("bias",which.stats)) out$bias[threshold] <- tmp$bias 
# 	if(is.element("ts",which.stats)) out$ts[threshold] <- tmp$ts 
# 	if(is.element("ets",which.stats)) out$ets[threshold] <- tmp$ets 
# 	if(is.element("pod",which.stats)) out$pod[threshold] <- tmp$pod 
# 	if(is.element("far",which.stats)) out$far[threshold] <- tmp$far 
# 	if(is.element("f",which.stats)) out$f[threshold] <- tmp$f 
# 	if(is.element("hk",which.stats)) out$hk[threshold] <- tmp$hk
#       } # end of if any vxstats besides mse stmts.
#       if(!is.null(climate)) {
# 	if(!is.null(thresholds)) {
#            X2 <- X
#            Y2 <- Y
# 	   Clim <- Climate
#            X2[X<thresholds[threshold,2]] <- 0
#            Y2[Y<thresholds[threshold,1]] <- 0
# 	   Clim[climate<thresholds[threshold,1]] <- 0
# 	   denom <- sqrt(sum(colSums((Y2-Clim)^2,na.rm=TRUE),na.rm=TRUE))*sqrt(sum(colSums((X2-Clim)^2,na.rm=TRUE),na.rm=TRUE))
#            numer <- sum(diag(t(Y2-Clim) %*% (X2-Clim)),na.rm=TRUE)
#         } else {
# 	   denom <- sqrt(sum(colSums((Y-Climate)^2,na.rm=TRUE),na.rm=TRUE))*sqrt(sum(colSums((X-Climate)^2,na.rm=TRUE),na.rm=TRUE))
# 	   numer <- sum(diag(t(Y-Climate) %*% (X-Climate)),na.rm=TRUE)
# 	}
# 	out$acc[threshold] <- numer/denom
#       } # end of if do acc stmts.
#    } # end of for 'threshold' loop.
# 
#    attr(out, "time.point") <- time.point
#    attr(out, "obs") <- obs
#    attr(out, "model") <- model   
#    class(out) <- "wavePurifyVx"
#    if(verbose) print(Sys.time() - begin.time)
#    return(out)
# } # end of 'wavePurifyVx' function.

# summary.wavePurifyVx <- function(object, ...) {
# 
#    object.attr <- attributes(object)
# 
#    if(!is.null(object.attr$thresholds)) {
# 	u <- Thresholds <- object.attr$thresholds
#    	q <- dim(u)[1]
#    	cat("\n", "Thresholds applied are:\n")
#    	print(Thresholds)
# 	   if(is.null(object.attr$qs)) {
# 		if(all(u[,1]==u[,2])) ulab <- as.character(u[,1])
# 		else ulab <- as.character(1:q)
# 	   } else ulab <- object.attr$qs
# 	   if(!is.null(object$bias)) {
# 		bias <- matrix(object$bias, nrow=1)
# 		colnames(bias) <- ulab
# 		cat("\n", "Frequency Bias: \n")
# 		print(bias)
# 	   }
# 	   if(!is.null(object$ts)) {
# 		ts <- matrix(object$ts, nrow=1)
# 		colnames(ts) <- ulab
# 		cat("\n", "Threat Score: \n")
# 		print(ts)
# 	   }
# 	   if(!is.null(object$ets))  {
# 		ets <- matrix(object$ets, nrow=1)
# 		colnames(ets) <- ulab
# 		cat("\n", "Gilbert Skill Score: \n")
# 		print(ets)
# 	   }
# 	   if(!is.null(object$pod))  {
# 		pod <- matrix(object$pod, nrow=1)
# 		colnames(pod) <- ulab
# 		cat("\n", "Probability of Detecting an Event (POD): \n")
# 		print(pod)
# 	   }
# 	   if(!is.null(object$far)) {
# 		far <- matrix(object$far, nrow=1)
# 		colnames(far) <- ulab
# 		cat("\n", "False Alarm Ratio: \n")
# 		print(far)
# 	   }
# 	   if(!is.null(object$f)) {
# 		f <- matrix(object$f, nrow=1)
# 		colnames(f) <- ulab
# 		cat("\n", "False Alarm Rate: \n")
# 		print(f)
# 	   }
# 	   if(!is.null(object$hk)) {
# 		hk <- matrix(object$hk, nrow=1)
# 		colnames(hk) <- ulab
# 		cat("\n", "Hanssen-Kuipers Score: \n")
# 		print(hk)
# 	   }
# 	   if(!is.null(object$mse)) {
# 		mse <- matrix(object$mse, nrow=1)
# 		colnames(mse) <- ulab
# 		cat("\n", "MSE: \n")
# 		print(mse)
# 	   }
# 	   if(!is.null(object$acc)) {
# 		acc <- matrix(object$acc, nrow=1)
# 		colnames(acc) <- ulab
# 		cat("\n", "ACC: \n")
# 		print(acc)
# 	   }
#    } else {
# 	if(!is.null(object$mse)) {
#                 mse <- matrix(object$mse, nrow=1)
# 		cat("\n", "MSE: \n")
#                 print(mse)
#            }
#            if(!is.null(object$acc)) {
#                 acc <- matrix(object$acc, nrow=1)
# 		cat("\n", "ACC: \n")
#                 print(acc)
#            } 
#    } # end of whether or not thresholds applied stmts.
#    invisible()
# } # end of 'summary.wavePurifyVx' function.

# plot.wavePurifyVx <- function(x, ..., col, zlim, horizontal=TRUE, type=c("stats", "fields")) {
# 
#    type <- tolower(type)
#    type <- match.arg(type)
# 
#    if(missing(col)) col <- c("gray", tim.colors(64))
# 
#    a <- attributes(x)
#    loc.byrow <- a$loc.byrow
#    nm <- names(x)
# 
#     op <- par() 
#     par( oma = c(0, 0, 2, 0) )
# 
#    denX <- !is.null(x$X2)
# 
#    if(denX && type=="fields") {
#         domap <- a$map
# 	if(is.null(domap)) domap <- FALSE
# 
# 	loc <- a$loc
# 
# 	if(domap) {
# 	    locr <- apply(loc,2,range,finite=TRUE)
# 	    ax <- list(x=pretty(round(loc[,1], digits=2)), y=pretty(round(loc[,2], digits=2)))
# 	}
# 
# 	X <- x$X
# 	Xhat <- x$Xhat
# 	X2 <- x$X2
# 	Y2 <- x$Y2
# 
# 	xd <- a$xdim
# 	if(is.null(xd)) xd <- dim(X)
# 	if(is.null(loc)) loc <- cbind(rep(1:xd[1],xd[2]), rep(1:xd[2],each=xd[1]))
# 
# 	ptitle <- a$data.name
# 	if(is.null(ptitle)) ptitle <- c("","")
# 	else {
# 	    if(length(ptitle) == a$nforecast + 2) ptitle <- ptitle[-1]
# 
# 	    if(a$nforecast > 1) {
# 	        if(!is.numeric(a$model)) {
# 		    modnum <- (1:length(ptitle))[ptitle==a$model] - 1
# 	        } else modnum <- a$model
# 	            ptitle <- c(ptitle[1], ptitle[modnum + 1])
# 	    }
# 
# 	    tiid <- a$time
# 	    if(length(tiid) > 1) ptitle <- paste(ptitle, " (t = ", tiid[a$time.point], ")", sep="")
# 	} # end of if 'ptitle' is NULL stmts.
# 
# 	if(!is.null(x$Climate)) {
# 		Clim <- x$Climate
# 		Clim2 <- x$Climate2
# 		if(missing(zlim)) zl <- range( c(c(X), c(Xhat), c(X2), c(Y2), c(Clim), c(Clim2)), finite=TRUE)
# 	} else if(missing(zlim)) zl <- range( c(c(X), c(Xhat), c(X2), c(Y2)), finite=TRUE)
# 
# 	proj <- a$projection
# 	if(is.null(proj)) proj <- FALSE
# 
# 	if(proj) {
# 	    if(domap) {
# 		map(xlim=locr[,1], ylim=locr[,2], type="n")
# 		axis(1, at=ax$x, labels=ax$x)
# 		axis(2, at=ax$y, labels=ax$y)
# 		poly.image(matrix(loc[,1], xd[1], xd[2], byrow=loc.byrow), matrix(loc[,2], xd[1], xd[2], byrow=loc.byrow),
# 		    X, main=ptitle[1], col=col, zlim=zl, add=TRUE)
# 		map(add=TRUE, lwd=1.5)
# 		map(add=TRUE, database="state")
# 
# 		map(xlim=locr[,1], ylim=locr[,2], type="n")
# 		axis(1, at=ax$x, labels=ax$x)
# 		axis(2, at=ax$y, labels=ax$y)
#                 poly.image(matrix(loc[,1], xd[1], xd[2], byrow=loc.byrow), matrix(loc[,2], xd[1], xd[2], byrow=loc.byrow),
#                     Xhat, main=ptitle[2], col=col, zlim=zl, add=TRUE)
#                 map(add=TRUE, lwd=1.5)
#                 map(add=TRUE, database="state")
# 
# 		if(!is.null(x$Climate)) {
#                     map(xlim=locr[,1], ylim=locr[,2], type="n")
# 		axis(1, at=ax$x, labels=ax$x)
# 		axis(2, at=ax$y, labels=ax$y)
#                     poly.image(matrix(loc[,1], xd[1], xd[2], byrow=loc.byrow), matrix(loc[,2], xd[1], xd[2],
#                         byrow=loc.byrow), Clim, main="Climate", col=col, zlim=zl, add=TRUE)
#                     map(add=TRUE, lwd=1.5)
#                     map(add=TRUE, database="state")
#                 }
# 
# 		map(xlim=locr[,1], ylim=locr[,2], type="n")
# 		axis(1, at=ax$x, labels=ax$x)
# 		axis(2, at=ax$y, labels=ax$y)
#                 poly.image(matrix(loc[,1], xd[1], xd[2], byrow=loc.byrow), matrix(loc[,2], xd[1], xd[2], byrow=loc.byrow),
#                     X2, main=paste(ptitle[1], " de-noised", sep=""), col=col, zlim=zl, add=TRUE)
#                 map(add=TRUE, lwd=1.5)
#                 map(add=TRUE, database="state")
# 
# 		map(xlim=locr[,1], ylim=locr[,2], type="n")
# 		axis(1, at=ax$x, labels=ax$x)
# 		axis(2, at=ax$y, labels=ax$y)
#                 poly.image(matrix(loc[,1], xd[1], xd[2], byrow=loc.byrow), matrix(loc[,2], xd[1], xd[2], byrow=loc.byrow),
#                     Y2, main=paste(ptitle[2], " de-noised", sep=""), col=col, zlim=zl, add=TRUE)
#                 map(add=TRUE, lwd=1.5)
#                 map(add=TRUE, database="state")
# 
# 		if(!is.null(x$Climate)) {
#                     map(xlim=locr[,1], ylim=locr[,2], type="n")
# 		axis(1, at=ax$x, labels=ax$x)
# 		axis(2, at=ax$y, labels=ax$y)
#                     poly.image(matrix(loc[,1], xd[1], xd[2], byrow=loc.byrow), matrix(loc[,2], xd[1], xd[2],
#                         byrow=loc.byrow), Clim2, main="Climate de-noised", col=col, zlim=zl, add=TRUE)
#                     map(add=TRUE, lwd=1.5)
#                     map(add=TRUE, database="state")
#                 }
# 
# 	    } else {
# 		poly.image(matrix(loc[,1], xd[1], xd[2], byrow=loc.byrow), matrix(loc[,2], xd[1], xd[2], byrow=loc.byrow),
#                     X, main=ptitle[1], col=col, zlim=zl)
# 		poly.image(matrix(loc[,1], xd[1], xd[2], byrow=loc.byrow), matrix(loc[,2], xd[1], xd[2], byrow=loc.byrow),
#                     Xhat, main=ptitle[2], col=col, zlim=zl)
# 		if(!is.null(x$Climate)) {
# 		    poly.image(matrix(loc[,1], xd[1], xd[2], byrow=loc.byrow),
# 				matrix(loc[,2], xd[1], xd[2], byrow=loc.byrow), Clim, main="Climate", col=col, zlim=zl)
# 		}
# 		poly.image(matrix(loc[,1], xd[1], xd[2], byrow=loc.byrow), matrix(loc[,2], xd[1], xd[2], byrow=loc.byrow),
#                     X2, main=paste(ptitle[1], " de-noised", sep=""), col=col, zlim=zl)
# 		poly.image(matrix(loc[,1], xd[1], xd[2], byrow=loc.byrow), matrix(loc[,2], xd[1], xd[2], byrow=loc.byrow),
#                     Y2, main=paste(ptitle[2], " de-noised", sep=""), col=col, zlim=zl)
# 		if(!is.null(x$Climate)) {
#                     poly.image(matrix(loc[,1], xd[1], xd[2], byrow=loc.byrow),
#                                 matrix(loc[,2], xd[1], xd[2], byrow=loc.byrow), Clim2, main="Climate de-noised",
# 				col=col, zlim=zl)
#                 }
# 	    }
# 	} else {
# 	    if(domap) {
#                 map(xlim=locr[,1], ylim=locr[,2], type="n")
# 		axis(1, at=ax$x, labels=ax$x)
# 		axis(2, at=ax$y, labels=ax$y)
#                 image(as.image(X, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main=ptitle[1], col=col, zlim=zl, add=TRUE)
#                 map(add=TRUE, lwd=1.5)
#                 map(add=TRUE, database="state")
# 
#                 map(xlim=locr[,1], ylim=locr[,2], type="n")
# 		axis(1, at=ax$x, labels=ax$x)
# 		axis(2, at=ax$y, labels=ax$y)
#                 image(as.image(Xhat, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main=ptitle[2], col=col, zlim=zl, add=TRUE)
#                 map(add=TRUE, lwd=1.5)
#                 map(add=TRUE, database="state")
# 
#                 if(!is.null(x$Climate)) {
#                     map(xlim=locr[,1], ylim=locr[,2], type="n")
# 		axis(1, at=ax$x, labels=ax$x)
# 		axis(2, at=ax$y, labels=ax$y)
#                     image(as.image(Clim, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main="Climate", col=col, zlim=zl, add=TRUE)
#                     map(add=TRUE, lwd=1.5)
#                     map(add=TRUE, database="state")
#                 }
# 
#                 map(xlim=locr[,1], ylim=locr[,2], type="n")
# 		axis(1, at=ax$x, labels=ax$x)
# 		axis(2, at=ax$y, labels=ax$y)
#                 image(as.image(X2, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main=paste(ptitle[1], " de-noised", sep=""),
# 		    col=col, zlim=zl, add=TRUE)
#                 map(add=TRUE, lwd=1.5)
#                 map(add=TRUE, database="state")
# 
#                 map(xlim=locr[,1], ylim=locr[,2], type="n")
# 		axis(1, at=ax$x, labels=ax$x)
# 		axis(2, at=ax$y, labels=ax$y)
#                 image(as.image(Y2, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main=paste(ptitle[2], " de-noised", sep=""),
# 		    col=col, zlim=zl, add=TRUE)
#                 map(add=TRUE, lwd=1.5)
#                 map(add=TRUE, database="state")
# 
#                 if(!is.null(x$Climate)) {
#                     map(xlim=locr[,1], ylim=locr[,2], type="n")
# 		axis(1, at=ax$x, labels=ax$x)
# 		axis(2, at=ax$y, labels=ax$y)
#                     image(as.image(Clim2, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main="Climate de-noised", col=col, zlim=zl, add=TRUE)
#                     map(add=TRUE, lwd=1.5)
#                     map(add=TRUE, database="state")
#                 }
# 
#             } else {
#                 image(as.image(X, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main=ptitle[1], col=col, zlim=zl)
#                 image(as.image(Xhat, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main=ptitle[2], col=col, zlim=zl)
#                 if(!is.null(x$Climate)) {
#                     image(as.image(Clim, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main="Climate", col=col, zlim=zl)
#                 }
#                 image(as.image(X2, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main=paste(ptitle[1], " de-noised", sep=""),
# 		    col=col, zlim=zl)
#                 image(as.image(Y2, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main=paste(ptitle[2], " de-noised", sep=""),
# 		    col=col, zlim=zl)
#                 if(!is.null(x$Climate)) {
#                     image(as.image(Clim2, nx=xd[1], ny=xd[2], x=loc, na.rm=TRUE), main="Climate de-noised",
#                                 col=col, zlim=zl)
#                 }
#             }
# 	} # end of if else 'projection' stmts.
#         image.plot(X, zlim=zl, col=col, horizontal=horizontal, legend.only=TRUE)	
# 	return(invisible())
#    } else if(!denX && type == "fields") stop("plot.wavePurifyVx: return.fields must be TRUE in call to wavePurifyVx to use type fields.")
# 
#     if(!is.null(a$thresholds)) {
# 
#       q <- dim( a$thresholds$X )[ 1 ]
# 
#       if(q  >= 2 ) {
# 	 if(is.null(a$units)) unt <- NULL
#          else unt <- paste("(", a$units, ")", sep="")
#          if(any(is.element(c("bias","ts","ets","pod","far","hk"),nm))) {
# 		if(is.element("hk",nm)) yl <- c(-1,1)
# 		else if(is.element("ets",nm)) yl <- c(-1/3,1)
# 		else yl <- c(0,1)
# 		plot(1:q, runif(q,-1,1), type="n", xaxt="n", ylim=yl, xlab=paste("Threshold ", unt, sep=""), ylab="")
# 		legnames <- c("TS", "GSS", "POD", "FAR", "F", "HK", "Bias")[is.element(c("ts", "ets", "pod", "far", "f", "hk", "bias"),nm)]
# 		nl <- length(legnames)
# 		if(nl==7) cl <- c(1:6,8)
# 		else cl <- 1:nl
# 		legend("topright", legend=legnames, lty=cl, col=cl, bty="n")
# 		i <- 1
# 		if(is.element("ts", nm)) {
# 		   lines(1:q, x$ts, col=i, lty=i, lwd=1.5)
# 		   i <- i+1
# 		}
# 		if(is.element("ets", nm)) {
# 		   lines(1:q, x$ets, col=i, lty=i, lwd=1.5)
#                    i <- i+1
# 		}
# 		if(is.element("pod", nm)) {
#                    lines(1:q, x$pod, col=i, lty=i, lwd=1.5)
#                    i <- i+1
#                 }
# 		if(is.element("far", nm)) {
#                    lines(1:q, x$far, col=i, lty=i, lwd=1.5)
#                    i <- i+1
#                 }
# 		if(is.element("f", nm)) {
#                    lines(1:q, x$f, col=i, lty=i, lwd=1.5)
#                    i <- i+1
#                 }
# 		if(is.element("hk", nm)) {
#                    lines(1:q, x$hk, col=i, lty=i, lwd=1.5)
#                    i <- i+1
#                 }
# 		if(is.element("bias", nm)) {
# 		   if(i==7) i <- 8
# 		   brg <- range(x$bias, finite=TRUE)
# 		   if(length(unique(brg))==1) brg <- c(brg[1]-0.01, brg[2]+0.01)
# 		   par(usr=c(1,q,brg))
#                    lines(1:q, x$bias, col=i, lty=i, lwd=1.5)
# 		   bax <- pretty(sort(unique(x$bias)))
# 		   axis(4, at=bax, labels=bax, col=i)
#                 }
#       	   } # end of if 'any c("bias", "ts", "ets", "pod", "far", "f", "hk")' to make in plot stmts.
# 	   if(is.element("mse",nm)) plot(1:q, x$mse, type="l", col="darkblue", xlab=paste("Threshold ", unt, sep=""), ylab="MSE")
# 	   if(is.element("acc",nm)) plot(1:q, x$acc, type="l", col="darkblue", xlab=paste("Threshold ", unt, sep=""), ylab="ACC")
# 	} # end of if more than one threshold stmt.
#    } # end of if 'thresholds' stmts.
# 
#    mtext( a$msg, line = 0.05, outer = TRUE )
# 
#     par( oma = op$oma )
# 
#    invisible()
# 
# } # end of 'plot.wavePurifyVx' function.

Try the SpatialVx package in your browser

Any scripts or data that you put into this service are public.

SpatialVx documentation built on March 28, 2021, 1:10 a.m.