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
|
#
#
# hardcore.S
#
# $Revision: 1.16 $ $Date: 2022/03/07 02:06:12 $
#
# The Hard core process
#
# Hardcore() create an instance of the Hard Core process
# [an object of class 'interact']
#
#
# -------------------------------------------------------------------
#
Hardcore <- local({
BlankHardcore <-
list(
name = "Hard core process",
creator = "Hardcore",
family = "pairwise.family", # evaluated later
pot = function(d, par) {
v <- 0 * d
v[ d <= par$hc ] <- (-Inf)
attr(v, "IsOffset") <- TRUE
v
},
par = list(hc = NULL), # filled in later
parnames = "hard core distance",
hasInf = TRUE,
selfstart = function(X, self) {
# self starter for Hardcore
nX <- npoints(X)
if(nX < 2) {
# not enough points to make any decisions
return(self)
}
md <- minnndist(X)
if(md == 0) {
warning(paste("Pattern contains duplicated points:",
"impossible under Hardcore model"))
return(self)
}
if(!is.na(hc <- self$par$hc)) {
# value fixed by user or previous invocation
# check it
if(md < hc)
warning(paste("Hard core distance is too large;",
"some data points will have zero probability"))
return(self)
}
# take hc = minimum interpoint distance * n/(n+1)
hcX <- md * nX/(nX+1)
Hardcore(hc = hcX)
},
init = function(self) {
hc <- self$par$hc
if(length(hc) != 1)
stop("hard core distance must be a single value")
if(!is.na(hc) && !(is.numeric(hc) && hc > 0))
stop("hard core distance hc must be a positive number, or NA")
},
update = NULL, # default OK
print = NULL, # default OK
interpret = function(coeffs, self) {
return(NULL)
},
valid = function(coeffs, self) {
return(TRUE)
},
project = function(coeffs, self) {
return(NULL)
},
irange = function(self, coeffs=NA, epsilon=0, ...) {
return(self$par$hc)
},
hardcore = function(self, coeffs=NA, epsilon=0, ...) {
return(self$par$hc)
},
version=NULL, # evaluated 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,
splitInf=FALSE, ...) {
# fast evaluator for Hardcore interaction
if(!all(correction %in% c("border", "none")))
return(NULL)
if(spatstat.options("fasteval") == "test")
message("Using fast eval for Hardcore")
hc <- potpars$hc
# call evaluator for Strauss process
counts <- strausscounts(U, X, hc, EqualPairs)
forbid <- (counts != 0)
if(!splitInf) {
## usual case
v <- matrix(ifelseAB(forbid, -Inf, 0), ncol=1L)
} else {
## separate hard core
v <- matrix(0, nrow=npoints(U), ncol=1L)
attr(v, "-Inf") <- forbid
}
attr(v, "IsOffset") <- TRUE
return(v)
},
Mayer=function(coeffs, self) {
# second Mayer cluster integral
hc <- self$par$hc
return(pi * hc^2)
},
Percy=function(d, coeffs, par, ...) {
## term used in Percus-Yevick type approximation
H <- par$hc
t <- abs(d/(2*H))
t <- pmin.int(t, 1)
y <- 2 * H^2 * (pi - acos(t) + t * sqrt(1 - t^2))
return(y)
}
)
class(BlankHardcore) <- "interact"
Hardcore <- function(hc=NA) {
instantiate.interact(BlankHardcore, list(hc=hc))
}
Hardcore <- intermaker(Hardcore, BlankHardcore)
Hardcore
})
|
