#
# quadclass.S
#
# Class 'quad' to define quadrature schemes
# in (rectangular) windows in two dimensions.
#
# $Revision: 4.15 $ $Date: 2010/07/10 10:22:09 $
#
# An object of class 'quad' contains the following entries:
#
# $data: an object of class 'ppp'
# defining the OBSERVATION window,
# giving the locations (& marks) of the data points.
#
# $dummy: object of class 'ppp'
# defining the QUADRATURE window,
# giving the locations (& marks) of the dummy points.
#
# $w: vector giving the nonnegative weights for the
# data and dummy points (data first, followed by dummy)
#
# w may also have an attribute attr(w, "zeroes")
# equivalent to (w == 0). If this is absent
# then all points are known to have positive weights.
#
# $param:
# parameters that were used to compute the weights
# and possibly to create the dummy points (see below).
#
# The combined (data+dummy) vectors of x, y coordinates of the points,
# and their weights, are extracted using standard functions
# x.quad(), y.quad(), w.quad() etc.
#
# ----------------------------------------------------------------------
# Note about parameters:
#
# If the quadrature scheme was created by quadscheme(),
# then $param contains
#
# $param$weight
# list containing the values of all parameters
# actually used to compute the weights.
#
# $param$dummy
# list containing the values of all parameters
# actually used to construct the dummy pattern
# via default.dummy();
# or NULL if the dummy pattern was provided externally
#
# If you constructed the quadrature scheme manually, this
# structure may not be present.
#
#-------------------------------------------------------------
quad <- function(data, dummy, w, param=NULL) {
data <- as.ppp(data)
dummy <- as.ppp(dummy)
n <- data$n + dummy$n
if(missing(w))
w <- rep(1, n)
else {
w <- as.vector(w)
if(length(w) != n)
stop("length of weights vector w is not equal to total number of points")
}
if(is.null(attr(w, "zeroes")) && any( w == 0))
attr(w, "zeroes") <- (w == 0)
Q <- list(data=data, dummy=dummy, w=w, param=param)
class(Q) <- "quad"
invisible(Q)
}
# ------------------ extractor functions ----------------------
x.quad <- function(Q) {
verifyclass(Q, "quad")
c(Q$data$x, Q$dummy$x)
}
y.quad <- function(Q) {
verifyclass(Q, "quad")
c(Q$data$y, Q$dummy$y)
}
w.quad <- function(Q) {
verifyclass(Q, "quad")
Q$w
}
param.quad <- function(Q) {
verifyclass(Q, "quad")
Q$param
}
n.quad <- function(Q) {
verifyclass(Q, "quad")
Q$data$n + Q$dummy$n
}
marks.quad <- function(x, dfok=FALSE, ...) {
verifyclass(x, "quad")
dat <- x$data
dum <- x$dummy
if(dfok) warning("ignored dfok = TRUE; not implemented")
mdat <- marks(dat, dfok=FALSE, ...)
mdum <- marks(dum, dfok=FALSE, ...)
if(is.null(mdat) && is.null(mdum))
return(NULL)
if(is.null(mdat))
mdat <- rep(NA, dat$n)
if(is.null(mdum))
mdum <- rep(NA, dum$n)
mall <- c(mdat, mdum)
if(is.factor(mdat) && is.factor(mdum) && all(levels(mdat) == levels(mdum))) {
mall <- factor(mall)
levels(mall) <- levels(mdat)
}
return(mall)
}
is.marked.quad <- function(X, na.action="warn", ...) {
marx <- marks(X, ...)
if(is.null(marx))
return(FALSE)
if(any(is.na(marx)))
switch(na.action,
warn = {
warning(paste("some mark values are NA in the point pattern",
deparse(substitute(X))))
},
fatal = {
return(FALSE)
},
ignore = {}
)
return(TRUE)
}
is.multitype.quad <- function(X, na.action="warn", ...) {
marx <- marks(X, ...)
if(is.null(marx))
return(FALSE)
if(any(is.na(marx)))
switch(na.action,
warn = {
warning(paste("some mark values are NA in the point pattern",
deparse(substitute(X))))
},
fatal = {
return(FALSE)
},
ignore = {}
)
return(!is.data.frame(marx) && is.factor(marx))
}
is.data <- function(Q) {
verifyclass(Q, "quad")
return(c(rep(TRUE, Q$data$n),
rep(FALSE, Q$dummy$n)))
}
equals.quad <- function(Q) {
# return matrix E such that E[i,j] = (X[i] == U[j])
# where X = Q$data and U = union.quad(Q)
n <- Q$data$n
m <- Q$dummy$n
E <- matrix(FALSE, nrow=n, ncol=n+m)
diag(E) <- TRUE
E
}
equalsfun.quad <- function(Q) {
stopifnot(inherits(Q, "quad"))
return(function(i,j) { i == j })
}
equalpairs.quad <- function(Q) {
# return two-column matrix E such that
# X[E[i,1]] == U[E[i,2]] for all i
# where X = Q$data and U = union.quad(Q)
n <- Q$data$n
return(matrix(rep(1:n,2), ncol=2))
}
union.quad <- function(Q) {
verifyclass(Q, "quad")
ppp(x= c(Q$data$x, Q$dummy$x),
y= c(Q$data$y, Q$dummy$y),
window=Q$dummy$window,
marks=marks.quad(Q),
check=FALSE)
}
#
# Plot a quadrature scheme
#
#
plot.quad <- function(x, ..., main=deparse(substitute(x)), dum=list()) {
verifyclass(x, "quad")
data <- x$data
dummy <- x$dummy
dummyplot <- function(x, ..., pch=".", add=TRUE) {
plot(x, pch=pch, add=add, ...)
}
if(!is.marked(data)) {
plot(data, main=main, ...)
do.call("dummyplot", append(list(
dummy,
main=paste(main, "\n dummy points")
),
dum))
} else if(is.multitype(data)) {
oldpar <- par(ask = interactive() &&
(.Device %in% c("X11", "GTK", "windows", "Macintosh")))
on.exit(par(oldpar))
data.marks <- marks(data)
dummy.marks <- marks(dummy)
types <- levels(data.marks)
for(k in types) {
maink <- paste(main, "\n mark = ", k, sep="")
plot(unmark(data[data.marks == k]), main=maink, ...)
do.call("dummyplot", append(list(unmark(dummy[dummy.marks == k])), dum))
}
} else {
plot(data, ..., main=main)
addplot <- function(x, ..., add=TRUE, main=deparse(substitute(x))) {
plot(x, ..., main=main, add=add)
}
do.call("addplot", append(list(dummy), dum))
}
invisible(NULL)
}
# subset operator
"[.quad" <- function(x, ...) {
U <- union.quad(x)
Z <- is.data(x)
w <- w.quad(x)
# determine serial numbers of points to be included
V <- U %mark% seq(U$n)
i <- marks(V[...])
# extract corresponding subsets of vectors
Z <- Z[i]
w <- w[i]
# take subset of points, using any type of subset index
U <- U[...]
# stick together
quad(U[Z], U[!Z], w)
}