https://github.com/cran/RandomFields
Tip revision: 683e381531c37e8e7224edd899422f119d926418 authored by Martin Schlather on 21 January 2014, 00:00:00 UTC
version 3.0.10
version 3.0.10
Tip revision: 683e381
rf.R
# source("rf.R")
# see getNset.R for the variable .methods
## sudo apt-get install tcl8.5-dev
## sudo apt-get install tk8.5-dev
## see R Installation and Administration
##./configure --with-tcl-config=/usr/lib/tcl8.5/tclConfig.sh --with-tk-config=/usr/lib/tk8.5/tkConfig.sh
# sudo tar xvf ~/TMP/bwidget-1.9.5.tar
# readline <- function(...) return("")
.onAttach <- function (lib, pkg) {
#packageStartupMessage("\n")
}
#.onUnload <- function(lib, pkg){
# }
#Implementierung von Cox & Isham's non-separable model
#individuelle Praeferenzliste:
# 1) erste praeferenz
# 2) if # pkte < 10 / 50 -> Gauss??
# 4) if nugget/variance > 0.01 -> CE (and not spectral, e.g.)
# 3) if C(maxh)/C(0) > 0.05 -> TBM else CE
.onUnload <- function(lib, pkg){
RFoptions(storing=FALSE) ## delete everything
}
RFdistr <- function(model, x, q, p, n, dim=1, ...) {
## note: * if x is a matrix, the points are expected to be given row-wise
## * if not anisotropic Covariance expects distances as values of x
##
## here, in contrast to Covariance, nonstatCovMatrix needs only x
RFoptOld <- internal.rfoptions(..., RELAX=isFormulaModel(model))
on.exit(RFoptions(LIST=RFoptOld[[1]]))
model<- list("Distr", PrepareModel2(model), dim=as.integer(dim));
if (!missing(x)) {
model$x <- if (is.matrix(x)) t(x) else x
}
if (!missing(q)) {
model$q <- if (is.matrix(q)) t(q) else q
}
if (!missing(p)) {
model$p <- if (is.matrix(p)) t(p) else p
}
if (!missing(n)) {
model$n <- n
}
rfInitSimulation(model=model, x=matrix(0, ncol=dim, nrow=1),
y=NULL, z=NULL, T=NULL, grid=TRUE,
reg = MODEL.USER,
dosimulate=FALSE)
return(.Call("EvaluateModel", double(0), as.integer(MODEL.USER),
PACKAGE="RandomFields"))
}
RFddistr <- function(model, x, dim=1,...) {
if (hasArg("q") || hasArg("p") || hasArg("n")) stop("unknown argument(s)");
RFdistr(model=model, x=x, dim=dim,...)
}
RFpdistr <- function(model, q, dim=1,...) {
if (hasArg("x") || hasArg("p") || hasArg("n")) stop("unknown argument(s)");
RFdistr(model=model, q=q, dim=dim,...)
}
RFqdistr <- function(model, p, dim=1,...){
if (hasArg("x") || hasArg("q") || hasArg("n")) stop("unknown argument(s)");
RFdistr(model=model, p=p, dim=dim,...)
}
RFrdistr <- function(model, n, dim=1,...) {
if (hasArg("x") || hasArg("q") || hasArg("p")) stop("unknown argument(s)");
RFdistr(model=model, n=n, dim=dim,...)
}
RFcov <- function(model, x, y = NULL, z = NULL, T=NULL, grid,
distances, dim,
# dim = ifelse(is.matrix(x), ncol(x), 1),
fctcall=c("Cov", "CovMatrix", "Variogram",
"Pseudovariogram", "Fctn"), ...) {
## note: * if x is a matrix, the points are expected to be given row-wise
## * if not anisotropic Covariance expects distances as values of x
##
## here, in contrast to Covariance, nonstatCovMatrix needs only x
RFoptOld <- internal.rfoptions(xyz_notation=2*(!is.null(y) && !is.matrix(y)),
..., RELAX=isFormulaModel(model))
on.exit(RFoptions(LIST=RFoptOld[[1]]))
fctcall <- match.arg(fctcall)
p <- list(fctcall, PrepareModel2(model));
rfInitSimulation(model=p,
x=x, y=y, z=z, T=T, grid=grid,
distances=distances, spdim=dim, reg = MODEL.USER,
dosimulate=FALSE)
result <- .Call("EvaluateModel", double(0), as.integer(MODEL.USER),
PACKAGE="RandomFields")
return(result)
}
RFcovmatrix <- function(...) {
RFcov(..., fctcall="CovMatrix")
}
RFvariogram <- function (model, x, y=NULL, ...)
RFcov(model=model, x=x, y, fctcall="Variogram", ...)
RFpseudovariogram <- function(model, x, y=NULL, ...)
RFcov(model=model, x=x, y, fctcall="Pseudovariogram", ...)
RFfctn <- function(...) {
RFcov(..., fctcall="Fctn")
}
######################################################################
######################################################################
prepare4RFspDataFrame <- function(model.orig=NULL,
info, x, y, z, T, grid, data, RFopt) {
vdim <- info$vdim
locinfo <- info$loc
## names
if (missing(data))
names <- rfPrepareNames(model=model.orig, locinfo=locinfo)
else
names <- rfPrepareData(model.orig, x, y, z, T, grid, data,
cathegories=NULL, names.only=TRUE)
if (!is.null(names$variab.names)) {
if (vdim != length(names$variab.names))
stop(paste("you passed a formula for 'model' with left-hand side '",
paste(names$variab.names, collapse=","),
"', but vdim of the model equals ", vdim, sep=""))
}
coord.names.incl.T <- names$coord.names
## coords or GridTopology
if (locinfo$grid) { ## grid=TRUE
coords <- NULL
xgr <- cbind(locinfo$xgr, locinfo$T)
colnames(xgr) <- coord.names.incl.T
xgr[is.na(xgr)] <- 0
gridTopology <- GridTopology(xgr[1, ], xgr[2, ], xgr[3, ])
} else { ## grid == FALSE
gridTopology <- NULL
# cbind of locations from x-matrix and T (if given)
coords <- apply(t(locinfo$x), 2, rep,
times=(locinfo$totpts/locinfo$spatialtotpts))
if (locinfo$Time) {
T <- locinfo$T
coord <- cbind(coord, rep(seq(T[1], by=T[2], len=T[3]),
each=locinfo$spatialtotpts))
}
colnames(coords) <- coord.names.incl.T
}
# Print(RFopt)
if (RFopt$general$printlevel>0 && RFopt$warn$newstyle) {
RFoptions(warn.newstyle = FALSE)
message("New output format of RFsimulate: object of class 'RFsp';\n",
"for conventional array format use 'RFoptions(spConform=FALSE)'.")
}
return(list(coords=coords, gridTopology=gridTopology, vdim=vdim, names=names))
}
rfDoSimulate <- function(n = 1, reg, spConform) {
stopifnot(length(n) == 1, n>0, is.finite(n))
RFopt <- RFoptions()
if (missing(spConform)) spConform <- RFopt$general$spConform
if (RFopt$gauss$paired && (n %% 2 != 0))
stop("if paired, then n must be an even number")
info <- RFgetModelInfo(RFopt$general$register)
# Print(info); str(info$loc); xxx
len <- info$loc$len
vdim <- info$vdim
total <- info$loc$totpts
if (is.null(total) || total <= 0)
stop("register ", RFopt$general$register, " does not look initialized")
error <- integer(1)
# result <- double(total * n * vdim)
# Print(info, total, vdim, len, info$simu$distr, n); xxx
result <- .Call("EvaluateModel", as.double(n), as.integer(reg), #userdefined,
PACKAGE="RandomFields", DUP=FALSE)
if (!spConform) return(result)
prep <- prepare4RFspDataFrame(model.orig=NULL, info=info, RFopt=RFopt)
attributes(result)$variab.names <- prep$names$variab.names
res2 <- conventional2RFspDataFrame(result, coords=prep$coords,
gridTopology=prep$gridTopology,
n=n, vdim=prep$vdim)
return(res2)
}
rfInitSimulation <- function(model, x, y = NULL, z = NULL, T=NULL, grid=FALSE,
distances, spdim,
reg, dosimulate=TRUE) {
stopifnot(xor(missing(x) || length(x)==0,
missing(distances) || length(distances)==0))
RFopt <- RFoptions()
if (!is.na(RFopt$general$seed)) set.seed(RFopt$general$seed)
neu <- CheckXT(x, y, z, T, grid=grid, distances=distances, spdim=spdim,
y.ok=!dosimulate)
ts.xdim <- as.integer(ncol(neu$x) + neu$Time)
if (missing(reg)) stop("rfInitSimulation may not be called by the user")
# userdefined <- GetParameterModelUser(model)
# Print(userdefined)
vdim <- .Call("Init",
as.integer(reg),
model,
if (neu$grid) neu$x else t(neu$x),
if (is.null(neu$y)) {
double(0)
} else if (neu$grid) neu$y else t(neu$y),
as.double(neu$T),
as.integer(neu$spacedim),
as.logical(neu$grid),
as.logical(neu$distances),
as.logical(neu$Time),
NAOK=TRUE,
PACKAGE="RandomFields")
}
RFsimulate <- function (model, x, y = NULL, z = NULL, T = NULL, grid,
data, distances, dim, err.model,
n = 1, ...) {
### preparations #########################################################
if (!missing(distances) && length(distances) > 0)
RFoptOld <- internal.rfoptions(xyz_notation=!is.null(y),
expected_number_simu=n, ...,
general.spConform = FALSE,
RELAX=isFormulaModel(model))
else
RFoptOld <- internal.rfoptions(xyz_notation=!is.null(y),
expected_number_simu=n, ...,
RELAX=isFormulaModel(model))
on.exit(RFoptions(LIST=RFoptOld[[1]]))
RFopt <- RFoptOld[[2]]
cond.simu <- !missing(data) && !is.null(data)
reg <- if (cond.simu) RFopt$general$condregister else RFopt$general$register
### simulate from stored register ########################################
mc <- as.list(match.call(expand.dots=FALSE))
if (length(mc)==1 ||
length(mc)==2 && !is.null(mc$n) ||
length(mc)==3 && !is.null(mc$n) && "..." %in% names(mc)) {
if (cond.simu) {
stop("repeated simulation of conditional simulation not programmed yet")
} else {
# userdefined <- GetParameterModelUser(PrepareModel2(model, ...))
return(rfDoSimulate(n=n, reg=reg, spConform=RFopt$general$spConform
#userdefined=userdefined
))
}
}
### preparations #########################################################
stopifnot(!missing(model) && !is.null(model))
model.orig <- model
model <- PrepareModel2(model, ...)
err.model <- if (missing(err.model)) NULL else PrepareModel2(err.model, ...)
### conditional simulation ###############################################
if (cond.simu) {
#Print(RFoptions()$general)
if (missing.x <- missing(x)) {
if (is(data, "RFsp")) {
if (!missing(dim))
stop("'dim' may not be given when 'data' is an RFsp object")
x.tmp <- coordinates(data)
} else { ## data not RFsp
start <- data.starting.col(data, halt=FALSE)
if (missing(dim)) {
if (ncol(data) == 2) dim <- 1
else if (start<0) stop("Please give 'dim' explicitely since x is missing and data is not an RFsp object.")
else dim <- start - 1
} else {
if (start < 0) colnames(data)[dim+1] <- "data"
else if (start - 1 != dim) stop(paste("value of 'dim' does not match the column name convention of the data: the name 'data' is found in position", start, "which does not equal 'dim'+1 (dim=", dim, ")"))
}
x.tmp <- data[, 1:dim]
}
} else {
x.tmp <- x
}
p <- list("Simulate", model);
rfInitSimulation(model=p, x=x.tmp, y=y, z=z, T=T, grid=grid,
distances=distances, spdim=dim,
reg=reg, dosimulate=-1)
info <- RFgetModelInfo(reg)
grid <- info$loc$grid
# xxxx
# Print(info$type, info)
res <- switch(info$role,
Gauss =
rfCondGauss(model=model.orig, x=x, y=y, z=z, T=T,
grid=grid, n=n, data=data,
err.model=err.model, ...),
BrownResnick = stop(paste("conditional simulation for ",
"BRprocesses not programmed yet.")),
Smith = stop(paste("conditional simulation for Smith ",
"processes not programmed yet.")),
Schlather = stop(paste("conditional simulation for ",
"Schlather processes not programmed yet.")),
Poisson = stop(paste("conditional simulation for Poisson ",
"processes not programmed yet.")),
stop(paste(info$role, "not programmed yet"))
)
## if (missing.x==TRUE), result can already be returned here
## but rfPrepareData has to be called again to get the full vector
## of given locations and the given data
if (missing.x) {
all <- rfPrepareData(model.orig=model.orig,
x, y, z, T, grid, data, cathegories=NULL)
coords <- all$x
colnames(coords) <- all$coord.names
tmp.res <- conventional2RFspDataFrame(data=res, coords=coords,
gridTopology=NULL,
n=n, vdim=info$vdim)
res <- tmp.res@data
tmp.res@coords <- all$fullgiven
tmp.res@data <- as.data.frame(all$fulldata)
tmp.res@data[all$data.na, ] <- res
if (is(tmp.res, "RFspatialPointsDataFrame"))
try(tmp.res <- as(tmp.res, "RFspatialGridDataFrame"), silent=TRUE)
if (is(tmp.res, "RFpointsDataFrame"))
try(tmp.res <- as(tmp.res, "RFgridDataFrame"), silent=TRUE)
return(if (RFopt$general$spConform) tmp.res else tmp.res@data)
}
## end of cond simu
} else {
## unconditional simulation #############################################
if(!is.null(err.model))
warning("model for measurement error is ignored in unconditional simulation")
stat <- (!missing(distances) && length(distances) > 0) || !is.matrix(x) || is.null(y)
p <- list("Simulate", model);
#Print(unix.time(
rfInitSimulation(model=p, x=x, y=y, z=z, T=T, grid=grid,
distances=distances, spdim=dim,
reg=reg)#)) #$userdefined
info <- RFgetModelInfo(reg)
#Print(info$simu)
#
# Print(info);xxxxxx
## perform simulation ###################################################
#Print(unix.time(
if (n < 1) return(NULL)
res <- rfDoSimulate(n=n, reg=reg, spConform=FALSE)#))
} # end of uncond simu
## output: either conventional or RFsp #################################
if (!RFopt$general$spConform)
return(res)
prep <- prepare4RFspDataFrame(model.orig, info, x, y, z, T,
grid, data, RFopt)
attributes(res)$variab.names <- prep$names$variab.names
# Print(res, coords=prep$coords,
# gridTopology=prep$gridTopology,
# n=n, vdim=prep$vdim)
res2 <- conventional2RFspDataFrame(res, coords=prep$coords,
gridTopology=prep$gridTopology,
n=n, vdim=prep$vdim)
if (is.raster(x)) {
res2 <- raster::raster(res2)
projection(res2) <- projection(x)
}
return(res2)
}
### ist keine Methode im engeren Sinne. Habe ich aus Methods-RFsp.R
### rausgenommen, da bei jeglicher Aenderung in Methods-RFsp.R ich
### komplett neu installieren muss. Bei rf.R muss ich es nicht.
conventional2RFspDataFrame <-
function(data, coords=NULL, gridTopology=NULL, n=1, vdim=1)
{
if (!xor(is.null(coords), is.null(gridTopology)))
stop("one and only one of 'coords' and 'gridTopology' must be NULL")
variab.names <- attributes(data)$variab.names
## may be NULL, if called from 'RFsimulate', the left hand side of model, if
## model is a formula, is passed to 'variab.names'
attributes(data)$variab.names <- NULL
## grid case
if (length(coords) == 0) {# war is.null(coords) -- erfasst coords=list() nicht
grid <- convert2GridTopology(gridTopology)
timespacedim <- length(grid@cells.dim)
## naechste Zeile eingefuegt !! und (Martin 30.6.13) wieder
## auskommentiert. s. Bsp in 'RFsimulate'
## if (!is.null(dim(data)) && all(dim(data)[-1]==1)) data <- as.vector(data)
if (is.null(dim(data))) {
stopifnot(1 == timespacedim + (n > 1) + (vdim > 1))
} else {
if (length(dim(data)) != timespacedim + (n>1) + (vdim > 1)){
stop(paste(length(dim(data)),
"= length(dim(data)) != timespacedim + (n>1) + (vdim>1) =",
timespacedim, '+', (n>1), '+', (vdim > 1)))
}
}
if (vdim>1){
## new order of dimensions: space-time-dims, vdim, n
perm <- c( 1+(1:timespacedim), 1, if (n>1) timespacedim+2 else NULL)
data <- aperm(data, perm=perm)
}
if (timespacedim==1)
call <- "RFgridDataFrame"
else {
data <- reflection(data, 2, drop=FALSE)
call <- "RFspatialGridDataFrame"
}
}
## coords case
if (is.null(gridTopology)){
if (vdim>1){
n.dims <- length(dim(data))
perm <- c(2:(n.dims - (n>1)), 1, if (n>1) n.dims else NULL)
data <- aperm(data, perm=perm)
}
if (is.null(dim(coords)) || ncol(coords)==1)
call <- "RFpointsDataFrame"
else call <- "RFspatialPointsDataFrame"
}
## in both cases:
dim(data) <- NULL
data <- as.data.frame(matrix(data, ncol=n*vdim))
if (is.null(variab.names))
variab.names <- paste("variable", 1:vdim, sep="")
if (length(variab.names) == n*vdim)
names(data) <- variab.names
else
if (length(variab.names) == vdim)
names(data) <- paste(rep(variab.names, times=n),
if (n>1) ".n", if (n>1) rep(1:n, each=vdim),sep="")
else names(data) <- NULL
## Print(variab.names, names(data), coords)
if (is.null(coords)){
do.call(call, args=list(data=data, grid=grid,
RFparams=list(n=n, vdim=vdim)))
} else {
do.call(call, args=list(data=data, coords=coords,
RFparams=list(n=n, vdim=vdim)))
}
}
list2RMmodelExt <- function(x, isRFsp=FALSE,
coords, gridTopology, data.RFparams) {
stopifnot(is.list(x),
all(c("model", "trend", "ml.value", "residuals") %in% names(x)))
if (isRFsp) {
stopifnot(!missing(coords) &&
!missing(gridTopology) &&
!missing(data.RFparams))
## convert residuals to RFsp class
err <-
try({
lres <- length(x$residuals)
if (lres > 0) {
for (i in 1:lres) {
gT <- if (length(gridTopology) < i) NULL else gridTopology[[i]]
co <- if (length(coords)<i) NULL else coords[[i]]
if (!is.null(x$residuals[[i]]))
x$residuals[[i]]<-
conventional2RFspDataFrame(data=x$residuals[[i]],
coords=co,
gridTopology=gT,
n=data.RFparams[[i]]$n,
vdim=data.RFparams[[i]]$vdim)
}
}
}
, silent=!TRUE)
if(class(err)=="try-error")
warning(paste("residuals could not be coerced to class'RFsp';",
err))
}
return(new(ZF_MODELEXT,
list2RMmodel(x$model),
trend = list2RMmodel(x$trend),
likelihood = x$ml.value,
AICc = x$aicc,
residuals = x$residuals))
}