https://github.com/cran/spatstat
Tip revision: ace26c246ee6feb8779515fa668bec59b24a1fcc authored by Adrian Baddeley on 12 March 2007, 13:35:27 UTC
version 1.11-2
version 1.11-2
Tip revision: ace26c2
rmh.default.R
#
# $Id: rmh.default.R,v 1.46 2006/12/12 08:57:33 adrian Exp adrian $
#
rmh.default <- function(model,start=NULL,control=NULL, verbose=TRUE, ...) {
#
# Function rmh. To simulate realizations of 2-dimensional point
# patterns, given the conditional intensity function of the
# underlying process, via the Metropolis-Hastings algorithm.
#
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
#
# V A L I D A T E
#
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
if(verbose)
cat("Checking arguments..")
# validate arguments and fill in the defaults
model <- rmhmodel(model)
start <- rmhstart(start)
control <- rmhcontrol(control)
#### Multitype models
# Decide whether the model is multitype; if so, find the types.
types <- rmhResolveTypes(model, start, control)
ntypes <- length(types)
mtype <- (ntypes > 1)
model$types <- types
# If the model is multitype, check that the model parameters agree with types
# and digest them
if(mtype && !is.null(model$check))
model$fortran.par <- model$check(model$par, types)
######## Check for illegal combinations of model, start and control ########
# No expansion can be done if we are using x.start
if(start$given == "x") {
if(control$force.exp)
stop("Cannot expand window when using x.start.\n")
else
control$expand <- 1
}
# Warn about a silly value of fixall:
if(control$fixall & ntypes==1) {
warning("control$fixall applies only to multitype processes. Ignored. \n")
control$fixall <- FALSE
if(control$conditioning == "n.each.type")
control$conditioning <- "n.total"
}
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
#
# M O D E L P A R A M E T E R S
#
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
####### Determine windows ################################
if(verbose)
cat("determining simulation windows...")
# these may be NULL
w.model <- model$w
x.start <- start$x.start
trend <- model$trend
trendy <- !is.null(trend)
# window implied by trend image, if any
w.trend <-
if(is.im(trend))
as.owin(trend)
else if(is.list(trend) && any(ok <- unlist(lapply(trend, is.im))))
as.owin((trend[ok])[[1]])
else NULL
## Clipping window (for final result)
w.clip <-
if(!is.null(model$w))
model$w
else if(control$expand == 1) {
if(start$given == "x" && is.ppp(x.start))
x.start$window
else if(is.owin(w.trend))
w.trend
} else NULL
if(!is.owin(w.clip))
stop("Unable to determine window for pattern")
## Simulation window
expand <- control$expand
w.sim <-
if(is.owin(expand))
expand
else if(is.numeric(expand))
expand.owin(w.clip, expand)
else
stop(paste("Internal error: unrecognised value of",
sQuote("expand")))
expanded <- (!is.numeric(expand) || (expand > 1))
## Check the fine print
if(expanded) {
if(control$conditioning != "none")
stop(paste("If we're conditioning on the number of points,",
"we cannot clip the result to another window.\n"))
if(!is.subset.owin(w.clip, w.sim))
stop("Expanded simulation window does not contain clipping window")
}
####### Trend ################################
# Check that the expanded window fits inside the window
# upon which the trend(s) live if there are trends and
# if any trend is given by an image.
if(expanded && !is.null(trend)) {
trends <- if(!is.list(trend)) list(trend) else trend
images <- unlist(lapply(trends, is.im))
if(any(images)) {
iwindows <- lapply(trends[images], as.owin)
nimages <- length(iwindows)
misfit <- !unlist(lapply(iwindows,
function(x,w) { is.subset.owin(w,x) },
w = w.sim))
nmisfit <- sum(misfit)
if(nmisfit > 1)
stop(paste("Expanded simulation window is not contained in",
"several of the trend windows.\n",
"Bailing out.\n"))
else if(nmisfit == 1) {
warning(paste("Expanded simulation window is not contained in",
if(nimages == 1) "the trend window.\n"
else "one of the trend windows.\n",
"Expanding to this trend window (only).\n"))
w.sim <- iwindows[misfit]
}
}
}
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
#
# S T A R T I N G S T A T E
#
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
###################### Starting state data ############################
switch(start$given,
none={
stop("No starting state given")
},
x = {
# x.start was given
# coerce it to a ppp object
if(!is.ppp(x.start))
x.start <- as.ppp(x.start, w.sim)
npts <- x.start$n
},
n = {
# n.start was given
n.start <- start$n.start
# Adjust the number of points in the starting state in accordance
# with the expansion that has occurred.
if(expanded)
n.start <- ceiling(n.start * area.owin(w.sim)/area.owin(w.clip))
#
npts <- sum(n.start) # The ``sum()'' is redundant if n.start
# is scalar; no harm, but.
},
stop("Internal error: start$given unrecognized"))
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
#
# C O N T R O L P A R A M E T E R S
#
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
################### Periodic boundary conditions #########################
# If periodic is TRUE we have to be simulating in a rectangular window.
periodic <- control$periodic
if(periodic && w.sim$type != "rectangle")
stop("Need rectangular window for periodic simulation.\n")
# parameter passed to Fortran:
period <-
if(periodic)
c(diff(w.sim$xrange), diff(w.sim$yrange))
else
c(-1,-1)
#### vector of proposal probabilities
if(!mtype)
ptypes <- 1
else {
ptypes <- control$ptypes
if(is.null(ptypes)) {
# default values
ptypes <- switch(start$given,
none = ,
n = rep(1/ntypes,ntypes),
x = table(marks(x.start, dfok=FALSE))/x.start$n
)
} else {
# Validate ptypes
if(length(ptypes) != ntypes | sum(ptypes) != 1)
stop("Argument ptypes is mis-specified.\n")
}
}
########################################################################
# Normalising constant for proposal density
#
# Integral of trend over the expanded window (or area of window):
# Iota == Integral Of Trend (or) Area.
if(trendy) {
if(verbose)
cat("Evaluating trend integral...")
tlist <- if(is.function(trend) || is.im(trend)) list(trend) else trend
tsummaries <- lapply(tlist,
function(x, w) {
tmp <- as.im(x, w)[w, drop=FALSE]
return(summary(tmp))
},
w=w.sim)
nbg <- unlist(lapply(tsummaries, function(x) { x$min < 0 }))
if(any(nbg))
stop("Trend has negative values")
iota <- unlist(lapply(tsummaries, function(x) { x$integral }))
tmax <- unlist(lapply(tsummaries, function(x) { x$max }))
} else {
iota <- area.owin(w.sim)
tmax <- NULL
}
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
#
# A.S. EMPTY PROCESS
#
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
a.s.empty <- FALSE
#
# Empty pattern, simulated conditional on n
#
if(npts == 0 && control$conditioning != "none") {
if(verbose)
cat("Initial pattern has 0 points, and simulation is conditional on the number of points - returning an empty pattern\n")
a.s.empty <- TRUE
}
#
# If beta = 0, the process is almost surely empty
#
if(all(model$par[["beta"]] < .Machine$double.eps)) {
if(control$conditioning == "none") {
# return empty pattern
if(verbose)
cat("beta = 0 implies an empty pattern\n")
a.s.empty <- TRUE
} else
stop("beta = 0 implies an empty pattern, but we are simulating conditional on a nonzero number of points")
}
if(a.s.empty) {
# create empty pattern, to be returned
empty <- ppp(numeric(0), numeric(0), window=w.clip)
if(mtype) {
vide <- factor(types[integer(0)], levels=types)
empty <- empty %mark% vide
}
}
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
#
# PACKAGE UP
#
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
######### Store decisions
Model <- model
Start <- start
Control <- control
Model$w <- w.clip
Model$types <- types
Control$expand <- if(expanded) w.sim else 1
Control$force.exp <- expanded
Control$force.noexp <- !expanded
Control$internal <- list(w.sim=w.sim,
ptypes=ptypes,
period=period)
Model$internal <- list(a.s.empty=a.s.empty,
empty=if(a.s.empty) empty else NULL,
mtype=mtype,
trendy=trendy,
iota=iota,
tmax=tmax)
Start$internal <- list(npts=npts)
InfoList <- list(model=Model, start=Start, control=Control)
class(InfoList) <- c("rmhInfoList", class(InfoList))
# go
rmhEngine(InfoList, verbose=verbose, reseed=FALSE, kitchensink=TRUE, ...)
}
#--------------- rmhEngine -------------------------------------------
#
# This is the interface to the Fortran code.
#
# InfoList is a list of pre-digested, validated arguments
# obtained from rmh.default.
#
# This function is called by rmh.default to generate one simulated
# realisation of the model.
# It's called repeatedly by ho.engine and qqplot.ppm to generate multiple
# realisations (saving time by not repeating the argument checking
# in rmh.default).
# arguments:
# reseed: whether to reset the random seed to a new, random value
# kitchensink: whether to tack InfoList on to the return value as an attribute
# preponly: whether to just return InfoList without simulating
#
# rmh.default digests arguments and calls rmhEngine with kitchensink=T
#
# qqplot.ppm first gets InfoList by calling rmh.default with preponly=T
# (which digests the model arguments and calls rmhEngine
# with preponly=T, returning InfoList),
# then repeatedly calls rmhEngine(InfoList) to simulate.
#
# -------------------------------------------------------
rmhEngine <- function(InfoList, ...,
verbose=FALSE, reseed=TRUE, kitchensink=FALSE,
preponly=FALSE) {
# Internal Use Only!
# This is the interface to the Fortran code.
if(!inherits(InfoList, "rmhInfoList"))
stop("data not in correct format for internal function rmhEngine")
if(preponly)
return(InfoList)
model <- InfoList$model
start <- InfoList$start
control <- InfoList$control
w.sim <- control$internal$w.sim
w.clip <- model$w
types <- model$types
ntypes <- length(types)
ptypes <- control$internal$ptypes
period <- control$internal$period
mtype <- model$internal$mtype
trend <- model$trend
trendy <- model$internal$trendy
iota <- model$internal$iota
tmax <- model$internal$tmax
npts <- start$internal$npts
n.start <- start$n.start
x.start <- start$x.start
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
#
# E M P T Y P A T T E R N
#
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
if(model$internal$a.s.empty) {
empty <- model$internal$empty
attr(empty, "info") <- InfoList
return(empty)
}
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
#
# S I M U L A T I O N
#
#==+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
#############################################
####
#### Random number initialisation
####
#############################################
if(reseed || is.null(start$seed)) {
# generate a (new) random seed
ss <- start$seed <- rmhseed()
if(kitchensink)
InfoList$start$seed <- start$seed
}
set.seed(start$seed$build.seed)
#############################################
####
#### Poisson case
####
#############################################
if(model$cif == 'poisson') {
intensity <- if(!trendy) model$par[["beta"]] else model$trend
Xsim <-
switch(control$conditioning,
none= {
# Poisson process
if(!mtype)
rpoispp(intensity, win=w.sim, ...)
else
rmpoispp(intensity, win=w.sim, types=types)
},
n.total = {
# Binomial/multinomial process with fixed total number of points
if(!mtype)
rpoint(npts, intensity, win=w.sim, verbose=verbose)
else
rmpoint(npts, intensity, win=w.sim, types=types,
verbose=verbose)
},
n.each.type = {
# Multinomial process with fixed number of points of each type
npts.each <-
switch(start$given,
n = n.start,
x = as.integer(table(marks(x.start, dfok=FALSE))),
stop("No starting state given; can't condition on fixed number of points"))
rmpoint(npts.each, intensity, win=w.sim, types=types,
verbose=verbose)
},
stop("Internal error: control$conditioning unrecognised")
)
Xclip <- Xsim[w.clip]
attr(Xclip, "info") <- InfoList
return(Xclip)
}
########################################################################
# M e t r o p o l i s H a s t i n g s s i m u l a t i o n
########################################################################
if(verbose)
cat("Starting simulation.\nInitial state...")
#### Build starting state
#### First the marks, if any.
#### The marks must be integers 1 to ntypes, for passing to Fortran
marks <-
if(!mtype)
0
else
switch(start$given,
n = {
# n.start given
if(control$conditioning=="n.each.type")
rep(1:ntypes,n.start)
else
sample(1:ntypes,npts,TRUE,ptypes)
},
x = {
# x.start given
as.integer(marks(x.start, dfok=FALSE))
},
stop("internal error: start$given unrecognised")
)
#
# First the x, y coordinates
#
switch(start$given,
x = {
x <- x.start$x
y <- x.start$y
},
n = {
xy <-
if(!trendy)
runifpoint(npts, w.sim, ...)
else
rpoint.multi(npts, trend, tmax,
factor(marks,levels=1:ntypes), w.sim, ...)
x <- xy$x
y <- xy$y
})
#######################################################################
# Start to call Fortran
######################################################################
# Get the parameters in Fortran-ese
par <- model$fortran.par
nmbr <- model$fortran.id
# Absorb the constants or vectors `iota' and 'ptypes' into the beta parameters
# This assumes that the beta parameters are the first 'ntypes' entries
# of the parameter vector passed to Fortran.
par[1:ntypes] <- (iota/ptypes) * par[1:ntypes]
# Algorithm control parameters
nrep <- control$nrep
cond <- control$cond
# If we are simulating a Geyer saturation process we need to set up some
# ``auxiliary information''.
need.aux <- model$need.aux
aux <-
if(!need.aux)
0
else
.Fortran(
"initaux",
nmbr=as.integer(nmbr),
par=as.double(par),
period=as.double(period),
x=as.double(x),
y=as.double(y),
npts=as.integer(npts),
aux=integer(npts),
PACKAGE="spatstat"
)$aux
# The vectors x and y (and perhaps marks) which hold the generated
# process may grow. We need to allow storage space for them to grow
# in. Unless we are conditioning on the number of points, we have no
# real idea how big they will grow. Hence we start off with storage
# space which has at least twice the length of the ``initial state'', and
# structure things so that the storage space may be incremented
# without losing the ``state'' which has already been generated.
if(cond == 1) {
nincr <- npad <- max(npts, 50)
padding <- numeric(npad)
x <- c(x, padding)
y <- c(y, padding)
if(ntypes>1) marks <- c(marks, padding)
if(need.aux) aux <- c(aux, padding)
} else {
nincr <- npad <- 0
}
npmax <- npts + npad
mrep <- 1
if(verbose)
cat("Proposal points...")
# If the pattern is multitype, generate the mark proposals.
mprop <- if(ntypes>1)
sample(1:ntypes,nrep,TRUE,prob=ptypes) else 0
# Generate the ``proposal points'' in the expanded window.
xy <-
if(trendy)
rpoint.multi(nrep,trend,tmax,factor(mprop, levels=1:ntypes),w.sim,...)
else
runifpoint(nrep,w.sim)
xprop <- xy$x
yprop <- xy$y
if(verbose)
cat("Start simulation.\n")
# Determine Fortran subroutine name (this is not very safe ...)
mhname <- paste("mh", nmbr, sep="")
# The repetition is to allow the storage space to be incremented if
# necessary.
repeat {
# Call the Metropolis-Hastings simulator:
rslt <- .Fortran(
mhname,
par=as.double(par),
period=as.double(period),
xprop=as.double(xprop),
yprop=as.double(yprop),
mprop=as.integer(mprop),
ntypes=as.integer(ntypes),
iseed=as.integer(start$seed$iseed),
nrep=as.integer(nrep),
mrep=as.integer(mrep),
p=as.double(control$p),
q=as.double(control$q),
npmax=as.integer(npmax),
nverb=as.integer(control$nverb),
x=as.double(x),
y=as.double(y),
marks=as.integer(marks),
aux=as.integer(aux),
npts=as.integer(npts),
fixall=as.logical(control$fixall),
PACKAGE="spatstat"
)
npts <- rslt$npts
mrep <- rslt$mrep
# If mrep > nrep we have completed the nrep Metropolis Hastings steps.
# Exit the loop.
if(mrep > nrep) break
# If not, then we came back early because we were about to run out
# of storage space. Increase the storage space and re-call the
# methas subroutine.
# Internal consistency check
if(npts < npmax)
stop(paste("Internal error; Fortran code exited unexpectedly\n",
"(without completing nrep steps and without reaching\n",
"the limit of storage capacity)\n"))
if(verbose) {
cat('Number of points equal to ',npmax,';\n',sep='')
cat('increasing storage space and continuing.\n')
}
npmax <- npmax + nincr
x <- c(rslt$x,numeric(nincr))
y <- c(rslt$y,numeric(nincr))
marks <- if(ntypes>1) c(rslt$marks,numeric(nincr)) else 0
aux <- if(need.aux) c(rslt$aux,numeric(nincr)) else 0
iseed <- rslt$iseed
}
###### END OF LOOP ###################
# Extract the point pattern returned from Fortran
npts <- rslt$npts
indices <- if(npts == 0) numeric(0) else (1:npts)
x <- rslt$x[indices]
y <- rslt$y[indices]
xxx <- ppp(x=x, y=y, window=w.sim, check=FALSE)
if(mtype) {
marx <- factor(rslt$marks[indices],levels=types)
xxx <- xxx %mark% marx
}
# Now clip the pattern to the ``clipping'' window:
xxx <- xxx[w.clip]
# Append to the result information about how it was generated.
if(kitchensink)
attr(xxx, "info") <- InfoList
return(xxx)
}