1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360
|
#
# interact.S
#
#
# $Revision: 1.32 $ $Date: 2022/03/07 04:00:24 $
#
# Class 'interact' representing the interpoint interaction
# of a point process model
# (e.g. Strauss process with a given threshold r)
#
# Class 'isf' representing a generic interaction structure
# (e.g. pairwise interactions)
#
# These do NOT specify the "trend" part of the model,
# only the "interaction" component.
#
# The analogy is:
#
# glm() ppm()
#
# model formula trend formula
#
# family interaction
#
# That is, the 'systematic' trend part of a point process
# model is specified by a 'trend' formula argument to ppm(),
# and the interpoint interaction is specified as an 'interact'
# object.
#
# You only need to know about these classes if you want to
# implement a new point process model.
#
# THE DISTINCTION:
# An object of class 'isf' describes an interaction structure
# e.g. pairwise interaction, triple interaction,
# pairwise-with-saturation, Dirichlet interaction.
# Think of it as determining the "order" of interaction
# but not the specific interaction potential function.
#
# An object of class 'interact' completely defines the interpoint
# interactions in a specific point process model, except for the
# regular parameters of the interaction, which are to be estimated
# by ppm() or otherwise. An 'interact' object specifies the values
# of all the 'nuisance' or 'irregular' parameters. An example
# is the Strauss process with a given, fixed threshold r
# but with the parameters beta and gamma undetermined.
#
# DETAILS:
#
# An object of class 'isf' contains the following:
#
# $name Name of the interaction structure
# e.g. "pairwise"
#
# $print How to 'print()' this object
# [A function; invoked by the 'print' method
# 'print.isf()']
#
# $eval A function which evaluates the canonical
# sufficient statistic for an interaction
# of this general class (e.g. any pairwise
# interaction.)
#
# If lambda(u,X) denotes the conditional intensity at a point u
# for the point pattern X, then we assume
# log lambda(u, X) = theta . S(u,X)
# where theta is the vector of regular parameters,
# and we call S(u,X) the sufficient statistic.
#
# A typical calling sequence for the $eval function is
#
# (f$eval)(X, U, E, potentials, potargs, correction)
#
# where X is the data point pattern, U is the list of points u
# at which the sufficient statistic S(u,X) is to be evaluated,
# E is a logical matrix equivalent to (X[i] == U[j]),
# $potentials defines the specific potential function(s) and
# $potargs contains any nuisance/irregular parameters of these
# potentials [the $potargs are passed to the $potentials without
# needing to be understood by $eval.]
# $correction is the name of the edge correction method.
#
#
# An object of class 'interact' contains the following:
#
#
# $name Name of the specific potential
# e.g. "Strauss"
#
# $family Object of class "isf" describing
# the interaction structure
#
# $pot The interaction potential function(s)
# -- usually a function or list of functions.
# (passed as an argument to $family$eval)
#
# $par list of any nuisance/irregular parameters
# (passed as an argument to $family$eval)
#
# $parnames vector of long names/descriptions
# of the parameters in 'par'
#
# $init() initialisation action
# or NULL indicating none required
#
# $update() A function to modify $par
# [Invoked by 'update.interact()']
# or NULL indicating a default action
#
# $print How to 'print()' this object
# [Invoked by 'print' method 'print.interact()']
# or NULL indicating a default action
#
# --------------------------------------------------------------------------
print.isf <- function(x, ...) {
if(is.null(x)) return(invisible(NULL))
verifyclass(x, "isf")
if(!is.null(x$print))
(x$print)(x)
invisible(NULL)
}
print.interact <- function(x, ..., family, brief=FALSE, banner=TRUE) {
verifyclass(x, "interact")
if(missing(family)) family <- waxlyrical('extras')
#' Print name of model
if(banner) {
if(family && !brief && !is.null(xf <- x$family))
print.isf(xf)
splat(if(!brief) "Interaction:" else NULL, x$name, sep="")
}
# Now print the parameters
if(!is.null(x$print)) {
(x$print)(x)
} else {
# default
# just print the parameter names and their values
pwords <- x$parnames
parval <- x$par
pwords <- paste(toupper(substring(pwords, 1, 1)),
substring(pwords, 2), sep="")
isnum <- sapply(parval, is.numeric)
parval[isnum] <- lapply(parval[isnum], signif,
digits=getOption("digits"))
splat(paste(paste0(pwords, ":\t", parval), collapse="\n"))
}
invisible(NULL)
}
is.interact <- function(x) { inherits(x, "interact") }
update.interact <- function(object, ...) {
verifyclass(object, "interact")
if(!is.null(object$update))
(object$update)(object, ...)
else {
# Default
# First update the version
if(outdated.interact(object))
object <- reincarnate.interact(object)
# just match the arguments in "..."
# with those in object$par and update them
want <- list(...)
if(length(want) > 0) {
m <- match(names(want),names(object$par))
nbg <- is.na(m)
if(any(nbg)) {
which <- paste((names(want))[nbg])
warning(paste("Arguments not matched: ", which))
}
m <- m[!nbg]
object$par[m] <- want
}
# call object's own initialisation routine
if(!is.null(object$init))
(object$init)(object)
object
}
}
is.poisson.interact <- function(x) {
verifyclass(x, "interact")
is.null(x$family)
}
parameters.interact <- function(model, ...) {
model$par
}
# Test whether interact object was made by an older version of spatstat
outdated.interact <- function(object) {
ver <- object$version
older <- is.null(ver) || (package_version(ver) < versionstring.spatstat())
return(older)
}
# Test whether the functions in the interaction object
# expect the coefficient vector to contain ALL coefficients,
# or only the interaction coefficients.
# This change was introduced in 1.11-0, at the same time
# as interact objects were given version numbers.
newstyle.coeff.handling <- function(object) {
stopifnot(inherits(object, "interact"))
ver <- object$version
old <- is.null(ver) || (package_version(ver) < "1.11")
return(!old)
}
# ######
#
# Re-create an interact object in the current version of spatstat
#
#
reincarnate.interact <- function(object) {
# re-creates an interact object in the current version of spatstat
if(!is.null(object$update)) {
newobject <- (object$update)(object)
return(newobject)
}
par <- object$par
# pot <- object$pot
name <- object$name
# get creator function
creator <- object$creator
if(is.null(creator)) {
# old version: look up list
creator <- .Spatstat.Old.InteractionList[[name]]
if(is.null(creator))
stop(paste("Don't know how to update", sQuote(name),
"to current version of spatstat"))
}
if(is.character(creator))
creator <- get(creator)
if(!is.function(creator) && !is.expression(creator))
stop("Internal error: creator is not a function or expression")
# call creator
if(is.expression(creator))
newobject <- eval(creator)
else {
# creator is a function
# It's assumed that the creator function's arguments are
# either identical to components of 'par' (the usual case)
# or to one of the components of the object itself (Ord, Saturated)
# or to printfun=object$print (Pairwise).
argnames <- names(formals(creator))
available <- append(par, object)
available <- append(available, list(printfun=object$print))
ok <- argnames %in% names(available)
if(!all(ok))
stop(paste("Internal error:",
ngettext(sum(!ok), "argument", "arguments"),
paste(sQuote(argnames[!ok]), collapse=", "),
"in creator function were not understood"))
newobject <- do.call(creator, available[argnames])
}
if(!inherits(newobject, "interact"))
stop("Internal error: creator did not return an object of class interact")
return(newobject)
}
# This list is necessary to deal with older formats of 'interact' objects
# which did not include the creator name
.Spatstat.Old.InteractionList <-
list("Diggle-Gratton process" = "DiggleGratton",
"Geyer saturation process" = "Geyer",
"Lennard-Jones potential" = "LennardJones",
"Multitype Strauss process" = "MultiStrauss",
"Multitype Strauss Hardcore process" = "MultiStraussHard",
"Ord process with threshold potential"="OrdThresh",
"Piecewise constant pairwise interaction process"="PairPiece",
"Poisson process" = "Poisson",
"Strauss process" = "Strauss",
"Strauss - hard core process" = "StraussHard",
"Soft core process" = "Softcore",
# weird ones:
"Ord process with user-defined potential" = expression(Ord(object$pot)),
"Saturated process with user-defined potential"
=expression(Saturated(object$pot)),
"user-defined pairwise interaction process"=
expression(
Pairwise(object$pot,
par=object$par,
parnames=object$parnames,
printfun=object$print))
)
as.interact <- function(object) {
UseMethod("as.interact")
}
as.interact.interact <- function(object) {
verifyclass(object, "interact")
return(object)
}
as.isf <- function(object) {
if(inherits(object, "isf")) return(object)
object <- as.interact(object)
return(object$family)
}
interactionfamilyname <- function(object) { as.isf(object)$name }
interactionorder <- function(object) {
UseMethod("interactionorder")
}
interactionorder.isf <- function(object) {
return(object$order %orifnull% Inf)
}
interactionorder.interact <- function(object) {
## order may be specified in the interaction object itself (e.g. in a hybrid, or Poisson)
## but is usually determined by the interaction family
object$order %orifnull% interactionorder(object$family)
}
interactionorder.fii <- interactionorder.ppm <- function(object) {
interactionorder(as.interact(object))
}
# Extract version string from interact object
versionstring.interact <- function(object) {
verifyclass(object, "interact")
v <- object$version
return(v) # NULL before 1.11-0
}
#### internal code for streamlining initialisation of interactions
#
# x should be a partially-completed 'interact' object
#
instantiate.interact <- function(x, par=NULL) {
if(is.character(x$family)) x$family <- get(x$family)
# set parameter values
x$par <- par
# validate parameters
if(!is.null(x$init)) x$init(x)
x$version <- versionstring.spatstat()
return(x)
}
|