Skip to content

GitLab

Commit 21408a35 authored by alge's avatar alge
Browse files

neue cvs root auf sfs

parent d7085292
Show whitespace changes
Inline Side-by-side
R/bk.grid.R
View file @ 21408a35
......@@ -300,7 +300,7 @@ bk.grid <- function(point.obj,
var.mod.obj,
xsw=NULL,ysw=NULL,xne=NULL,yne=NULL,
dx=NULL,dy=NULL,nx=NULL,ny=NULL,
angle=NULL,
angle=0,
maxdist = NULL,
extrap = F,
border=NULL,
......@@ -360,8 +360,8 @@ bk.grid <- function(point.obj,
method <-"gqr"
}
method<-switch(method,direct=2,gqr=1,ols=0)
lwork <- glsfit.workquery(n,ntrend,method)
# lwork <- glsfit.workquery(n,ntrend,method)
lwork <- 3000
if(prior$ntr!=ntrend)
stop("model order of priors does not match!")
npr<-prior$n
......
makefile
View file @ 21408a35
all:
(cd .. ; R INSTALL -l $(R_LIBS) --create-lib-so baykrig)
(cd ~/ ; R INSTALL -l $(R_LIBS) baykrig)
clean:
rm -f src/*.o src/*/*.o src/*/*/*.o src/*.so
......
man/bk.grid.Rd
View file @ 21408a35
......@@ -57,33 +57,32 @@ bk.grid(point.obj, at, prior, var.mod.obj, xsw=NULL, ysw=NULL, xne=NULL, yne=NUL
\examples{
\testonly{
# prepare variables from other example pages if they are not already there:
require(rgeostat)
if(length(ls(pat="leman.78"))==0)data(leman)
if(length(ls(pat="leman.bank"))==0)data(leman.bank)
if(length(ls(pat="leman.pt"))==0) leman.pt<-point(leman.88)
if(length(ls(pat="leman.pr"))==0) leman.pr<-pair(leman.pt,maxdist=30)
if(length(ls(pat="leman.ev"))==0)
leman.ev<-est.variogram(leman.pt,leman.pr,"cadpbm")
if(length(ls(pat="leman.sph"))==0)
leman.sph<-fit.variogram("spherical",leman.ev,0.1,0.2,20,iter=50,plot.it=T)
if(length(ls(pat="leman.prior"))==0) example(empirical.prior)
if(length(ls(pat="leman.prior.1"))==0) example(empirical.prior)
}
leman.bk <-
bk.grid(point = leman.pt, at = "cadpbm",
prior=leman.prior,var.mod.obj = leman.sph,
leman.bk1 <-
bk.grid(point = leman.88.pt, at = "cadpbm",
prior=leman.prior.1,var.mod.obj = leman.88.vmsph,
xsw=min(leman.bank$x),ysw=min(leman.bank$y),
xne=max(leman.bank$x), yne=max(leman.bank$y),
nx=100, ny=100, trend=1, rsearch = 10,
extrap = F,border=leman.bank, duplicate="mean")
plot(leman.bk)
# compare with
require(rgeostat)
leman.krige <-
krige.grid(point = leman.pt, at = "cadpbm",
var.mod.obj = leman.sph,
plot(leman.bk1)
leman.bk0 <-
bk.grid(point = leman.88.pt, at = "cadpbm",
prior=leman.prior.0,var.mod.obj = leman.88.vmsph,
xsw=min(leman.bank$x),ysw=min(leman.bank$y),
xne=max(leman.bank$x), yne=max(leman.bank$y),
nx=100, ny=100, trend=1, rsearch = 10,
extrap = F, border=leman.bank)
nx=100, ny=100, trend=0, rsearch = 10,
extrap = F,border=leman.bank, duplicate="mean")
plot(leman.bk0)
# compare with
require(rgeostat)
leman.krg<-krige.grid.c(min(leman.bank$x),min(leman.bank$y),
max(leman.bank$x), max(leman.bank$y), nx=100, ny=100, point =
leman.88.pt, at = "cadpbm", var.mod.obj = leman.88.vmsph,
rsearch = 7.5,extrap = F,border=leman.bank)
plot(leman.krg,show.snb=T)
}
\keyword{ baykrig }
man/empirical.prior.Rd
View file @ 21408a35
......@@ -37,23 +37,49 @@ empirical.prior(x, formula=NULL, var.mod=NULL, prior=NULL, method="gqr", namx=NU
\examples{
\testonly{
# prepare variables from other example pages if they are not already there:
}
# need some preparation steps using sgeostat functions
require(rgeostat)
# the data set
if(length(ls(pat="leman.78"))==0)data(leman)
if(length(ls(pat="leman.bank"))==0) data(leman.bank)
if(length(ls(pat="leman.pt"))==0) leman.pt<-point(leman.88)
if(length(ls(pat="leman.pr"))==0) leman.pr<-pair(leman.pt,maxdist=30)
if(length(ls(pat="leman.ev"))==0)
leman.ev<-est.variogram(leman.pt,leman.pr,"cadpbm")
if(length(ls(pat="leman.sph"))==0)
leman.sph<-fit.variogram("spherical",leman.ev,0.1,0.2,20,iter=50,plot.it=T)
}
if(length(ls(pat="leman.bank"))==0)data(leman.bank)
# data(leman) contains three measurement grids:
# 1978, 1983, 1988
# generate the point objects:
if(length(ls(pat="leman.78.pt"))==0) leman.78.pt<-point(leman.78)
if(length(ls(pat="leman.83.pt"))==0) leman.83.pt<-point(leman.83)
if(length(ls(pat="leman.88.pt"))==0) leman.88.pt<-point(leman.88)
# use 1988 data for variogram estimation, star with pairs
if(length(ls(pat="leman.88.pr.30"))==0)
leman.88.pr.30<-pair(leman.88.pt,maxdist=30)
# do the estimation
if(length(ls(pat="leman.88.ev30"))==0)
leman.88.ev30<-est.variogram(leman.88.pt,leman.88.pr.30,"cadpbm")
# fit a sperical variogram
if(length(ls(pat="leman.88.vmsph"))==0)
leman.88.vmsph<-fit.variogram("spherical",leman.88.ev30,0.1,0.2,25,plot.it=T)
plot(leman.88.ev30,var.mod.obj=leman.88.vmsph)
# now prepare the prior guesses:
# simple linear trend
leman.prior.1<-empirical.prior(leman.83,cadpbm~x+y,leman.88.vmsph,duplicate="mean")
leman.prior.1<-empirical.prior(leman.78,cadpbm~x+y,leman.88.vmsph,prior=leman.prior.1,duplicate="mean")
# or constant mean model
leman.prior.0<-empirical.prior(leman.83,cadpbm~1,leman.88.vmsph,duplicate="mean",namx="x",namy="y")
leman.prior.0<-empirical.prior(leman.78,cadpbm~1,leman.88.vmsph,prior=leman.prior.0,duplicate="mean",namx="x",namy="y")
# (some outliers have to be removed:)
leman.prior<-empirical.prior(x = leman.83[-c(69,209),],
formula = cadpbm ~ x + y + 1,
var.mod=leman.sph, duplicate="mean")
leman.prior<-empirical.prior(x = leman.78[-c(69,208),],
formula = cadpbm ~ x + y + 1,
var.mod=leman.sph, prior=leman.prior,
duplicate="mean")
#leman.prior<-empirical.prior(x = leman.83[-c(69,209),],
# formula = cadpbm ~ x + y + 1,
# var.mod=leman.sph, duplicate="mean")
#leman.prior<-empirical.prior(x = leman.78[-c(69,208),],
# formula = cadpbm ~ x + y + 1,
# var.mod=leman.sph, prior=leman.prior,
# duplicate="mean")
}
\keyword{ baykrig }
src/Makevars
View file @ 21408a35
OBJS=\
linpack/dgefa.o \
linpack/dgedi.o \
bkwrp.o \
lsfit.o \
lapack/util/ilaenv.o \
BLAS=blas/dger.o \
blas/dtrsv.o \
blas/dtpmv.o \
blas/dgemm.o \
blas/dtrsm.o \
blas/dspr.o \
blas/dtrmv.o \
blas/dsymv.o \
blas/dtrmm.o \
blas/dsyrk.o \
blas/dgemv.o \
blas/dsyr.o \
blas/dsyr2.o \
blas/dsyr2k.o \
LAPACK=lapack/util/ilaenv.o \
lapack/util/ieeeck.o \
lapack/util/lsame.o \
lapack/double/dpotrf.o \
......@@ -76,12 +87,12 @@ lapack/double/dsterf.o \
lapack/double/dsyev.o \
lapack/double/dsytd2.o \
lapack/double/dsytrd.o \
design.o \
bk.o \
bktile.o \
dsysvq.o \
dgelse.o \
slatec/src/j4save.o \
preload/xerbla.o
LINPACK= linpack/dgefa.o \
linpack/dgedi.o
SLATEC=slatec/src/j4save.o \
slatec/src/fdump.o \
slatec/src/i1mach.o \
slatec/src/xercnt.o \
......@@ -90,30 +101,25 @@ slatec/src/dpsort.o \
slatec/src/xerhlt.o \
slatec/src/xersve.o \
slatec/src/xerprn.o \
slatec/src/xermsg.o \
slatec/src/xermsg.o
OBJS=$(LINPACK) $(BLAS) $(LAPACK) $(SLATEC) \
bkwrp.o \
lsfit.o \
design.o \
bk.o \
dsysvq.o \
dgelse.o \
bkgrid.o \
preload/xerbla.o \
errmsg.o \
srchnb.o \
dgggle.o \
tools.o \
bkpts.o \
blas/dger.o \
blas/dtrsv.o \
blas/dtpmv.o \
blas/dgemm.o \
blas/dtrsm.o \
blas/dspr.o \
blas/dtrmv.o \
blas/dsymv.o \
blas/dtrmm.o \
blas/dsyrk.o \
blas/dgemv.o \
blas/dsyr.o \
blas/dsyr2.o \
blas/dsyr2k.o \
glsfit.o \
covfn.o \
matpr.o \
bk_grid.o
# bktile.o \
src/bk.f
View file @ 21408a35
......@@ -42,7 +42,7 @@ c . LDPRIV,
. NSMIN,
. NSMAX,
. FWORK,
. FWORK2,
. FWRK2,
. LDFWRK,
. F0WORK,
. LDF0WK,
......@@ -90,7 +90,7 @@ c . LDPRIV,
DOUBLE PRECISION LAT0(*),LON0(*),LAT(*),LON(*),Z(*),
. COVMAT(LDCOV,*),C0VEC(LDC0,*),COV0,COVPAR(3),
. FWORK(LDFWRK,*),F0WORK(LDF0WK,*),MU(LDMU,*),
. Z0(*),LAMBDA(LDLMBD,*),VAR0(*),FWORK2(LDFWRK,*),
. Z0(*),LAMBDA(LDLMBD,*),VAR0(*),FWRK2(LDFWRK,*),
. RSEARCH,KWORK(LDKWRK,*),RHSWORK(LDKWRK,*),
. FPWORK(LDFWRK,*),FPFWORK(LDFWRK,*),
. FPF0WRK(LDFWRK,*),DIST(*),MUPR(LDMPR,*),
......@@ -255,8 +255,8 @@ c determine search neighbourhood
c prepare the design matrix
CALL DESIGN(LON,LAT,N,INDSNB,NS,LON0,LAT0,N0,INDDO,NDO,
. FWORK,LDFWRK,F0WORK,LDF0WK,NTREND,TREND,IERR)
c make a copy in fwork2
CALL DSUBMM(FWORK,NS,NTREND,1,1,NS,NTREND,LDFWRK,FWORK2,LDFWRK,1)
c make a copy in fwrk2
CALL DSUBMM(FWORK,NS,NTREND,1,1,NS,NTREND,LDFWRK,FWRK2,LDFWRK,1)
c extract appropriate parts of covariance matrix and Z
DO 10 I=1,NS
ZSRNB(I)=Z(INDSNB(I))
......@@ -270,7 +270,7 @@ c do (generalised) least squares fit in search neighbourhood
c CALL GLSFIT(FWORK,N,NTREND,LDFWRK,ZSRNB,CVSRNB,LDCOV,BETA,U,
c . COVBTA,CHLUP,CMINV,WORK1,WORK,LWORK,IPVT,IWORK,IERR)
IF (GLSMTH.EQ.0) THEN
CALL LSFIT(FWORK,FWORK2,NS,NTREND,LDFWRK,ZSRNB,BETA,
CALL LSFIT(FWORK,FWRK2,NS,NTREND,LDFWRK,ZSRNB,BETA,
. ERRBTA,DEV,ERRDEV,
. COVBTA,LDCVBT,SGSQR,
. CMINV,LDCINV,
......@@ -284,7 +284,7 @@ c . COVBTA,CHLUP,CMINV,WORK1,WORK,LWORK,IPVT,IWORK,IERR)
RETURN
END IF
ELSE
CALL GLSFIT(FWORK,FWORK2,NS,NTREND,LDFWRK,ZSRNB,CVSRNB,LDCOV,
CALL GLSFIT(FWORK,FWRK2,NS,NTREND,LDFWRK,ZSRNB,CVSRNB,LDCOV,
. BETA,ERRBTA,DEV,ERRDEV,
. COVBTA,LDCVBT,SGSQR,
. CHLUP,LDCLUP,CMINV,LDCINV,
......
src/bk_grid.c
View file @ 21408a35
......@@ -46,14 +46,14 @@ void bk_grid(double *xsw,
ldclup=(*n),ldcinv=(*n),ldzg=(*nx);
double *covmat, *c0vec, *muwrk, *phiwrk, *beta, errbta,
*dev, errdev, *covbta, *cvsrnb, *zsrnb,
*fwork, *fwork2, *f0work, *dist, *kwork,
*fwork, *fwrk2, *f0work, *dist, *kwork,
*rhswork, *fpwork, *fpfwork, *fpf0wrk, *chlup,
*cminv, *work, *ferr, *berr,
*mu, *lambda, cov0, lambd0;
int *indsnb, *indsnw, *indsrt, *ipiv, *ipvt, *iwork;
#ifdef TRANSIENT
#ifndef TRANSIENT
covmat =Calloc((size_t)(*n)*(*n),double);
c0vec =Calloc((size_t)(*n),double);
muwrk =Calloc((size_t)((*ntrend)*(*npr)),double);
......@@ -64,7 +64,7 @@ void bk_grid(double *xsw,
cvsrnb =Calloc((size_t)(*n)*(*n),double);
zsrnb =Calloc((size_t)(*n),double);
fwork =Calloc((size_t)(*n)*(*ntrend),double);
fwork2 =Calloc((size_t)(*n)*(*ntrend),double);
fwrk2 =Calloc((size_t)(*n)*(*ntrend),double);
f0work =Calloc((size_t)(*ntrend),double);
dist =Calloc((size_t)(*n),double);
indsnb =Calloc((size_t)(*n),int);
......@@ -96,7 +96,7 @@ void bk_grid(double *xsw,
cvsrnb =(double *) R_alloc((*n)*(*n),sizeof(double));
zsrnb =(double *) R_alloc((*n),sizeof(double));
fwork =(double *) R_alloc((*n)*(*ntrend),sizeof(double));
fwork2 =(double *) R_alloc((*n)*(*ntrend),sizeof(double));
fwrk2 =(double *) R_alloc((*n)*(*ntrend),sizeof(double));
f0work =(double *) R_alloc((*ntrend),sizeof(double));
dist =(double *) R_alloc((*n),sizeof(double));
indsnb =(int *) R_alloc((*n),sizeof(int));
......@@ -172,7 +172,7 @@ void bk_grid(double *xsw,
nsmin,
nsmax,
fwork,
fwork2,
fwrk2,
&ldfwrk,
f0work,
dist,
......@@ -204,8 +204,7 @@ void bk_grid(double *xsw,
ierr,
glsmth);
#ifdef TRANSIENT
#ifndef TRANSIENT
Free(lambda);
Free(mu);
Free(iwork);
......@@ -226,17 +225,55 @@ void bk_grid(double *xsw,
Free(indsnb);
Free(dist);
Free(f0work);
Free(fwork2);
// Free(fwrk2);
Free(fwork);
Free(zsrnb);
Free(cvsrnb);
Free(dev);
Free(beta);
Free(covbta);
Free(covbta); // crash bei trend=1
Free(phiwrk);
Free(muwrk);
Free(c0vec);
Free(covmat);
Free(fwrk2); // crash bei trend=0
#else
/*
free(lambda);
free(mu);
free(iwork);
free(berr);
free(ferr);
free(ipiv);
free(ipvt);
free(work);
free(cminv);
free(chlup);
free(fpf0wrk);
free(fpfwork);
free(fpwork);
free(rhswork);
free(kwork);
free(indsrt);
free(indsnw);
free(indsnb);
free(dist);
free(f0work);
free(fwrk2);
free(fwork);
free(zsrnb);
free(cvsrnb);
free(dev);
free(beta);
free(covbta);
free(phiwrk);
free(muwrk);
free(c0vec);
free(covmat);
*/
#endif
}
......
src/bk_grid.h
View file @ 21408a35
......@@ -97,7 +97,7 @@ void F77_NAME(bkgrid)(double *xsw,
int *nsmin,
int *nsmax,
double *fwork,
double *fwork2,
double *fwrk2,
int *ldfwrk,
double *f0work,
double *dist,
......
src/bkgrid.f
View file @ 21408a35
......@@ -25,7 +25,7 @@
. CVSRNB,
. RSEARCH,
. FWORK,
. FWORK2,
. FWRK2,
. LDFWRK,
. F0WORK,
. KWORK,
......@@ -52,7 +52,7 @@
. COVMAT(LDCOV,*),C0VEC(*),COV0,
. RSEARCH,FWORK(LDFWRK,*),F0WORK(*),
. KWORK(LDKWRK,*),RHSWORK(*),MU(*),
. LAMBDA(*),FWORK2(LDFWRK,*),
. LAMBDA(*),FWRK2(LDFWRK,*),
. COVPAR(*),
. FPWORK(LDFWRK,*),FPFWORK(LDFWRK,*),
. FPF0WRK(LDFWRK,*),MUPR(LDMPR,*),
......@@ -143,7 +143,7 @@ c wrapped parameters orig. params(.
. INTVEC(11+INTVEC(9)), NSMIN
. INTVEC(12+INTVEC(9)), NSMAX
. FWORK,
. FWORK2,
. FWRK2,
. LDFWRK,
. F0WORK,
. DBLVEC(INTVEC(1)+INTVEC(2)+2*INTVEC(8)+2*INTVEC(5)+3),DIST
......@@ -233,7 +233,7 @@ c wrapped parameters orig. params(.
. NSMIN,
. NSMAX,
. FWORK,
. FWORK2,
. FWRK2,
. LDFWRK,
. F0WORK,
. DIST,
......@@ -275,7 +275,7 @@ c wrapped parameters orig. params(.
. LON(*),LAT(*),Z(*),COVMAT(LDCOV,*),C0VEC(*),COV0,
. RSEARCH,FWORK(LDFWRK,*),F0WORK(*),
. KWORK(LDKWRK,*),RHSWORK(*),MU(*),
. LAMBDA(*),FWORK2(LDFWRK,*),
. LAMBDA(*),FWRK2(LDFWRK,*),
. COVPAR(*),
. FPWORK(LDFWRK,*),FPFWORK(LDFWRK,*),
. FPF0WRK(LDFWRK,*),DIST(*),MUPR(LDMPR,*),
......@@ -487,7 +487,7 @@ c the main work is now done by BK:
. NSMIN,
. NSMAX,
. FWORK,
. FWORK2,
. FWRK2,
. LDFWRK,
. F0WORK,
. NTREND,
......@@ -523,7 +523,6 @@ c the main work is now done by BK:
. VAR0,
. IERR,
. GLSMTH)
c name="xg\0"
c call matpr(name,xg,nx,1,nx,dbglvl)
c name="yg\0"
......
src/bktile.f
View file @ 21408a35
......@@ -27,7 +27,7 @@ c . LDPRIV,
. CVSRNB,
. RSEARCH,
. FWORK,
. FWORK2,
. FWRK2,
. LDFWRK,
. F0WORK,
. KWORK,
......@@ -54,7 +54,7 @@ c . LDPRIV,
. COVMAT(LDCOV,*),C0VEC(*),COV0,
. RSEARCH,FWORK(LDFWRK,*),F0WORK(*),
. KWORK(LDKWRK,*),RHSWORK(*),MU(*),
. LAMBDA(*),FWORK2(LDFWRK,*),
. LAMBDA(*),FWRK2(LDFWRK,*),
. COVPAR(*),
. FPWORK(LDFWRK,*),FPFWORK(LDFWRK,*),
. FPF0WRK(LDFWRK,*),MUPR(LDMPR,*),
......@@ -146,7 +146,7 @@ c . LDPRIV,
. INTVEC(11+INTVEC(9)), NSMIN
. INTVEC(12+INTVEC(9)), NSMAX
. FWORK,
. FWORK2,
. FWRK2,
. LDFWRK,
. F0WORK,
. DBLVEC(INTVEC(1)+INTVEC(2)+2*INTVEC(8)+2*INTVEC(5)+3),DIST
......@@ -223,7 +223,7 @@ c . LDPRIV,
. NSMIN,
. NSMAX,
. FWORK,
. FWORK2,
. FWRK2,
. F0WORK,
. DIST,
. INDSNB,
......@@ -256,7 +256,7 @@ c . LDPRIV,
. KWORK(NKWORK,*),RHSWORK(NKWORK,*),MU(NTREND,*),
. LAMBDA(N,*),X0(*),Y0(*),Z0(*),VAR0(*),
. XGWORK(NX,*),YGWORK(NX,*),COVPAR(*),
. FWORK2(N,*)
. FWRK2(N,*)
c subroutine for kriging prediction on tiles of a grid
c
......@@ -434,7 +434,7 @@ c fill work arrays with current tile data:
c the main work now is done be KRIGE:
CALL BK(X0,Y0,DO0,INDDO,LIPT,LON,LAT,Z,N,
. COVTYPE,COVPAR,COVMAT,C0,COV0,TREND,
. NTREND,RSEARCH,NSEARCH,NSMIN,NSMAX,FWORK,FWORK2,
. NTREND,RSEARCH,NSEARCH,NSMIN,NSMAX,FWORK,FWRK2,
. F0WORK,
. DIST,INDSNB,INDSNA,INDSRT,KWORK,NKWORK,RHSWORK,
. IPIV,MODE,MU,Z0,LAMBDA,VAR0,IERR)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment