#
#
#    pairsat.family.S
#
#    $Revision: 1.15 $	$Date: 2004/01/08 15:02:42 $
#
#    The saturated pairwise interaction family of point process models
#
#    (an extension of Geyer's saturation process to all pairwise interactions)
#
#    pairsat.family:         object of class 'isf'
#                     defining saturated pairwise interaction
#	
#
# -------------------------------------------------------------------
#	

pairsat.family <-
  list(
         name  = "saturated pairwise",
         print = function(self) {
                      cat("Saturated pairwise interaction family\n")
         },
         eval  = function(X,U,Equal,pairpot,potpars,correction) {
  #
  # This auxiliary function is not meant to be called by the user.
  # It computes the distances between points,
  # evaluates the pair potential and applies edge corrections.
  #
  # Arguments:
  #   X           data point pattern                      'ppp' object
  #   U           points at which to evaluate potential   list(x,y) suffices
  #   Equal       logical matrix X[i] == U[j]             matrix or NULL
  #   pairpot     potential function (see above)          function(d, p)
  #   potpars     auxiliary parameters for pairpot        list(......)
  #   correction  edge correction type                    (string)
  #
  #   Note the Geyer saturation threshold must be given in 'potpars$saturate'
  #
  # Value:
  #    matrix of values of the total pair potential
  #    induced by the pattern X at each location given in U.
  #    The rows of this matrix correspond to the rows of U (the sample points);
  #    the k columns are the coordinates of the k-dimensional potential.
  #
  # Note:
  # The pair potential function 'pairpot' will be called as
  #    pairpot(M, potpars)   where M is a matrix of interpoint distances.
  # It must return a matrix with the same dimensions as M
  # or an array with its first two dimensions the same as the dimensions of M.
  ##########################################################################

# coercion should be unnecessary, but this is useful for debugging
X <- as.ppp(X)
U <- as.ppp(U, X$window)   # i.e. X$window is DEFAULT window

# saturation parameter
saturate <- potpars$saturate

if(is.null(saturate)) {
  # pairwise interaction 
  V <- pairwise.family$eval(X, U, Equal, pairpot, potpars, correction)
  return(V)
}

# first check that all data points are included in the quadrature points
missingdata <- !as.vector(matrowany(Equal))
somemissing <- any(missingdata)
if(somemissing) {
  # add the missing data points
  U <- superimpose(U, X[missingdata])
  # extend the equality matrix to maintain Equal[i,j] = (X[i] == U[j])
  originalcolumns <- seq(ncol(Equal))
  sn <- seq(X$n)
  id <- outer(sn, sn[missingdata], "==")
  Equal <- cbind(Equal, id)
}

# compute the pair potentials POT and the unsaturated potential sums V

V <- pairwise.family$eval(X, U, Equal, pairpot, potpars, correction)
POT <- attr(V, "POT")

#################################################################
################## saturation part ##############################
#################################################################

#
# (a) compute SATURATED potential sums
V.sat <- array(pmin(saturate, V), dim=dim(V))

#
# (b) compute effect of addition/deletion of dummy/data point j
# on the UNSATURATED potential sum of each data point i
#
# Identify data points
     # is.data <- apply(Equal, 2, any)
is.data <- as.vector(matcolany(Equal))  # logical vector corresp. to rows of V

# Extract potential sums for data points only
V.data <- V[is.data, , drop=FALSE]

# replicate them so that V.dat.rep[i,j,k] = V.data[i, k]
V.dat.rep <- aperm(array(V.data, dim=c(dim(V.data), U$n)), c(1,3,2))

# make a logical array   col.is.data[i,j,k] = is.data[j]
dip <- dim(POT)
mat <- matrix(, nrow=dip[1], ncol=dip[2])
izdat <- is.data[col(mat)] 
col.is.data <- array(izdat, dim=dip) # automatically replicates

# compute value of unsaturated potential sum for each data point i
# obtained after addition/deletion of each dummy/data point j
                                  
V.after <- V.dat.rep + ifelse(col.is.data, -POT, POT)
#
#
# (c) difference of SATURATED potential sums for each data point i
# before & after increment/decrement of each dummy/data point j
#
# saturated values after increment/decrement
V.after.sat <- array(pmin(saturate, V.after), dim=dim(V.after))
# saturated values before
V.dat.rep.sat <- array(pmin(saturate, V.dat.rep), dim=dim(V.dat.rep))
# difference
V.delta <- V.after.sat - V.dat.rep.sat
V.delta <- ifelse(col.is.data, -V.delta, V.delta)
#
# (d) Sum (c) over all data points i
V.delta.sum <- apply(V.delta, c(2,3), sum)
#
# (e) Result
V <- V.sat + V.delta.sum

##########################################
# remove any columns that were added
if(somemissing)
      V <- V[originalcolumns, , drop=FALSE]

return(V)

}     ######### end of function $eval                            

)     ######### end of list

class(pairsat.family) <- "isf"