R/RFempvario-Methods-plots.R

In RandomFields: Simulation and Analysis of Random Fields

## Authors 
## Martin Schlather, schlather@math.uni-mannheim.de
##
##
## Copyright (C) 2012 -- 2014 Alexander Malinowski & Martin Schlather
##               2015 -- 2017 Martin Schlather
##
## This program is free software; you can redistribute it and/or
## modify it under the terms of the GNU General Public License
## as published by the Free Software Foundation; either version 3
## of the License, or (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.  



## Methods for classes 'RFempVario' and 'RFfit'  #######################

summary_RFempVariog <- function(OP, object, ...) {
  variogram <- do.call(OP, list(object, "empirical"))
  if (is.array(variogram)) {
    dims <- dim(variogram)
    dims <- dims[dims > 1]
    dim(variogram) <- dims
  }
  l <- list(centers=do.call(OP, list(object, "centers")),
            variogram=variogram)
  obj <- do.call(OP, list(object, "sd"))
  if (length(obj) > 0) {
    l$sd <- obj
    if (is.array(variogram)) dim(l$sd) <- dims
  }

  obj <- do.call(OP, list(object, "phi.centers"))
  if (length(obj) > 1) l$phi <- obj
  obj <- do.call(OP, list(object, "theta.centers"))
  if (length(obj) > 1) l$theta <- obj
  obj <- do.call(OP, list(object, "T"))
  if (length(obj)>0) l$T <- obj
  obj <- do.call(OP, list(object, "vdim"))
  if (obj > 1) l$vdim <- obj
  class(l) <- "summary.RFempVariog"
  l
}

summary.RFempVariog <- function(object, ...) summary_RFempVariog("@",object,...)
summary.RF_empVariog <- function(object,...) summary_RFempVariog("$",object,...)



print.summary.RFempVariog <- function(x, ...) {
  cat("Object of class 'RFempVariog'\n")
  str(x, no.list=TRUE) #
  invisible(x)
}

print.RFempVariog <- function(x, ...) {
  print.summary.RFempVariog(summary.RFempVariog(x, ...))
}
print.RF_empVariog <- function(x, ...) {
  print.summary.RFempVariog(summary.RF_empVariog(x, ...))
}


setMethod(f="show", signature="RFempVariog",
          definition=function(object) print.RFempVariog(object)) #

## coercion methods
setAs("RFempVariog", "list",
      def=function(from) {
        li <- list()
        for (name in slotNames(from)) {
          li[[name]] <- eval(parse(text=paste("from@", name, sep="")))
        }
        return(li)
      })

setAs("RFfit", "RFempVariog", def=function(from) from@ev)

list2RFempVariog <- function(li) {
  RFopt.coords <- RFoptions()$coords
  return(new("RFempVariog",
             centers=li$centers,
             empirical=li$empirical,
             var=var,
             sd=li$sd,
             n.bin=li$n.bin,
             phi.centers=li$phi.centers,
             theta.centers=li$theta.centers,
             T=li$Tbins,
             vdim = li$vdim,
             coordunits = RFopt.coords$coordunits,
             varunits = RFopt.coords$varunits,
             call=li$call,
	     method=li$method))
}


## plot method


