diff --git a/devel/fresh/.Rhistory b/devel/fresh/.Rhistory index 5466b4ffe190941de39cf5fd2ab1f962ccc574f4..2b6f7c8e44010b5b6c886a5d56996cbce7e190b3 100644 --- a/devel/fresh/.Rhistory +++ b/devel/fresh/.Rhistory @@ -1,32 +1,32 @@ +gc() +ls() +leman.bk.0 <- 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=0, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean") +library(baykrig) +leman.bk.0 <- 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=0, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean") +plot(leman.bk.0) +leman.prior.0 +leman.bk.0 <- 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=50, ny=50, trend=0, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean") +plot(leman.bk.0) +source("leman-example.R") +library(baykrig) +leman.bk.0 <- 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=50, ny=50, trend=0, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean") +leman.bk.0 <- 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=0, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean") +plot(leman.bk.0) +dev.print("leman-0.ps") +dev.print(file="leman-0.ps") +library(baykrig) library(baykrig) library(baykrig) library(baykrig) 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=10, ny=10, 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=10, ny=10, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean") -plot(leman.bk1) -leman.bk1$z -leman.bk1$snb -leman.bk1$ -q() -source("leman-example.R") -leman.prior.0<-empirical.prior(leman.83,cadpbm~1,leman.88.vmsph,duplicate="mean",name.x="x",name.y="y") -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,duplicate="mean",namx="x",namy="y")) -leman.prior.0<-empirical.prior(leman.78,cadpbm~1,leman.88.vmsph,prior=leman.prior,duplicate="mean",namx="x",namy="y") -save.image() -leman.bk.0 <- 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=50, ny=50, trend=0, rsearch = 10, extrap = T,border=leman.bank, duplicate="mean") -leman.prior.0 -rm(leman.prior.0) -leman.prior.0<-empirical.prior(leman.83,cadpbm~1,leman.88.vmsph,duplicate="mean",namx="x",namy="y") -leman.prior.0 -leman.prior.0<-empirical.prior(leman.78,cadpbm~1,leman.88.vmsph,prior=leman.prior.0,duplicate="mean",namx="x",namy="y") -save.image() -leman.prior.0 -leman.bk.0 <- 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=50, ny=50, trend=0, rsearch = 10, extrap = T,border=leman.bank, duplicate="mean") -str(leman.bk.0) -save.image() -leman.bk.0$z -gc() q() +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=100, ny=100, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean") +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=100, ny=100, 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=100, ny=100, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean",method="ols") +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=100, ny=100, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean") diff --git a/src/bk.f b/src/bk.f index 13ab47c1468ef83018c70f72d44b8cb19c418538..3a12c3a00ab14dccd10b36496f1547a0ae145f17 100644 --- a/src/bk.f +++ b/src/bk.f @@ -246,17 +246,21 @@ c call imatpr(do0,n0,1,n0,name,dbglvl) IF (NDO.EQ.0) THEN RETURN END IF + c determine search neighbourhood CALL SRCHNB(LON0,LAT0,DO0,INDDO,N0,NDO,LON,LAT,N,NS, . RSEARCH,NSEARCH,NSMIN,NSMAX, . INDSRT,INDSNB,INDSNW,DIST) + 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 fwrk2 CALL DESIGN(LON,LAT,N,INDSNB,NS,LON0,LAT0,N0,INDDO,NDO, . FWRK2,LDFWRK,F0WORK,LDF0WK,NTREND,TREND,IERR) + c 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 @@ -267,6 +271,7 @@ c extract appropriate parts of covariance matrix and Z 10 CONTINUE + 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) @@ -286,7 +291,7 @@ c . COVBTA,CHLUP,CMINV,WORK1,WORK,LWORK,IPVT,IWORK,IERR) END IF ELSE c CALL GLSFIT(FWORK,FWRK2,NS,NTREND,LDFWRK,ZSRNB,CVSRNB,LDCOV, - CALL GLSFIT(FWORK,FWork,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, @@ -304,7 +309,10 @@ c call matpr(name,beta,ntrend,1,ntrend,1) c name="covbta\0" c call matpr(name,covbta,ntrend,ntrend,ntrend,1) END IF -cccc goto 123 + + goto 123 + + c merge priors with search neighbourhood: c average all prior guesses with estimation in search neighbourhood @@ -519,7 +527,7 @@ c name="z0" c call matpr(z0,n0,1,n0,name,dbglvl) write(*,*)"z: ",z0(1) IWORK(1)=NS -cccc 123 continue + 123 continue RETURN END