[Rd] sweep sanity checking?

[Petr Savicky]

When I was preparing the patch of sweep submitted on July 25, I was
unaware of the code by Heather Turner. She suggested a very elegant
solution, if STATS is a vector and we want to use meaningful recycling
in full generality. I would like to suggest a combined solution, which
uses Heather Turner's algorithm if check.margin=FALSE (default) and STATS
is a vector and my previous algorithm, if check.margin=TRUE or STATS is
an array. The suggestion is

  # combined from the original code of sweep without warnings and from
  # https://stat.ethz.ch/pipermail/r-help/2005-June/073989.html by Robin Hankin
  # https://stat.ethz.ch/pipermail/r-help/2005-June/074001.html by Heather Turner
  # https://stat.ethz.ch/pipermail/r-devel/2007-June/046217.html by Ben Bolker
  # with some further modifications by Petr Savicky
  sweep <- function(x, MARGIN, STATS, FUN = "-", check.margin=FALSE, ...)
  {
      FUN <- match.fun(FUN)
      dims <- dim(x)
      dimmargin <- dims[MARGIN]
      if (is.null(dim(STATS))) {
          dimstats <- length(STATS)
      } else {
          dimstats <- dim(STATS)
          check.margin <- TRUE
      }
      s <- length(STATS)
      if (s > prod(dimmargin)) {
          warning("length of STATS greater than the extent of dim(x)[MARGIN]")
      } else if (check.margin) {
          dimmargin <- dimmargin[dimmargin > 1]
          dimstats <- dimstats[dimstats > 1]
          if (length(dimstats) > length(dimmargin) ||
              any(dimstats != dimmargin[seq(along.with=dimstats)]))
              warning("length(STATS) or dim(STATS) do not match dim(x)[MARGIN]")
      } else {
          cumDim <- c(1, cumprod(dimmargin))
          upper <- min(cumDim[cumDim >= s])
          lower <- max(cumDim[cumDim <= s])
          if (upper %% s != 0 || s %% lower != 0)
              warning("STATS does not recycle exactly across MARGIN")
      }
      perm <- c(MARGIN, (1:length(dims))[ - MARGIN])
      FUN(x, aperm(array(STATS, dims[perm]), order(perm)), ...)
  }

Heather presented four examples testing her code:
  sweep(array(1:24, dim = c(4,3,2)), 1, 1:2)    # no warning
  sweep(array(1:24, dim = c(4,3,2)), 1, 1:12)   # no warning
  sweep(array(1:24, dim = c(4,3,2)), 1, 1:24)   # no warning
  sweep(array(1:24, dim = c(4,3,2)), 1:2, 1:3)  # warning
The second and third example are not really correct, since STATS extends
also to dimensions not included in MARGIN. The problem is better visible
for example in
  sweep(array(1:24, dim = c(4,4,3,3,2,2)), c(1,3), 1:12)
where MARGIN clearly has to contain two dimensions explicitly.
So, I use the examples with a larger margin corresponding to STATS
as follows
  sweep(array(1:24, dim = c(4,3,2)), 1, 1:2)    # no warning
  sweep(array(1:24, dim = c(4,3,2)), 1:2, 1:12) # no warning
  sweep(array(1:24, dim = c(4,3,2)), 1:3, 1:24) # no warning
  sweep(array(1:24, dim = c(4,3,2)), 1:2, 1:3)  # warning
The current proposal for sweep indeed gives no warning in the first
three examples and gives a warning in the last one.

I did not use the suggestion to call the option warn with default
warn = getOption("warn"). The reason is that there are two
different decisions:
 (1) whether to generate a warning
 (2) what to do with the warning, if it is generated.
The warn option influences (2): the warning may be suppressed,
printed after the return to the top level, printed immediately or
it may be converted to an error. I think that the option
controling (2) should not be mixed with an option which
controls (1). If R has an option controling to which extent
recycling is allowed, then this could be used, but warn
has a different purpose.

I appreciate feedback.

Petr.