# file spatial/kr.q copyright (C) 1994-9 W. N. Venables and B. D. Ripley
#
spfmat <- function(x, y, np)
{
    n <- length(x)
    npar <- ((np + 1) * (np + 2))/2
    .C("VR_fmat",
       f = double(n * npar),
       as.double(x),
       as.double(y),
       as.integer(n),
       as.integer(np))$f
}

surf.ls <- function(np, x, y, z)
{
    if(!is.loaded(symbol.C("VR_frset")))
        stop("Compiled code has not been dynamically loaded")
    if (np > 6) stop("np exceeds 6")
    if(is.data.frame(x)) {
        if(missing(y)) y <- x$y
        if(missing(z)) z <- x$z
        x <- x$x
    }
    rx <- range(x)
    ry <- range(y)
    .C("VR_frset",
       as.double(rx[1]),
       as.double(rx[2]),
       as.double(ry[1]),
       as.double(ry[2])
       )
    n <- length(x)
    npar <- ((np + 1) * (np + 2))/2
    f <- spfmat(x, y, np)
    Z <- .C("VR_ls",
            as.double(x),
            as.double(y),
            as.double(z),
            as.integer(n),
            as.integer(np),
            as.integer(npar),
            f = f,
            r = double((npar * (npar + 1))/2),
            beta = double(npar),
            wz = double(n),
            ifail = as.integer(0))
    res <- list(x=x, y=y, z=z, np=np, f=f, r=Z$r, beta=Z$beta,
                wz=Z$wz, rx=rx, ry=ry, call=match.call())
    class(res) <- "trls"
    res
}

surf.gls <- function(np, covmod, x, y, z, nx=1000, ...)
{
    if (np > 6) stop("np exceeds 6")
    if(!is.loaded(symbol.C("VR_frset")))
        stop("Compiled code has not been dynamically loaded")
    if(is.data.frame(x)) {
        if(missing(y)) y <- x$y
        if(missing(z)) z <- x$z
        x <- x$x
    }
    rx <- range(x)
    ry <- range(y)
    .C("VR_frset",
       as.double(rx[1]),
       as.double(rx[2]),
       as.double(ry[1]),
       as.double(ry[2])
       )
    covmod <- covmod
    # replaces
    #  args <- list(...)
    #  if(length(args)) {
    #    pm <- pmatch(names(args), names(covmod), nomatch=0)
    #    if(any(pm == 0)) warning(paste("some of ... do not match"))
    #    covmod[pm] <- unlist(args)
    #  }

    arguments <- list(...)
    if (length(arguments)) {
        onames <- names(formals(covmod))
        pm <- pmatch(names(arguments), onames, nomatch = 0)
        if (any(pm == 0)) warning(paste("some of ... do not match"))
        names(arguments[pm > 0]) <- onames[pm]
        oargs <- formals(covmod)
        oargs[pm] <- arguments[pm > 0]
        formals(covmod) <- oargs
    }
    mm <- 1.5*sqrt((rx[2]-rx[1])^2 + (ry[2]-ry[1])^2)
    alph <- c(mm/nx, covmod(seq(0, mm, mm/nx)))
    .C("VR_alset",
       as.double(alph),
       as.integer(length(alph))
       )
    n <- length(x)
    npar <- ((np + 1) * (np + 2))/2
    f <- spfmat(x, y, np)
    Z <- .C("VR_gls",
            as.double(x),
            as.double(y),
            as.double(z),
            as.integer(n),
            as.integer(np),
            as.integer(npar),
            as.double(f),
            l = double((n * (n + 1))/2),
            r = double((npar * (npar + 1))/2),
            beta = double(npar),
            wz = double(n),
            yy = double(n),
            W = double(n),
            ifail = as.integer(0),
            l1f = double(n * npar)
            )
    if(Z$ifail > 0) stop("Rank failure in Choleski decomposition")
    if(nx > 1000) alph <- alph[1]
    res <- list(x=x, y=y, z=z, np=np, f=f, alph=alph, l=Z$l, r=Z$r,
                beta=Z$beta, wz=Z$wz, yy=Z$yy, W=Z$W, l1f=Z$l1f, rx=rx, ry=ry,
                covmod=covmod, call=match.call())
    class(res) <- c("trgls", "trls")
    res
}

