# source("RFtest.Covariance.R") if (file.exists("source.R")) source("source.R") epsilon <- 0.000000001 tf <- function(x, l, i=1) { x <- as.matrix(x) if (!is.null(l[[i]]$s)) { return(sqrt(apply(x^2, 1, sum))/l[[i]]$s) } else { return(sqrt(apply((x %*% matrix(l[[i]]$a, nrow=ncol(x)))^2, 1, sum))) } } tfst <- function(x, m) { if (is.matrix(x)) { if (is.null(m$a)) { stop("tfst -- no $aniso"); ERR } else { x <- x %*% matrix(m$a, nrow=ncol(x)) return(cbind(sqrt(apply(x[,-ncol(x),drop=FALSE]^2, 1, sum)), abs(x[,ncol(x)]))) } } else { stop("tfst -- iso"); ERR } } st.fct <- function(x, m) { eval(parse(text=paste(letters[1:6], "<- m$k[", 1:6,"];"))) y <- abs(tfst(x, m)) switch(e, psi <- sqrt(y[,ncol(y)]^c +1)^d, psi <- sqrt((d^(-1) * y[,ncol(y)]^c + 1) / (y[,ncol(y)]^c + 1)), psi <- sqrt(-log(y[,ncol(y)]^c + 1/d)/log(d)) ) if (ncol(y)==2) { cov <- c("stable", "whittle", "cauchy")[b] phi <- CovarianceFct(y[,1] / psi, model=cov, param=c(0, m$v, 0, 1, a)) } else { stopifnot(ncol(y)==1); ERR phi <- 1 } return(psi^(-f) * phi) } aniso <- function(dim) matrix(runif(dim^2), ncol=dim) id <- function(dim) diag(dim) x0 <- list(c(1:4), cbind(1:4, c(3,2,5,6))) x1 <- list(cbind(0:4, c(0,3,2,5,6)), cbind(0:4, c(0,3,2,5,6), c(0,2:5))) models <- list( list(list(list(model="exp", var=3, scale =5)), function(x, l) { exp(-tf(x,l,1)) * l[[1]]$v}, x0), list(list(list(model="exp", var=7, scale =5),"+", list(model="gencau", var=3, scale =5, k=c(2,3))), function(x, l) { exp(-tf(x,l,1)) * l[[1]]$v + (1+tf(x,l,3)^l[[3]]$k[1])^(-l[[3]]$k[2]/l[[3]]$k[1]) * l[[3]]$v }, x0), list(list(list(model="exp", var=7, scale =5),"*", list(model="gencau", var=3, scale =5, k=c(2,3))), function(x, l) { exp(-tf(x,l,1)) * l[[1]]$v * (1+tf(x,l,3)^l[[3]]$k[1])^(-l[[3]]$k[2]/l[[3]]$k[1]) * l[[3]]$v }, x0), list(list(list(model="exp", var=7, aniso=aniso),"*", list(model="gencau", var=3, aniso=aniso, k=c(2,3))), function(x, l) { exp(-tf(x,l,1)) * l[[1]]$v * (1+tf(x,l,3)^l[[3]]$k[1])^(-l[[3]]$k[2]/l[[3]]$k[1]) * l[[3]]$v }, x0), list(list(list(model="nsst", var=7, aniso=aniso, kappa=c(1.1, 1, 1.5,0.7, 1, 3))), function(x,l) st.fct(x, l[[1]]), x1), list(list(list(model="nsst", var=7, aniso=aniso, kappa=c(1.1, 2, 1.5,0.7, 2, 3))), function(x,l) st.fct(x, l[[1]]), x1), list(list(list(model="nsst", var=7, aniso=aniso, kappa=c(1.1, 3, 1.5,0.7, 3, 3))), function(x,l) st.fct(x, l[[1]]), x1), list(list(list(model="nsst", var=7, aniso=aniso, kappa=c(1.3, 3, 1.2,0.5, 2, 4)),"+", list(model="nsst", var=5, aniso=aniso, kappa=c(1.2, 1, 1.5,0.7, 2, 3))), function(x,l) st.fct(x, l[[1]]) + st.fct(x, l[[3]]), x1), list(list(list(model="nsst", var=7, aniso=aniso, kappa=c(1.3, 3, 1.2,0.5, 2, 4)),"*", list(model="nsst", var=5, aniso=aniso, kappa=c(1.2, 1, 1.5,0.7, 2, 3)), "+", list(model="nsst", var=5, aniso=aniso, kappa=c(1.2, 1, 1.5,0.7, 1, 3))), function(x,l)st.fct(x,l[[1]]) * st.fct(x,l[[3]]) + st.fct(x,l[[5]]), x1), ) { for (model in models) { cat("\n",paste(c(lapply(model[[1]], function(l) l[[1]]), recursive=TRUE)),"..") for (x in model[[3]]) { f <- model[[2]] M <- model[[1]] for (i in 1:length(model[[1]])) if (is.function(M[[i]]$aniso)) M[[i]]$aniso <- M[[i]]$aniso(ncol(as.matrix(x))) #print(M) if ( any(abs(CovarianceFct(as.matrix(x), M) - f(x, M)) > f(x, M) * epsilon)){ #print(model) print("X") print(x) print("C code") print(CovarianceFct(x, M)) print("f eval") print(f(x, M)) print("-----------") print(abs(CovarianceFct(x, M) - f(x, M))) print("X-----------") print(f(x, M) * epsilon) stop("differences"); ERR } else cat(". OK") } } cat("\n") }