Nothing
#
#
# hybrid.R
#
# $Revision: 1.12 $ $Date: 2022/03/07 03:35:48 $
#
# Hybrid of several interactions
#
# Hybrid() create a hybrid of several interactions
# [an object of class 'interact']
#
#
# -------------------------------------------------------------------
#
Hybrid <- local({
Hybrid <- function(...) {
interlist <- list(...)
n <- length(interlist)
if(n == 0)
stop("No arguments given")
#' arguments may be interaction objects or ppm objects
isinter <- unlist(lapply(interlist, is.interact))
isppm <- unlist(lapply(interlist, is.ppm))
if(any(nbg <- !(isinter | isppm)))
stop(paste(ngettext(sum(nbg), "Argument", "Arguments"),
paste(which(nbg), collapse=", "),
ngettext(sum(nbg), "is not an interaction",
"are not interactions")))
#' ensure the list contains only interaction objects
if(any(isppm))
interlist[isppm] <- lapply(interlist[isppm], as.interact)
#' recursively expand any components that are themselves hybrids
while(any(ishybrid <- unlist(lapply(interlist, is.hybrid)))) {
i <- min(which(ishybrid))
n <- length(interlist)
expandi <- interlist[[i]]$par
interlist <- c(if(i > 1) interlist[1:(i-1L)] else NULL,
expandi,
if(i < n) interlist[(i+1L):n] else NULL)
}
#'
ncomponents <- length(interlist)
if(ncomponents == 1) {
#' single interaction - return it
return(interlist[[1L]])
}
#' ensure all components have names
names(interlist) <- good.names(names(interlist),
"HybridComponent", 1:ncomponents)
#' check for infinite potential values
haveInf <- lapply(interlist, getElement, name="hasInf")
haveInf <- !sapply(haveInf, identical, y=FALSE)
hasInf <- any(haveInf)
#' determine interaction order
maxorder <- max(sapply(interlist, interactionorder))
#' build object
out <-
list(
name = "Hybrid interaction",
creator = "Hybrid",
family = hybrid.family,
order = maxorder, # Overrides family$order
pot = NULL,
par = interlist,
parnames = names(interlist),
hasInf = hasInf,
selfstart = function(X, self) {
ilist <- self$par
sslist <- lapply(ilist, getElement, name="selfstart")
has.ss <- sapply(sslist, is.function)
if(any(has.ss)) {
ilist[has.ss] <- lapply(ilist[has.ss], invokeSelfStart, Y=X)
self$par <- ilist
}
return(self)
},
init = NULL,
update = NULL, # default OK
print = function(self, ..., family=FALSE, brief=FALSE) {
if(family)
print.isf(self$family)
ncomponents <- length(self$par)
clabs <- self$parnames
splat("Hybrid of", ncomponents, "components:",
commasep(sQuote(clabs)))
for(i in 1:ncomponents) {
splat(paste0(clabs[i], ":"))
print(self$par[[i]], ..., family=family, brief=brief)
}
parbreak()
return(invisible(NULL))
},
interpret = function(coeffs, self) {
interlist <- self$par
result <- list(param=list(),
inames=character(0),
printable=list())
for(i in 1:length(interlist)) {
interI <- interlist[[i]]
nameI <- names(interlist)[[i]]
nameI. <- paste(nameI, ".", sep="")
#' find coefficients with prefix that exactly matches nameI.
Cname <- names(coeffs)
prefixlength <- nchar(nameI.)
Cprefix <- substr(Cname, 1, prefixlength)
relevant <- (Cprefix == nameI.)
#' extract them
if(any(relevant)) {
Crelevant <- coeffs[relevant]
names(Crelevant) <-
substr(Cname[relevant], prefixlength+1, max(nchar(Cname)))
#' invoke the self-interpretation of interI
interpretI <- interI$interpret
if(is.function(interpretI)) {
resultI <- interpretI(Crelevant, interI)
paramI <- resultI$param
prinI <- resultI$printable
inamesI <- resultI$inames
inamesI <- paste(nameI, inamesI)
if(length(prinI) > 0) {
result$param <- append(result$param, paramI)
result$printable <- append(result$printable, list(prinI))
result$inames <- c(result$inames, inamesI)
}
}
}
}
return(result)
},
valid = function(coeffs, self) {
#' check validity via mechanism used for 'rmhmodel'
siminfo <- .Spatstat.Rmhinfo[["Hybrid interaction"]]
Z <- siminfo(coeffs, self)
cifs <- Z$cif
pars <- Z$par
ntypes <- Z$ntypes
if((Ncif <- length(cifs)) == 1) {
#' single cif
pars <- append(pars, list(beta=rep.int(1, ntypes)))
} else {
for(i in 1:Ncif)
pars[[i]] <- append(pars[[i]], list(beta=rep.int(1, ntypes[i])))
}
RM <- rmhmodel(cif=cifs, par=pars, types=1:max(ntypes),
stopinvalid=FALSE)
return(RM$integrable)
},
project = function(coeffs, self) {
if((self$valid)(coeffs, self)) return(NULL)
#' separate into components
spl <- splitHybridInteraction(coeffs, self)
interlist <- spl$interlist
coeflist <- spl$coeflist
#' compute projection for each component interaction
Ncif <- length(interlist)
projlist <- vector(mode="list", length=Ncif)
nproj <- integer(Ncif)
for(i in 1:Ncif) {
coefsI <- coeflist[[i]]
interI <- interlist[[i]]
if(!is.interact(interI))
stop("Internal error: interlist entry is not an interaction")
projI <- interI$project
if(is.null(projI))
stop(paste("Projection is not yet implemented for a",
interI$name))
p <- projI(coefsI, interI)
#' p can be NULL (indicating no projection required for interI)
#' or a single interaction or a list of interactions.
if(is.null(p)) {
if(Ncif == 1) return(NULL) # no projection required
p <- list(NULL)
nproj[i] <- 0
} else if(is.interact(p)) {
p <- list(p)
nproj[i] <- 1L
} else if(is.list(p) && all(unlist(lapply(p, is.interact)))) {
nproj[i] <- length(p)
} else
stop("Internal error: result of projection had wrong format")
projlist[[i]] <- p
}
#' for interaction i there are nproj[i] **new** interactions to try.
if(all(nproj == 0))
return(NULL)
if(spatstat.options("project.fast")) {
#' Single interaction required.
#' Extract first entry from each list
#' (there should be only one entry, but...)
qlist <- lapply(projlist, "[[", i=1L)
#' replace NULL entries by corresponding original interactions
isnul <- unlist(lapply(qlist, is.null))
if(all(isnul))
return(NULL)
if(any(isnul))
qlist[isnul] <- interlist[isnul]
names(qlist) <- names(interlist)
#' build hybrid and return
result <- do.call(Hybrid, qlist)
return(result)
}
#' Full case
result <- list()
for(i in which(nproj > 0)) {
ntry <- nproj[i]
tries <- projlist[[i]]
for(j in 1:ntry) {
#' assemble list of component interactions for hybrid
qlist <- interlist
qlist[[i]] <- tries[[j]]
#' eliminate Poisson
ispois <- unlist(lapply(qlist, is.poisson))
if(all(ispois)) {
#' collapse to single Poisson
h <- Poisson()
} else {
if(any(ispois)) qlist <- qlist[!ispois]
h <- do.call(Hybrid, qlist)
}
result <- append(result, list(h))
}
}
#' 'result' is a list of interactions, each a hybrid
if(length(result) == 1)
result <- result[[1L]]
return(result)
},
irange = function(self, coeffs=NA, epsilon=0, ...) {
interlist <- self$par
answer <- 0
for(i in 1:length(interlist)) {
interI <- interlist[[i]]
nameI <- names(interlist)[[i]]
nameI. <- paste(nameI, ".", sep="")
#' find coefficients with prefix that exactly matches nameI.
if(all(is.na(coeffs)))
Crelevant <- NA
else {
Cname <- names(coeffs)
prefixlength <- nchar(nameI.)
Cprefix <- substr(Cname, 1, prefixlength)
relevant <- (Cprefix == nameI.)
#' extract them
Crelevant <- coeffs[relevant]
names(Crelevant) <-
substr(Cname[relevant], prefixlength+1, max(nchar(Cname)))
}
#' compute reach
reachI <- interI$irange
if(is.function(reachI)) {
resultI <- reachI(interI,
coeffs=Crelevant, epsilon=epsilon, ...)
answer <- max(answer, resultI)
}
}
return(answer)
},
hardcore = function(self, coeffs=NA, epsilon=0, ...) {
interlist <- self$par
n <- length(interlist)
results <- vector(mode="list", length=n)
for(i in seq_len(n)) {
interI <- interlist[[i]]
nameI <- names(interlist)[[i]]
nameI. <- paste(nameI, ".", sep="")
#' find coefficients with prefix that exactly matches nameI.
if(all(is.na(coeffs)))
Crelevant <- NA
else {
Cname <- names(coeffs)
prefixlength <- nchar(nameI.)
Cprefix <- substr(Cname, 1, prefixlength)
relevant <- (Cprefix == nameI.)
#' extract them
Crelevant <- coeffs[relevant]
names(Crelevant) <-
substr(Cname[relevant], prefixlength+1, max(nchar(Cname)))
}
#' compute hard core for component i
hardI <- interI$hardcore
if(is.function(hardI))
results[[i]] <- hardI(interI,
coeffs=Crelevant, epsilon=epsilon, ...)
}
## collate answers
results <- results[!sapply(results, is.null)]
if(length(results) == 0)
return(0)
values <- Reduce(function(...) pmax(..., na.rm=TRUE), results)
dims <- lapply(results, dim)
dims <- dims[!sapply(dims, is.null)]
if(length(dims) == 0)
return(values)
dims <- unique(dims)
if(length(dims) > 1)
stop("Incompatible matrix dimensions of hardcore distance matrices in hybrid")
d <- dims[[1]]
dn <- unique(lapply(results, dimnames))[[1]]
answer <- matrix(values, d[1], d[2], dimnames=dn)
return(answer)
},
version=versionstring.spatstat()
)
class(out) <- "interact"
return(out)
}
invokeSelfStart <- function(inte, Y) {
ss <- inte$selfstart
if(!is.function(ss)) return(inte)
return(ss(Y, inte))
}
Hybrid
})
is.hybrid <- function(x) { UseMethod("is.hybrid") }
is.hybrid.interact <- function(x) {
return(is.interact(x) && (x$name == "Hybrid interaction"))
}
is.hybrid.ppm <- function(x) {
return(is.hybrid(as.interact(x)))
}
splitHybridInteraction <- function(coeffs, inte) {
# For hybrids, $par is a list of the component interactions,
# but coeffs is a numeric vector.
# Split the coefficient vector into the relevant coeffs for each interaction
interlist <- inte$par
N <- length(interlist)
coeflist <- vector(mode="list", length=N)
for(i in 1:N) {
interI <- interlist[[i]]
# forbid hybrids-of-hybrids - these should not occur anyway
if(interI$name == "Hybrid interaction")
stop("A hybrid-of-hybrid interactions is not implemented")
# nameI is the tag that identifies I-th component in hybrid
nameI <- names(interlist)[[i]]
nameI. <- paste(nameI, ".", sep="")
# find coefficients with prefix that exactly matches nameI.
Cname <- names(coeffs)
prefixlength <- nchar(nameI.)
Cprefix <- substr(Cname, 1, prefixlength)
relevant <- (Cprefix == nameI.)
# extract coefficients
# (there may be none, if this interaction is Poisson or an 'offset')
coeffsI <- coeffs[relevant]
# remove the prefix so the coefficients are recognisable to interaction
if(any(relevant))
names(coeffsI) <-
substr(Cname[relevant], prefixlength+1, max(nchar(Cname)))
# store
coeflist[[i]] <- coeffsI
}
names(coeflist) <- names(interlist)
return(list(coeflist=coeflist, interlist=interlist))
}
Hybrid <- intermaker(Hybrid, list(creator="Hybrid",
name="general hybrid Gibbs process",
par=list("..."=42),
parnames=list("any list of interactions")))
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.