badgey.R
#
#
# badgey.S
#
# $Revision: 1.16 $ $Date: 2016/04/25 02:34:40 $
#
# Hybrid Geyer process
#
# BadGey() create an instance of the process
# [an object of class 'interact']
#
#
# -------------------------------------------------------------------
#
BadGey <- local({
# ........... auxiliary functions ..............
delBG <- function(i, r, sat) {
r <- r[-i]
if(length(r) == length(sat)) {
r <- r[-i]
sat <- sat[-i]
} else if(length(sat) == 1) {
r <- r[-i]
} else stop("Mismatch in dimensions of arguments r and sat")
nr <- length(r)
if(nr == 0) return(Poisson())
if(nr == 1) return(Geyer(r, sat))
return(BadGey(r, sat))
}
# .............. template ....................
BlankBG <-
list(
name = "hybrid Geyer process",
creator = "BadGey",
family = "pairsat.family", # will be evaluated later
pot = function(d, par) {
r <- par$r
nr <- length(r)
out <- array(FALSE, dim=c(dim(d), nr))
for(i in 1:nr)
out[,,i] <- (d <= r[i])
out
},
par = list(r = NULL, sat=NULL), # to fill in later
parnames = c("interaction radii", "saturation parameters"),
init = function(self) {
r <- self$par$r
sat <- self$par$sat
if(!is.numeric(r) || !all(r > 0))
stop("interaction radii r must be positive numbers")
if(length(r) > 1 && !all(diff(r) > 0))
stop("interaction radii r must be strictly increasing")
if(!is.numeric(sat) || any(sat < 0))
stop("saturation parameters must be nonnegative numbers")
if(length(sat) != length(r) && length(sat) != 1)
stop("vectors r and sat must have equal length")
},
update = NULL, # default OK
print = NULL, # default OK
interpret = function(coeffs, self) {
r <- self$par$r
npiece <- length(r)
# extract coefficients
gammas <- exp(as.numeric(coeffs))
# name them
gn <- gammas
names(gn) <- paste("[0,", r, ")", sep="")
#
return(list(param=list(gammas=gammas),
inames="interaction parameters gamma_i",
printable=dround(gn)))
},
valid = function(coeffs, self) {
# interaction parameters gamma must be
# non-NA
# finite, if sat > 0
# less than 1, if sat = Inf
gamma <- (self$interpret)(coeffs, self)$param$gammas
sat <- self$par$sat
if(anyNA(gamma))
return(FALSE)
return(all((is.finite(gamma) | sat == 0)
& (gamma <= 1 | sat != Inf)))
},
project = function(coeffs, self){
loggammas <- as.numeric(coeffs)
sat <- self$par$sat
r <- self$par$r
good <- is.finite(loggammas) & (is.finite(sat) | loggammas <= 0)
if(all(good))
return(NULL)
if(!any(good))
return(Poisson())
bad <- !good
if(spatstat.options("project.fast") || sum(bad) == 1) {
# remove smallest threshold with an unidentifiable parameter
firstbad <- min(which(bad))
return(delBG(firstbad, r, sat))
} else {
# consider all candidate submodels
subs <- lapply(which(bad), delBG, r=r, sat=sat)
return(subs)
}
},
irange = function(self, coeffs=NA, epsilon=0, ...) {
r <- self$par$r
sat <- self$par$sat
if(all(is.na(coeffs)))
return(2 * max(r))
gamma <- (self$interpret)(coeffs, self)$param$gammas
gamma[is.na(gamma)] <- 1
active <- (abs(log(gamma)) > epsilon) & (sat > 0)
if(!any(active))
return(0)
else return(2 * max(r[active]))
},
version=NULL, # to be added later
# fast evaluation is available for the border correction only
can.do.fast=function(X,correction,par) {
return(all(correction %in% c("border", "none")))
},
fasteval=function(X,U,EqualPairs,pairpot,potpars,correction,
..., halfway=FALSE) {
# fast evaluator for BadGey interaction
if(!all(correction %in% c("border", "none")))
return(NULL)
if(spatstat.options("fasteval") == "test")
message("Using fast eval for BadGey")
r <- potpars$r
sat <- potpars$sat
# ensure r and sat have equal length
if(length(r) != length(sat)) {
if(length(r) == 1)
r <- rep.int(r, length(sat))
else if(length(sat) == 1)
sat <- rep.int(sat, length(r))
else stop("lengths of r and sat do not match")
}
# first ensure all data points are in U
nX <- npoints(X)
nU <- npoints(U)
Xseq <- seq_len(nX)
if(length(EqualPairs) == 0) {
# no data points currently included
missingdata <- rep.int(TRUE, nX)
} else {
Xused <- EqualPairs[,1]
missingdata <- !(Xseq %in% Xused)
}
somemissing <- any(missingdata)
if(somemissing) {
# add the missing data points
nmiss <- sum(missingdata)
U <- superimpose(U, X[missingdata], W=X$window)
# correspondingly augment the list of equal pairs
originalrows <- seq_len(nU)
newXindex <- Xseq[missingdata]
newUindex <- nU + seq_len(nmiss)
EqualPairs <- rbind(EqualPairs, cbind(newXindex, newUindex))
nU <- nU + nmiss
}
nterms <- length(r)
answer <- matrix(, nrow=nU, ncol=nterms)
for(k in 1:nterms) {
# first determine saturated pair counts
counts <- strausscounts(U, X, r[k], EqualPairs)
satcounts <- pmin.int(sat[k], counts)
# trapdoor used by suffstat()
if(halfway)
answer[,k] <- satcounts
else if(sat[k] == Inf)
answer[,k] <- 2 * satcounts
else {
# extract counts for data points
Uindex <- EqualPairs[,2]
Xindex <- EqualPairs[,1]
Xcounts <- integer(npoints(X))
Xcounts[Xindex] <- counts[Uindex]
# evaluate change in saturated counts of other data points
change <- geyercounts(U, X, r[k], sat[k], Xcounts, EqualPairs)
answer[,k] <- satcounts + change
}
}
if(somemissing)
answer <- answer[originalrows, , drop=FALSE]
return(answer)
}
)
class(BlankBG) <- "interact"
BadGey <- function(r, sat) {
instantiate.interact(BlankBG, list(r=r, sat=sat))
}
BadGey <- intermaker(BadGey, BlankBG)
BadGey
})