plotRFempVariogUnbinned <- function(x, coordunits, varunits, varnames,
                                    legend, ..., plotmethod="plot") {
  stop("currently not used ") # !!
  graphics <- RFoptions()$graphics
  dots = mergeWithGlobal(list(...))
  dotnames <- names(dots)
  coords <- GridTopology2gridVectors(cbind(x@centers$x, x@centers$T))  
  if (length(coords)>1) {    
    coords[[1]] <- sort(unique((coords[[1]] - min(coords[[1]])) *
			       rep(c(-1, 1), each=length(coords[[1]]))))
    lab.names <- dimnames(x@centers$x)[[2]]
    
    if (!(x@centers$spatialdim==1 && x@centers$has.time.comp) &&
	!(x@centers$x[3,2]==dim(x@empirical)[2]))
      coords[[2]] <- sort(unique((coords[[2]] - min(coords[[2]])) *
				 rep(c(-1, 1), each=length(coords[[2]]))))
    else
      lab.names[2] <- "T"
  } else {
    x@centers <- coords[[1]]-min(coords[[1]])
    do.call(graphics::plot, args=c(dots, list(x=x, plot.nbin=FALSE)))
    return()
  }
  
  if (!("main" %in% dotnames)) {
    main <- "Variogram image plot"
    if (length(varnames)>0) main <- paste(main, "for", varnames)
    dots$main <- main
  }
  lab.names <- paste(lab.names, "-distance", sep="")
  
  idx <- lab.names != "T-distance"
  if (any(idx) && all(coordunits[idx] != ""))
    lab.names[idx] <-
  paste(lab.names[idx], " [", coordunits[idx], "]", sep="")
  if (!all(idx) && all(coordunits[!idx] != ""))
    lab.names[!idx] <-
  paste(lab.names[!idx], " [", coordunits[!idx], "]", sep="")
  
  if (!("xlab" %in% dotnames)) dots$xlab <- lab.names[1]
  if (!("ylab" %in% dotnames)) dots$ylab <- lab.names[2]
  if (!("xlim" %in% dotnames)) dots$xlim <- range(coords[[1]])
  if (!("ylim" %in% dotnames)) dots$ylim <- range(coords[[2]])
  
  idx1 <- coords[[1]] >= dots$xlim[1] & coords[[1]] <= dots$xlim[2]
  idx2 <- coords[[2]] >= dots$ylim[1] & coords[[2]] <= dots$ylim[2]
  coords[[1]] <- coords[[1]][idx1]
  coords[[2]] <- coords[[2]][idx2]
  
  dims <- dim(x@empirical)
  ev.plot <- matrix(x@empirical, nrow=dims[1], ncol=dims[2])
  ev.plot <- ev.plot[idx1, idx2]
  
  range.ev <- range(ev.plot)  #x@empirical
  col <- (if ("col" %in% dotnames) dots$col else
          default.image.par(NULL, NULL)$default.col)

  if (graphics$split.screen && legend) {
    scr <- split.screen( rbind(c(0,.85,0,1), c(.85,1,0,1)))
    screen(scr[1])
    par(mar=c(4,4,3,1))
  } else {
    scr <- NULL
    par(mar=c(4,4,3,1))
  }

  dots$type <- NULL
  do.call(plotmethod,
	  args=c(dots, list(x=coords[[1]], y=coords[[2]], z=ev.plot)))
					#      xlab=lab.names[1],
					#      ylab=lab.names[2],
					#      main="Variogram image plot"
  
  if (legend) {
    stop("das ist doch quark?!")
    
    
    if (graphics$split.screen) {
      screen(scr[2])
      par(mar=c(4,0,3,3))
      z.legend <- seq(range.ev[1], range.ev[2], length=length(col))
      image(y=z.legend, x=c(0,1), z=rbind(z.legend, z.legend), axes=F, col=col,
	    xlab="")
      axis(4)
      box()
    } else {
      my.legend(min(coords[[1]]), max(coords[[2]]), range(ev.plot),
		col=col, bg="white")
    }
  }
  
  if (graphics$close_screen) {
    close.screen(scr)
    scr <- NULL
  }
  return(invisible(scr))
}


