Nothing
################################################################################
### Spatial interaction functions for twinstim's epidemic component.
### Specific implementations are in separate files (e.g.: Gaussian, power law).
###
### Copyright (C) 2009-2015,2017 Sebastian Meyer
###
### This file is part of the R package "surveillance",
### free software under the terms of the GNU General Public License, version 2,
### a copy of which is available at https://www.R-project.org/Licenses/.
################################################################################
#####################
### "Constructor" ###
#####################
siaf <- function (f, F, Fcircle, effRange, deriv, Deriv, simulate, npars,
validpars = NULL)
{
npars <- as.integer(npars)
if (length(npars) != 1 || npars < 0L) {
stop("'siaf$npars' must be a single nonnegative number")
}
f <- .checknargs3(f, "siaf$f")
F <- if (missing(F) || is.null(F)) siaf.fallback.F else {
F <- match.fun(F)
if (length(formals(F)) < 4L)
stop("siaf$F() must accept >=4 arguments ",
"(polydomain, f, pars, type)")
F
}
haspars <- npars > 0L
if (!haspars || missing(deriv)) deriv <- NULL
if (!is.null(deriv)) deriv <- .checknargs3(deriv, "siaf$deriv")
if (missing(effRange)) effRange <- NULL
if (missing(Fcircle) || is.null(Fcircle)) {
Fcircle <- NULL
if (!is.null(effRange)) {
message("'siaf$effRange' only works in conjunction with 'siaf$Fcircle'")
effRange <- NULL
}
}
if (!is.null(Fcircle)) Fcircle <- .checknargs3(Fcircle, "siaf$Fcircle")
if (!is.null(effRange)) {
effRange <- match.fun(effRange)
if (length(formals(effRange)) < 1L) {
stop("the 'siaf$effRange' function must accept a parameter vector")
}
}
Deriv <- if (is.null(deriv)) NULL else if (missing(Deriv) || is.null(Deriv))
siaf.fallback.Deriv else {
Deriv <- match.fun(Deriv)
if (length(formals(Deriv)) < 4L)
stop("siaf$Deriv() must accept >=4 arguments ",
"(polydomain, deriv, pars, type)")
Deriv
}
## Check if simulation function has proper format
if (missing(simulate)) simulate <- NULL
if (!is.null(simulate)) {
simulate <- .checknargs3(simulate, "siaf$simulate")
if (length(formals(simulate)) == 3L)
formals(simulate) <- c(formals(simulate), alist(ub=))
}
## Check if the validpars are of correct form
validpars <- if (!haspars || is.null(validpars))
NULL else match.fun(validpars)
## Done, return result.
list(f = f, F = F, Fcircle = Fcircle, effRange = effRange,
deriv = deriv, Deriv = Deriv,
simulate = simulate,
npars = npars, validpars = validpars)
}
##########################################
### Constant spatial interaction/dispersal
##########################################
siaf.constant <- function ()
{
res <- list(
## use explicit quote()ing to prevent notes from codetools::checkUsage
f = as.function(c(alist(s=, pars=NULL, types=NULL),
quote(rep.int(1, length(s)/2))),
##<- nrow() would take extra time in standardGeneric()
envir = .GlobalEnv),
## integration over polydomains (F) is handled specially in twinstim
Fcircle = as.function(c(alist(r=, pars=NULL, type=NULL),
quote(pi*r^2)),
envir = .GlobalEnv),
simulate = as.function(c(alist(n=, pars=NULL, type=NULL, ub=),
quote(runifdisc(n, ub))),
envir = getNamespace("surveillance")),
npars = 0L
)
attr(res, "constant") <- TRUE
res
}
##########################################
### Naive defaults for the siaf primitives
##########################################
## numerical integration of f over a polygonal domain (single "owin" and type)
siaf.fallback.F <- function (polydomain, f, pars, type, method = "SV", ...)
{
if (identical(method,"SV")) {
polyCub.SV(polyregion = polydomain, f = f, pars, type,
alpha = 0, ...) # since max at origin
} else {
polyCub(polyregion = polydomain, f = f, method = method,
pars, type, ...)
}
}
## numerical integration of f over a circular domain
getFcircle <- function (siaf, control.F = list()) {
if (is.null(siaf$Fcircle)) {
function (r, pars, type) {
disc <- discpoly(c(0,0), r, npoly = 64, class = "owin")
do.call(siaf$F, c(alist(disc, siaf$f, pars, type), control.F))
}
} else {
siaf$Fcircle
}
}
## numerical integration of deriv over a polygonal domain
siaf.fallback.Deriv <- function (polydomain, deriv, pars, type,
method = "SV", ...)
{
deriv1 <- function (s, paridx)
deriv(s, pars, type)[,paridx,drop=TRUE]
intderiv1 <- function (paridx)
polyCub(polyregion = polydomain, f = deriv1, method = method,
paridx = paridx, ...)
vapply(X = seq_along(pars), FUN = intderiv1,
FUN.VALUE = 0, USE.NAMES = FALSE)
}
####################################
### Simulation via polar coordinates (used, e.g., for siaf.powerlaw)
####################################
## Simulate from an isotropic spatial interaction function
## f_{2D}(s) \propto f(||s||), ||s|| <= ub.
## within a maximum distance 'ub' via polar coordinates and the inverse
## transformation method:
## p_{2D}(r,theta) = r * f_{2D}(x,y) \propto r*f(r)
## => angle theta ~ U(0,2*pi) and sample r according to r*f(r)
siaf.simulatePC <- function (intrfr) # e.g., intrfr.powerlaw
{
as.function(c(alist(n=, siafpars=, type=, ub=), substitute({
## Note: in simEpidataCS, simulation is always bounded to eps.s and to
## the largest extend of W, thus, 'ub' is finite
stopifnot(is.finite(ub))
## Normalizing constant of r*f(r) on [0;ub]
normconst <- intrfr(ub, siafpars, type)
## => cumulative distribution function
CDF <- function (q) intrfr(q, siafpars, type) / normconst
## For inversion sampling, we need the quantile function CDF^-1
## However, this is not available in closed form, so we use uniroot
## (which requires a finite upper bound)
QF <- function (p) uniroot(function(q) CDF(q)-p, lower=0, upper=ub)$root
## Now sample r as QF(U), where U ~ U(0,1)
r <- vapply(X=runif(n), FUN=QF, FUN.VALUE=0, USE.NAMES=FALSE)
## Check simulation of r via kernel estimate:
## plot(density(r, from=0, to=ub)); curve(p(x)/normconst,add=TRUE,col=2)
## now rotate each point by a random angle to cover all directions
theta <- runif(n, 0, 2*pi)
r * cbind(cos(theta), sin(theta))
})), envir=parent.frame())
}
################################################
### Check F, Fcircle, deriv, Deriv, and simulate
################################################
checksiaf <- function (siaf, pargrid, type = 1, tolerance = 1e-5,
method = "SV", ...)
{
stopifnot(is.list(siaf), is.numeric(pargrid), !is.na(pargrid),
length(pargrid) > 0)
pargrid <- as.matrix(pargrid)
stopifnot(siaf$npars == ncol(pargrid))
## Check 'F'
if (!is.null(siaf$F)) {
cat("'F' vs. cubature using method = \"", method ,"\" ... ", sep="")
comp.F <- checksiaf.F(siaf$F, siaf$f, pargrid, type=type,
method=method, ...)
cat(attr(comp.F, "all.equal") <-
all.equal(comp.F[,1], comp.F[,2],
check.attributes=FALSE, tolerance=tolerance),
"\n")
}
## Check 'Fcircle'
if (!is.null(siaf$Fcircle)) {
cat("'Fcircle' vs. cubature using method = \"",method,"\" ... ", sep="")
comp.Fcircle <- checksiaf.Fcircle(siaf$Fcircle, siaf$f, pargrid,
type=type, method=method, ...)
cat(attr(comp.Fcircle, "all.equal") <-
all.equal(comp.Fcircle[,1], comp.Fcircle[,2],
check.attributes=FALSE, tolerance=tolerance),
"\n")
}
## Check 'deriv'
if (!is.null(siaf$deriv)) {
cat("'deriv' vs. numerical derivative ... ")
if (requireNamespace("maxLik", quietly=TRUE)) {
maxRelDiffs.deriv <- checksiaf.deriv(siaf$deriv, siaf$f, pargrid,
type=type)
cat(attr(maxRelDiffs.deriv, "all.equal") <-
if (any(maxRelDiffs.deriv > tolerance))
paste("maxRelDiff =", max(maxRelDiffs.deriv)) else TRUE,
"\n")
} else cat("Failed: need package", sQuote("maxLik"), "\n")
}
## Check 'Deriv'
if (!is.null(siaf$Deriv)) {
cat("'Deriv' vs. cubature using method = \"", method ,"\" ... ", sep="")
comp.Deriv <- checksiaf.Deriv(siaf$Deriv, siaf$deriv, pargrid,
type=type, method=method, ...)
if (siaf$npars > 1) cat("\n")
attr(comp.Deriv, "all.equal") <-
sapply(seq_len(siaf$npars), function (j) {
if (siaf$npars > 1) cat("\tsiaf parameter ", j, ": ", sep="")
ae <- all.equal(comp.Deriv[,j], comp.Deriv[,siaf$npars+j],
check.attributes=FALSE, tolerance=tolerance)
cat(ae, "\n")
ae
})
}
## Check 'simulate'
if (interactive() && !is.null(siaf$simulate)) {
cat("Simulating ... ")
checksiaf.simulate(siaf$simulate, siaf$f, pargrid[1,], type=type)
cat("(-> check the plot)\n")
}
## invisibly return check results
invisible(mget(c("comp.F", "comp.Fcircle",
"maxRelDiffs.deriv", "comp.Deriv"),
ifnotfound=list(NULL), inherits=FALSE))
}
checksiaf.F <- function (F, f, pargrid, type=1, method="SV", ...)
{
res <- t(apply(pargrid, 1, function (pars) {
given <- F(LETTERR, f, pars, type)
num <- siaf.fallback.F(polydomain = LETTERR, f = f, pars = pars,
type = type, method = method, ...)
c(given, num)
}))
colnames(res) <- c("F", method)
res
}
checksiaf.Fcircle <- function (Fcircle, f, pargrid, type=1,
rs=c(1,5,10,50,100), method="SV", ...)
{
pargrid <- pargrid[rep(1:nrow(pargrid), each=length(rs)),,drop=FALSE]
rpargrid <- cbind(rs, pargrid, deparse.level=0)
res <- t(apply(rpargrid, 1, function (x) {
disc <- discpoly(c(0,0), x[1L], npoly = 128, class = "owin")
c(ana = Fcircle(x[1L], x[-1L], type),
num = siaf.fallback.F(polydomain = disc, f = f, pars = x[-1L],
type = type, method = method, ...))
}))
res
}
checksiaf.deriv <- function (deriv, f, pargrid, type=1, rmax=100)
{
rgrid <- seq(-rmax,rmax,length.out=21) / sqrt(2)
rgrid <- rgrid[rgrid != 0] # some siafs are always 1 at (0,0) (deriv=0)
sgrid <- cbind(rgrid, rgrid)
apply(pargrid, 1, function (pars) {
maxLik::compareDerivatives(f, deriv, t0=pars, s=sgrid,
print=FALSE)$maxRelDiffGrad
## Note: numDeriv::grad() would only allow one location s at a time
})
}
checksiaf.Deriv <- function (Deriv, deriv, pargrid, type=1, method="SV", ...)
{
res <- t(apply(pargrid, 1, function (pars) {
given <- Deriv(LETTERR, deriv, pars, type)
num <- siaf.fallback.Deriv(polydomain = LETTERR, deriv = deriv, pars = pars,
type = type, method = method, ...)
c(given, num)
}))
paridxs <- seq_len(ncol(pargrid))
colnames(res) <- c(paste("Deriv",paridxs,sep="."),
paste(method,paridxs,sep="."))
res
}
checksiaf.simulate <- function (simulate, f, pars, type=1, B=3000, ub=10,
plot=interactive())
{
## Simulate B points on the disc with radius 'ub'
simpoints <- simulate(B, pars, type=type, ub=ub)
if (plot) {
## Graphical check in 2D
opar <- par(mfrow=c(2,1), mar=c(4,3,2,1)); on.exit(par(opar))
plot(as.im.function(function(x,y,...) f(cbind(x,y), pars, type),
W=discpoly(c(0,0), ub, class="owin")),
axes=TRUE, main="Simulation from the spatial kernel")
points(simpoints, cex=0.2)
kdens <- kde2d(simpoints[,1], simpoints[,2], n=100)
contour(kdens, add=TRUE, col=2, lwd=2,
labcex=1.5, vfont=c("sans serif", "bold"))
##x11(); image(kdens, add=TRUE)
## Graphical check of distance distribution
truehist(sqrt(rowSums(simpoints^2)), xlab="Distance")
rfr <- function (r) r*f(cbind(r,0), pars, type)
rfrnorm <- integrate(rfr, 0, ub)$value
do.call("curve", list(quote(rfr(x)/rfrnorm), add=TRUE, col=2, lwd=2))
##<- use do.call-construct to prevent codetools::checkUsage from noting "x"
}
## invisibly return simulated points
invisible(simpoints)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.