https://github.com/cran/spatstat
Tip revision: 489e4596bd2a4db09fabcc2b72c6831e6b662a94 authored by Adrian Baddeley on 12 May 2020, 16:10:02 UTC
version 1.64-1
version 1.64-1
Tip revision: 489e459
rmhmodel.R
#
#
# rmhmodel.R
#
# $Revision: 1.77 $ $Date: 2020/01/08 01:21:59 $
#
#
rmhmodel <- function(...) {
UseMethod("rmhmodel")
}
rmhmodel.rmhmodel <- function(model, ...) {
# Check for outdated internal format
# C.par was replaced by C.beta and C.ipar in spatstat 1.22-3
if(outdated <- !is.null(model$C.par))
warning("Outdated internal format of rmhmodel object; rebuilding it")
if(outdated || (length(list(...)) > 0))
model <- rmhmodel.list(unclass(model), ...)
return(model)
}
rmhmodel.list <- function(model, ...) {
argnames <- c("cif","par","w","trend","types")
ok <- argnames %in% names(model)
do.call(rmhmodel.default,
resolve.defaults(list(...), model[argnames[ok]]))
}
rmhmodel.default <- local({
rmhmodel.default <- function(...,
cif=NULL, par=NULL, w=NULL, trend=NULL, types=NULL) {
rmhmodelDefault(..., cif=cif, par=par, w=w, trend=trend, types=types)
}
rmhmodelDefault <- function(...,
cif=NULL, par=NULL, w=NULL, trend=NULL, types=NULL,
stopinvalid=TRUE) {
if(length(list(...)) > 0)
stop(paste("rmhmodel.default: syntax should be",
"rmhmodel(cif, par, w, trend, types)",
"with arguments given by name if they are present"),
call. = FALSE)
## Validate parameters
if(is.null(cif)) stop("cif is missing or NULL")
if(is.null(par)) stop("par is missing or NULL")
if(!is.null(w))
w <- as.owin(w)
if(!is.character(cif))
stop("cif should be a character string")
betamultiplier <- 1
Ncif <- length(cif)
if(Ncif > 1) {
## hybrid
## check for Poisson components
ispois <- (cif == 'poisson')
if(any(ispois)) {
## validate Poisson components
Npois <- sum(ispois)
poismodels <- vector(mode="list", length=Npois)
parpois <- par[ispois]
for(i in 1:Npois)
poismodels[[i]] <- rmhmodel(cif='poisson', par=parpois[[i]],
w=w, trend=NULL, types=types,
stopinvalid=FALSE)
## consolidate Poisson intensity parameters
poisbetalist <- lapply(poismodels, getElement, name="C.beta")
poisbeta <- Reduce("*", poisbetalist)
if(all(ispois)) {
## model collapses to a Poisson process
cif <- 'poisson'
Ncif <- 1
par <- list(beta=poisbeta)
betamultiplier <- 1
} else {
## remove Poisson components
cif <- cif[!ispois]
Ncif <- sum(!ispois)
par <- par[!ispois]
if(Ncif == 1) # revert to single-cif format
par <- par[[1]]
## absorb beta parameters
betamultiplier <- poisbeta
}
}
}
if(Ncif > 1) {
## genuine hybrid
models <- vector(mode="list", length=Ncif)
check <- vector(mode="list", length=Ncif)
for(i in 1:Ncif)
models[[i]] <- rmhmodel(cif=cif[i], par=par[[i]],
w=w, trend=NULL, types=types,
stopinvalid=FALSE)
C.id <- unlist(lapply(models, getElement, name="C.id"))
C.betalist <- lapply(models, getElement, name="C.beta")
C.iparlist <- lapply(models, getElement, name="C.ipar")
## absorb beta multiplier into beta parameter of first component
C.betalist[[1]] <- C.betalist[[1]] * betamultiplier
## concatenate for use in C
C.beta <- unlist(C.betalist)
C.ipar <- unlist(C.iparlist)
check <- lapply(models, getElement, name="check")
maxr <- max(unlist(lapply(models, getElement, name="reach")))
ismulti <- unlist(lapply(models, getElement, name="multitype.interact"))
multi <- any(ismulti)
## determine whether model exists
integ <- unlist(lapply(models, getElement, name="integrable"))
stabi <- unlist(lapply(models, getElement, name="stabilising"))
integrable <- all(integ) || any(stabi)
stabilising <- any(stabi)
## string explanations of conditions for validity
expl <- lapply(models, getElement, name="explainvalid")
integ.ex <- unlist(lapply(expl, getElement, name="integrable"))
stabi.ex <- unlist(lapply(expl, getElement, name="stabilising"))
stabi.oper <- !(stabi.ex %in% c("TRUE", "FALSE"))
integ.oper <- !(integ.ex %in% c("TRUE", "FALSE"))
compnames <- if(!anyDuplicated(C.id)) paste("cif", sQuote(C.id)) else
paste("component", 1:Ncif, paren(sQuote(C.id)))
if(!integrable && stopinvalid) {
## model is not integrable: explain why
ifail <- !integ & integ.oper
ireason <- paste(compnames[ifail], "should satisfy",
paren(integ.ex[ifail], "{"))
ireason <- verbalogic(ireason, "and")
if(sum(ifail) <= 1) {
## There's only one offending cif, so stability is redundant
sreason <- "FALSE"
} else {
sfail <- !stabi & stabi.oper
sreason <- paste(compnames[sfail], "should satisfy",
paren(stabi.ex[sfail], "{"))
sreason <- verbalogic(sreason, "or")
}
reason <- verbalogic(c(ireason, sreason), "or")
stop(paste("rmhmodel: hybrid model is not integrable; ", reason),
call.=FALSE)
} else {
## construct strings summarising conditions for validity
if(!any(integ.oper))
ireason <- as.character(integrable)
else {
ireason <- paste(compnames[integ.oper], "should satisfy",
paren(integ.ex[integ.oper], "{"))
ireason <- verbalogic(ireason, "and")
}
if(!any(stabi.oper))
sreason <- as.character(stabilising)
else {
sreason <- paste(compnames[stabi.oper], "should satisfy",
paren(stabi.ex[stabi.oper], "{"))
sreason <- verbalogic(sreason, "or")
}
ireason <- verbalogic(c(ireason, sreason), "or")
explainvalid <- list(integrable=ireason,
stabilising=sreason)
}
out <- list(cif=cif,
par=par,
w=w,
trend=trend,
types=types,
C.id=C.id,
C.beta=C.beta,
C.ipar=C.ipar,
C.betalist=C.betalist,
C.iparlist=C.iparlist,
check=check,
multitype.interact=multi,
integrable=integrable,
stabilising=stabilising,
explainvalid=explainvalid,
reach=maxr)
class(out) <- c("rmhmodel", class(out))
return(out)
}
## non-hybrid
## Check that this is a recognised model
## and look up the rules for this model
rules <- spatstatRmhInfo(cif)
## Map the name of the cif from R to C
## (the names are normally identical in R and C,
## except "poisson" -> NA)
C.id <- rules$C.id
## Check that the C name is recognised in C
if(!is.na(C.id)) {
z <- .C("knownCif",
cifname=as.character(C.id),
answer=as.integer(0),
PACKAGE = "spatstat")
ok <- as.logical(z$answer)
if(!ok)
stop(paste("Internal error: the cif", sQuote(C.id),
"is not recognised in the C code"))
}
## Validate the model parameters and reformat them
check <- rules$parhandler
checkedpar <-
if(!rules$multitype)
check(par)
else if(!is.null(types))
check(par, types)
else
## types vector not given - defer checking
NULL
if(!is.null(checkedpar)) {
stopifnot(is.list(checkedpar))
stopifnot(!is.null(names(checkedpar)) && all(nzchar(names(checkedpar))))
stopifnot(names(checkedpar)[[1]] == "beta")
C.beta <- unlist(checkedpar[[1]])
C.beta <- C.beta * betamultiplier
C.ipar <- as.numeric(unlist(checkedpar[-1]))
} else {
C.beta <- C.ipar <- NULL
}
## Determine whether model is integrable
integrable <- rules$validity(par, "integrable")
explainvalid <- rules$explainvalid
if(!integrable && stopinvalid)
stop(paste("rmhmodel: the model is not integrable; it should satisfy",
explainvalid$integrable),
call.=FALSE)
## Determine whether cif is stabilising
## (i.e. any hybrid including this cif will be integrable)
stabilising <- rules$validity(par, "stabilising")
## Calculate reach of model
mreach <- rules$reach(par)
###################################################################
## return augmented list
out <- list(cif=cif,
par=par,
w=w,
trend=trend,
types=types,
C.id=C.id,
C.beta=C.beta,
C.ipar=C.ipar,
check= if(is.null(C.ipar)) check else NULL,
multitype.interact=rules$multitype,
integrable=integrable,
stabilising=stabilising,
explainvalid=explainvalid,
reach=mreach
)
class(out) <- c("rmhmodel", class(out))
return(out)
}
rmhmodel.default
})
print.rmhmodel <- function(x, ...) {
verifyclass(x, "rmhmodel")
splat("Metropolis-Hastings algorithm, model parameters\n")
Ncif <- length(x$cif)
splat("Conditional intensity:",
if(Ncif == 1) "cif=" else "hybrid of cifs",
commasep(sQuote(x$cif)))
if(!is.null(x$types)) {
if(length(x$types) == 1)
splat("Univariate process.")
else {
cat("Multitype process with types =\n")
print(x$types)
if(!x$multitype.interact)
splat("Interaction does not depend on type")
}
} else if(x$multitype.interact) {
splat("Multitype process, types not yet specified.")
} else {
typ <- try(rmhResolveTypes(x, rmhstart(), rmhcontrol()))
if(!inherits(typ, "try-error")) {
ntyp <- length(typ)
if(ntyp > 1) {
splat("Data imply a multitype process with", ntyp, "types of points.")
splat("Interaction does not depend on type.")
}
}
}
cat("\nNumerical parameters: par =\n")
print(x$par)
if(is.null(x$C.ipar))
splat("Parameters have not yet been checked for compatibility with types.")
if(is.owin(x$w)) print(x$w) else splat("Window: not specified.")
cat("\nTrend: ")
tren <- x$trend
if(is.null(tren)) {
cat("none.\n")
} else {
if(is.list(tren)) cat(paste0("List of ", length(tren), ":\n"))
print(tren)
}
if(!is.null(x$integrable) && !x$integrable)
cat("\n*Warning: model is not integrable and cannot be simulated*\n")
return(invisible(NULL))
}
reach.rmhmodel <- function(x, ...) {
if(length(list(...)) == 0)
return(x$reach)
# reach must be recomputed
cif <- x$cif
Ncif <- length(cif)
pars <- if(Ncif == 1) list(x$par) else x$par
maxr <- 0
for(i in seq_len(Ncif)) {
cif.i <- cif[i]
par.i <- pars[[i]]
rules <- spatstatRmhInfo(cif.i)
rchfun <- rules$reach
if(!is.function(rchfun))
stop(paste("Internal error: reach is unknown for cif=", sQuote(cif.i)),
call.=FALSE)
r.i <- rchfun(par.i, ...)
maxr <- max(maxr, r.i, na.rm=TRUE)
}
return(maxr)
}
is.poisson.rmhmodel <- function(x) {
verifyclass(x, "rmhmodel")
identical(x$cif, 'poisson')
}
is.stationary.rmhmodel <- function(x) {
verifyclass(x, "rmhmodel")
tren <- x$trend
return(is.null(tren) || is.numeric(tren))
}
as.owin.rmhmodel <- function(W, ..., fatal=FALSE) {
# W is the rmhmodel object. It contains a window w
ans <- W$w
if(is.owin(ans)) return(ans)
if(fatal) stop("rmhmodel object does not contain a window")
return(NULL)
}
domain.rmhmodel <- Window.rmhmodel <- function(X, ...) { as.owin(X) }
is.expandable.rmhmodel <- local({
ok <- function(z) { is.null(z) || is.numeric(z) || is.function(z) }
is.expandable.rmhmodel <- function(x) {
tren <- x$tren
ans <- if(!is.list(tren)) ok(tren) else all(sapply(tren, ok))
return(ans)
}
is.expandable.rmhmodel
})
##### Table of rules for handling rmh models ##################
spatstatRmhInfo <- function(cifname) {
rules <- .Spatstat.RmhTable[[cifname]]
if(is.null(rules))
stop(paste("Unrecognised cif:", sQuote(cifname)), call.=FALSE)
return(rules)
}
.Spatstat.RmhTable <-
list(
#
# 0. Poisson (special case)
#
'poisson'=
list(
C.id=NA,
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the Poisson process"
with(par, forbidNA(beta, ctxt))
par <- check.named.list(par, "beta", ctxt)
with(par, explain.ifnot(all(beta >= 0), ctxt))
return(par)
},
validity=function(par, kind) {
switch(kind,
integrable=TRUE,
stabilising=FALSE)
},
explainvalid=list(integrable="TRUE",stabilising="FALSE"),
reach = function(par, ...) { return(0) },
hardcore = function(par, ...) { return(0) },
temper = function(par, invtemp) { return(par^invtemp) }
),
#
# 1. Strauss.
#
'strauss'=
list(
C.id="strauss",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the strauss cif"
par <- check.named.list(par, c("beta","gamma","r"), ctxt)
# treat r=NA as absence of interaction
par <- within(par, if(is.na(r)) { r <- 0; gamma <- 1 })
with(par, check.finite(beta, ctxt))
with(par, check.finite(gamma, ctxt))
with(par, check.finite(r, ctxt))
with(par, check.1.real(gamma, ctxt))
with(par, check.1.real(r, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
with(par, explain.ifnot(gamma >= 0, ctxt))
with(par, explain.ifnot(r >= 0, ctxt))
return(par)
},
validity=function(par, kind) {
gamma <- par$gamma
switch(kind,
integrable=(gamma <= 1),
stabilising=(gamma == 0)
)
},
explainvalid=list(
integrable="gamma <= 1",
stabilising="gamma == 0"),
reach = function(par, ...) {
r <- par[["r"]]
g <- par[["gamma"]]
return(if(g == 1) 0 else r)
},
hardcore = function(par, ..., epsilon=0) {
r <- par[["r"]]
g <- par[["gamma"]]
return(if(g <= epsilon) r else 0)
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
gamma <- gamma^invtemp
})
}
),
#
# 2. Strauss with hardcore.
#
'straush' =
list(
C.id="straush",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the straush cif"
par <- check.named.list(par, c("beta","gamma","r","hc"), ctxt)
# treat hc=NA as absence of hard core
par <- within(par, if(is.na(hc)) { hc <- 0 } )
# treat r=NA as absence of interaction
par <- within(par, if(is.na(r)) { r <- hc; gamma <- 1 } )
with(par, check.finite(beta, ctxt))
with(par, check.finite(gamma, ctxt))
with(par, check.finite(r, ctxt))
with(par, check.finite(hc, ctxt))
with(par, check.1.real(gamma, ctxt))
with(par, check.1.real(r, ctxt))
with(par, check.1.real(hc, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
with(par, explain.ifnot(gamma >= 0, ctxt))
with(par, explain.ifnot(r >= 0, ctxt))
with(par, explain.ifnot(hc >= 0, ctxt))
with(par, explain.ifnot(hc <= r, ctxt))
return(par)
},
validity=function(par, kind) {
hc <- par$hc
gamma <- par$gamma
switch(kind,
integrable=(hc > 0 || gamma <= 1),
stabilising=(hc > 0)
)
},
explainvalid=list(
integrable="hc > 0 or gamma <= 1",
stabilising="hc > 0"),
reach = function(par, ...) {
h <- par[["hc"]]
r <- par[["r"]]
g <- par[["gamma"]]
return(if(g == 1) h else r)
},
hardcore = function(par, ..., epsilon=0) {
h <- par[["hc"]]
r <- par[["r"]]
g <- par[["gamma"]]
return(if(g <= epsilon) r else h)
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
gamma <- gamma^invtemp
})
}
),
#
# 3. Softcore.
#
'sftcr' =
list(
C.id="sftcr",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the sftcr cif"
par <- check.named.list(par, c("beta","sigma","kappa"), ctxt)
with(par, check.finite(beta, ctxt))
with(par, check.finite(sigma, ctxt))
with(par, check.finite(kappa, ctxt))
with(par, check.1.real(sigma, ctxt))
with(par, check.1.real(kappa, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
with(par, explain.ifnot(sigma >= 0, ctxt))
with(par, explain.ifnot(kappa >= 0 && kappa <= 1, ctxt))
return(par)
},
validity=function(par, kind) {
switch(kind,
integrable=TRUE,
stabilising=FALSE)
},
explainvalid=list(integrable="TRUE",stabilising="FALSE"),
reach = function(par, ..., epsilon=0) {
if(epsilon==0)
return(Inf)
kappa <- par[["kappa"]]
sigma <- par[["sigma"]]
return(sigma/(epsilon^(kappa/2)))
},
hardcore = function(par, ..., epsilon=0) {
if(epsilon==0)
return(0)
kappa <- par[["kappa"]]
sigma <- par[["sigma"]]
return(sigma/((-log(epsilon))^(kappa/2)))
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
sigma <- sigma * (invtemp^(kappa/2))
})
}
),
#
# 4. Multitype Strauss.
#
'straussm' =
list(
C.id="straussm",
multitype=TRUE,
parhandler=function(par, types) {
ctxt <- "For the straussm cif"
par <- check.named.list(par, c("beta","gamma","radii"), ctxt)
beta <- par$beta
gamma <- par$gamma
r <- par$radii
ntypes <- length(types)
check.finite(beta, ctxt)
check.nvector(beta, ntypes, TRUE, "types")
MultiPair.checkmatrix(gamma, ntypes, "par$gamma")
gamma[is.na(gamma)] <- 1
check.finite(gamma, ctxt)
MultiPair.checkmatrix(r, ntypes, "par$radii")
if(any(nar <- is.na(r))) {
r[nar] <- 0
gamma[nar] <- 1
}
check.finite(r, ctxt)
explain.ifnot(all(beta >= 0), ctxt)
explain.ifnot(all(gamma >= 0), ctxt)
explain.ifnot(all(r >= 0), ctxt)
par <- list(beta=beta, gamma=gamma, r=r)
return(par)
},
validity=function(par, kind) {
gamma <- par$gamma
radii <- par$radii
dg <- diag(gamma)
dr <- diag(radii)
hard <-!is.na(dg) & (dg == 0) & !is.na(dr) & (dr > 0)
operative <- !is.na(gamma) & !is.na(radii) & (radii > 0)
switch(kind,
stabilising=all(hard),
integrable=all(hard) || all(gamma[operative] <= 1))
},
explainvalid=list(
integrable=paste(
"gamma[i,j] <= 1 for all i and j,",
"or gamma[i,i] = 0 for all i"),
stabilising="gamma[i,i] = 0 for all i"),
reach = function(par, ...) {
r <- par$radii
g <- par$gamma
operative <- ! (is.na(r) | (g == 1))
return(max(0, r[operative]))
},
hardcore = function(par, ..., epsilon=0) {
r <- par$radii
g <- par$gamma
return(max(0, r[!is.na(r) & g <= epsilon]))
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
gamma <- gamma^invtemp
})
}
),
#
# 5. Multitype Strauss with hardcore.
#
'straushm' =
list(
C.id="straushm",
multitype=TRUE,
parhandler=function(par, types) {
ctxt="For the straushm cif"
par <- check.named.list(par,
c("beta","gamma","iradii","hradii"),
ctxt)
beta <- par$beta
gamma <- par$gamma
iradii <- par$iradii
hradii <- par$hradii
ntypes <- length(types)
check.nvector(beta, ntypes, TRUE, "types")
check.finite(beta, ctxt)
MultiPair.checkmatrix(gamma, ntypes, "par$gamma")
gamma[is.na(gamma)] <- 1
check.finite(gamma, ctxt)
MultiPair.checkmatrix(iradii, ntypes, "par$iradii")
if(any(nar <- is.na(iradii))) {
iradii[nar] <- 0
gamma[nar] <- 1
}
check.finite(iradii, ctxt)
MultiPair.checkmatrix(hradii, ntypes, "par$hradii")
nah <- is.na(hradii)
hradii[nah] <- 0
check.finite(hradii, ctxt)
explain.ifnot(all(beta >= 0), ctxt)
explain.ifnot(all(gamma >= 0), ctxt)
explain.ifnot(all(iradii >= 0), ctxt)
explain.ifnot(all(hradii >= 0), ctxt)
comparable <- !nar & !nah
explain.ifnot(all((iradii >= hradii)[comparable]), ctxt)
par <- list(beta=beta,gamma=gamma,iradii=iradii,hradii=hradii)
return(par)
},
validity=function(par, kind) {
gamma <- par$gamma
iradii <- par$iradii
hradii <- par$hradii
dh <- diag(hradii)
dg <- diag(gamma)
dr <- diag(iradii)
hhard <- !is.na(dh) & (dh > 0)
ihard <- !is.na(dr) & (dr > 0) & !is.na(dg) & (dg == 0)
hard <- hhard | ihard
operative <- !is.na(gamma) & !is.na(iradii) & (iradii > 0)
switch(kind,
stabilising=all(hard),
integrable={
all(hard) || all(gamma[operative] <= 1)
})
},
explainvalid=list(
integrable=paste(
"hradii[i,i] > 0 or gamma[i,i] = 0 for all i, or",
"gamma[i,j] <= 1 for all i and j"),
stabilising="hradii[i,i] > 0 or gamma[i,i] = 0 for all i"),
reach=function(par, ...) {
r <- par$iradii
h <- par$hradii
g <- par$gamma
roperative <- ! (is.na(r) | (g == 1))
hoperative <- ! is.na(h)
return(max(0, r[roperative], h[hoperative]))
},
hardcore = function(par, ..., epsilon=0) {
r <- par$radii
h <- par$hradii
g <- par$gamma
return(max(h[!is.na(h)],
r[!is.na(r) & g <= epsilon]))
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
gamma <- gamma^invtemp
})
}
),
#
# 6. Diggle-Gates-Stibbard interaction
# (function number 1 from Diggle, Gates, and Stibbard)
'dgs' =
list(
C.id="dgs",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the dgs cif"
par <- check.named.list(par, c("beta","rho"), ctxt)
with(par, check.finite(beta, ctxt))
with(par, check.finite(rho, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
with(par, check.1.real(rho, ctxt))
with(par, explain.ifnot(rho >= 0, ctxt))
return(par)
},
validity=function(par, kind) {
switch(kind,
integrable=TRUE,
stabilising=FALSE)
},
explainvalid=list(integrable="TRUE", stabilising="FALSE"),
reach=function(par, ...) {
return(par[["rho"]])
},
hardcore=function(par, ..., epsilon=0) {
if(epsilon == 0) return(0)
return(par[["rho"]] * (2/pi) * asin(sqrt(epsilon)))
},
temper = NULL # not a loglinear model
),
#
# 7. Diggle-Gratton interaction
# (function number 2 from Diggle, Gates, and Stibbard).
'diggra' =
list(
C.id="diggra",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the diggra cif"
par <- check.named.list(par, c("beta","kappa","delta","rho"),
ctxt)
with(par, check.finite(beta, ctxt))
with(par, check.finite(kappa, ctxt))
with(par, check.finite(delta, ctxt))
with(par, check.finite(rho, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
with(par, check.1.real(kappa, ctxt))
with(par, check.1.real(delta, ctxt))
with(par, check.1.real(rho, ctxt))
with(par, explain.ifnot(kappa >= 0, ctxt))
with(par, explain.ifnot(delta >= 0, ctxt))
with(par, explain.ifnot(rho >= 0, ctxt))
with(par, explain.ifnot(delta < rho, ctxt))
return(par)
},
validity=function(par, kind) {
switch(kind,
integrable=TRUE,
stabilising=FALSE)
},
explainvalid=list(integrable="TRUE",stabilising="FALSE"),
reach=function(par, ...) {
return(par[["rho"]])
},
hardcore=function(par, ..., epsilon=0) {
return(par[["delta"]])
},
temper = function(par, invtemp) {
within(par, {
kappa <- kappa * invtemp
})
}),
#
# 8. Geyer saturation model
#
'geyer' =
list(
C.id="geyer",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the geyer cif"
par <- check.named.list(par, c("beta","gamma","r","sat"), ctxt)
with(par, check.1.real(gamma, ctxt))
with(par, check.1.real(r, ctxt))
with(par, check.1.real(sat, ctxt))
par <- within(par, sat <- min(sat, .Machine$integer.max-100))
par <- within(par, if(is.na(gamma)) { r <- 0; gamma <- 1 })
with(par, check.finite(beta, ctxt))
with(par, check.finite(gamma, ctxt))
with(par, check.finite(r, ctxt))
with(par, check.finite(sat, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
return(par)
},
validity=function(par, kind) {
switch(kind,
integrable=TRUE,
stabilising=FALSE)
},
explainvalid=list(integrable="TRUE", stabilising="FALSE"),
reach = function(par, ...) {
r <- par[["r"]]
g <- par[["gamma"]]
return(if(g == 1) 0 else 2 * r)
},
hardcore = function(par, ..., epsilon=0) {
r <- par[["r"]]
g <- par[["gamma"]]
return(if(g <= epsilon) r else 0)
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
gamma <- gamma^invtemp
})
}
),
#
# 9. The ``lookup'' device. This permits simulating, at least
# approximately, ANY pairwise interaction function model
# with isotropic pair interaction (i.e. depending only on distance).
# The pair interaction function is provided as a vector of
# distances and corresponding function values which are used
# as a ``lookup table'' by the C code.
#
'lookup' =
list(
C.id="lookup",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the lookup cif"
par <- check.named.list(par, c("beta","h"), ctxt, "r")
with(par, check.finite(beta, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
beta <- par[["beta"]]
h.init <- par[["h"]]
r <- par[["r"]]
if(is.null(r)) {
if(!is.stepfun(h.init))
stop(paste("For cif=lookup, if component r of",
"par is absent then component h must",
"be a stepfun object."))
if(!is.cadlag(h.init))
stop(paste("The lookup pairwise interaction step",
"function must be right continuous,\n",
"i.e. built using the default values of the",
sQuote("f"), "and", sQuote("right"),
"arguments for stepfun."))
r <- knots(h.init)
h0 <- get("yleft",envir=environment(h.init))
h <- h.init(r)
nlook <- length(r)
if(!isTRUE(all.equal(h[nlook],1)))
stop(paste("The lookup interaction step function",
"must be equal to 1 for", dQuote("large"),
"distances."))
if(r[1] <= 0)
stop(paste("The first jump point (knot) of the lookup",
"interaction step function must be",
"strictly positive."))
h <- c(h0,h)
} else {
h <- h.init
nlook <- length(r)
if(length(h) != nlook)
stop("Mismatch of lengths of h and r lookup vectors.")
if(anyNA(r))
stop("Missing values not allowed in r lookup vector.")
if(is.unsorted(r))
stop("The r lookup vector must be in increasing order.")
if(r[1] <= 0)
stop(paste("The first entry of the lookup vector r",
"should be strictly positive."))
h <- c(h,1)
}
if(any(h < 0))
stop(paste("Negative values in the lookup",
"pairwise interaction function."))
if(h[1] > 0 & any(h > 1))
stop(paste("Lookup pairwise interaction function does",
"not define a valid point process."))
rmax <- r[nlook]
r <- c(0,r)
nlook <- nlook+1
deltar <- mean(diff(r))
if(isTRUE(all.equal(diff(r),rep.int(deltar,nlook-1)))) {
par <- list(beta=beta,nlook=nlook,
equisp=1,
deltar=deltar,rmax=rmax, h=h)
} else {
par <- list(beta=beta,nlook=nlook,
equisp=0,
deltar=deltar,rmax=rmax, h=h,
r=r)
}
return(par)
},
validity=function(par, kind) {
h <- par$h
if(is.stepfun(h))
h <- eval(expression(c(yleft,y)),envir=environment(h))
switch(kind,
integrable={
(h[1] == 0) || all(h <= 1)
},
stabilising={ h[1] == 0 })
},
explainvalid=list(
integrable="h[1] == 0 or h[i] <= 1 for all i",
stabilising="h[1] == 0"),
reach = function(par, ...) {
r <- par[["r"]]
h <- par[["h"]]
if(is.null(r))
r <- knots(h)
return(max(r))
},
hardcore = function(par, ..., epsilon=0) {
r <- par[["r"]]
h <- par[["h"]]
if(is.null(r))
r <- knots(h)
return(max(0, r[h <= epsilon]))
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
h <- h^invtemp
})
}
),
#
# 10. Area interaction
#
'areaint'=
list(
C.id="areaint",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the areaint cif"
par <- check.named.list(par, c("beta","eta","r"), ctxt)
par <- within(par, if(is.na(r)) { r <- 0 })
with(par, check.finite(beta, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
with(par, check.1.real(eta, ctxt))
with(par, check.1.real(r, ctxt))
with(par, check.finite(eta, ctxt))
with(par, check.finite(r, ctxt))
with(par, explain.ifnot(eta >= 0, ctxt))
with(par, explain.ifnot(r >= 0, ctxt))
return(par)
},
validity=function(par, kind) {
switch(kind,
integrable=TRUE,
stabilising=FALSE)
},
explainvalid=list(integrable="TRUE", stabilising="FALSE"),
reach = function(par, ...) {
r <- par[["r"]]
eta <- par[["eta"]]
return(if(eta == 1) 0 else (2 * r))
},
hardcore = function(par, ..., epsilon=0) {
r <- par[["r"]]
eta <- par[["eta"]]
if(eta > epsilon) return(0)
if(eta == 0) return(2 * r)
# linear approximation
return(2 * r * eta/epsilon)
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
eta <- eta^invtemp
})
}
),
#
# 11. The ``badgey'' (Baddeley-Geyer) model.
#
'badgey' =
list(
C.id="badgey",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the badgey cif"
par <- check.named.list(par, c("beta","gamma","r","sat"), ctxt)
par <- within(par, sat <- pmin(sat, .Machine$integer.max-100))
par <- within(par, gamma[is.na(gamma) | is.na(r)] <- 1)
par <- within(par, r[is.na(r)] <- 0)
with(par, check.finite(beta, ctxt))
with(par, check.finite(gamma, ctxt))
with(par, check.finite(r, ctxt))
with(par, check.finite(sat, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
with(par, explain.ifnot(all(gamma >= 0), ctxt))
with(par, explain.ifnot(all(r >= 0), ctxt))
with(par, explain.ifnot(all(sat >= 0), ctxt))
with(par, explain.ifnot(length(gamma) == length(r), ctxt))
gamma <- par[["gamma"]]
r <- par[["r"]]
sat <- par[["sat"]]
if(length(sat)==1) sat <- rep.int(sat,length(gamma))
else explain.ifnot(length(sat) == length(gamma), ctxt)
mmm <- cbind(gamma,r,sat)
mmm <- mmm[fave.order(r),]
ndisc <- length(r)
par <- list(beta=par$beta,ndisc=ndisc,parms=as.vector(t(mmm)))
return(par)
},
validity=function(par, kind) {
switch(kind,
integrable=TRUE,
stabilising=FALSE)
},
explainvalid=list(integrable="TRUE", stabilising="FALSE"),
reach = function(par, ...) {
r <- par[["r"]]
gamma <- par[["gamma"]]
operative <- (gamma != 1)
return(if(!any(operative)) 0 else (2 * max(r[operative])))
},
hardcore = function(par, ..., epsilon=0) {
r <- par[["r"]]
gamma <- par[["gamma"]]
return(max(0, r[gamma <= epsilon]))
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
gamma <- gamma^invtemp
})
}
),
#
# 12. The hard core process
'hardcore' =
list(
C.id="hardcore",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the hardcore cif"
par <- check.named.list(par, c("beta", "hc"), ctxt)
par <- within(par, if(is.na(hc)) { hc <- 0 })
with(par, check.finite(beta, ctxt))
with(par, check.finite(hc, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
with(par, check.1.real(hc, ctxt))
return(par)
},
validity=function(par, kind) {
hc <- par$hc
switch(kind,
integrable=TRUE,
stabilising=(hc > 0))
},
explainvalid=list(integrable="TRUE", stabilising="hc > 0"),
reach = function(par, ...) {
hc <- par[["hc"]]
return(hc)
},
hardcore = function(par, ...) {
hc <- par[["hc"]]
return(hc)
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
})
}
),
#
# Lucky 13. Fiksel process
'fiksel' =
list(
C.id="fiksel",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the Fiksel cif"
par <- check.named.list(par,
c("beta", "r", "hc", "kappa", "a"),
ctxt)
with(par, check.finite(beta, ctxt))
with(par, check.finite(r, ctxt))
with(par, check.finite(hc, ctxt))
with(par, check.finite(kappa, ctxt))
with(par, check.finite(a, ctxt))
with(par, check.1.real(r, ctxt))
with(par, check.1.real(hc, ctxt))
with(par, check.1.real(kappa, ctxt))
with(par, check.1.real(a, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
with(par, explain.ifnot(hc >= 0, ctxt))
with(par, explain.ifnot(r > hc, ctxt))
return(par)
},
validity=function(par, kind) {
hc <- par$hc
a <- par$a
switch(kind,
integrable=(hc > 0 || a <= 0),
stabilising=(hc > 0))
},
explainvalid=list(
integrable="hc > 0 or a <= 0",
stabilising="hc > 0"),
reach = function(par, ...) {
r <- par[["r"]]
hc <- par[["hc"]]
a <- par[["a"]]
return(if(a != 0) r else hc)
},
hardcore = function(par, ...) {
return(par[["hc"]])
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
a <- a * invtemp
})
}
),
#
# 14. Lennard-Jones
'lennard' =
list(
C.id="lennard",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the Lennard-Jones cif"
par <- check.named.list(par,
c("beta", "sigma", "epsilon"),
ctxt)
with(par, check.finite(beta, ctxt))
with(par, check.finite(sigma, ctxt))
with(par, check.finite(epsilon, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
with(par, check.1.real(sigma, ctxt))
with(par, check.1.real(epsilon, ctxt))
with(par, explain.ifnot(sigma > 0, ctxt))
with(par, explain.ifnot(epsilon > 0, ctxt))
return(par)
},
validity=function(par, kind) {
switch(kind,
integrable=(par$sigma > 0),
stabilising=FALSE)
},
explainvalid=list(
integrable="sigma > 0",
stabilising="FALSE"),
reach = function(par, ...) {
sigma <- par[["sigma"]]
return(2.5 * sigma)
},
hardcore = function(par, ...) {
sigma <- par[["sigma"]]
return(sigma/2.5)
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
epsilon <- epsilon * invtemp
})
}
),
#
# 15. Multitype hardcore.
#
'multihard' =
list(
C.id="multihard",
multitype=TRUE,
parhandler=function(par, types) {
ctxt="For the multihard cif"
par <- check.named.list(par,
c("beta","hradii"),
ctxt)
beta <- par$beta
hradii <- par$hradii
ntypes <- length(types)
check.nvector(beta, ntypes, TRUE, "types")
check.finite(beta, ctxt)
MultiPair.checkmatrix(hradii, ntypes, "par$hradii")
hradii[is.na(hradii)] <- 0
check.finite(hradii, ctxt)
explain.ifnot(all(beta >= 0), ctxt)
explain.ifnot(all(hradii >= 0), ctxt)
par <- list(beta=beta,hradii=hradii)
return(par)
},
validity=function(par, kind) {
switch(kind,
integrable=return(TRUE),
stabilising={
hself <- diag(par$hradii)
repel <- !is.na(hself) & (hself > 0)
return(all(repel))
})
},
explainvalid=list(
integrable="TRUE",
stabilising="hradii[i,i] > 0 for all i"),
reach=function(par, ...) {
return(max(0, par$hradii, na.rm=TRUE))
},
hardcore=function(par, ..., epsilon=0) {
return(max(0, par$hradii, na.rm=TRUE))
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
})
}
),
#
# 16. Triplets.
#
'triplets'=
list(
C.id="triplets",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the triplets cif"
par <- check.named.list(par, c("beta","gamma","r"), ctxt)
# treat r=NA as absence of interaction
par <- within(par, if(is.na(r)) { r <- 0; gamma <- 1 })
with(par, check.finite(beta, ctxt))
with(par, check.finite(gamma, ctxt))
with(par, check.finite(r, ctxt))
with(par, check.1.real(gamma, ctxt))
with(par, check.1.real(r, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
with(par, explain.ifnot(gamma >= 0, ctxt))
with(par, explain.ifnot(r >= 0, ctxt))
return(par)
},
validity=function(par, kind) {
gamma <- par$gamma
switch(kind,
integrable=(gamma <= 1),
stabilising=(gamma == 0)
)
},
explainvalid=list(
integrable="gamma <= 1",
stabilising="gamma == 0"),
reach = function(par, ...) {
r <- par[["r"]]
g <- par[["gamma"]]
return(if(g == 1) 0 else r)
},
hardcore = function(par, ...) {
return(0)
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
gamma <- gamma^invtemp
})
}
),
#
# 17. Penttinen.
#
'penttinen'=
list(
C.id="penttinen",
multitype=FALSE,
parhandler=function(par, ...) {
ctxt <- "For the penttinen cif"
par <- check.named.list(par, c("beta", "gamma", "r"), ctxt)
# treat r=NA as absence of interaction
par <- within(par, if(is.na(r)) { r <- 0; gamma <- 1 })
with(par, check.finite(beta, ctxt))
with(par, check.finite(gamma, ctxt))
with(par, check.finite(r, ctxt))
with(par, check.1.real(gamma, ctxt))
with(par, check.1.real(r, ctxt))
with(par, explain.ifnot(all(beta >= 0), ctxt))
with(par, explain.ifnot(gamma >= 0, ctxt))
with(par, explain.ifnot(r > 0, ctxt))
return(par)
},
validity=function(par, kind) {
gamma <- par$gamma
switch(kind,
integrable=(gamma <= 1),
stabilising=(gamma == 0)
)
},
explainvalid=list(
integrable="gamma <= 1",
stabilising="gamma == 0"),
reach = function(par, ...) {
r <- par[["r"]]
g <- par[["gamma"]]
return(if(g == 1) 0 else (2 * r))
},
hardcore = function(par, ..., epsilon=0) {
r <- par[["r"]]
g <- par[["gamma"]]
return(if(g <= epsilon) (2 * r) else 0)
},
temper = function(par, invtemp) {
within(par, {
beta <- beta^invtemp
gamma <- gamma^invtemp
})
}
)
# end of list '.Spatstat.RmhTable'
)