## Authors ## Martin Schlather, schlather@math.uni-mannheim.de ## ## ## Copyright (C) 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. ## source("modelling.R") predictGauss <- function(Reg, Reg.predict= if (missing(model)) Reg else RFopt$register$predict_register, model, x, y = NULL, z = NULL, T=NULL, grid=NULL, data, distances, dim, kriging_variance, ...) { relax <- !missing(model) && !is.null(model) && is(model, "formula") RFoptOld <- internal.rfoptions(..., RELAX=relax) on.exit(RFoptions(LIST=RFoptOld[[1]])) RFopt <- RFoptOld[[2]] if (missing(model) || is.null(model)) { ## no measurement errors if (Reg != Reg.predict) stop("'Reg.predict' may not be given.") model <- list("predict", register=Reg) } else { ## here, 'Reg' contains a measurement structure, and model is without it if (Reg == Reg.predict) stop("'Reg.predict' must be different from 'Reg'") model <- list("predict", PrepareModel2(model), register=Reg) ## to do : warum nachfolgendes? hat doch keine Wirkung!? splittingC <- function(model, preceding, factor) { const <- sapply(model[-1], function(m) (is.numeric(m) && !is.na(m)) || (m[[1]] == R_CONST && !is.na(m[[2]])) ) if (all(const)) { model <- c(model[1], if (preceding > 0) rep(0, preceding), model[-1]) return(list(SYMBOL_MULT, model, if (!missing(factor)) list(factor))) } for (i in 2:length(model)) { vdim <- preceding + (if (i==2) 0 else GetDimension(model[[i-1]])) m <- ReplaceC(model[[i]]) L <- GetDimension(m) model[[i]] <- setvector(m, preceding = vdim, len = L, factor=factor) } model[[1]] <- SYMBOL_PLUS names(model) <- NULL L <- GetDimension(model[[length(model)]]) model <- SetDimension(model, L) } } rfInit(model=model, x=x, y=y, z=z, T=T, grid=grid, distances=distances, dim=dim, reg = Reg.predict, dosimulate=FALSE) .Call(C_EvaluateModel, double(0), as.integer(Reg.predict)) } GetModelEffects <- function(Z) { .Call(C_SetAndGetModelLikelihood, MODEL_AUX, list("RFloglikelihood", data = Z$data, Z$model), trafo.to.C_UnifyXT(Z$coord), original)$effect } ModelAbbreviations <- function(model){ L <- length(model) N <- character(L) ##Print(model) for (mm in 1:L) { m <- unlist(model[[mm]]) dollar <- m==DOLLAR[1] | m==DOLLAR[2] idx <- which( names(m) != "" & !dollar) ## names(m):names of parameters ## Print(m, mm, dollar, idx, names(m)) for (s in which(dollar)) { ## Print(s, which(idx > s)) z <- idx[min(which(idx > s))] ## to which model belongs the dollar param? m[s] <- m[z] ## put the name of subsequent model in position of dollar pos m <- m[-z] ## delete subseqeuent model } idx <- names(m) != "" m[idx] <- paste(names(m)[idx], m[idx], sep="=") N[mm] <- paste(abbreviate(m), collapse =".") ## Print(idx, m, N) } return(N) } ModelParts <- function(model, effects, complete) { ## model immer schon aufbrtt model <- PrepareModel2(model) ##, params=params, ...) model <- if ((model[[1]] %in% SYMBOL_PLUS)) model[-1] else list(model) if (complete) { return(model[effects >= RandomEffect]) } else { err <- model[effects == ErrorEffect] if (length(err) == 1) err <- err[[1]] else err <- c(SYMBOL_PLUS, err) m <- model[effects == RandomEffect] if (length(m) == 1) m <- m[[1]] else m <- c(SYMBOL_PLUS, m) return(list(model=m, err.model=err)) } } xRFranef <- function(fit, method="ml", OP='@') { ## see ranef random effect #Z <- do.call(OP, list(object, "Z")) # model <- fit[method] # parts <- ModelParts(model, effects=GetModelEffects(Z), complete=TRUE) # ans <- RFinterpolate(model, x = Z$coord, data = Z$data, given = Z$coord, # err.model=parts) # return(ans) } FinImputIntern <- function(data, simu, coords, coordnames, data.col, vdim, spConform, fillall=FALSE) { n <- length(data) / (vdim * coords$restotal) #Print(data, all, tail(all$simu), spConform); if (is(data, "RFsp")) { if (spConform) { data@data[ , ] <- as.vector(simu) return(data) } else { values <- as.matrix(data@data) values[is.na(values) | fillall] <- simu return(cbind(coordinates(data), values)) } } else { ## not RFsp #Print("for testing") if (coords$grid) { ## to do stop("not programmed yet") } else { ## coords <- all$x colnames(coords$x) <- coordnames values <- data[, data.col] values[is.na(values) | fillall] <- simu if (!spConform) return(cbind(coords$x, values)) tmp.all <- conventional2RFspDataFrame(data=values, coords=coords$x, gridTopology=NULL, n=n, vdim=vdim, vdim_close_together=FALSE) if (is(tmp.all, "RFspatialPointsDataFrame")) try(tmp.all <- as(tmp.all, "RFspatialGridDataFrame"), silent=TRUE) if (is(tmp.all, "RFpointsDataFrame")) try(tmp.all <- as(tmp.all, "RFgridDataFrame"), silent=TRUE) } return(tmp.all) } } FinishImputing <- function(data, simu, Z, spConform, fillall) { ## to do: grid if (is.list(data)) { for (i in 1:length(data)) data[[i]] <- FinImputIntern(data=data[[i]], simu=simu[[i]], coords=Z$coord[[i]], coordnames=Z$coordnames, data.col=Z$data.col, vdim=Z$vdim, spConform = spConform, fillall=fillall) return(data) } return(FinImputIntern(data=data[[1]], simu=simu, coords=Z$coord[[1]], coordnames=Z$coordnames, data.col=Z$data.col, vdim=Z$vdim, spConform=spConform, fillall=fillall)) } ExpandGrid <- function(x) { #### ACHTUNG! ZWINGENDE REIHENFOLGE if (x$grid) { # 0 x$x <- as.matrix(do.call(expand.grid, lapply(apply(cbind(x$x, x$T), 2, function(x) list(seq(x[1],by=x[2],length.out=x[3]))), function(x) x[[1]]))) } else if (x$has.time.comp) { dim.x <- if (is.vector(x$x)) c(length(x$x), 1) else dim(x$x) x$x <- cbind(matrix(rep(t(x$x), times=x$T[3]), ncol=dim.x[2], byrow=FALSE), rep(seq(x$T[1], by=x$T[2], length.out=x$T[3]), each=dim.x[1])) } if (length(x$y) > 0) stop("no expansion within a kernel definition") # x$y <- double(0) #1 x$T <- double(0) #2 x$grid <- FALSE #3 # x$spatialdim <- ncol(x$x) #4 x$has.time.comp <- FALSE #5 # x$di st.given <- FALSE #6 x$restotal <- nrow(x$x) #7 x$l <- x$restotal #8 return(x) } rfPrepareData <- function(model, x, y=NULL, z=NULL, T=NULL, distances=NULL, dim, grid, data, given=NULL, params, RFopt, reg, err.model = NULL, err.params, ...) { ## NOTE: if err.model is a list of list, the err.model is ## not added to the krige model, but both are taken as they are ## internal behaviour to get the random effects with the ## existing algorithm!! if (is(model, "RFfit")) { message("To continue with the output of 'RFfit' better use 'predict' or give the components separately.") model <- model["ml"] } if (!missing(distances) && length(distances)>0) stop("option distances not programmed yet.") # Print(model=model, data=data, given=given, T, ...) missing.x <- missing(x) || length(x) == 0 imputing <- missing.x && length(distances) == 0 if (length(given) == 0) { ## so either within the data or the same the x-values Z <- UnifyData(model=model, data=data, RFopt=RFopt, params=params, ...) if (Z$matrix.indep.of.x.assumed) { if (missing.x) stop("coordinates cannot be detected") Z <- UnifyData(model=model, x=x, y=y, z=z, T=T, RFopt=RFopt, distances=distances, dim=dim, grid=grid, data=data, params=params, ...) } } else { Z <- UnifyData(model=model, data=data, x=given, RFopt=RFopt, params=params,...) } model <- krige <- Z$model if (length(Z$data) != 1) stop("exactly one data set must be given.") dim_data <- base::dim(Z$data[[1]]) Z$data[[1]] <- as.double(Z$data[[1]]) repet <- Z$repetitions new.dim_data <- c(prod(dim_data) / repet, repet) base::dim(Z$data[[1]]) <- new.dim_data data.na <- is.na(Z$data[[1]]) data.na.var <- rowSums(data.na) base::dim(Z$data[[1]]) <- dim_data base::dim(data.na.var) <- c(length(data.na.var) / Z$vdim , Z$vdim) data.na.loc <- rowSums(data.na.var > 0) > 0 any.data.na <- any(data.na.loc) split <- any(data.na.var > 0 & data.na.var != repet) if (any.data.na && Z$coord[[1]]$dist.given) stop("missing values not programmed yet for given distancs") if (imputing) { if (Z$vdim > 1) stop("imputing does not work in the multivariate case") if (repet == 1) { if (RFopt$krige$fillall || !any.data.na) { data.na <- rep(TRUE, length(data.na)) ## nur ## um Daten im Endergebnis einzutragen new <- Z$coord[[1]] } else { new <- ExpandGrid(Z$coord[[1]]) new$x <- new$x[data.na.loc, , drop=FALSE] } } else new <- NULL } else { ## needed in soilRd, in condsimu!! new <- UnifyXT(x, y, z, T, grid=grid, distances=distances, dim=dim) ## new <- Z$coord[[1]] ## why was this set??? ## Print(new, UnifyXT(x, y, z, T, grid=grid, distances=distances, dim=dim)); kkkff if (Z$tsdim != new$spatialdim + new$has.time.comp) stop("coodinate dimension of locations with and without data, ", "respectively, do not match.") } na_rm_lines <- FALSE if (any.data.na) { na_rm_lines <- RFopt$general$na_rm_lines if (na_rm_lines && (!imputing || repet==1)) { Z$data[[1]] <- Z$data[[1]][!data.na.loc, , drop=FALSE] Z$coord[[1]] <- ExpandGrid(Z$coord[[1]]) Z$coord[[1]]$x <- Z$coord[[1]]$x[!data.na.loc, , drop=FALSE] } else if (split) { data <- list() for (i in 1:repet) { dim(Z$data[[1]]) <- c(length(Z$data)[[1]] / (Z$vdim*repet), Z$vdim, repet) data[[i]] <- Z$data[[1]][ , , i, drop=FALSE] } Z$data <- data } } Z$na_rm_lines <- na_rm_lines ## krige enthaelt err.model ## model ohne err.model if (length(err.model) == 0) { model <- list(model) } else if (is(err.model, "RMmodel") || (is.list(err.model) && !is.list(err.model[[1]]))) { linpart <- RFlinearpart(model=err.model, params=err.params, new$x, set=1, ...) if (length(linpart$X) > 0 || any(linpart$Y != 0)) stop("a trend is not allowed for the error model.") krige <- list(SYMBOL_PLUS, PrepareModel2(err.model, params=err.params, ...), model) model <- list(model) } else if (is.list(err.model) && !is.character(err.model[[1]])) { model <- list() if (missing(err.params)) err.params <- list() for (i in 1:length(err.model)) { linpart <- RFlinearpart(model=err.model[[i]], new$x, set=1) if (length(linpart$X) > 0 || any(linpart$Y != 0)) stop("a trend is not allowed for the error model.") model[[i]] <- PrepareModel2(err.model[[i]], params=err.params[[i]], ...) } } else if (is.vector(err.model) && length(err.model)==1 && is.na(err.model) && missing(err.params)) { model <- list(ModelParts(model, effects=GetModelEffects(Z), complete = FALSE)$model) } else stop("The argument 'err.model' cannot be interpreted. or 'err.params' is given where it should not") return(list(Z=Z, X=new, krige = krige, # model for the matrix to be inverted in SK model=model, # model for the vector in SK imputing=imputing, data.na = if (imputing) data.na)) } RFinterpolate <- function(model, x, y=NULL, z=NULL, T=NULL, grid=NULL, distances, dim, data, given=NULL, params, err.model=NULL, err.params, ignore.trend=FALSE, ...) { if (is(model, "RFfit")) { message("To continue with the output of 'RFfit' use 'predict' or give the components separately") model <- model["ml"] } if (!missing(distances) && length(distances) > 0) stop("'distances' not programmed yet.") opt <- list(...) i <- pmatch(names(opt), c("MARGIN")) opt <- opt[is.na(i)] RFoptOld <- do.call("internal.rfoptions", c(opt, RELAX=is(model, "formula"))) on.exit(RFoptions(LIST=RFoptOld[[1]])) RFopt <- RFoptOld[[2]] boxcox <- .Call(C_get_boxcox) ## eingabe wird anstonsten auch als vdim_close erwartet -- ## dies ist nocht nicht programmiert! Ausgabe ist schon programmiert ## CondSimu ist auch noch nicht programmiert if (RFopt$general$vdim_close_together) stop("'vdim_close_together' must be FALSE") reg <- MODEL_KRIGE return.variance <- RFopt$krige$return_variance all <- rfPrepareData(model=model, x=x, y=y, z=z, T=T, distances=distances, dim=dim, grid=grid, data=data, given=given, params=params, RFopt=RFopt, reg=reg, err.model = err.model, err.params=err.params, ...) #Print(all); imputing <- all$imputing tsdim <- as.integer(all$Z$tsdim) repet <- as.integer(all$Z$repetitions) vdim <- all$Z$vdim if (!imputing) { coordnames.incl.T <- c(if (!is.null(all$Z$coordnames)) all$Z$coordnames else paste(COORD_NAMES_GENERAL[1], 1:all$Z$spatialdim, sep=""), if (all$Z$has.time.comp) COORD_NAMES_GENERAL[2] else NULL) if (all$X$grid) { coords <- list(x=NULL, T=NULL) xgr <- cbind(all$X$x, all$X$T) colnames(xgr) <- coordnames.incl.T gridTopology <- sp::GridTopology(xgr[1, ], xgr[2, ], xgr[3, ]) ## bis 3.0.70 hier eine alternative } else { coords <- list(x=all$X$x, T=all$X$T) ## wenn bei gegeben unklar was zu tun ist. Ansonsten if (length(coords$T) == 0) colnames(coords$x) <- coordnames.incl.T gridTopology <- NULL } } nx <- all$X$restotal dimension <- if (all$X$grid) c(if (length(all$X$x) > 0) all$X$x[3, ], if (length(all$X$T) > 0) all$X$T[3]) else nx # to do:grid newdim <- c(dimension, if (vdim>1) vdim, if (repet>1) repet) if (imputing && return.variance) { return.variance <- FALSE warning("with imputing, currently the variance cannot be returned") } L <- length(all$model) if (length(all$Z$data) > 1) { Res <- array(dim=c(nx, vdim, repet)) for (i in 1:length(all$Z$data)) { Res[ , , i] <- RFinterpolate(model=model, x=x, y=y, z=z, T=T, grid=grid, distances=distances, dim=dim, data = all$Z$data[[i]], given = all$Z$coord, err.model=err.model, ..., spConform = FALSE, return.variance=FALSE) } dim(Res) <- c(nx * vdim, repet) } else { ## length(all$Z$data) == 1 exact <- RFopt$general$exact maxn <- RFopt$krige$locmaxn ngiven <- as.integer(all$Z$coord[[1]]$restotal) ## number of given points split <- RFopt$krige$locsplitn[1] split <- ngiven > maxn || (!is.na(exact) && !exact && ngiven > split) data <- RFboxcox(all$Z$data[[1]], vdim=vdim) .Call(C_set_boxcox, c(Inf, 0)) if (imputing) { Res <- rep(list(data), L) } else { Res <- rep(list(matrix(nrow=nx, ncol=repet * vdim)), L) } if (return.variance) sigma2 <- rep(list(NULL), L) ## currently just a dummy names(Res) <- ModelAbbreviations(all$model) if (split) { ## to do: all$X <- ExpandGrid(all$X) ## to do all$Z$coord[[1]] <- ExpandGrid(all$Z$coord[[1]]) ## to do ## neighbourhood kriging ! if (!is.na(exact) && exact) stop("number of conditioning locations too large for an exact result.") if (ngiven > maxn && is.na(exact) && RFopt$basic$printlevel>=PL_IMPORTANT) message("performing neighbourhood kriging") stop("neighbourhood kriging currently not programmed") ## calculate the boxes for the locations where we will interpolate idx <- GetNeighbourhoods(Z=Z, X=all$X, ## given locations; to do: grid splitfactor=RFopt$krige$locsplitfactor, maxn=RFopt$krige$locmaxn, split_vec = RFopt$krige$locsplitn, ) totalparts <- length(idx[[2]]) if (totalparts > 1) RFoptions(general.pch="") pr <- totalparts > 1 && RFopt$general$pch != "" &&RFopt$general$pch != " " for (p in 1:totalparts) { stopifnot((Nx <- as.integer(length(idx[[3]][[p]]))) > 0) if (pr && p %% 5==0) cat(RFopt$general$pch) givenidx <- unlist(idx[[1]][idx[[2]][[p]]]) if (ignore.trend) initRFlikelihood(all$krige,# including error structure Reg=reg, grid=FALSE, x=all$Z$coord[[1]]$x[givenidx, , drop=FALSE], data=data[givenidx, , drop=FALSE], ignore.trend = ignore.trend) else RFlikelihood(all$krige, # including error structure Reg=reg, grid=FALSE, x=all$Z$coord[[1]]$x[givenidx, , drop=FALSE], data=data[givenidx, , drop=FALSE], likelihood_register = reg) for (m in i:L) { res <- predictGauss(Reg=reg, model=all$model[[m]], # without error structure x = all$X[idx[[3]][[p]], ], grid = FALSE, kriging_variance=FALSE) if (imputing) { ## TO DO : idx[[3]] passt nicht, da sowohl fuer Daten ## als auch coordinaten verwendet wird. Bei repet > 1 ## ist da ein Problem -- ueberpruefen ob repet=1 where <- all$data.na[idx[[3]][[p]]] ## to do:grid isNA <- is.na(Res[[m]][where, ]) Res[[m]][where, ][isNA] <- res[isNA] } else { Res[[m]][idx[[3]][[p]], ] <- res } } ## for p in totalparts if (pr) cat("\n") } } else { ## no splits if (ignore.trend) initRFlikelihood(all$krige, Reg=reg, x=all$Z$coord, data=data, ignore.trend = ignore.trend) else RFlikelihood(all$krige, x=all$Z$coord, data=data, likelihood_register = reg) for (m in 1:L) { Res[[m]] <- predictGauss(Reg=reg, model=all$model[[m]], x=all$X, kriging_variance=FALSE) } } } ## !is.list(Z) Z <- all$Z ## achtung! oben kann sich noch all$Z aendern! X <- all$X spConform <- RFopt$general$spConform for (m in 1:L) { if (FALSE) {## jonas if (return.variance && length(newdim <- c(dimension, if (vdim>1) vdim)) > 1) base::dim(sigma2[[m]]) <- newdim } if (length(newdim)>1) base::dim(Res[[m]]) <- newdim else Res[[m]] <- as.vector(Res[[m]]) if (!is.null(Z$varnames)) attributes(Res[[m]])$varnames <- Z$varnames Res[[m]] <- RFboxcox(data=Res[[m]], boxcox = boxcox, inverse=TRUE,vdim=vdim) } if (!spConform && !imputing) { if (vdim > 1 && RFopt$general$vdim_close_together) { Resperm <- c(length(dimension)+1, 1:length(dimension), if(repet>1) length(dimension)+2) for (m in 1:L) { Res[[m]] <- aperm(Res[[m]], perm=Resperm) if (return.variance) sigma2[[m]] <- aperm(sigma2[[m]], perm=Resperm[1:(length(dimension)+1)]) } } if (return.variance) for (m in 1:L) Res[[m]] <- list(estim = Res[[m]], var = sigma2[[m]]) return( if (L == 1) Res[[1]] else Res) } if (imputing) { for (m in 1:L) { Res[[m]] <- FinishImputing(data=data, simu=Res[[m]], Z=Z, spConform=spConform, fillall = RFopt$krige$fillall) ## to do : grid if (return.variance){ var <- FinishImputing(data=data, simu=sigma2[[m]], Z=Z, spConform=spConform, fillall = RFopt$krige$fillall)# to do : grid if (spConform) { names(var@data) <- paste("var.", names(var@data), sep="") Res[[m]]@.RFparams$has.variance <- TRUE Res[[m]] <- cbind(Res[[m]], var) } else Res[[m]] <- list(Res[[m]], var) } } return( if (L == 1) Res[[1]] else Res) } else { for (m in 1:L) { Res[[m]] <- conventional2RFspDataFrame(Res[[m]], coords=coords$x, gridTopology=gridTopology, n=repet, vdim=vdim, T = coords$T, vdim_close_together = RFopt$general$vdim_close_together) if (return.variance){ var <- conventional2RFspDataFrame(sigma2[[m]], coords=coords$x, gridTopology=gridTopology, n=1, vdim=vdim, T = coords$T, vdim_close_together = RFopt$general$vdim_close_together) names(var@data) <- paste("var.", names(var@data), sep="") Res[[m]]@.RFparams$has.variance <- TRUE Res[[m]] <- cbind(Res[[m]], var) } } } # Res@.RFparams$krige.method <- # c("Simple Kriging", "Ordinary Kriging", "Kriging the Mean", # "Universal Kriging", "Intrinsic Kriging")[krige.meth.nr] ## Res@.RFparams$var <- sigma2 ## sehr unelegant. ## * plot(Res) sollte zwei Bilder zeigen ## * var(Res) sollte sigma2 zurueckliefern ## * summary(Res) auch summary der varianz, falls vorhanden ## * summary(Res) auch die Kriging methode if (is.raster(x)) { for (m in 1:L) { Res[[m]] <- raster::raster(Res[[m]]) raster::projection(Res[[m]]) <- raster::projection(x) } } return( if (L == 1) Res[[1]] else Res) } rfCondGauss <- function(model, x, y=NULL, z=NULL, T=NULL, grid, n=1, data, # first coordinates, then data given=NULL, ## alternative for coordinates of data params, err.model=NULL, err.params, ...) { # ... wegen der Variablen dots <- list(...) if ("spConform" %in% names(dots)) dots$spConform <- NULL RFoptOld <- internal.rfoptions(..., RELAX=is(model, "formula")) on.exit(RFoptions(LIST=RFoptOld[[1]])) RFopt <- RFoptOld[[2]] boxcox <- .Call(C_get_boxcox) cond.reg <- RFopt$registers$register all <- rfPrepareData(model=model, x=x, y=y, z=z, T=T, grid=grid, data=data, given=given, params=params, RFopt=RFopt, reg=MODEL_KRIGE, err.model = err.model, err.params=err.params,...) Z <- all$Z X <- all$X simu.grid <- X$grid tsdim <- Z$tsdim vdim <- Z$vdim data <- RFboxcox(Z$data[[1]], vdim=vdim) .Call(C_set_boxcox, c(Inf, 0)) if (all$Z$repetitions != 1) stop("conditional simulation allows only for a single data set") txt <- "kriging in space time dimensions>3 where not all the point ly on a grid is not possible yet" ## if 4 dimensional then the last coordinates should ly on a grid ## now check whether and if so, which of the given points belong to the ## points where conditional simulation takes place coord <- ExpandGrid(Z$coord[[1]]) simu <- NULL if (simu.grid) { xgr <- cbind(X$x, X$T) l <- ncol(xgr) ind <- 1 + (t(coord$x) - xgr[1, ]) / xgr[2, ] index <- round(ind) outside.grid <- apply((abs(ind-index)>RFopt$general$gridtolerance) | (index<1) | (index > 1 + xgr[3, ]), 2, any) if (any(outside.grid)) { ## at least some data points are not on the grid: ## simulate as if there is no grid simu.grid <- FALSE if (l>3) stop(txt) xx <- if (l==1) ## dim x locations matrix(seq(from=xgr[1], by=xgr[2], len=xgr[3]), nrow=1) else eval(parse(text=paste("t(expand.grid(", paste("seq(from=xgr[1,", 1:l, "], by=xgr[2,", 1:l, "], len=xgr[3,", 1:l, "])", collapse=","), "))"))) ll <- eval(parse(text=paste("c(", paste("length(seq(from=xgr[1,", 1:l, "], by=xgr[2,", 1:l, "], len=xgr[3,", 1:l, "]))", collapse=","), ")"))) new.index <- rep(0,ncol(index)) ## data points that are on the grid, must be registered, ## so that they can be used as conditioning points of the grid if (!all(outside.grid)) { new.index[!outside.grid] <- 1 + colSums((index[, !outside.grid, drop=FALSE]-1) * cumprod(c(1, ll[-length(ll)]))) } index <- new.index new.index <- NULL } else { ## data points are all lying on the grid simu <- do.call(RFsimulate, args=c(list(model=all$krige, x=X$x, # y=y, z=z, T=X$T, grid=X$grid, n=n, register=cond.reg, seed = NA), dots, list(spConform=FALSE))) ## for all the other cases of simulation see, below ## index <- t(index) index <- 1 + t(index - 1) %*% c(1, xgr[3, -ncol(xgr)]) if (is.vector(simu)) dim(simu) <- c(length(simu), 1) else if (!is.matrix(simu)) { ## 3.1.26 is programmed differently nvdim <- (vdim > 1) + (n > 1) if (nvdim > 0) { d <- dim(simu) last <- length(d) + 1 - (1 : nvdim) dim(simu) <- c(prod(d[-last]), prod(d[last])) } else dim(simu) <- c(length(simu), 1) } } } else { ## not simu.grid xx <- t(X$x) ## dim x locations ## the next step can be pretty time consuming!!! ## to.do: programme it in C ## ## identification of the points that are given twice, as points to ## be simulated and as data points (up to a tolerance distance !) ## this is important in case of nugget effect, since otherwise ## the user will be surprised not to get the value of the data at ## that point one2ncol.xx <- 1:ncol(xx) index <- apply(coord$x, 1, function(u){ i <- one2ncol.xx[colSums(abs(xx - u)) < RFopt$general$gridtolerance] if (length(i)==0) return(0) if (length(i)==1) return(i) return(NA) }) } if (!simu.grid) { ## otherwise the simulation has already been performed (see above) tol <- RFopt$general$gridtolerance * nrow(xx) if (any(is.na(index))) stop("identification of the given data points is not unique - `tol' too large?") if (any(notfound <- (index==0))) { index[notfound] <- (ncol(xx) + 1) : (ncol(xx) + sum(notfound)) } xx <- rbind(t(xx), coord$x[notfound, , drop=FALSE]) simu <- do.call(RFsimulate, args=c(list(model=all$krige, x=xx, grid=FALSE, n=n, register = cond.reg, seed = NA), dots, spConform=FALSE, examples_reduced = FALSE)) if (is.vector(simu)) dim(simu) <- c(length(simu), 1) rm("xx") } if (is.null(simu)) stop("random field simulation failed") simu.given <- simu[index, ] simu <- as.vector(simu[1:X$restotal, ]) # as.vector is necessary !! Otherwis ## is not recognized as a vector # Print(simu, simu.given, simu.grid, index, as.vector(Z$data[[1]]), simu.given, X) ## to do: als Naeherung bei UK, OK: ## kriging(data, method="A") + simu - kriging(simu, method="O") ! stopifnot(length(X$y)==0, length(X$z)==0) simu <- simu + RFinterpolate(x=X, model=model, err.model = err.model, register=MODEL_KRIGE, given = coord, data = as.vector(data) - simu.given, spConform=FALSE, ignore.trend = TRUE) simu <- RFboxcox(data=simu, boxcox = boxcox, inverse=TRUE, vdim=vdim) if (all$imputing) { return(FinishImputing(data=Z$data[[1]], simu=simu, Z=Z, spConform=RFopt$general$spConform, fillall = RFopt$krige$fillall)) } attributes(simu)$varnames <- Z$varnames attributes(simu)$coordnames <- Z$coordnames return(simu) }