trmat <- function(obj, xl, xu, yl, yu, n)
{
    if(!inherits(obj, "trls")) stop("object not a fitted trend surface")
    .C("VR_frset",
       as.double(obj$rx[1]),
       as.double(obj$rx[2]),
       as.double(obj$ry[1]),
       as.double(obj$ry[2])
       )
    dx <- (xu - xl)/n
    dy <- (yu - yl)/n
    x <- seq(xl, xu, dx)
    y <- seq(yl, yu, dy)
    z <- matrix(nrow = length(x), ncol = length(y))
    for(i in seq(along = y))
        z[, i] <- trval(obj, x, rep(y[i], length(x)))
    invisible(list(x = x, y = y, z = z))
}

trval <- function(obj, x, y)
{
    n <- length(x)
    .C("VR_valn",
       z = double(n),
       as.double(x),
       as.double(y),
       as.integer(n),
       as.double(obj$beta),
       as.integer(obj$np))$z
}

prmat <- function(obj, xl, xu, yl, yu, n)
{
    if(!inherits(obj, "trgls")) stop("object not from kriging")
    if(n > 999) stop("n is too large")
    .C("VR_frset",
       as.double(obj$rx[1]),
       as.double(obj$rx[2]),
       as.double(obj$ry[1]),
       as.double(obj$ry[2])
       )
    alph <- obj$alph
    if(length(alph) <= 1) {
        mm <- 1.5*sqrt((obj$rx[2]-obj$rx[1])^2 + (obj$ry[2]-obj$ry[1])^2)
        alph <- c(alph[1], obj$covmod(seq(0, mm, alph[1])))
    }
    .C("VR_alset",
       as.double(alph),
       as.integer(length(alph))
       )
    dx <- (xu - xl)/n
    dy <- (yu - yl)/n
    xs <- seq(xl, xu, dx)
    ys <- seq(yl, yu, dy)
    z <- matrix(nrow = length(xs), ncol = length(ys))
    for(i in seq(along = ys))
        z[, i] <- trval(obj, xs, rep(ys[i], length(xs))) +
            predval(obj, xs, rep(ys[i], length(xs)))
    invisible(list(x = xs, y = ys, z = z))
}

predval <- function(obj, xp, yp)
{
    npt <- length(xp)
    .C("VR_krpred",
       z = double(npt),
       as.double(xp),
       as.double(yp),
       as.double(obj$x),
       as.double(obj$y),
       as.integer(npt),
       as.integer(length(obj$x)),
       as.double(obj$yy)
       )$z
}

semat <- function(obj, xl, xu, yl, yu, n, se)
{
    if(!inherits(obj, "trgls")) stop("object not from kriging")
    if(n > 999) stop("n is too large")
    .C("VR_frset",
       as.double(obj$rx[1]),
       as.double(obj$rx[2]),
       as.double(obj$ry[1]),
       as.double(obj$ry[2])
       )
    alph <- obj$alph
    if(length(alph) <= 1) {
        mm <- 1.5*sqrt((obj$rx[2]-obj$rx[1])^2 + (obj$ry[2]-obj$ry[1])^2)
        alph <- c(alph[1], obj$covmod(seq(0, mm, alph[1])))
    }
    .C("VR_alset",
       as.double(alph),
       as.integer(length(alph))
       )
    dx <- (xu - xl)/n
    dy <- (yu - yl)/n
    xs <- seq(xl, xu, dx)
    ys <- seq(yl, yu, dy)
    z <- matrix(nrow = length(xs), ncol = length(ys))
    np <- obj$np
    npar <- ((np + 1) * (np + 2))/2
    if(missing(se))
        se <- sqrt(sum(obj$W^2)/(length(obj$x) - npar))
    for(i in seq(along = ys))
        z[, i] <- se * sqrt(seval(obj, xs, rep(ys[i], length(xs))))
    invisible(list(x = xs, y = ys, z = z))
}

seval <- function(obj, xp, yp)
{
    npt <- length(xp)
    np <- obj$np
    npar <- ((np + 1) * (np + 2))/2
    .C("VR_prvar",
       z = double(npt),
       as.double(xp),
       as.double(yp),
       as.integer(npt),
       as.double(obj$x),
       as.double(obj$y),
       as.double(obj$l),
       as.double(obj$r),
       as.integer(length(obj$x)),
       as.integer(np),
       as.integer(npar),
       as.double(obj$l1f))$z
}

