#
#	quadclass.S
#
#	Class 'quad' to define quadrature schemes
#	in (rectangular) windows in two dimensions.
#
#	$Revision: 4.14 $	$Date: 2009/07/28 00:55:25 $
#
# 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 })
}
  
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)
}