Kriging <- function(krige.method, x, y=NULL, z=NULL, T=NULL, grid, gridtriple=FALSE, model, param, given, data, trend, pch=".") { krige.methodlist <- c("S", "O") if (is.na(krige.method.nr <- pmatch(krige.method,krige.methodlist))) stop(paste("kriging method not identifiable from the list.", "Possible values are", paste(krige.methodlist, collapse=","))) x <- CheckXT(x=x, y=y, z=z, T=T, grid=grid, gridtriple=gridtriple) y <- z <- NULL pm <- PrepareModel(model=model, param=param, timespacedim=x$spacedim + x$Time, trend=trend) is.matrix.data <- is.matrix(data) && (ncol(data)>1) data <- as.matrix(data) given <- as.matrix(given) xdim <- ncol(given) if (pm$timespacedim!= xdim) stop("dimensions of the kriged points and the given points do not match") stopifnot(pm$timespacedim == ncol(x$x) + !is.null(x$T)) ## program wrong nd <- nrow(given) pos <- integer(nd) ## lexicographical ordering of vectors --- necessary to check ## whether any location is given twice, but with different value of ## the data .C("Ordering", as.double(t(given)), as.integer(nd), as.integer(xdim), pos, PACKAGE="RandomFields", DUP=FALSE) pos <- pos + 1 given <- given[pos, , drop=FALSE] data <- data[pos, , drop=FALSE] ## are locations given twice with different data values? dup <- c(FALSE, apply(abs(given[-nd, , drop=FALSE] - given[-1, , drop=FALSE]), 1, sum)) if (any(dup <- c(FALSE, apply(abs(given[-1, , drop=FALSE] - given[-nd, , drop=FALSE]), 1, sum)==0))) { if (any(data[dup, ] != data[c(dup[-1], FALSE), ])) stop("duplicated conditioning points with different data values") given <- given[!dup, , drop=FALSE] data <- data[!dup, , drop=FALSE] } tgiven <- t(given) nd <- nrow(given) ## the next passage might be replaced by a direct call of ## covarianceMatrix... covmatrix <- diag(nd) * (CovarianceFct(if (pm$anisotropy) t(rep(0, pm$timespacedim)) else 0, model, param) / 2) low.tri <- lower.tri(covmatrix, diag=FALSE) covmatrix[low.tri] <- CovarianceFct(if (pm$anisotropy) matrix(.C("vectordist", as.double(given), as.integer(dim(given)), vd=double(xdim * nd * (nd - 1)/2), as.integer(FALSE), PACKAGE="RandomFields")$vd, ncol=xdim) else as.matrix(dist(given))[low.tri], model, param) covmatrix <- covmatrix + t(covmatrix) given <- NULL if (grid) { zz <- cbind(x$x, x$T) eval(parse(text=paste("dimension <- c(", paste(paste("length(seq(",zz[1,],",", zz[2,],",", zz[3,],"))"), collapse=","), ")"))) ## `x' will be used in apply within kriging if ( (l <- ncol(zz))==1 ) x <- matrix(seq(zz[1],zz[2],zz[3]),ncol=1) else { text <- paste("x <- expand.grid(", paste("seq(zz[1,",1:l, "],zz[2,",1:l, "],zz[3,",1:l,"])",collapse=","), ")") eval(parse(text=text)) } } else x <- x$x env <- environment() nd.step <- ceiling(nrow(x) / 79) ## for the user's entertainment assign("enu", 0, envir=env) ## ="= nn <- as.integer(ncol(tgiven)) error <- integer(1) .C("InitUncheckedCovFct", as.integer(pm$covnr), as.double(pm$param), as.integer(length(pm$param)), as.integer(pm$timespacedim), as.integer(xdim), as.integer(length(pm$covnr)), as.integer(pm$anisotropy), as.integer(pm$op), as.integer(RFparameters()$PracticalRange), error, PACKAGE="RandomFields", DUP=FALSE); if (error) stop(" Error in definition of covariance function") switch(krige.method.nr, { ## simple kriging# stopifnot(is.null(pm$trend)) if (!(is.matrix(try(data <- as.matrix(solve(covmatrix, data-pm$mean)))))) { res <- x * NA } else { ## time consuming "apply" -- to do : replace by c call res <- apply(x, 1, function(z){ if ((enu %% nd.step) == 0) cat(pch); assign("enu", enu+1, envir=env); .C("UncheckedCovFct", tgiven - z, nn, cov=double(nn), PACKAGE="RandomFields", DUP=FALSE)$cov %*% data} ) + pm$mean } }, { ## ordinary kriging stopifnot(is.null(pm$trend)) covmatrix <- rbind(cbind(covmatrix,1), c(rep(1,nd),0)) if (!(is.matrix(try(data <- solve(covmatrix, rbind(as.matrix(data),0)))))) { res <- x * NA } else { ## time consuming "apply" -- to do : replace by c call res <- apply(x, 1, function(z){ if ((enu %% nd.step) == 0) cat(pch); assign("enu", enu+1, envir=env); c(.C("UncheckedCovFct", tgiven-z, nn, cov=double(nn), PACKAGE="RandomFields", DUP=FALSE)$cov, 1) %*% data }) } } , { ## universal kriging stopifnot(!is.null(trend)) stop("not programmed yet") } ) # switch if (pch!="") cat("\n") ncol.data <- ncol(data) x <- data <- NULL if (is.matrix.data) { if (grid) return(array(t(res), dim=c(dimension, ncol.data))) else return(t(res)) } else { if (grid) return(array(res, dim=dimension)) else return(res) } } CondSimu <- function(krige.method, x, y=NULL, z=NULL, T=NULL, grid, gridtriple=FALSE, model, param, method=NULL, given, data, trend, n=1, register=0, err.model=NULL, err.param=NULL, err.method=NULL, err.register=1, tol=1E-5, pch=".", na.rm=FALSE ) { op.list <- c("+","*") if (is.character(method) && (!is.na(pmatch(method, c("S","O"))))) stop("Sorry. The parameters of the function `CondSimu' as been redefined. Use `krige.method' instead of `method'. See help(CondSimu) for details") x <- CheckXT(x, y, z, T, grid, gridtriple) pm <- PrepareModel(model, param, x$spacedim + x$Time, trend, method=method) y <- z <- NULL total <- x$total krige.mean <- 0 ## pm$mean + err$mean ?? krige.trend <- pm$trend if (!is.null(krige.trend)) stop("not programmed yet") if (!is.null(err.model)) { err <- PrepareModel(err.model, err.param,x$spacedim+x$Time,method=err.method) if (xor(pm$anisotropy, err$anisotropy)) stop("both, the data model and the error model must be either istropic or anisotropic") if (xor(is.null(pm$method), is.null(err$method))) stop("the simulation method must be defined for the data model and the error model or for none of them") if (((length(err$cov)>1 && err$cov!=.C("GetModelNr", "nugget", 1, nr=integer(1), PACKAGE="RandomFields")$nr ) || err$timespacedim>1) && (err$timespacedim != pm$timespacedim)) stop("time-space dimensions do not match") if (!is.null(err$trend)) stop("trend for the error term not allowed") krige <- convert.to.readable(list(covnr=c(pm$covnr, err$covnr), param=c(pm$param, err$param), anisotropy = pm$anisotropy, op = c(pm$op, pmatch("+", op.list) - 1, err$op), mean = krige.mean, trend = krige.trend, method = c(pm$method, err$method), timespacedim = pm$timespacedim ) ) } else krige <-convert.to.readable(list(covnr=pm$covnr, param=pm$param, anisotropy = pm$anisotropy, op = pm$op, mean = krige.mean, trend = krige.trend, method = pm$method, timespacedim = pm$timespacedim ) ) simu.grid <- grid given <- as.matrix(given) if (nrow(given)!=length(data)) { cat("dimension of 'given' does not match 'data'") return(NA) } if (na.rm && any(data.idx <- is.na(data))) { data <- data[data.idx] given <- given[data.idx, ] } 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 if (grid) { zz <- cbind(x$x, x$T) ind <- 1 + (t(given) - zz[1,]) / zz[3,] index <- round(ind) endpoints <- 1 + floor((zz[2,]-zz[1,])/zz[3,]) outside.grid <- apply((abs(ind-index)>tol) | (index<1) | (index>endpoints), 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 l <- ncol(zz) if (l>3) stop(txt) if (l==1) xx <- matrix(seq(x$x[1],x$x[2],x$x[3]),nrow=1) else eval(parse(text=paste("xx <- t(as.matrix(expand.grid(", paste("seq(zz[1,",1:l, "],zz[2,",1:l, "],zz[3,",1:l,"])",collapse=","), ")))"))) eval(parse(text=paste("ll <- c(", paste("length(seq(zz[1,",1:l,"],zz[2,",1:l,"],zz[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 + apply((index[,!outside.grid,drop=FALSE]-1) * cumprod(c(1,ll[-length(ll)])), 2, sum) } index <- new.index new.index <- NULL } else { ## data points are all lying on the grid z <- GaussRF(x=x$x, T=x$T, grid=TRUE, model=model, param=param, method=method, n=n, register=register, gridtriple=TRUE) ## for all the other cases of simulation see, below index <- t(index) } } else { xx <- t(as.matrix(x$x)) ## not a grid if (!is.null(x$T)) { if (ncol(xx) > 2) stop(txt) T <- seq(x$T[1], x$T[2], x$T[3]) ## multiple the xx structure by length of T, ## then add first component of T to first part ... last component of T ## to last part xx <- rbind(matrix(xx, nrow=nrow(xx), ncol=ncol(xx) * length(T)), as.vector(t(matrix(T, nrow=length(T),ncol(xx))))) } ## 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(as.matrix(given),1,function(z){ i <- one2ncol.xx[apply(abs(xx - z),2,sum) < tol] 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 <- tol * 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), given[notfound,,drop=FALSE]) z <- GaussRF(x=xx, grid=FALSE, model=model, param=param, method=method, n=n, register=register) xx <- NULL } if (is.null(z)) stop("random fields simulation failed") if (n==1) { ## z values at the `given' locations zgiven <- z[index] z <- z[1:total] } else { ## this is a bit more complicated since index is either a vector or ## a matrix of dimension dim(z)-1 zgiven <- apply(z, length(dim(z)), function(m) m[index]) z <- as.vector(apply(z,length(dim(z)),function(m) m[1:total])) } if (!is.null(err.model)) { error <- GaussRF(given, grid=FALSE, model=err.model, param=err.param, method=err.method, n=n, register=err.register) if (is.null(error)) stop("error field simulation failed") zgiven <- zgiven + as.vector(error) error <- NULL } z + Kriging(krige.method=krige.method, x=x$x, grid=grid, model=krige, given=given,data=as.vector(data)-zgiven, gridtriple=TRUE, pch=pch) }