Revision a88578bf75299850a5f3ff3efba8d520c28d6a81 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 b2f6358
RFtest.precision.R
#q()
# source("RFtest.precision.R")
## file.remove("xx.dat")
if (file.exists("source.R")) {
source("source.R")
if (EXTENDED.TESTING) {
if (FALSE) {
n <- .C("getCov", integer(1), double(0))[[1]]
C <- matrix(.C("getCov", as.integer(n), double(n*n))[[2]], ncol=n)
}
if (FALSE) {
z <- NULL
a <- 1.1
repeat {
z <- c(z, (1-runif(100000))^{-1/a} - a/(a-1))
print(c(length(z),mean(z)))
}
}
if (FALSE) {
library(RandomFields, lib="~/TMP")
model <- "exp"
param <- c(0,1,0,10)
x <- seq(0,5,0.1)
##plot(x, Variogram(x, model=model, param=param))
z <- GaussRF(x,x, model=model, param=param, grid=TRUE, me="TBM2", n=1000)
unix.time(e <- EmpiricalVariogram(x,x,data=z, bin=c(0.05,seq(0.1,3,0.1)),
grid=TRUE))
plot (e$c, e$e, pch=16)
lines(x, Variogram(x, model=model, param=param))
}
##runif(1);save(file="randomseed", .Random.seed)
##load("randomseed")
.path <- "/home/schlather/article/R/NEW.RF/RandomFields/R/"
if (!file.exists(paste(.path, "rf.R", sep=""))) q()
if (file.exists("source.R")) source("source.R") else
if (file.exists("../source.R")) source("../source.R")
x <- seq(0, 10, 0.1)
y <- .C("I0ML0", y=as.double(x), as.integer(length(x)))$y
print(cbind(x,y))
models <- NULL ## if NULL then all models are considered
##model <- NULL
##models <- c("exp", "sp")
##models <- "nnst"
##models <- "exp"
## ****** NOTE falls cauchy in jump.models, dann auch cauchytbm ignoriert !!
jump.models <- c("cone", "bessel", "cauchy",
"cauchytbm", "circular", #nochmal!
"cubic", "dampedcosine",
## "exponent", #: method=TBM2 T=TRUE; funktioniert nicht!!
"nsst", "nsst2", ## beide nicht ueberprueft!!
##"cauchy" ## ci-emb funktioniert nicht (gut);
## run again with grids <- TRUE
"fractalB"
)
##models <- "spherical"
# tbm2:
jump.models <- c("cone", "bessel",
##"cauchy", ## numerical difficulties
##"cauchytbm", ## numerical diffuclties!!!
"circular", "cubic", "dampedcosine", "exponential", "FD",
"fractgauss", "gauss",
"gencauchy", #"gengneiting"
)
##jump.models <- c("cone")
jump.methods <- c("none", "circ", "local", "TBM2", "TBM3", "spectral",
"direct", "nugget", "add.MPP", "hyper", "other")
jump.methods <- "none"
jump.methods <- c("none", "circ", "local", "TBM3", "spectral", "direct",
"nugget", "add.MPP", "hyper", "other")
models <- "cauchytbm"
jump.methods <- c("none", "circ", "local", "TBM2", "spectral", "direct",
"nugget", "add.MPP", "hyper", "other")
if (exists("time.list")) rm("time.list")
##time.list <- list(NULL) ## is set (and then to list(NULL), no temporal
## components are considered
T <- runif(1,1,5); T <- c(T/2, 3*T/2, T);
##time.list <- list(T)
##time.list <- list(NULL)
p <- 0 # percent of randomly skipped tests
## p <- 90
dimensions <- 1:3
grids <- c(FALSE, TRUE)
##grids <- TRUE
max.unix.time <- 1 ## fuer 10 Durchlaeufe.
## 1 entspricht 5min Rechenzeit bei 3000 rep
max.kappas <- 20 ## maximum number of considered combination of kappa
## if the value is passed a random subset is taken
repetitions <- 3000 ## number of simulation the estimation is based on
rep.factor <- 10 ## if deviation then simulation is redone with increased
## number of simulations
##repetitions <- 10
tol = 0.05 ## tolerate deviation (in %)
##repetitions <- 15; tol <- 100
setparameters <- function(n.sf, pch="!") {
sf <- c(5, 20) # linesimufactor
prnt <- c(1, 6) # >6 gives too much output
trials <- c(3, 6)
RFparameters(TBM2.linesimufactor=sf[n.sf],
TBM3D2.linesimufactor=sf[n.sf],
TBMCE.force=TRUE,
TBMCE.trials=trials[n.sf],
TBM2.lines=120,#120,
TBM2.num=TRUE,
CE.trials=4, #trials[n.sf],
CE.force=FALSE,
TBM3D3.linesimufactor=sf[n.sf], pch=pch,
Print=prnt[n.sf], Storing=TRUE, #TBM.method="di",
direct.method=0)
}
ENVIR = environment()
refined.simulation <- function(tol, xx, tt, grid, model, repetitions, method,
v, rep.factor, pch="!") {
DeleteAllRegisters()
assign("z", NULL, envir=ENVIR)
zaehler <- 0
dev <- Inf
setparameters(2, pch)
while (abs(dev) > tol) {
zaehler <- zaehler + 1
cat("\nlarge deviance ")
assign("z", cbind(matrix(GaussRF(xx, T=tt, grid=grid, gridtriple=TRUE,
model=model, n=repetitions, method=method)
, ncol=repetitions), z), envir=ENVIR)
RFparameters(Print=1)
e <- 0.5 * mean((z[1,] - z[nrow(z), ])^2)
dev <- if (abs(v)>1e-10) (e-v)/v else e * 100000
cat("; simu=", e, " (", formatC(dev * 100, dig=2),
"%)",
##formatC((0.5*(diff(z[c(1,4), 1]))^2-v)/v * 100, dig=2),
##"; ", ncol(z),
##paste(formatC(z[c(1,4), 1], dig=3), collapse=", "),
sep="")
if (zaehler > rep.factor && abs(dev) > tol) {
str(model)
cat("too large deviance:", abs(e-v), ">", tol * v,"\n")
DeleteAllRegisters()
return(FALSE)
}
} # while dev > tol
DeleteAllRegisters()
return(TRUE)
}
allmodels <- dimnames(GetModelList())[[1]]
allmethods <- dimnames(GetModelList(FALSE))[[2]]
if (is.null(models)) models <- allmodels
models <- pmatch(models, allmodels, dup=TRUE)
methods <- GetModelList()[models, , drop=FALSE]
models <- allmodels[models]
variance <- runif(1,1,10)
scale <- runif(1,0.1,1)
T <- runif(1,1,5)
T <- c(T/2, 3*T/2, T)
if (!exists("time.list")) time.list <- list(NULL, T)
if (file.exists("xx.dat")) {
e <- NULL
load("xx.dat")
##repetitions <- 1 ## number of simulation the estimation is based on
##rep.factor <- 100000 ## if deviation then simulation
## is redone with increased
## tol <- 0
vv <- v
ee <- e
v <- Variogram(cbind(xx[2,,drop=FALSE]- xx[1,,drop=FALSE], tt[3]),
model=model)
cat("\nRECHECK: \n", model$model, "grid=",grid, "; d=", d,
"; T=", !is.null(tt),
"; ani=", anisotropy, "; k=", paste(nk, collapse=","),
"; vario=", v, " ", method, sep="")
str(model); str(xx); str(tt)
print(ncol(xx))
if (FALSE) {
print(as.matrix(if (grid) {
switch(ncol(xx),
expand.grid(xx[1:2], tt[1:2]),
expand.grid(xx[1:2,1], xx[1:2,2], tt[1:2]),
expand.grid(xx[1:2,1], xx[1:2,2], xx[1:2,3], tt[1:2]),
)
} else {
switch(ncol(xx),
expand.grid(xx, tt[1:2]),
expand.grid(xx[,1], xx[,2], tt[1:2]),
expand.grid(xx[,1], xx[,2], xx[,3], tt[1:2]),
)
}) %*% model[[1]]$aniso)
}
v <- Variogram(cbind(xx[2,,drop=FALSE]- xx[1,,drop=FALSE], tt[3]),
model=model)
if (refined.simulation(tol, xx, tt, grid, model, repetitions, method,
v, rep.factor, pch="!")) {
file.remove("xx.dat")
stop("OK")
}
else {
n <- .C("getCov", integer(1), double(0))[[1]]
C <- .C("getCov", integer(n), double(n*n))[[2]]
stop("recheck failed")
}
}
for (mi in 1:length(models)) {
m <- models[mi]
if (!all(is.na(pmatch(jump.models, m, dup=TRUE)))) next;
cat("\n\n\n\n ************ ", m, "var=", variance, "scale=", scale,
" ************")
gridx <- runif(3, 0.1, 1)
gridx <- rbind(gridx/2, 3*gridx/2, gridx)
x <- list(rbind(runif(3), runif(3)), gridx)
for (grid in grids) {
for (d in dimensions) {
xx <- x[[grid+1]][, 1:d, drop=FALSE]
for (tt in time.list) {
spt <- d + !is.null(tt)
kappa.range <- parameter.range(m, spt)
if (!is.null(kappa.range)) {
if (any(is.na(unlist(kappa.range)))) next;
kappas <- NULL
f <- 2
variable.kappas <- kappa.range$th[[1]][2,] != kappa.range$th[[1]][1,]
for (k in 1:length(kappa.range$th)) {
txt <- paste("expand.grid(",
paste("seq(",kappa.range$pr[[k]][1,],
",", kappa.range$pr[[k]][2,],
",len=", 1 + f * variable.kappas, ")", sep="",
collapse=", "), ")")
kappas <- rbind(kappas, as.matrix(eval(parse(text=txt))))
}
dimnames(kappas) <- NULL
if (nrow(kappas)> max.kappas)
kappas <- kappas[ order(runif(nrow(kappas)))[1:max.kappas], ]
kappas <- lapply(apply(kappas, 1, function(x) list(x)),
function(x) x[[1]])
} else kappas <- list(NULL)
aniso <- matrix(c(-0.5,0.5,1, 0,1.2,0.7, 0.3,0.1,0.9), ncol=3) *
runif(1,1,10)
## print(kappas); yyyy
for (anisotropy in (!is.null(tt)):1) {
if (anisotropy) {
ani <- aniso[1:d, 1:d]
if (!is.null(tt))
ani <- rbind(cbind(ani,0), runif(d + 1, -1.5, 1.5))
model <- list(list(model=m, var=variance, aniso=ani))
} else
model <- list(list(model=m, var=variance, scale=scale))
## print(kappas); xxxx
for (nk in kappas) {
model[[1]]$kappa <- nk
RFparameters(Print=2)
v <- Variogram(cbind(xx[2,,drop=FALSE]- xx[1,,drop=FALSE], tt[3]),
model=model)
# print(v)
# print(cbind(xx[2,,drop=FALSE]- xx[1,,drop=FALSE]))
# print(tt)
# str(model)
stopifnot(!is.na(v))
cat("\n\n", m, ": grid=",grid, "; d=", d, "; T=", !is.null(tt),
"; ani=", anisotropy, "; k=", paste(nk, collapse=","),
"; vario=", v, sep="")
last.failed <- FALSE
for (me in 1:length(allmethods)) {
if (!all(is.na(pmatch(jump.methods, allmethods[me], dup=TRUE))))
next;
if (methods[mi, me] && (grid || allmethods[me]!="circu")) {
setparameters(1)
method <- allmethods[me]
if (runif(1,0,100) < p) {cat("%"); next}
if (!last.failed) cat("\n")
cat("method=", allmethods[me], sep="");
DeleteRegister(0)
save(file="xx.dat", v, xx, tt, model, grid, repetitions,
method, d, anisotropy, nk)
# RFparameters(Print=5)
InitGaussRF(xx, T=tt, grid=grid,gridtriple=TRUE,
model=model, method=method)
ut <- unix.time(z <-GaussRF(xx, T=tt, grid=grid,gridtriple=TRUE,
model=model, n=10, method=method))[1]
if (last.failed <- is.null(z)) {
# cat(" failed")
if (FALSE)
{
if ((grid || me!=1) && (is.null(tt) || me!=5)
&& (me!=5 || d<=2) && (me!=3 || d<=2)
)
readline(" -- what now?")
}
## RFparameters(Print=20)
next
}
if (ut>max.unix.time) {
cat(" too time consuming:", ut / 10 , "s per simulation")
next
} else cat(" ")
z <-
matrix(c(as.double(z),
as.double(GaussRF(xx, T=tt, grid=grid,
gridtriple=TRUE, model=model,
n=repetitions-10, method=method))),
ncol=repetitions)
e <- 0.5 * mean((z[1,] - z[nrow(z), ])^2)
dev <- if (abs(v)>1e-10) (e-v)/v else e * 100000
cat("; simu=", e, " (", formatC(dev * 100, dig=2),"%)", sep="")
if (abs(dev) > tol &&
!refined.simulation(tol, xx, tt, grid, model, repetitions,
method, v, rep.factor)) {
save(file="xx.dat", v, xx, tt, model, grid, repetitions,
method, e, d, anisotropy, nk)
stop("check failed")
} else file.remove("xx.dat")
}
} # for method
} # for kappas
} # for anisotropy
} # for time
} # for dimensions
} # for grid yes/no
} # for models
if (FALSE) {
.C("printKEY", as.integer(0))
}
} # extended testing
} # exists source
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...