correlogram <- function(krig, nint, plotit=TRUE, ...)
{
    if(!is.loaded(symbol.C("VR_correlogram")))
        stop("Compiled code has not been dynamically loaded")
    z <- .C("VR_correlogram",
            xp = double(nint),
            yp = double(nint),
            nint = as.integer(nint),
            as.double(krig$x),
            as.double(krig$y),
            if(krig$np > 0) as.double(krig$wz) else as.double(krig$z),
            as.integer(length(krig$x)),
            cnt = integer(nint)
            )
    xp <- z$xp[1:z$nint]
    yp <- z$yp[1:z$nint]
    z <- list(x = xp, y = yp, cnt = z$cnt[1:z$nint])
    if(plotit)
        if(exists(".Device")) {
            plot(xp, yp, type = "p", ylim = c(-1, 1), ...)
            abline(0, 0)
            invisible(z)
        }
        else {
            warning("Device not active")
            return(z)
        }
    else z
}

variogram <- function(krig, nint, plotit=TRUE, ...)
{
    if(!is.loaded(symbol.C("VR_variogram")))
        stop("Compiled code has not been dynamically loaded")
    z <- .C("VR_variogram",
            xp = double(nint),
            yp = double(nint),
            nint = as.integer(nint),
            as.double(krig$x),
            as.double(krig$y),
            if(krig$np > 0) as.double(krig$wz) else as.double(krig$z),
            as.integer(length(krig$x)),
            cnt = integer(nint)
            )
    xp <- z$xp[1:z$nint]
    yp <- z$yp[1:z$nint]
    if(xp[1] > 0) {xp <- c(0, xp); yp <- c(0, yp)}
    z <- list(x = xp, y = yp, cnt = z$cnt[1:z$nint])
    if(plotit)
        if(exists(".Device")) {
            plot(xp, yp, type = "p", ...)
            invisible(z)
        }
        else {
            warning("Device not active")
            return(z)
        }
    else z
}

expcov <- function(r, d, alpha=0, se=1)
{
    se^2*(alpha*(r < d/10000) + (1-alpha)*exp(-r/d))
}

gaucov <- function(r, d, alpha=0, se=1)
{
    se^2*(alpha*(r < d/10000) + (1-alpha)*exp(-(r/d)^2))
}

sphercov <- function(r, d, alpha=0, se=1, D=2)
{
    r <- r/d
    if(D==2) {
        t <- 1 - (2/pi)*(r*sqrt(1-r^2) + asin(r))
    } else {
        t <- 1 - 1.5*r + r^3/2
    }
    se^2*(alpha*(r < 1/10000) + (1-alpha)*t*(r < 1))
}
# file spatial/pp.q copyright (C) 1994-9 W. N. Venables and B. D. Ripley
#
ppinit <- function(file)
{
    tfile <- file
    if(exists(".sp.lib.name", envir=.GlobalEnv)) {
        t1file <- paste(get(".sp.lib.name", envir=.GlobalEnv), file, sep="/")
        if(file.exists(t1file)) tfile <- t1file
    }
    h <- scan(tfile, list(xl = 0, xu = 0, yl = 0, yu = 0, fac = 0),
              n = 5, skip = 2, quiet = TRUE)
    pp <- scan(tfile, list(x = 0, y = 0), skip = 3, quiet = TRUE)
    pp$x <- pp$x/h$fac
    pp$y <- pp$y/h$fac
    pp$area <- c(xl=h$xl/h$fac, xu=h$xu/h$fac, yl=h$yl/h$fac, yu=h$yu/h$fac)
    if(!is.loaded(symbol.C("VR_ppset")))
        stop("Compiled code has not been dynamically loaded")
    ppregion(pp)
    invisible(pp)
}

Kfn <- function(pp, fs, k = 100)
{
    zz <- (c(range(pp$x), range(pp$y)) - ppgetregion())*c(1,-1,1,-1)
    if(any(zz < 0)) stop("some points outside region")
    z <- .C("VR_sp_pp2",
            as.double(pp$x),
            as.double(pp$y),
            as.integer(length(pp$x)),
            k1 = as.integer(k),
            h = double(k),
            dmin = double(1),
            lm = double(1),
            as.double(fs))
    list(y = z$h[1:z$k1], x = (seq(1:z$k1) * fs)/k, k = k,
         dmin = z$dmin, lm = max(z$dmin, z$lm),
         call=match.call())
}

