Commit e0e9e178 authored by alge's avatar alge
Browse files

single point prediction, not working

parent c65f212e
matlist.cbind<-function(lst){
# cbind a list of matrices
nmat<-length(lst); ret<-NULL;
# mdim<-dim(lst[[1]])
for(i in 1:nmat) {
# if(dim(lst[[i]])[[1]]!=mdim[[1]] && dim(lst[[i]])[[2]]!=mdim[[2]])
# stop('matrix dimensions not conform in matlist.cbind')
ret<-cbind(ret,lst[[i]])
}
ret
}
check.gridparams<-function(angle,xsw,xne,ysw,yne,
dx,dy,nx,ny){
if(is.null(angle)) angle <- 0
if(is.null(xsw)) xsw <- min(point.obj$x)
if(is.null(xne)) xne <- max(point.obj$x)
if(is.null(ysw)) ysw <- min(point.obj$y)
if(is.null(yne)) yne <- max(point.obj$y)
dgx <- xne-xsw
dgy <- yne-ysw
if(is.null(nx)){
if(is.null(dx)) dx <- dgx/20
nx <- ceiling(dgx/dx)+1
} else {
if(!is.null(dx))
stop("either \"dx\" or \"nx\" should be given!")
if(nx==0)
stop("\"dx\" cannot be 0!")
dx <- dgx/nx
}
if(is.null(ny)){
if(is.null(dy)) dy <- dgy/20
ny <- ceiling(dgy/dy)+1
} else {
if(!is.null(dy))
stop("either \"dy\" or \"ny\" should be given!")
if(ny==0)
stop("\"dy\" cannot be 0!")
dy <- dgy/ny
}
nz <- nx * ny
list(nz=nz,nx=nx,ny=ny,dx=dx,dy=dy,
xsw=xsw,xne=xne,ysw=ysw,yne=yne)
}
remove.duplicates <- function(point.obj,at,
duplicate,dupfun) {
# eliminate duplicates:
xy <- paste(point.obj$x, point.obj$y, sep =",")
idup <- match(xy, xy)
if(duplicate=="user" && !is.function(dupfun))
stop("duplicate=\"user\" requires dupfun to be set to a function")
if(duplicate!="error")
{
centre <- function(x) {
switch(duplicate,
mean = mean(x),
median = median(x),
user = dupfun(x))
}
if(duplicate!="strip"){
point.obj[[match(at, names(point.obj))]] <- unlist(lapply(split(point.obj[[match(at, names(point.obj))]],idup), centre))
ord <- !duplicated(xy)
point.obj$x <- point.obj$x[ord]
point.obj$y <- point.obj$y[ord]
n <- length(point.obj$x)
}
else{
ord <- (hist(idup,plot=F,freq=T,breaks=seq(0.5,max(idup)+0.5,1))$counts==1)
point.obj$x <- point.obj$x[ord]
point.obj$y <- point.obj$y[ord]
point.obj[[match(at, names(point.obj))]] <- point.obj[[match(at, names(point.obj))]][ord]
n <- length(point.obj$x)
}
}
else
if(any(duplicated(xy)))
stop("duplicate data points")
point.obj
}
check.border<-function(extrap,xsw,xne,ysw,yne,nx,ny,border,point.obj){
dog <- matrix(1, nx, ny)
if (!extrap) {
tmpgrd <- cbind(rep(seq(xsw,xne,length=nx),ny),sort(rep(seq(ysw,yne,length=ny),nx)))
if(is.null(border))
dog <- in.chull(tmpgrd[,1], tmpgrd[,2],point.obj$x,point.obj$y)
else {
if(is.null(border$x) | is.null(border$y) | length(border$x)!=length(border$y))
stop("border argument wrong!")
dog <- in.polygon(tmpgrd[,1], tmpgrd[,2], border$x,border$y)
}
# workaround for int <-> unsigned int problem on alpha platform:
dog <- abs(as.numeric(dog))
dog <- matrix(dog, nx, ny,byrow=F)
}
dog
}
check.searchparams <- function(maxdist,rsearch,nsearch,nsmin,nsmax,n){
if(!is.null(maxdist) && is.null(rsaerch))
rsearch<-maxdist
if(rsearch>0 & nsearch>0)
stop("specify only one of rsearch and nsearch\n")
if(nsmin>nsmax)
stop("nsmin>nsmax\n")
if(nsmin<0)
nsmin<-0
if(nsmax<0)
nsmax<-n
if(nsearch<0)
stop("nsearch negative!\n")
list(rsearch=rsearch,nsmin=nsmin,nsmax=nsmax)
}
check.krigedata <- function(point.obj,at,var.mod.obj,mode){
at <- point.obj[[match(at, names(point.obj))]]
n <- length(point.obj$x)
if (!inherits(point.obj, "point"))
stop("point.obj must be of class, \"point\".\n")
if (!inherits(var.mod.obj, "variogram.model"))
stop("var.mod.obj must be of class, \"variogram.model\".\n")
if(mode==1 && length(at)!=n) stop("length of x and z differ\n")
list(n=n)
}
bk.grid <- function(point.obj,
......@@ -169,7 +33,7 @@ bk.grid <- function(point.obj,
ny <- tmp$ny
nz <- tmp$nz
dog <- check.border(extrap,xsw,xne,ysw,yne,nx,ny,border,point.obj)
dog <- check.border.grid(extrap,xsw,xne,ysw,yne,nx,ny,border,point.obj)
point.obj <- remove.duplicates(point.obj,at,duplicate,dupfun)
......@@ -204,6 +68,7 @@ bk.grid <- function(point.obj,
method<-switch(method,direct=2,gqr=1,ols=0)
# lwork <- glsfit.workquery(n,ntrend,method)
lwork <- 3000
if(prior$ntr!=ntrend)
stop("model order of priors does not match!")
npr<-prior$n
......
......@@ -17,34 +17,20 @@ bk.points <- function(point.obj,
method = "gqr")
{
lcbind<-function(lst){
nmat<-length(lst); ret<-NULL;
for(i in 1:nmat) ret<-cbind(ret,lst[[i]])
ret
}
point.obj <- remove.duplicates(point.obj,at,duplicate,dupfun)
tmp <- check.krigedata(point.obj,at,var.mod.obj,mode)
n <- tmp$n
nx<-length(x)
ny<-length(y)
if(nx!=ny)stop("length of x and y differ!")
npts<-nx
n<-length(point.obj$x)
at <- point.obj[[match(at, names(point.obj))]]
if (!inherits(point.obj, "point"))
stop("point.obj must be of class, \"point\".\n")
if (!inherits(var.mod.obj, "variogram.model"))
stop("var.mod.obj must be of class, \"variogram.model\".\n")
if(mode==1 && length(at)!=n) stop("length of x and z differ\n")
if(rsearch>0 & nsearch>0)
stop("specify only one of rsearch and nsearch\n")
if(nsmin>nsmax)
stop("nsmin>nsmax\n")
if(!is.null(maxdist))rsearch<-maxdist
# if(rsearch>0){
# if(nsmin==0) nsmin<-ceiling(n*0.1)
# if(nsmax==0) nsmax<-ceiling(n*0.9)
# }
npts <- length(x)
if(length(y)!=npts)
stop("x and y have different lenght!")
tmp <- check.searchparams(maxdist,rsearch,nsearch,nsmin,nsmax,n)
rsearch <- tmp$rsearch
nsmin <- tmp$nsmin
nsmax <- tmp$nsmax
if(trend==0) ntrend<-1
if(trend==1) ntrend<-3
if(trend==2) ntrend<-6
......@@ -55,161 +41,71 @@ bk.points <- function(point.obj,
spherical=3,
linear=4,
0)
c0<-0
cov0<-0
covmat<-matrix(0,n,n)
snbbit<-rep(0,1+n*npts)
snbbit[1]<-1
# eliminate duplicates:
xy <- paste(point.obj$x, point.obj$y, sep =",")
idup <- match(xy, xy)
if(duplicate=="user" && !is.function(dupfun))
stop("duplicate=\"user\" requires dupfun to be set to a function")
if(duplicate!="error")
{
centre <- function(x) {
switch(duplicate,
mean = mean(x),
median = median(x),
user = dupfun(x))
}
if(duplicate!="strip"){
at <- unlist(lapply(split(at,idup), centre))
ord <- !duplicated(xy)
point.obj$x <- point.obj$x[ord]
point.obj$y <- point.obj$y[ord]
n <- length(point.obj$x)
}
else{
ord <- (hist(idup,plot=F,freq=T,breaks=seq(0.5,max(idup)+0.5,1))$counts==1)
point.obj$x <- point.obj$x[ord]
point.obj$y <- point.obj$y[ord]
at <- at[ord]
n <- length(point.obj$x)
}
}
else
if(any(duplicated(xy)))
stop("duplicate data points")
dopts<-rep(1,length=npts)
if (!extrap) {
if(is.null(border))
dopts <- in.chull(x,y,point.obj$x,point.obj$y)
else {
if(is.null(border$x) | is.null(border$y) | length(border$x)!=length(border$y))
stop("border argument wrong!")
dopts <- in.polygon(x,y, border$x,border$y)
}
# workaround for int <-> unsigned int problem on alpha platform:
dopts <- abs(as.numeric(dopts))
}
# determine optimum array sizes:
# determine optimum array sizes:
if(!is.null(method)){
if(method!="gqr" && method!="direct" && method!="ols")
stop("method (used for glsfit) should be one of \"gqr\", \"ols\" or \"direct\"!")
if(method!="gqr" && method!="direct" && method!="ols")
stop("method (used for glsfit) should be one of \"gqr\", \"ols\" or \"direct\"!")
} else {
method <-"gqr"
method <-"gqr"
}
method<-switch(method,direct=2,gqr=1,ols=0)
method<-switch(method,direct=2,gqr=1,ols=0)
# lwork <- glsfit.workquery(n,ntrend,method)
lwork <- 3000
lwork <- glsfit.workquery(n,ntrend,method)
if(prior$ntr!=ntrend)
stop("model order of priors does not match!")
npr<-prior$n
typpr<-prior$info
typpr[prior$type=="subjective"]<-typpr[prior$type=="subjective"]*(-1)
# in FORTRAN code:
# WRKVEC= muwrk beta errbeta dev errdev zsrnb dist work ferr berr
# length=ntrend+ntrend +1 +n +1 +n +n +lwork +n +n
wrkvec<-double(ntrend+ntrend+1+n+1+n+n+lwork+n+n)
# DATALL= XPTS + YPTS + ZPTS + VARPTS + LON + LAT + Z
# length= NPTS NPTS NPTS NPTS N N N
datall<-c(x,y,double(npts),double(npts),point.obj$x,point.obj$y,at)
snbbit<-rep(0,1+n*npts)
snbbit[1]<-1
dopts<-rep(1,length=npts)
if (!extrap) {
if(is.null(border))
dopts <- in.chull(x,y,point.obj$x,point.obj$y)
else {
if(is.null(border$x) | is.null(border$y) | length(border$x)!=length(border$y))
stop("border argument wrong!")
dopts <- in.polygon(x,y, border$x,border$y)
}
# workaround for int <-> unsigned int problem on alpha platform:
dopts <- abs(as.numeric(dopts))
}
dopts<-check.border.points(extrap,x,y,npts,border,point.obj)
ans<-.Fortran("bkptwr",
n1=as.integer(n),
datall=as.double(datall),
n2=as.integer(n),
wrkvec=double(ntrend+ntrend +1 +n +1 +n +n +lwork +n +n),
napts=integer(npts),
npts=as.integer(npts),
dopts=as.integer(dopts),
n=as.integer(n),
covtype=as.integer(covtype),
covpar=as.double(var.mod.obj$parameters),
covmat=as.double(covmat),
ldcov=as.integer(n),
c0vec=double(n),
cov0=as.double(cov0),
trend=as.integer(trend),
ntrend=as.integer(ntrend),
mupr=as.double(lcbind(prior$mu)),,
ldmpr=as.integer(ntrend),,
phipr = as.double(lcbind(prior$phi)),
ldphpr = as.integer(ntrend),
prinv = as.double(lcbind(prior$phiinv)),
ldpriv = as.integer(ntrend),
phiwrk = double(ntrend*ntrend),
ldphwk = as.integer(ntrend),
lonpr = as.double(prior$lon),
latpr = as.double(prior$lat),
covbta = double(ntrend*ntrend),
ldcvbt = as.integer(ntrend),
cvsrnb = double(n*n),
npr = as.integer(prior$n),
typpr = as.integer(typpr),
rsearch=as.double(rsearch),
nsearch=as.integer(nsearch),
nsmin=as.integer(nsmin),
nsmax=as.integer(nsmax),
fwork=double(n*ntrend),
ldfwrk=as.integer(n),
f0work=double(ntrend),
dist=double(n),
indsnb=integer(n),
indsnw=integer(n),
indsrt=integer(n),
kwork=double(n*n),
ldkwrk=as.integer(n),
rhswork=double(n),
fpwork = double(n*ntrend),
fpfwork = double(n*n),
fpf0wrk = double(n),
chlup = double(n*n),
ldclup = as.integer(n),
cminv = double(n*n),
ldcmnv = as.integer(n),
lwork = as.integer(lwork),
ipvt = integer(ntrend),
ipiv=integer(n+ntrend),
iwork=integer(3*n),
mode=as.integer(mode),
mu=double(ntrend),
lambda=double(n),
lambd0=double(1),
bits=as.integer(c(integer(npts),snbbit)),
ierr=integer(1),
glsmth=as.integer(method),
.Package="baykrig")
ans<-.C("bk_point",
xp = as.double(x),
yp = as.double(y),
zp = double(npts),
varp = double(npts),
dop = as.integer(dopts),
np = as.integer(npts),
lon = as.double(point.obj$x),
lat = as.double(point.obj$y),
z = as.double(point.obj[[match(at, names(point.obj))]]),
extrap = as.integer(extrap),
n = as.integer(n),
covtype = as.integer(covtype),
covpar = as.double(var.mod.obj$parameters),
trend = as.integer(trend),
ntrend = as.integer(ntrend),
mupr = as.double(matlist.cbind(prior$mu)),
ldmpr = as.integer(ntrend),
phipr = as.double(matlist.cbind(prior$phi)),
ldphpr = as.integer(ntrend),
lonpr = as.double(prior$lon),
latpr = as.double(prior$lat),
npr = as.integer(npr),
typpr = as.integer(typpr),
rsearch = as.double(rsearch),
nsearch = as.integer(nsearch),
nsmin = as.integer(nsmin),
nsmax = as.integer(nsmax),
lwork = as.integer(lwork),
mode = as.integer(mode),
lambda = double(n*npts),
lambd0 = double(npts),
bits = as.integer(c(integer(npts),snbbit)),
ierr = integer(1),
glsmth = as.integer(method),
# ,.Package= "baykrig"
)
# ans<-krige.solve(s$x,s$y,point.obj$x,point.obj$y,
......@@ -222,5 +118,19 @@ bk.points <- function(point.obj,
retval$z[ans$napts==1]<-NA
retval$var[ans$napts==1]<-NA
retval
retval<-list(x=ans$xp,
y=ans$yp,
z=ans$zp,
var=ans$varp,
done=ans$dop,
snb=matrix(ans$bits[(npts+2):(nz+n*npts+1)],nrow=n,ncol=npts,byrow=F),
lambda=matrix(ans$lambda,nrow=n,ncol=npts,byrow=FALSE),
lambda0=ans$lambd0
)
retval$z[retval$done<=0] <- NA
retval$var[retval$done<=0] <- NA
retval$lambda0[retval$done<=0] <- NA
retval$data<-point.obj
retval$at<-at
}
check.searchparams <- function(maxdist,rsearch,nsearch,nsmin,nsmax,n){
if(!is.null(maxdist) && is.null(rsaerch))
rsearch<-maxdist
if(rsearch>0 & nsearch>0)
stop("specify only one of rsearch and nsearch\n")
if(nsmin>nsmax)
stop("nsmin>nsmax\n")
if(nsmin<0)
nsmin<-0
if(nsmax<0)
nsmax<-n
if(nsearch<0)
stop("nsearch negative!\n")
list(rsearch=rsearch,nsmin=nsmin,nsmax=nsmax)
}
check.krigedata <- function(point.obj,at,var.mod.obj,mode){
at <- point.obj[[match(at, names(point.obj))]]
n <- length(point.obj$x)
if (!inherits(point.obj, "point"))
stop("point.obj must be of class, \"point\".\n")
if (!inherits(var.mod.obj, "variogram.model"))
stop("var.mod.obj must be of class, \"variogram.model\".\n")
if(mode==1 && length(at)!=n) stop("length of x and z differ\n")
list(n=n)
}
remove.duplicates <- function(point.obj,at,
duplicate,dupfun) {
# eliminate duplicates:
xy <- paste(point.obj$x, point.obj$y, sep =",")
idup <- match(xy, xy)
if(duplicate=="user" && !is.function(dupfun))
stop("duplicate=\"user\" requires dupfun to be set to a function")
if(duplicate!="error")
{
centre <- function(x) {
switch(duplicate,
mean = mean(x),
median = median(x),
user = dupfun(x))
}
if(duplicate!="strip"){
point.obj[[match(at, names(point.obj))]] <- unlist(lapply(split(point.obj[[match(at, names(point.obj))]],idup), centre))
ord <- !duplicated(xy)
point.obj$x <- point.obj$x[ord]
point.obj$y <- point.obj$y[ord]
n <- length(point.obj$x)
}
else{
ord <- (hist(idup,plot=F,freq=T,breaks=seq(0.5,max(idup)+0.5,1))$counts==1)
point.obj$x <- point.obj$x[ord]
point.obj$y <- point.obj$y[ord]
point.obj[[match(at, names(point.obj))]] <- point.obj[[match(at, names(point.obj))]][ord]
n <- length(point.obj$x)
}
}
else
if(any(duplicated(xy)))
stop("duplicate data points")
point.obj
}
matlist.cbind<-function(lst){
# cbind a list of matrices
nmat<-length(lst); ret<-NULL;
# mdim<-dim(lst[[1]])
for(i in 1:nmat) {
# if(dim(lst[[i]])[[1]]!=mdim[[1]] && dim(lst[[i]])[[2]]!=mdim[[2]])
# stop('matrix dimensions not conform in matlist.cbind')
ret<-cbind(ret,lst[[i]])
}
ret
}
check.gridparams<-function(angle,xsw,xne,ysw,yne,
dx,dy,nx,ny){
if(is.null(angle)) angle <- 0
if(is.null(xsw)) xsw <- min(point.obj$x)
if(is.null(xne)) xne <- max(point.obj$x)
if(is.null(ysw)) ysw <- min(point.obj$y)
if(is.null(yne)) yne <- max(point.obj$y)
dgx <- xne-xsw
dgy <- yne-ysw
if(is.null(nx)){
if(is.null(dx)) dx <- dgx/20
nx <- ceiling(dgx/dx)+1
} else {
if(!is.null(dx))
stop("either \"dx\" or \"nx\" should be given!")
if(nx==0)
stop("\"dx\" cannot be 0!")
dx <- dgx/nx
}
if(is.null(ny)){
if(is.null(dy)) dy <- dgy/20
ny <- ceiling(dgy/dy)+1
} else {
if(!is.null(dy))
stop("either \"dy\" or \"ny\" should be given!")
if(ny==0)
stop("\"dy\" cannot be 0!")
dy <- dgy/ny
}
nz <- nx * ny
list(nz=nz,nx=nx,ny=ny,dx=dx,dy=dy,
xsw=xsw,xne=xne,ysw=ysw,yne=yne)
}
check.border.grid<-function(extrap,xsw,xne,ysw,yne,nx,ny,border,point.obj){
dog <- matrix(1, nx, ny)
if (!extrap) {
tmpgrd <- cbind(rep(seq(xsw,xne,length=nx),ny),sort(rep(seq(ysw,yne,length=ny),nx)))
if(is.null(border))
dog <- in.chull(tmpgrd[,1], tmpgrd[,2],point.obj$x,point.obj$y)
else {
if(is.null(border$x) | is.null(border$y) | length(border$x)!=length(border$y))
stop("border argument wrong!")
dog <- in.polygon(tmpgrd[,1], tmpgrd[,2], border$x,border$y)
}
# workaround for int <-> unsigned int problem on alpha platform:
dog <- abs(as.numeric(dog))
dog <- matrix(dog, nx, ny,byrow=F)
}
dog
}
check.border.points<-function(extrap,x,y,npts,border,point.obj){
dopts<-rep(1,length=npts)
if (!extrap) {
if(is.null(border))
dopts <- in.chull(x,y,point.obj$x,point.obj$y)
else {
if(is.null(border$x) | is.null(border$y) | length(border$x)!=length(border$y))
stop("border argument wrong!")
dopts <- in.polygon(x,y, border$x,border$y)
}
# workaround for int <-> unsigned int problem on alpha platform:
dopts <- abs(as.numeric(dopts))
}
dopts
}