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
|
#
#
# badgey.S
#
# $Revision: 1.17 $ $Date: 2018/03/15 07:37:41 $
#
# 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"),
hasInf = FALSE,
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[,1L]
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[,2L]
Xindex <- EqualPairs[,1L]
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
})
|
