Commit e4510cc1 authored by agebhard's avatar agebhard
Browse files

return of lambda ... optional (get.lm). crashes if get.lm==T.

parent a897e6ab
...@@ -18,7 +18,8 @@ bk.grid <- function(point.obj, ...@@ -18,7 +18,8 @@ bk.grid <- function(point.obj,
mode=3, mode=3,
duplicate = "error", duplicate = "error",
dupfun = NULL, dupfun = NULL,
method="gqr") method="gqr",
get.lm=FALSE)
{ {
tmp <- check.gridparams(angle,xsw,xne,ysw,yne, tmp <- check.gridparams(angle,xsw,xne,ysw,yne,
dx,dy,nx,ny) dx,dy,nx,ny)
...@@ -78,6 +79,15 @@ bk.grid <- function(point.obj, ...@@ -78,6 +79,15 @@ bk.grid <- function(point.obj,
snbbit<-rep(0,1+n*nz) snbbit<-rep(0,1+n*nz)
snbbit[1]<-1 snbbit[1]<-1
if(get.lm){
mu <- double(ntrend*nz)
lambda <- double(n*nz)
lambd0 <- double(nz)
} else {
mu <- double(ntrend)
lambda <- double(n)
lambd0 <- double(1)
}
ans<-.C("bk_grid", ans<-.C("bk_grid",
...@@ -121,10 +131,12 @@ bk.grid <- function(point.obj, ...@@ -121,10 +131,12 @@ bk.grid <- function(point.obj,
nsmax = as.integer(nsmax), nsmax = as.integer(nsmax),
lwork = as.integer(lwork), lwork = as.integer(lwork),
mode = as.integer(mode), mode = as.integer(mode),
lambda = double(n*nz), mu = as.double(mu),
lambd0 = double(nz), lambda = as.double(lambda),
lambd0 = double(lambd0),
bits = as.integer(c(integer(nz),snbbit)), bits = as.integer(c(integer(nz),snbbit)),
ierr = integer(1), ierr = integer(1),
get.lm = as.integer(get.lm),
glsmth = as.integer(method), glsmth = as.integer(method),
# ,.Package= "baykrig" # ,.Package= "baykrig"
) )
...@@ -134,10 +146,13 @@ bk.grid <- function(point.obj, ...@@ -134,10 +146,13 @@ bk.grid <- function(point.obj,
z=matrix(ans$zg,nx,ny), z=matrix(ans$zg,nx,ny),
var=matrix(ans$varg,nx,ny), var=matrix(ans$varg,nx,ny),
done=matrix(ans$dog, nx, ny), done=matrix(ans$dog, nx, ny),
snb=matrix(ans$bits[(nz+2):(nz+n*nz+1)],nrow=n,ncol=nz,byrow=F), snb=matrix(ans$bits[(nz+2):(nz+n*nz+1)],nrow=n,ncol=nz,byrow=F)
lambda=matrix(ans$lambda,nrow=n,ncol=nz,byrow=FALSE),
lambda0=matrix(ans$lambd0,nx,ny)
) )
if(get.lm){
retval$lambda <- matrix(ans$lambda,nrow=n,ncol=nz,byrow=FALSE)
retval$lambda0 <- matrix(ans$lambd0,nx,ny)
retval$mu <- matrix(ans$mu,nrow=n,ncol=ntrend,byrow=FALSE)
}
retval$z[retval$done<=0] <- NA retval$z[retval$done<=0] <- NA
retval$var[retval$done<=0] <- NA retval$var[retval$done<=0] <- NA
retval$lambda0[retval$done<=0] <- NA retval$lambda0[retval$done<=0] <- NA
......
leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,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",method="ols") ls()
library(baykrig) library(baykrig)
library(baykrig) library(baykrig)
leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,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=5, ny=5, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
library(baykrig) library(baykrig)
leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,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=5, ny=5, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,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=50, ny=50, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,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=150, ny=50, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
library(baykrig) library(baykrig)
leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,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=10, ny=10, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
plot(leman.bk)
leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,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)
leman.bk.2<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,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=200, ny=200, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
library(baykrig) library(baykrig)
leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,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=10, ny=10, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean") leman.bk.pt <- bk.points(point = leman.88.pt , at = "cadpbm", prior=leman.prior, var.mod.obj = leman.88.vmsph, x=c(530,550),y=c(142,147), trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
plot(leman.bk)
library(baykrig) library(baykrig)
leman.bk<- bk.grid(point = leman.88.pt , at = "cadpbm", prior=leman.prior,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=50, ny=50, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean") leman.bk.pt <- bk.points(point = leman.88.pt , at = "cadpbm", prior=leman.prior.1, var.mod.obj = leman.88.vmsph, x=c(530,550),y=c(142,147), trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
plot(leman.bk) example(prior)
plot(leman.bk,show.snb=T) example(emp.prior)
plot(leman.bk,show.snb=T,bias=T) ls("package:baykrig")
dev.print(file="bias.ps") example(empirical.prior)
leman.bk.pt <- bk.points(point = leman.88.pt , at = "cadpbm", prior=leman.prior.1, var.mod.obj = leman.88.vmsph, x=c(530),y=c(142), trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean")
leman.bk.pt
str(leman.bk.pt)
q() q()
library(baykrig) library(baykrig)
plot(leman.bk,show.snb=T,bias=T) example(bk.grid)
ls()
library(baykrig) library(baykrig)
leman.bk.pt <- bk.points(point = leman.88.pt , at = "cadpbm", prior=leman.prior, var.mod.obj = leman.88.vmsph, x=c(530,550),y=c(142,147), trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean") example(empirical.prior)
ls() save.image()
help(bk.grid)
library(baykrig) library(baykrig)
leman.bk.pt <- bk.points(point = leman.88.pt , at = "cadpbm", prior=leman.prior, var.mod.obj = leman.88.vmsph, x=c(530,550),y=c(142,147), trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean") library(baykrig)
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",get.lm=T)
library(baykrig)
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",get.lm=F)
str(leman.bk1)
str(leman.bk1$var)
(leman.bk1$var)
q()
#include "bk_grid.h" #include "bk_grid.h"
void F77_NAME(bk_grid)(double *xsw, void bk_grid__(double *xsw,
double *ysw, double *ysw,
double *xne, double *xne,
double *yne, double *yne,
...@@ -40,10 +40,12 @@ void F77_NAME(bk_grid)(double *xsw, ...@@ -40,10 +40,12 @@ void F77_NAME(bk_grid)(double *xsw,
int *nsmax, int *nsmax,
int *lwork, int *lwork,
int *mode, int *mode,
double *mu,
double *lambda, double *lambda,
double *lambd0, double *lambd0,
int *searchnb, int *searchnb,
int *ierr, int *ierr,
int *retlm,
int *glsmth){ int *glsmth){
/* simple Fortran wrapper */ /* simple Fortran wrapper */
bk_grid(xsw, bk_grid(xsw,
...@@ -86,10 +88,12 @@ bk_grid(xsw, ...@@ -86,10 +88,12 @@ bk_grid(xsw,
nsmax, nsmax,
lwork, lwork,
mode, mode,
mu,
lambda, lambda,
lambd0, lambd0,
searchnb, searchnb,
ierr, ierr,
retlm,
glsmth); glsmth);
} }
...@@ -133,10 +137,12 @@ void bk_grid(double *xsw, ...@@ -133,10 +137,12 @@ void bk_grid(double *xsw,
int *nsmax, int *nsmax,
int *lwork, int *lwork,
int *mode, int *mode,
double *mu,
double *lambda, double *lambda,
double *lambd0, double *lambd0,
int *searchnb, int *searchnb,
int *ierr, int *ierr,
int *retlm,
int *glsmth){ int *glsmth){
int nz=(*nx)*(*ny), ldc0=(*n),ldphwk=(*ntrend), int nz=(*nx)*(*ny), ldc0=(*n),ldphwk=(*ntrend),
...@@ -148,7 +154,6 @@ void bk_grid(double *xsw, ...@@ -148,7 +154,6 @@ void bk_grid(double *xsw,
*fwork, *fwrk2, *f0work, *dist, *kwork, *fwork, *fwrk2, *f0work, *dist, *kwork,
*rhswork, *fpwork, *fpfwork, *fpf0wrk, *chlup, *rhswork, *fpwork, *fpfwork, *fpf0wrk, *chlup,
*cminv, *work, *ferr, *berr, *cminv, *work, *ferr, *berr,
*mu,
cov0; cov0;
int *indsnb, *indsnw, *indsrt, *ipiv, *ipvt, *iwork; int *indsnb, *indsnw, *indsrt, *ipiv, *ipvt, *iwork;
...@@ -183,7 +188,6 @@ void bk_grid(double *xsw, ...@@ -183,7 +188,6 @@ void bk_grid(double *xsw,
ferr =Calloc((size_t)(*n),double); ferr =Calloc((size_t)(*n),double);
berr =Calloc((size_t)(*n),double); berr =Calloc((size_t)(*n),double);
iwork =Calloc((size_t)(3*(*n)),int); iwork =Calloc((size_t)(3*(*n)),int);
mu =Calloc((size_t)(*ntrend),double);
#else #else
c0vec =(double *) R_alloc((*n),sizeof(double)); c0vec =(double *) R_alloc((*n),sizeof(double));
muwrk =(double *) R_alloc((*ntrend),sizeof(double)); muwrk =(double *) R_alloc((*ntrend),sizeof(double));
...@@ -213,7 +217,6 @@ void bk_grid(double *xsw, ...@@ -213,7 +217,6 @@ void bk_grid(double *xsw,
ferr =(double *) R_alloc((*n),sizeof(double)); ferr =(double *) R_alloc((*n),sizeof(double));
berr =(double *) R_alloc((*n),sizeof(double)); berr =(double *) R_alloc((*n),sizeof(double));
iwork =(int *) R_alloc(3*(*n),sizeof(int)); iwork =(int *) R_alloc(3*(*n),sizeof(int));
mu =(double *) R_alloc((*ntrend),sizeof(double));
#endif #endif
F77_CALL(bkgrid)(xsw, F77_CALL(bkgrid)(xsw,
...@@ -301,10 +304,10 @@ void bk_grid(double *xsw, ...@@ -301,10 +304,10 @@ void bk_grid(double *xsw,
lambd0, lambd0,
searchnb, searchnb,
ierr, ierr,
retlm,
glsmth); glsmth);
#ifndef TRANSIENT #ifndef TRANSIENT
Free(mu);
Free(iwork); Free(iwork);
Free(berr); Free(berr);
Free(ferr); Free(ferr);
...@@ -338,7 +341,6 @@ void bk_grid(double *xsw, ...@@ -338,7 +341,6 @@ void bk_grid(double *xsw,
#else #else
/* /*
free(mu);
free(iwork); free(iwork);
free(berr); free(berr);
free(ferr); free(ferr);
......
...@@ -43,13 +43,15 @@ void bk_grid(double *xsw, ...@@ -43,13 +43,15 @@ void bk_grid(double *xsw,
int *nsmax, int *nsmax,
int *lwork, int *lwork,
int *mode, int *mode,
double *mu,
double *lambda, double *lambda,
double *lambd0, double *lambd0,
int *searchnb, int *searchnb,
int *ierr, int *ierr,
int *retlm,
int *glsmth); int *glsmth);
void F77_NAME(bk_grid)(double *xsw, void bk_grid__(double *xsw,
double *ysw, double *ysw,
double *xne, double *xne,
double *yne, double *yne,
...@@ -89,10 +91,12 @@ void F77_NAME(bk_grid)(double *xsw, ...@@ -89,10 +91,12 @@ void F77_NAME(bk_grid)(double *xsw,
int *nsmax, int *nsmax,
int *lwork, int *lwork,
int *mode, int *mode,
double *mu,
double *lambda, double *lambda,
double *lambd0, double *lambd0,
int *searchnb, int *searchnb,
int *ierr, int *ierr,
int *retlm,
int *glsmth); int *glsmth);
void F77_NAME(bkgrid)(double *xsw, void F77_NAME(bkgrid)(double *xsw,
...@@ -180,6 +184,7 @@ void F77_NAME(bkgrid)(double *xsw, ...@@ -180,6 +184,7 @@ void F77_NAME(bkgrid)(double *xsw,
double *lambd0, double *lambd0,
int *searchnb, int *searchnb,
int *ierr, int *ierr,
int *retlm,
int *glsmth); int *glsmth);
...@@ -83,9 +83,10 @@ ...@@ -83,9 +83,10 @@
. LAMBD0, . LAMBD0,
. BITS, . BITS,
. IERR, . IERR,
. RETLM,
. GLSMTH) . GLSMTH)
IMPLICIT NONE IMPLICIT NONE
INTEGER NX,NY,NZ,N,COVTYPE,TREND,NTREND,LDZG, INTEGER NX,NY,NZ,N,COVTYPE,TREND,NTREND,LDZG,RETLM,
. NSEARCH,NSMIN,NSMAX,MODE,IERR,INDSNB(*),INDSNW(*), . NSEARCH,NSMIN,NSMAX,MODE,IERR,INDSNB(*),INDSNW(*),
. INDSRT(*),IPIV(*),IPT,EXTRAP,DOG(NX,*),LDKWRK,EXTCOV, . INDSRT(*),IPIV(*),IPT,EXTRAP,DOG(NX,*),LDKWRK,EXTCOV,
. LDCOV,LDFWRK,LDMPR,LDPHPR,LDPHWK,LDCVBT,LWORK,LDLMBD, . LDCOV,LDFWRK,LDMPR,LDPHPR,LDPHWK,LDCVBT,LWORK,LDLMBD,
...@@ -94,7 +95,7 @@ ...@@ -94,7 +95,7 @@
. XG(*),YG(*),ZG(LDZG,*),VARG(LDZG,*), . XG(*),YG(*),ZG(LDZG,*),VARG(LDZG,*),
. LON(*),LAT(*),Z(*),COVMAT(LDCOV,*),C0VEC(*),COV0, . LON(*),LAT(*),Z(*),COVMAT(LDCOV,*),C0VEC(*),COV0,
. RSEARCH,FWORK(LDFWRK,*),F0WORK(*), . RSEARCH,FWORK(LDFWRK,*),F0WORK(*),
. KWORK(LDKWRK,*),RHSWORK(*),MU(*), . KWORK(LDKWRK,*),RHSWORK(*),MU(NTREND,*),
. LAMBDA(LDLMBD,*),FWRK2(LDFWRK,*), . LAMBDA(LDLMBD,*),FWRK2(LDFWRK,*),
. COVPAR(*), . COVPAR(*),
. FPWORK(LDFWRK,*),FPFWORK(LDFWRK,*), . FPWORK(LDFWRK,*),FPFWORK(LDFWRK,*),
...@@ -255,10 +256,16 @@ c loop over all points and pass them to KRIGE: ...@@ -255,10 +256,16 @@ c loop over all points and pass them to KRIGE:
if (bits(1+nz).ne.0) then if (bits(1+nz).ne.0) then
usesbbt=1 usesbbt=1
end if end if
PCNT=0 if (retlm.eq.1) then
PCNT=0
else
PCNT=1
end if
DO 20 I=1,NX DO 20 I=1,NX
DO 10 J=1,NY DO 10 J=1,NY
PCNT=PCNT+1 if (retlm.eq.1) then
PCNT=PCNT+1
end if
BITS(I+NY*(J-1)) = PCNT BITS(I+NY*(J-1)) = PCNT
X0=XG(I) X0=XG(I)
Y0=YG(J) Y0=YG(J)
...@@ -335,7 +342,7 @@ c the main work is now done by BK: ...@@ -335,7 +342,7 @@ c the main work is now done by BK:
. IPIV, . IPIV,
. IWORK, . IWORK,
. MODE, . MODE,
. MU, . MU(1,pcnt),
. NTREND, . NTREND,
. Z0, . Z0,
. NA0, . NA0,
...@@ -359,12 +366,16 @@ c extract results for this tile ...@@ -359,12 +366,16 @@ c extract results for this tile
DOG(I,J)=-1 DOG(I,J)=-1
ZG(I,J)=0 ZG(I,J)=0
VARG(I,J)=0 VARG(I,J)=0
LAMBD0(I,J)=0 if (retlm.eq.1) then
LAMBD0(I,J)=0
end if
c write(*,*)"x" c write(*,*)"x"
ELSE ELSE
ZG(I,J)=Z0 ZG(I,J)=Z0
VARG(I,J)=VAR0 VARG(I,J)=VAR0
LAMBD0(I,J)=L0 if (retlm.eq.1) then
LAMBD0(I,J)=L0
end if
c write(*,*)"o" c write(*,*)"o"
c VARG(I,J)=IWORK(1)*1.0D0 c VARG(I,J)=IWORK(1)*1.0D0
c VARG(I,J)=lambd0 c VARG(I,J)=lambd0
......
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