Kenvl <- function(fs, nsim, ...)
{
    dot.expression <- as.expression(substitute(...))
    h <- Kfn(pp = eval(dot.expression), fs)
    hx <- h$x
    hu <- h$y
    hl <- h$y
    ha <- h$y^2
    for(i in 2:nsim) {
        h <- Kfn(pp = eval(dot.expression), fs)$y
        hu <- pmax(hu, h)
        hl <- pmin(hl, h)
        ha <- ha + h^2
    }
    list(x = hx, lower = hl, upper = hu, aver = sqrt(ha/nsim),
         call=match.call())
}

Kaver <- function(fs, nsim, ...)
{
    dot.expression <- as.expression(substitute(...))
    h <- Kfn(pp = eval(dot.expression), fs)
    hx <- h$x
    ha <- h$y^2
    for(i in 2:nsim) {
        h <- Kfn(pp = eval(dot.expression), fs)$y
        ha <- ha + h^2
    }
    list(x = hx, y = sqrt(ha/nsim), call=match.call())
}

ppregion <- function(xl = 0, xu = 1, yl = 0, yu = 1)
{
    if(is.list(xl)) {
        if(length(xl$area)) .C("VR_ppset", as.double(xl$area))
        else .C("VR_ppset", as.double(c(xl$xl, xl$xu, xl$yl, xl$yu)))
    } else .C("VR_ppset", as.double(c(xl, xu, yl, yu)))
    invisible()
}

ppgetregion <- function()
{
    xx <- .C("VR_ppget", z=double(4))$z
    names(xx) <- c("xl", "xu", "yl", "yu")
    xx
}

Psim <- function(n)
{
    z <- .C("VR_pdata",
            as.integer(n),
            x = double(n),
            y = double(n))
    invisible(list(x = z$x, y = z$y, call=match.call()))
}

Strauss <- function(n, c = 0, r)
{
    init <-  0
    if(!exists(".ppx")) {
        init <-  1
        z <- .C("VR_pdata",
                as.integer(n),
                x = double(n),
                y = double(n))
        assign(".ppx", z$x)
        assign(".ppy", z$y)
    }
    z <- .C("VR_simpat",
            as.integer(n),
            x = as.double(.ppx),
            y = as.double(.ppy),
            as.double(c),
            as.double(r),
            as.integer(init))
    assign(".ppx", z$x)
    assign(".ppy", z$y)
    invisible(list(x = z$x, y = z$y, call=match.call()))
}

SSI <- function(n, r)
{
    z <- .C("VR_simmat",
            as.integer(n),
            x = double(n),
            y = double(n),
            as.double(r))
    invisible(list(x = z$x, y = z$y, call=match.call()))
}

pplikfn <- function(cc, R, n, x, y, ng, target, trace=FALSE)
{
    z <- .C("VR_plike",
            as.double(x),
            as.double(y),
            as.integer(n),
            as.double(cc),
            as.double(R),
            as.integer(ng),
            as.double(target),
            res=double(1)
            )
    if(trace) print(c(cc, z$res))
    z$res
}

pplik <- function(pp, R, ng=50, trace=FALSE)
{
    n <- length(pp$x)
    ar <- pp$area
    target <- n * (Kfn(pp, R,1)$y)^2 * pi /
        ((ar["xu"] - ar["xl"]) * (ar["yu"] - ar["yl"]))
    if(target == 0) return(0)
    tmp <- pplikfn(1, R, n, pp$x, pp$y, ng, target, FALSE)
    if(tmp <= 0) return(1)
    uniroot(pplikfn, c(0,1),
            R=R, n=n, x=pp$x, y=pp$y, ng=ng, target=target,
            trace=trace)$root
}

.First.lib <- function(lib, pkg)
{
    if(version$major==0 && version$minor < 64)
        stop("This version for R 0.64 or later")
    library.dynam("spatial", pkg, lib)
    assign(".sp.lib.name", paste(lib, pkg, "data", sep="/"), envir=.GlobalEnv)
}

