Revision b2f6358d2679f8a415bd4f01a706fc38b04bc0c5 authored by Martin Schlather on 26 April 2005, 00:00:00 UTC, committed by Gabor Csardi on 26 April 2005, 00:00:00 UTC
1 parent dec8043
modelling.R.old
Kriging <- function(krige.method, x, y=NULL, z=NULL, T=NULL,
grid, gridtriple=FALSE,
model, param, given, data, trend,
pch=".", return.variance=FALSE, unchecked=FALSE) {
krige.methodlist <- c("S", "O")
krige.method.nr <- pmatch(krige.method, krige.methodlist)
is.matrix.data <- is.matrix(data) && (ncol(data)>1)
if (unchecked) {
stopifnot(!grid)
stopifnot(missing(param))
pm <- model
xdim <- ncol(x)
tgiven <- t(given)
nd <- nrow(given)
} else {
if (is.na(krige.method.nr))
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)
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)
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 <- as.matrix(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
} # !unchecked
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")
## the next passage might be replaced by a direct call of
## covarianceMatrix...
covmatrix <-
diag(rep(0.5 *
.C("UncheckedCovFct",
as.double(if (pm$anisotropy) t(rep(0, pm$timespacedim)) else 0),
as.integer(1), var=double(1), PACKAGE="RandomFields")$var, nd))
low.tri <- lower.tri(covmatrix, diag=FALSE)
cat("\nanisotropy", pm$anisotropy, "\n")
print(if (pm$anisotropy)
t(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.double(as.matrix(dist(given))[low.tri]))
if (nd>1) {
low.tri <- lower.tri(covmatrix, diag=FALSE)
covmatrix[low.tri] <-
.C("UncheckedCovFct", if (pm$anisotropy)
t(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.double(as.matrix(dist(given))[low.tri]),
as.integer(nd * (nd - 1) / 2), cov=double(nd * (nd - 1) / 2),
PACKAGE="RandomFields")$cov
}
covmatrix <- covmatrix + t(covmatrix)
given <- NULL
print(covmatrix)
stopifnot(!any(is.na(covmatrix)))
print(c(krige.method.nr, return.variance))
# print(data)
switch(krige.method.nr,
{ ## simple kriging#
### time consuming "apply" -- to do : replace by c call
### res <- apply(x, 1,
### function(z){
### .C("UncheckedCovFct", as.double(tgiven - z), nn,
### cov=double(nn),
### PACKAGE="RandomFields", DUP=FALSE)$cov %*% data
### }
### ) + pm$mean
stopifnot(is.null(pm$trend))
res <- double(nrow(x) * ncol(data))
if (return.variance) {
if (!(is.matrix(try(invcov <- solve(covmatrix)))))
stop("covmatrix is singular")
sigma2 <- double(nrow(x))
.C("simpleKriging2", as.double(tgiven), as.double(x),
as.double(data-pm$mean), as.double(invcov),
as.integer(nrow(x)), nn, as.integer(ncol(x)),
as.integer(ncol(data)), as.double(pm$mean),
res, sigma2,
PACKAGE="RandomFields", DUP=FALSE)
} else {
if (!(is.matrix(try(invcov <- solve(covmatrix, data-pm$mean)))))
stop("covmatrix is singular")
.C("simpleKriging", as.double(tgiven), as.double(x),
as.double(invcov), as.integer(nrow(x)), nn, as.integer(ncol(x)),
as.integer(ncol(data)), as.double(pm$mean), res,
PACKAGE="RandomFields", DUP=FALSE)
}
res <- matrix(res, nrow=ncol(data))
}, {
## ordinary kriging
### time consuming "apply" -- to do : replace by c call
### covmatrix <- rbind(cbind(covmatrix,1), c(rep(1,nd),0))
### if (!(is.matrix(try(invcov <- solve(covmatrix,
### rbind(as.matrix(data), 0))))))
### stop("covmatrix is singular")
### res <-
### apply(x, 1,
### function(z){
### if ((enu %% nd.step) == 0) cat(pch);
### assign("enu", enu+1, envir=env);
### y <- c(.C("UncheckedCovFct", as.double(tgiven-z), nn,
### cov=double(nn),
### PACKAGE="RandomFields", DUP=FALSE)$cov,
### 1)
### y %*% invcov
### })
### RES0 <- res
stopifnot(is.null(pm$trend))
covmatrix <- rbind(cbind(covmatrix,1), c(rep(1,nd),0))
res <- double(nrow(x) * ncol(data))
print(tgiven)
print(x)
print(covmatrix)
print(dim(x))
print(data)
if (return.variance) {
if (!(is.matrix(try(invcov <- solve(covmatrix)))))
stop("covmatrix is singular")
print(invcov)
sigma2 <- double(nrow(x))
.C("ordinaryKriging2", as.double(tgiven), as.double(x),
as.double(data), as.double(invcov),
as.integer(nrow(x)), nn, as.integer(ncol(x)),
as.integer(ncol(data)), res, sigma2,
PACKAGE="RandomFields", DUP=FALSE)
} else {
if (!(is.matrix(try(invcov <- solve(covmatrix,
rbind(as.matrix(data), 0))))))
stop("covmatrix is singular")
res <- double(nrow(x) * ncol(data))
.C("ordinaryKriging", as.double(tgiven), as.double(x),
as.double(invcov),
as.integer(nrow(x)), nn, as.integer(ncol(x)),
as.integer(ncol(data)), res,
PACKAGE="RandomFields", DUP=FALSE)
}
res <- matrix(res, nrow=ncol(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) res <- array(t(res), dim=c(dimension, ncol.data))
else res <- t(res)
} else {
if (grid) res <- array(res, dim=dimension)
## else res <- res
}
return(if (return.variance)
list(estim=res, var=if (grid) array(sigma2, dim=dimension) else sigma2)
else 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=".",
paired=FALSE,
na.rm=FALSE
) {
op.list <- c("+","*")
if (paired && n==1) stop("if paired, then n must be greater than 1")
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, , drop=FALSE]
}
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 <- matrix(apply(z, length(dim(z)), function(m) m[index]), ncol=n)
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
}
# zgiven is matrix
if (paired)
c(z, -z) + Kriging(krige.method=krige.method,
x=x$x, grid=grid,
model=krige,
given=given,
data=as.vector(data) - cbind(zgiven, -zgiven),
gridtriple=TRUE, pch=pch, ret=FALSE)
else
z + Kriging(krige.method=krige.method,
x=x$x, grid=grid,
model=krige,
given=given, data=as.vector(data) - zgiven,
gridtriple=TRUE, pch=pch, ret=FALSE)
}
Computing file changes ...