RFplotEmpVariogram <- function(x, model = NULL, nmax.phi = NA, nmax.theta = NA,
			       nmax.T = NA,
			       plot.nbin = TRUE, plot.sd=FALSE, method = "ml",
			       variogram=TRUE,
			       boundaries = TRUE,
			       ...) {
  
  if (!variogram)
    stop("plot of estimated covariance functions not programmed yet.")
  
  graphics <- RFoptions()$graphics
  newx <- list()
  methodnames <- double()
  if (is(x, "RFfit") || is(x, "RF_fit")) {
    OP <- c("$", "@")[1 + is(x, "RFfit")]
    if(length(do.call(OP, list(x, "ev")))==0) stop("The fit does not contain an empirical variogram.")
    newx$autostart <- do.call(OP, list(x, "autostart"))
    newx$self <- do.call(OP, list(x, "self"))
    newx$plain <- do.call(OP, list(x, "plain"))
    newx$sqrt.nr <- do.call(OP, list(x, "sqrt.nr"))
    newx$sd.inv <- do.call(OP, list(x, "sd.inv"))
    newx$internal1 <- do.call(OP, list(x, "internal1"))
    newx$internal2 <- do.call(OP, list(x, "internal2"))
    newx$internal3 <- do.call(OP, list(x, "internal3"))
    newx$ml <- do.call(OP, list(x, "ml"))
    
    x <- do.call(OP, list(x, "ev"))
    if(is.null(method)) method <- if (length(newx$ml@name)>0) "ml" else "plain"
    
    methodidx<-names(newx[unlist(lapply(newx, function(x) is(x, CLASS_FIT)))])
    methodidx <- method %in% methodidx
    methodnames <- method[methodidx]
    nomethodnames <- method[!methodidx]
    
    if(!all(methodidx))
      warning( paste("The following method does not exist: ", nomethodnames))
  } else {
    if (!is(x, "RFempVariog") && !is(x, "RF_empVariog"))
      stop("method only for objects of class 'RFempVariog' or 'RFfit'")
    OP <- c("$", "@")[1 + is(x, "RFempVariog")]
  }
 
  ev <- do.call(OP, list(x, "empirical"))
  sd <- do.call(OP, list(x, "sd"))
  centers <- do.call(OP, list(x, "centers"))
  phi.centers <- do.call(OP, list(x, "phi.centers"))
  theta.centers <- do.call(OP, list(x, "theta.centers"))
  Time <- do.call(OP, list(x, "T"))
  n.bin <- do.call(OP, list(x, "n.bin"))
  coordunits <- do.call(OP, list(x, "coordunits"))
  empmethod <- do.call(OP, list(x, "method"))
  
  varnames <- if (is.matrix(ev)) dimnames(ev)[[2]][1] else names(ev)[1]
  if(empmethod > 2 && !is.null(model)) 
    warning("theoretical madogram model is not available at the moment")
  
  dots = list(...)
  dotnames <- names(dots)
  if (!("type" %in% dotnames)) dots$type <- "b" 
  cex <- if ("cex" %in% dotnames) dots$cex else .8
  par(cex=cex, xaxs="i")
  dots$cex <- NULL
  if (!("pch" %in% dotnames)) dots$pch <- 19
  
  if(!("xlim" %in% dotnames)) dots$xlim <- range(centers)
  ylim.not.in.dotnames <- !("ylim" %in% dotnames)
  xlab <- if ("xlab" %in% dotnames) dots$xlab else "distance"
  dots$xlab <- NULL
  ylab.ev <- if ("ylab" %in% dotnames) {
	       dots$ylab
	     } else if(empmethod == 2) {
	       "covariance"
	     } else if(empmethod < 2) {
	       "semivariance"
	     } else {
	       "madogram"
	     }
  dots$ylab <- NULL
  
  main0 <- if ("main" %in% dotnames) {
	     dots$main
	   } else if(empmethod == 2) {
	     "Covariance plot"
	   } else if(empmethod < 2) {
	     "Variogram plot"
	   } else {
	     "Madogram plot"
	   }
  dots$main <- NULL
  if (!is.null(main0)) oma.top <- 2 else oma.top <- 0
  
  
  has.sd <- !is.null(sd)
  
  if(!is.null(model)) {
    if (!is.list(model)) model <- list(model)
    if (!all(unlist(lapply(model, FUN=function(x) is(x, CLASS_CLIST)))))
      stop("model must be (a list of elements) of class 'CLASS_CLIST'")
    modelnames <-
    if(length(names(model)) > 0) names(model)
    else paste("model", 1:length(model))
    methodnames <- c(methodnames, modelnames)
    names(model) <- modelnames
    newx <- c(newx, model)
  }
  n.methods <- length(methodnames)
  
  
  n.phi <- min(nmax.phi, l.phi <- max(1,length(phi.centers)), na.rm=TRUE)
  n.theta <- min(nmax.theta, l.theta <- max(1,length(theta.centers)),
		 na.rm=TRUE)
  n.T <- min(nmax.T, l.T <- max(1,length(Time)), na.rm=TRUE)
  vdim <- dim(ev)
  vdim <- vdim[length(vdim)]
  
  if (n.phi > 6 || n.theta > 3 || n.T > 3)
    message("'If you feel the picture is overloaded, set the parameters 'nmax.phi', 'nmax.theta' and 'nmax.T'")
  
  halfphi.sector <- pi/(2*l.phi)
  halftheta.sector <- pi/(2*l.theta)
  phi.angles <- c(halfphi.sector, 0, -halfphi.sector) 
  theta.angles <- seq(-halftheta.sector, halftheta.sector, len=5) # len=odd!
  if (n.phi>1 && boundaries) {
    phi.angles <- phi.angles * 0.96
    theta.angles <- theta.angles * 0.96
  } 
  
  TandV <- (n.T > 1 && vdim > 1) && graphics$split_screen
  
  if (vdim>1 && length(varnames)==0)
    varnames <- paste("v", 1:vdim, sep="")
  
  range.nbin <- range(c(0, n.bin), na.rm=TRUE)
  ylim.nbin <- range.nbin * c(1,1.2)
  
  col.v <- col <-
    if ("col" %in% dotnames) rep(dots$col, len=n.phi) else 1:max(n.phi)
  dots$col <- NULL
  
  if (n.methods > 0){
    dotsRFfit <- dots
    dotsRFfit$type <- "l"
    dotsRFfit$lwd <- 2
    ltyRFfit.v <- 1:n.methods
    dotsRFfit$lty <- NULL
  }
  
  dir2vario <- function(dir.vec, x.eval, x.time, method.model, v1, v2){
    x.space <- as.matrix(x.eval) %*% t(dir.vec)            
    for (j in (1:4)) { ## orig 20
      if(empmethod == 2) { # covariance plot
	vario.vals <- try(RFcov(x = cbind(x.space, x.time),
				model = method.model,
				grid = FALSE,
				internal.examples_reduced=FALSE),
			  silent = TRUE)
      } else if(empmethod < 2) {
	vario.vals <- try(RFvariogram(x = cbind(x.space, x.time),
				      model = method.model,
				      grid = FALSE,
				      internal.examples_reduced=FALSE),
			  silent = TRUE)
      } else {
	break;
      }
      if(!is(vario.vals, "try-error")) {
	if (is.array(vario.vals)) {
	  return(vario.vals[, v1, v2])
	} else {
	  return(vario.vals)
	}
      }
      x.space <- cbind(x.space, 0)    
    }     
    return(NA)
  }
  
  oma.left <- 6
  Screens <- if (TandV) c(n.T, n.theta) else c(n.T * vdim * vdim, n.theta)
  ArrangeDevice(graphics, Screens)
  
  all.scr <- scr <- (if (!graphics$split_screen) NULL else
		    split.screen(if (TandV) c(vdim, vdim) else Screens))
  for (v1 in 1:vdim) {
    for (v2 in 1:vdim) {
      if (TandV) {
	scr <- split.screen(Screens, all.scr[v1 + (v2 - 1) * vdim])
	all.scr <- c(all.scr, scr)          
      }
      if (vdim == 1) {
        main <-
	if (!is.null(main0) && length(varnames)>0)
	  paste(main0, "for", varnames) else main0
      } else {
        main <-
	if (!TandV) main0
	else paste(main0, "for", varnames[v1], "vs.",  varnames[v2])
      }      
      if (ylim.not.in.dotnames)
        dots$ylim <- range(ev[,,,, v1, v2], na.rm=TRUE)
      for (iT in 1:n.T) {
        for (ith in 1:n.theta) {
          ## plot n.bin
          if (plot.nbin) {
            screen(scr[1])
            par(oma=c(4,oma.left,oma.top,0))
            scr2 <- split.screen(rbind(c(0,1,.2,1), c(0,1,0,.2)), screen=scr[1])
            all.scr <- c(all.scr, scr2)
            screen(scr2[2])
            par(mar=c(0,.5,0,.5))
            for (iph in 1:n.phi) {
              if (n.phi > 1) col <- col.v[iph]
              if (iph==1) {
                lab <- xylabs("bin centers", NULL, units=coordunits)
                plot(centers,
                     n.bin[ ,iph, ith, iT, v1, v2],
                     xlim=dots$xlim, ylim=ylim.nbin,
                     type=if (n.phi>1) "p" else "h",
                     col =if (n.phi>1) col else "darkgray", lwd=4,
                     pch=16, axes=FALSE, ylab="", xlab = lab$x)
                box()
                at <- seq(range.nbin[1], range.nbin[2], len=3)
                if (ith==1)
                  axis(2, las=1, at=at, labels=format(at, scient=-1, dig=2),
                       outer=TRUE)
                else  axis(2, las=1, at=at, labels=FALSE)
                
                if (iT==n.T && (n.T > 1 || (v1==vdim && v2==vdim))) axis(1)
                if (ith==1) title(ylab="n.bin", line=5, outer=TRUE, adj=0)
                
              } else {
                points(centers, n.bin[ ,iph, ith, iT, v1, v2],
                       type="p", col=col, pch=16)
              }
            }
            screen(scr2[1])
          } else {
            screen(scr[1])
            par(oma=c(4,oma.left,oma.top,0))
          }
          
          ## plot empirical
          ##if (ith==1) par(mar=c(0,6,1,1)) else par(mar=c(0,.5,1,.5))
          
          plotted.meth <- NULL  # needed for legend
          
          for (iph in 1:n.phi) {
            if (n.phi>1) col <- col.v[iph]
            if (iph==1) {
              do.call(graphics::plot,
                      c(dots, list(
				x=centers,
				y=ev[,iph,ith,iT,v1,v2],
                                        #ylim=ylim.ev, type=type, pch=19,
				col=col, axes=FALSE, ylab=""))
                      )
              box()
              axis(2, las=1, labels=(ith==1), outer=(ith==1))
              if (!plot.nbin) axis(1)          
              if (l.theta > 1 || l.T > 1 || vdim > 1) {
                L <- character(3)
                if (!TandV && vdim > 1) L[1] <- paste(varnames[c(v1,v2)],
                                                      collapse=":")
                if (l.T>1) L[2] <-
		paste("T=",signif(Time[iT], 3), " ", sep="")
                if (l.theta>1) L[3] <-
		paste(sep="", "theta=", signif(theta.centers[ith],3))
                legend("topleft", legend=paste(L[L!=""], collapse=","))
              }
              if (ith == 1) title(ylab=ylab.ev, line=5, outer=TRUE)
              if (has.sd && plot.sd)
                legend("topright", bty="n", #col=NA, lwd=NA,
                       legend="arrows represent +-1 sd intervals")
              
            } else {  ## iph > 1
              do.call(graphics::points,
                      c(dots, list(x=centers,
                                   y=ev[,iph,ith,iT,v1,v2],
                                   col=col))
                      )   #type=type,  pch=19
	    } ## for iph 
	    
	    k <- 1  
	    if (n.methods > 0) {          
	      for(i in 1:n.methods) {        
		method.model <- newx[[methodnames[i]]]
		
		if(length(method.model@name) == 0){
		  warning("The method '", methodnames[i], "' cannot be fitted.")
                  next
		}
		if (!is.null(phi.centers)) {
		  x.radial <- cbind(cos(phi.centers[iph]+phi.angles),
                                    sin(phi.centers[iph]+phi.angles))
		  if (!is.null(phi.centers))
		    x.radial <-
		  cbind(x.radial,
			cos(rep(theta.angles, each = length(phi.angles))
			    + theta.centers[ith]))
		  
		} else x.radial <- matrix(1, nrow=1, ncol=1)
		
		x.time <- NULL
		if (!is.null(Time)) x.time <- Time[iT]
		x.eval <- seq(from = max(dotsRFfit$xlim[1],1e-3),
			      to = dotsRFfit$xlim[2], len = 150) 
		
		## sehr genau Abschaetzung, indem mehrere (3)
		## Winkel angeschaut werden und dann der mittlere
		## Wert der Variogramme angeschaut wird, da ja auch
		## das emp. Variogramm ueber ein Winkelintervall gemittelt wird
		dummy.vals <- as.matrix(apply(x.radial, 1, FUN = dir2vario,
					      x.eval=x.eval, x.time=x.time,
					      method.model=method.model,
					      v1=v1, v2=v2))
		
		
		## Print(i, x.eval, dummy.vals);print.RMmodelFit(method.model)
                
		if(length(dummy.vals) > 0 && !all(is.na(dummy.vals))){
					# Print(n.phi, boundaries)
		  if (n.phi>1 && boundaries) {
		    do.call(graphics::matplot,
                            args=c(dotsRFfit, list(x=x.eval,
				   y= if (ncol(dummy.vals) == 1) dummy.vals
				   else t(apply(dummy.vals, 1, range)),
				   add = TRUE, col=col, lty = 3)))
		  } else {
		    do.call(graphics::points,
			    args=c(dotsRFfit, list( x=x.eval,
					 y = rowMeans(dummy.vals)  , 
					 col=col, lty = ltyRFfit.v[k])))
		  }             
		  k <- k+1
		  if(iph == 1) plotted.meth <- c(plotted.meth, methodnames[i])
		}	      
	      } ## for i in n.methods
	    } ## if n.methods > 0	  
	  
	    if (has.sd && plot.sd) {
	      sdnot0 <-  sd[ ,iph, ith, iT] != 0
	      arrows(centers[sdnot0],
		     ev[sdnot0 ,iph, ith,iT] - sd[sdnot0,iph,ith,iT],
		     centers[sdnot0],
		     ev[sdnot0 ,iph, ith,iT] + sd[sdnot0,iph,ith,iT],
		     code=2, angle=90, length=0.05, col=col)
	    }
	    
	  } # iph in nphi
	  
	
	  pifactor <- signif((phi.centers[1:n.phi]) / pi, 2)
	  
	  len.mnames <- length(plotted.meth)
	  string.emp <- "empirical"
	  if(len.mnames > 0) {
	    labels <-
	    if (n.phi>1) paste("\"phi=\",", rep(pifactor, each = len.mnames+1),
			       ", pi, ", "\", ") else "\""
	    labels <- paste("c(",
			    paste("expression(paste(",labels,
				  rep(c(string.emp, plotted.meth), l.phi),
				  "\"", "))", collapse=","), ")")
	    labels <- eval(parse(text=labels))
	  } else {
	    labels <- if (n.phi > 1) 
			eval(parse(text=paste(
			  "c(",
			  paste("expression(paste(\"phi=\",", pifactor,", pi))",
				collapse=","),
			  ")")))
	  }
	  
	  if (l.phi > 1 || len.mnames > 0)# && iT==1 && ith==1 # auskommentiert auf Sebs wunsch        
	legend("bottomright", col=rep(col.v, each = len.mnames+1),
	       lwd=1, pch=rep(c(19,rep(NA, len.mnames)), l.phi), bty="n",
	       legend=labels, lty =
	       rep(c(1, if(len.mnames==0) NULL else ltyRFfit.v[1:len.mnames]),
		   l.phi))
	
	  scr <- scr[-1]
	} # n.theta
      } # T
    } # vdim 1
  } # vdim 2
  dots$type <- NULL
  
  if (!is.null(main))
    do.call(graphics::title, args=c(dots, main=main, outer=TRUE))
  if (!is.null(xlab))
    do.call(graphics::title, args=c(dots, xlab=xlab, outer=TRUE))
  
  if (graphics$close_screen) {
    close.screen(all.scr)
    all.scr <- NULL
  }
  
  return(invisible(all.scr))
}

setMethod(f="plot", signature(x="RFempVariog", y="missing"),
          function(x, y, ...) RFplotEmpVariogram(x, ...))