"CovarianceFct" <- function (x, model, param, dim = 1) { if (!is.finite(param[1])) { param[1] <- 0 } stopifnot(length(x)!=0, all(is.finite(x)), is.character(model), all(is.finite(param)), length(dim) == 1, is.finite(dim) ) covnr <- as.integer(.C("GetModelNr", as.character(model), nr = integer(1))$nr) if (covnr < 0) { .C("PrintModelList") stop("given model cannot be (uniquely) identified from the above list") } kappas <- as.integer(.C("GetNrParameters", covnr, k = integer(1),DUP=FALSE)$k) if (length(param) < 4 + kappas) stop("not enough parameters for the covariance function") storage.mode(x) <- "double" storage.mode(param) <- "double" storage.mode(dim) <- "integer" return(.C("Covariance", x, as.integer(length(x)), covnr, param, as.integer(length(param)), dim, res=double(length(x)), DUP = FALSE)$res ) } "DeleteAllRegisters" <- function () { .C("DeleteAllKeys") return(NULL) } "DeleteRegister" <- function (nr = 0) { stopifnot( length(nr) == 1, is.finite(nr) ) .C("DeleteKey", as.integer(nr)) return(NULL) } "DoSimulateRF" <- function (n = 1, register = 0) { stopifnot(length(n) == 1, is.finite(n), length(register) == 1, is.finite(register) ) DoNameList <- c("DoSimulateRF","DoSimulateRF","DoMaxStableRF") MAXDIM <- 3 dimensions <- .C("GetKeyInfo", as.integer(register), total = integer(1), len = integer(MAXDIM), dim = integer(1), grid = integer(1), distr = integer(1), DUP=FALSE) param <- RFparameters() if (dimensions$total <= 0) { stop(paste("register", register, "does not look initialized")) } if (n == 1) { x <- .C(DoNameList[1+dimensions$distr], as.integer(register), result=double(dimensions$total), error=integer(1),DUP=FALSE) if (x$error) return(NULL) else if (dimensions$grid) return(array(x$result, dim = dimensions$len[1:dimensions$dim])) else return(x$result) } else { if (!param$Storing) stop("Use option `RFparameters(Storing=TRUE)' if n>1") if (dimensions$grid) { result <- array(NA, dim = c(dimensions$len[1:dimensions$dim], n)) s <- paste(rep(",",dimensions$dim),collapse="") } else { result <- matrix(nrow = dimensions$total, ncol = n) s <- "," } s <- paste("result[", s, "i] <- dummy$res") pch <- param$pch for (i in 1:n) { cat(pch) dummy <- .C(DoNameList[1+dimensions$distr], as.integer(register), res=double(dimensions$total), error=integer(1), DUP = FALSE, NAOK = TRUE) if (dummy$error) break eval(parse(text = s)) } if (pch!="") cat("\n") if (dummy$error) return(NULL) else return(result) } } "EmpiricalVariogram" <- function (x, y = NULL, z = NULL, data, grid, bin, gridtriple = FALSE) { stopifnot(all(is.finite(data)), length(bin)>=2, all(is.finite(bin)), bin[1] <= 0 ) new <- CheckXYZ(x, y, z, grid, gridtriple) repet <- as.integer(length(data)/new$total) if (length(data) != new$total * repet) stop("number of data does not match coordinates") centers <- pmax(0,(bin[-1] + bin[-length(bin)])/2) n.bins <- length(bin) - 1 emp.vario <- double(n.bins) n.bin <- integer(n.bins) .C("empiricalvariogram", as.double(new$x), as.double(new$y), as.double(new$z), as.integer(new$dim), as.integer(new$l), as.double(data), as.integer(repet), as.integer(grid), as.double(bin), as.integer(n.bins), as.integer(0), emp.vario, n.bin, DUP = FALSE) emp.vario[is.na(emp.vario) & (centers==0)] <- 0 return(list(centers=centers, emp.vario=emp.vario, n.bin=n.bin)) } ".First.lib" <- function (lib, pkg) { library.dynam("RandomFields", pkg, lib) library(mva) } "InitGaussRF" <- function(x, y = NULL, z = NULL, grid, model, param, method = NULL, register = 0, gridtriple = FALSE) { return(InitSimulateRF(x=x, y=y, z=z, grid=grid, model=model, param=param, method=method, register=register, gridtriple=gridtriple, distribution="Gauss")) } CheckXYZ <- function(x, y, z, grid, gridtriple){ stopifnot(length(x) != 0, is.logical(grid), is.logical(gridtriple) ) if (is.data.frame(x)) { if (ncol(x)==1) x <- as.vector(x) else x <- as.matrix(x) } stopifnot(all(is.finite(x))) if (is.matrix(x)) { if (!is.null(y) || !is.null(z)) { stop("If x is a matrix, then y and z may not be given") } dim <- ncol(x) if (dim > 3) { dim <- 3 warning("only the first three columns are considered as coordinates") } if (dim > 1) { y <- x[, 2] if (dim > 2) z <- x[, 3] } x <- x[, 1] l <- length(x) } else { l <- length(x) storage.mode(x) <- "double" if (is.null(y)) { dim <- 1 if (!is.null(z)) { stop("y is not given, but z") } } else { if ((l != length(y)) && (!grid || gridtriple)) { stop("x and y differ in length") } stopifnot(all(is.finite(y))) storage.mode(y) <- "double" if (is.null(z)) { dim <- 2 } else { dim <- 3 if (l != length(z) && (!grid || gridtriple)) { stop("x and z differ in length") } stopifnot(all(is.finite(z))) storage.mode(z) <- "double" } } } if (l == 1) stop("Use grid=FALSE if only a single point is simulated") else if (grid) { if (gridtriple) { if (l != 3) { stop("In case of simulating a grid with option gridtriple, exactly 3 numbers are needed for each direction") } lx <- length(seq(x[1],x[2],x[3])) ly <- lz <- 1 x[2] <- x[1] + (lx - 0.999) * x[3] if (dim > 1) { ly <- length(seq(y[1],y[2],y[3])) y[2] <- y[1] + (ly - 0.999) * y[3] if (dim > 2) { lz <- length(seq(z[1],z[2],z[3])) z[2] <- z[1] + (lz - 0.999) * z[3] } } total <- lx * ly * lz if (total==0) stop("incorrect grid specification (one or no points)") } else { eqdist <- function(x) { step <- diff(x) if (sum(abs(step - step[1]))/step[1] > 1e-10) { stop("Grid must have equal distances in each direction.") } return(c(x[1], x[length(x)]+0.001*step[1], step[1])) } total <- length(x) x <- eqdist(x) if (dim > 1) { total <- total * length(y) y <- eqdist(y) if (dim > 2) { total <- total * length(z) z <- eqdist(z) } } l <- 3 } } else { total <- length(x) } return(list(x=x, y=y, z=z, total=total, l=l, dim=dim)) } "InitSimulateRF" <- function (x, y = NULL, z = NULL, grid, model, param, method = NULL, register = 0, gridtriple = FALSE,distribution=NA) { distributionList <- c("Gauss","Poisson","MaxStable") InitNameList <- c("InitSimulateRF","InitSimulateRF","InitMaxStableRF") if (is.na(distribution)) { stop("This function is an internal function.\nUse `GaussRF', `InitGaussRF', `MaxStableRF', etc., instead.\n") } distrNr <- pmatch(distribution,distributionList) if (is.na(distrNr)) stop("Unknown distribution -- do not use this internal function directly") else { if ((distrNr==2) && !exists("PoissonRF")) stop("Sorry. Not programmed yet.") InitName <- InitNameList[distrNr] distrNr <- distrNr - 1; } stopifnot(is.character(model)) covnr <- as.integer(.C("GetModelNr", as.character(model), nr = integer(1))$nr) if (covnr < 0) { .C("PrintModelList") stop("given model cannot be (uniquely) identified from the above list") } if (is.null(method)) { method <- -1 } else { if (!is.character(method)) { stop("method must be NULL or a string") } method <- as.integer(.C("GetMethodNr", method, nr = integer(1))$nr) if (method < 0) { .C("PrintMethods") stop("given method cannot be (uniquely) identified from the above list") } } kappas <- as.integer(.C("GetNrParameters", covnr, k = integer(1))$k) if (length(param) < 4 + kappas) stop("not enough parameters for the covariance function") storage.mode(param) <- "double" new <- CheckXYZ(x, y, z, grid, gridtriple) return( max(0,.C(InitName, as.double(new$x), as.double(new$y), as.double(new$z), as.integer(new$dim), as.integer(new$l), as.integer(grid), covnr, param, as.integer(length(param)), as.integer(method), as.integer(distrNr), as.integer(register), error=integer(1), DUP = FALSE)$error) ) } "PrintModelList" <- function () { .C("PrintModelList") return(NULL) } "PrintMethodList" <- function () { .C("PrintMethods") return(NULL) } "RFparameters" <-function (...) { RFparameters.default(...) } "RFparameters.default" <- function (Storing = storing, PrintLevel = printlevel, PracticalRange = practicalrange, CE.force = ce.force, CE.mmin = ce.mmin, CE.tolRe = ce.tolRe, CE.tolIm = ce.tolIm, CE.trials = ce.trials, direct.checkprecision = directcheckprecision, direct.maxvariables = directmaxvariables, direct.method = directmethod, direct.requiredprecision = directrequiredprecision, spectral.lines = spectrallines, spectral.grid = spectralgrid, TBMCE.force = tbmceforce, TBMCE.mmin = tbmcemmin, TBMCE.tolRe = tbmcetolre, TBMCE.tolIm = tbmcetolim, TBMCE.trials = tbmcetrials, TBM2.lines = tbm2lines, TBM2.linesimufactor = tbm2linesimufactor, TBM2.linesimustep = tbm2linesimustep, TBM3D2.lines = tbm3D2lines, TBM3D2.linesimufactor = tbm3D2linesimufactor, TBM3D2.linesimustep = tbm3D2linesimustep, TBM3D3.lines = tbm3D3lines, TBM3D3.linesimufactor = tbm3D3linesimufactor, TBM3D3.linesimustep = tbm3D3linesimustep, MPP.approxzero=mppapproxzero, add.MPP.realisations=addmpprealisations, MPP.radius=mppradius, maxstable.maxGauss=maxstablemaxGauss, pch = pchx ) { x <- .C("SetParam", as.integer(1), storing = integer(1), printlevel = integer(1), practicalrange = integer(1), pch=" ") storing <- x$storing printlevel <- x$printlevel practicalrange <- as.logical(x$practicalrange) if (!is.finite(Storing + PrintLevel + PracticalRange)) stop("some parameters are not finite") pchx <- x$pch if (!is.character(pch)) stop("pch is not a character") ## do not allow integer values for users! PracticalRange <- as.logical(PracticalRange) .C("SetParam", as.integer(0), as.integer(Storing), as.integer(PrintLevel), as.integer(PracticalRange), pch) x <- .C("SetParamCircEmbed", as.integer(1), force = integer(1), tolRe = double(1), tolIm = double(1), trials = integer(1), mmin = integer(1)) ce.force <- x$force ce.tolRe <- x$tolRe ce.tolIm <- x$tolIm ce.trials <- x$trials ce.mmin <- x$mmin if (!is.finite(CE.force + CE.mmin + CE.tolRe + CE.tolIm + CE.trials)) stop("some parameters are not finite") .C("SetParamCircEmbed", as.integer(0), as.integer(CE.force), as.double(CE.tolRe), as.double(CE.tolIm), as.integer(CE.trials), as.integer(CE.mmin)) x <- .C("SetParamTBMCE", as.integer(1), tbmceforce = integer(1), tbmcetolre = double(1), tbmcetolim = double(1), tbmcetrials = integer(1), tbmcemmin = integer(1)) tbmceforce <- x$tbmceforce tbmcetolre <- x$tbmcetolre tbmcetolim <- x$tbmcetolim tbmcetrials <- x$tbmcetrials tbmcemmin <- x$tbmcemmin if (!is.finite(TBMCE.force + TBMCE.mmin + TBMCE.tolRe + TBMCE.tolIm + TBMCE.trials)) stop("some parameters are not finite") .C("SetParamTBMCE", as.integer(0), as.integer(TBMCE.force), as.double(TBMCE.tolRe), as.double(TBMCE.tolIm), as.integer(TBMCE.trials), as.integer(TBMCE.mmin)) x <- .C("SetParamTBM2", as.integer(1), tbm2lines = integer(1), tbm2linesimufactor = double(1), tbm2linesimustep = double(1)) tbm2lines <- x$tbm2lines tbm2linesimufactor <- x$tbm2linesimufactor tbm2linesimustep <- x$tbm2linesimustep if (!is.finite(TBM2.lines + TBM2.linesimufactor + TBM2.linesimustep)) stop("some parameters are not finite") .C("SetParamTBM2", as.integer(0), as.integer(TBM2.lines), as.double(TBM2.linesimufactor), as.double(TBM2.linesimustep)) x <- .C("SetParamTBM3D2", as.integer(1), tbm3D2lines = integer(1), tbm3D2linesimufactor = double(1), tbm3D2linesimustep = double(1)) tbm3D2lines <- x$tbm3D2lines tbm3D2linesimufactor <- x$tbm3D2linesimufactor tbm3D2linesimustep <- x$tbm3D2linesimustep if (!is.finite(TBM3D2.lines + TBM3D2.linesimufactor + TBM3D2.linesimustep)) stop("some parameters are not finite") .C("SetParamTBM3D2", as.integer(0), as.integer(TBM3D2.lines), as.double(TBM3D2.linesimufactor), as.double(TBM3D2.linesimustep)) x <- .C("SetParamTBM3D3", as.integer(1), tbm3D3lines = integer(1), tbm3D3linesimufactor = double(1), tbm3D3linesimustep = double(1)) tbm3D3lines <- x$tbm3D3lines tbm3D3linesimufactor <- x$tbm3D3linesimufactor tbm3D3linesimustep <- x$tbm3D3linesimustep if (!is.finite(TBM3D3.lines + TBM3D3.linesimufactor + TBM3D3.linesimustep)) stop("some parameters are not finite") .C("SetParamTBM3D3", as.integer(0), as.integer(TBM3D3.lines), as.double(TBM3D3.linesimufactor), as.double(TBM3D3.linesimustep)) x <- .C("SetParamSpectral", as.integer(1), spectrallines = integer(1), spectralgrid = integer(1)) spectrallines <- x$spectrallines spectralgrid <- x$spectralgrid if (!is.finite(spectral.lines + spectral.grid)) stop("some parameters are not finite") .C("SetParamSpectral", as.integer(0), as.integer(spectral.lines), as.integer(spectral.grid)) x <- .C("SetParamDirectGauss", as.integer(1), directmethod = integer(1), directcheckprecision = integer(1), directrequiredprecision = double(1), directmaxvariables = integer(1)) directmethod <- x$directmethod directcheckprecision <- x$directcheckprecision directrequiredprecision <- x$directrequiredprecision directmaxvariables <- x$directmaxvariables if (!is.finite(direct.checkprecision + direct.maxvariables + direct.method + direct.requiredprecision)) stop("some parameters are not finite") .C("SetParamDirectGauss", as.integer(0), as.integer(direct.method), as.integer(direct.checkprecision), as.double(direct.requiredprecision), as.integer(direct.maxvariables)) x <- .C("SetMPP", as.integer(1), mppapproxzero=double(1), addmpprealisations=double(1), mppradius=double(1)) mppapproxzero <- x$mppapproxzero; addmpprealisations <- x$addmpprealisations; mppradius <- x$mppradius; if (!is.finite(MPP.approxzero+add.MPP.realisations+MPP.radius)) stop("some parameters are not finite") .C("SetMPP", as.integer(0), as.double(MPP.approxzero),as.double(add.MPP.realisations), as.double(MPP.radius)) x <- .C("SetExtremes", as.integer(1), maxstablemaxGauss=double(1)) maxstablemaxGauss <- x$maxstablemaxGauss; if (!is.finite(maxstable.maxGauss)) stop("some parameters are not finite") .C("SetExtremes", as.integer(0), as.double(maxstable.maxGauss)) if (length(as.list(match.call()))>1) return(NULL) else return(list(Storing = as.logical(Storing), PrintLevel = PrintLevel, PracticalRange = as.logical(PracticalRange), CE.force = as.logical(CE.force), CE.mmin = CE.mmin, CE.tolRe = CE.tolRe, CE.tolIm = CE.tolIm, CE.trials = CE.trials, direct.checkprecision = as.logical(direct.checkprecision), direct.maxvariables = direct.maxvariables, direct.method = direct.method, direct.requiredprecision = direct.requiredprecision, spectral.lines = spectral.lines, spectral.grid = as.logical(spectral.grid), TBMCE.force = as.logical(TBMCE.force), TBMCE.mmin = TBMCE.mmin, TBMCE.tolRe = TBMCE.tolRe, TBMCE.tolIm = TBMCE.tolIm, TBMCE.trials = TBMCE.trials, TBM2.lines = TBM2.lines, TBM2.linesimufactor = TBM2.linesimufactor, TBM2.linesimustep = TBM2.linesimustep, TBM3D2.lines = TBM3D2.lines, TBM3D2.linesimufactor = TBM3D2.linesimufactor, TBM3D2.linesimustep = TBM3D2.linesimustep, TBM3D3.lines = TBM3D3.lines, TBM3D3.linesimufactor = TBM3D3.linesimufactor, TBM3D3.linesimustep = TBM3D3.linesimustep, MPP.approxzero=MPP.approxzero, add.MPP.realisations=add.MPP.realisations, MPP.radius=MPP.radius, maxstable.maxGauss=maxstable.maxGauss, pch=pch )) } "GaussRF" <- function (x, y = NULL, z = NULL, grid, model, param, method = NULL, n = 1, register = 0, gridtriple = FALSE) { if (InitSimulateRF(x=x, y=y, z=z, grid=grid, model=model, param=param, method=method, register=register, gridtriple=gridtriple, distribution="Gauss") <= 0) { return(DoSimulateRF(n=n, reg=register)) } else { return(NULL) } } "SimulateRF" <- function (x, y = NULL, z = NULL, grid, model, param, method = NULL, n = 1, register = 0, gridtriple = FALSE, distribution=NA) { if (InitSimulateRF(x=x, y=y, z=z, grid=grid, model=model, param=param, method=method, reg=register, gridtriple=gridtriple, distr=distribution) <= 0) { return(DoSimulateRF(n=n, reg=register)) } else { return(NULL) } } "Variogram" <- function (x, model, param, dim = 1) { if (!is.finite(param[1])) param[1] <- 0 stopifnot(length(x) != 0, all(is.finite(x)), is.character(model), all(is.finite(param)), length(dim) == 1, is.finite(dim)) covnr <- as.integer(.C("GetModelNr", as.character(model), nr = integer(1))$nr) if (covnr < 0) { .C("PrintModelList") stop("given model cannot be (uniquely) identified from the above list") } kappas <- as.integer(.C("GetNrParameters", covnr, k = integer(1))$k) if (length(param) < 4 + kappas) stop("not enough parameters for the variogram model") storage.mode(x) <- "double" storage.mode(param) <- "double" storage.mode(dim) <- "integer" return( .C("Variogram", x, as.integer(length(x)), covnr, param, as.integer(length(param)), dim, res=double(length(x)), DUP = FALSE)$res ) } GetModelNames <- function() { p <- .C("GetrfParameters",covmaxchar=integer(1),methodmaxchar=integer(1), covnr=integer(1),methodnr=integer(1)) l <- character(p$covnr) for (i in 1:p$covnr) { l[i] <- .C("GetModelName",as.integer(i-1), n=paste(rep(" ",p$covmaxchar),collapse=""))$n } return (l) } GetMethodNames <- function() { p <- .C("GetrfParameters",covmaxchar=integer(1),methodmaxchar=integer(1), covnr=integer(1),methodnr=integer(1)) l <- character(p$methodnr) for (i in 1:p$methodnr) { l[i] <- .C("GetMethodName",as.integer(i-1), n=paste(rep(" ",p$methodmaxchar),collapse=""))$n } return(l) } GetPracticalRange <- function(model,kappas=NULL) { covnr <- as.integer(.C("GetModelNr", as.character(model), nr = integer(1))$nr) if (covnr < 0) { .C("PrintModelList") stop("given model cannot be (uniquely) identified from the above list") } if (length(kappas)!=.C("GetNrParameters", covnr, k = integer(1),DUP=FALSE)$k) stop("incorrect number of parameters!") nat.scl <- double(1) error <- integer(1) .C("GetNaturalScaling", as.integer(covnr), as.double(c(0,1,0,1,kappas)), ## not stable w.r.t. to changings !! as.integer(4+ length(kappas)), ## dito as.integer(11), nat.scl, error,DUP=FALSE) if (error) stop("natural scaling could not be obtained") return(1.0 / nat.scl) }