#'
#' sparse3Darray.R
#'
#' Sparse 3D arrays represented as list(i,j,k,x)
#' 
#' Copyright (c) Adrian Baddeley, Ege Rubak and Rolf Turner 2016-2020
#' GNU Public Licence >= 2.0
#'
#' $Revision: 1.45 $  $Date: 2023/06/23 02:26:26 $
#'

sparse3Darray <- function(i=integer(0), j=integer(0), k=integer(0),
                          x=numeric(0),
                          dims=c(max(i),max(j),max(k)),
                          dimnames=NULL, strict=FALSE, nonzero=FALSE) {
  dat <- data.frame(i=as.integer(i),
                    j=as.integer(j),
                    k=as.integer(k),
                    x=x)
  if(typeof(x) == "complex") 
    warning(paste(
      "complex-valued sparse 3D arrays are supported in spatstat,",
      "but complex-valued sparse matrices",
      "are not yet supported by the Matrix package"),
      call.=FALSE)
  stopifnot(length(dims) == 3)
  dims <- as.integer(dims)
  if(!all(i >= 1 & i <= dims[1])) stop("indices i are outside range")
  if(!all(j >= 1 & j <= dims[2])) stop("indices j are outside range")
  if(!all(k >= 1 & k <= dims[3])) stop("indices k are outside range")
  if(!is.null(dimnames)) {
    stopifnot(is.list(dimnames))
    stopifnot(length(dimnames) == 3)
    notnull <- !sapply(dimnames, is.null)
    dimnames[notnull] <- lapply(dimnames[notnull], as.character)
  }
  if(nonzero || strict) {
    #' drop zeroes
    ok <- (x != RelevantZero(x))
    dat <- dat[ok, , drop=FALSE]
  }
  if(strict) {
    #' arrange in 'R order'
    dat <- dat[with(dat, order(k,j,i)), , drop=FALSE]
    #' duplicates will be adjacent
    dup <- with(dat, c(FALSE, diff(i) == 0 & diff(j) == 0 & diff(k) == 0))
    if(any(dup)) {
      #' accumulate values at the same array location
      retain <- !dup
      newrow <- cumsum(retain)
      newx <- as(tapply(dat$x, newrow, sum), typeof(dat$x))
      newdat <- dat[retain,,drop=FALSE]
      newdat$x <- newx
      dat <- newdat
    }
  }
  result <- append(as.list(dat),
                   list(dim=dims, dimnames=dimnames))
  class(result) <- "sparse3Darray"
  return(result)
}

as.sparse3Darray <- function(x, ...) {
  if(inherits(x, "sparse3Darray")) {
    y <- x
  } else if(inherits(x, c("matrix", "sparseMatrix"))) {
    z <- as(x, Class="TsparseMatrix")
    dn <- dimnames(x)
    dn <- if(is.null(dn)) NULL else c(dn, list(NULL))
    one <- if(length(z@i) > 0) 1L else integer(0)
    y <- sparse3Darray(i=z@i + 1L, j=z@j + 1L, k=one, x=z@x,
                       dims=c(dim(x), 1L), dimnames=dn)
  } else if(is.array(x)) {
    stopifnot(length(dim(x)) == 3)
    dimx <- dim(x)
    if(prod(dimx) == 0) {
      y <- sparse3Darray(, dims=dimx, dimnames=dimnames(x))
    } else {
      ijk <- which(x != RelevantZero(x), arr.ind=TRUE)
      ijk <- cbind(as.data.frame(ijk), x[ijk])
      y <- sparse3Darray(i=ijk[,1L], j=ijk[,2L], k=ijk[,3L], x=ijk[,4L],
                         dims=dimx, dimnames=dimnames(x))
    }
  } else if(inherits(x, "sparseVector")) {
    one <- if(length(x@i) > 0) 1L else integer(0)
    y <- sparse3Darray(i=x@i, j=one, k=one, x=x@x,
                       dims=c(x@length, 1L, 1L))
  } else if(is.null(dim(x)) && is.atomic(x)) {
    n <- length(x)
    dn <- names(x)
    if(!is.null(dn)) dn <- list(dn, NULL, NULL)
    one <- if(n > 0) 1L else integer(0)
    y <- sparse3Darray(i=seq_len(n), j=one, k=one, x=x,
                       dims=c(n, 1L, 1L), dimnames=dn)
    
  } else if(is.list(x) && length(x) > 0) {
    n <- length(x)
    if(all(sapply(x, is.matrix))) {
      z <- Reduce(abind, x)
      y <- as.sparse3Darray(z)
    } else if(all(sapply(x, inherits, what="sparseMatrix"))) {
      dimlist <- unique(lapply(x, dim))
      if(length(dimlist) > 1) stop("Dimensions of matrices do not match")
      dimx <- c(dimlist[[1L]], n)
      dnlist <- lapply(x, dimnames)
      isnul <- sapply(dnlist, is.null)
      dnlist <- unique(dnlist[!isnul])
      if(length(dnlist) > 1) stop("Dimnames of matrices do not match")
      dn <- if(length(dnlist) == 0) NULL else c(dnlist[[1L]], list(NULL))
      for(k in seq_len(n)) {
        mk <- as(x[[k]], "TsparseMatrix")
        kvalue <- if(length(mk@i) > 0) k else integer(0)
        dfk <- data.frame(i=mk@i + 1L, j=mk@j + 1L, k=kvalue, x=mk@x)
        df <- if(k == 1) dfk else rbind(df, dfk)
      }
      y <- sparse3Darray(i=df$i, j=df$j, k=df$k, x=df$x,
                         dims=dimx, dimnames=dn)
    } else {
      warning("I don't know how to convert a list to a sparse array")
      return(NULL)
    }
  } else {          
    warning("I don't know how to convert x to a sparse array")
    return(NULL)
  }
  return(y)
}

dim.sparse3Darray <- function(x) { x$dim }

"dim<-.sparse3Darray" <- function(x, value) {
  stopifnot(length(value) == 3)
  if(!all(inside.range(x$i, c(1, value[1]))))
    stop("indices i are outside new range")
  if(!all(inside.range(x$j, c(1, value[2]))))
    stop("indices j are outside new range")
  if(!all(inside.range(x$k, c(1, value[3]))))
    stop("indices k are outside new range")
  dimx <- dim(x)
  x$dim <- value
  if(!is.null(dimnames(x))) {
    dn <- dimnames(x)
    for(n in 1:3) {
      if(value[n] < dimx[n]) dn[[n]] <- dn[[n]][1:value[n]] else
      if(value[n] > dimx[n]) dn[n] <- list(NULL)
    }
    dimnames(x) <- dn
  }
  return(x)
}

dimnames.sparse3Darray <- function(x) { x$dimnames }

"dimnames<-.sparse3Darray" <- function(x, value) {
  if(!is.list(value)) value <- list(value)
  if(length(value) == 1) value <- rep(value, 3)
  x$dimnames <- value
  return(x)
}

print.sparse3Darray <- function(x, ...) {
  dimx <- dim(x)
  cat("Sparse 3D array of dimensions", paste(dimx, collapse="x"), fill=TRUE)
  if(prod(dimx) == 0)
    return(invisible(NULL))
  dn <- dimnames(x) %orifnull% rep(list(NULL), 3)
  d3 <- dimx[3]
  dn3 <- dn[[3]] %orifnull% as.character(seq_len(d3))
  df <- data.frame(i=x$i, j=x$j, k=x$k, x=x$x)
  pieces <- split(df, factor(df$k, levels=1:d3))
  dim2 <- dimx[1:2]
  dn2 <- dn[1:2]
  if(typeof(x$x) == "complex") {
    splat("Complex-valued")
    splat("\t\tReal component:")
    stuff <- capture.output(eval(Re(x)))
    cat(stuff[-1],sep="\n")
    splat("\n\n\t\tImaginary component:")
    stuff <- capture.output(eval(Im(x)))
    cat(stuff[-1],sep="\n")
  } else {
    for(k in seq_along(pieces)) {
      cat(paste0("\n\t[ , , ", dn3[k], "]\n\n"))
      Mi <- with(pieces[[k]],
                 sparseMatrix(i=i, j=j, x=x, dims=dim2, dimnames=dn2))
      stuff <- capture.output(eval(Mi))
      #' Remove 'sparse Matrix' header blurb
      stuff <- stuff[-1]
      if(is.blank(stuff[1]))
        stuff <- stuff[-1]
      cat(stuff, sep="\n")
    }
  }
  return(invisible(NULL))
}

aperm.sparse3Darray <- function(a, perm=NULL, resize=TRUE, ...) {
  if(is.null(perm)) return(a)
  stopifnot(length(perm) == 3)
  a <- unclass(a)
  a[c("i", "j", "k")] <- a[c("i", "j", "k")][perm]
  if(resize) {
    a$dim <- a$dim[perm]
    if(length(a$dimnames)==3) a$dimnames <- a$dimnames[perm]
  }
  class(a) <- c("sparse3Darray", class(a))
  return(a)
}

as.array.sparse3Darray <- function(x, ...) {
  zerovalue <- vector(mode=typeof(x$x), length=1L)
  z <- array(zerovalue, dim=dim(x), dimnames=dimnames(x))
  z[cbind(x$i,x$j,x$k)] <- x$x
  return(z)
}

"[.sparse3Darray" <- local({

  Extract <- function(x, i,j,k, drop=TRUE, ...) {
    dimx <- dim(x)
    dn <- dimnames(x) %orifnull% rep(list(NULL), 3)
    if(!missing(i) && length(dim(i)) == 2) {
      ## matrix index
      i <- as.matrix(i)
      if(!(missing(j) && missing(k)))
        stop("If i is a matrix, j and k should not be given", call.=FALSE)
      if(ncol(i) != 3)
        stop("If i is a matrix, it should have 3 columns", call.=FALSE)
      ## start with vector of 'zero' answers of the correct type
      answer <- sparseVector(x=RelevantEmpty(x$x),
                             i=integer(0),
                             length=nrow(i))
      ## values outside array return NA
      if(anybad <- !all(good <- inside3Darray(dim(x), i))) {
        bad <- !good
        answer[bad] <- NA
      }
      ## if entire array is zero, there is nothing to match
      if(length(x$x) == 0)
        return(answer)
      ## restrict attention to entries inside array
      igood <- if(anybad) i[good, , drop=FALSE] else i
      ## match desired indices to sparse entries
      varies <- (dimx > 1)
      nvary <- sum(varies)
      varying <- which(varies)
      if(nvary == 3) {
        ## ---- older code -----
        ## convert triples of integers to character codes
        #### icode <- apply(i, 1, paste, collapse=",") << is too slow >>
        ## icode <- paste(i[,1], i[,2], i[,3], sep=",")
        ## dcode <- paste(x$i, x$j, x$k, sep=",")
	## ------------------
	mgood <- matchIntegerDataFrames(igood, cbind(x$i, x$j, x$k))
      } else if(nvary == 2) {
        ## effectively a sparse matrix
        ## ---- older code -----
        ## icode <- paste(i[,varying[1]], i[,varying[2]], sep=",")
        ## ijk <- cbind(x$i, x$j, x$k)
        ## dcode <- paste(ijk[,varying[1]], ijk[,varying[2]], sep=",")
	## ------------------
	ijk <- cbind(x$i, x$j, x$k)
	mgood <- matchIntegerDataFrames(igood[,varying,drop=FALSE],
                                        ijk[,varying,drop=FALSE])
      } else if(nvary == 1) {
        ## effectively a sparse vector
        ## ---- older code -----
        ## icode <- i[,varying]
        ## dcode <- switch(varying, x$i, x$j, x$k)
	## ------------------
	mgood <- match(igood[,varying], switch(varying, x$i, x$j, x$k))
      } else {
        ## effectively a single value
        ## ---- older code -----
        ## icode <- rep(1, nrow(i))
        ## dcode <- 1  # since we know length(x$x) > 0
	mgood <- 1
      }
      ## insert any found elements
      found <- logical(nrow(i))
      found[good] <- foundgood <- !is.na(mgood)
      answer[found] <- x$x[mgood[foundgood]]
      return(answer)
    }
    if(!(missing(i) && missing(j) && missing(k))) {
      I <- grokIndexVector(if(missing(i)) NULL else i, dimx[1], dn[[1]])
      J <- grokIndexVector(if(missing(j)) NULL else j, dimx[2], dn[[2]])
      K <- grokIndexVector(if(missing(k)) NULL else k, dimx[3], dn[[3]])
      IJK <- list(I,J,K)
      if(!all(sapply(lapply(IJK, getElement, name="full"), is.null))) {
        ## some indices exceed array bounds;
        ## result is a full array containing NA's
        fullindices   <- lapply(IJK, fullIndexSequence)
        strictindices <- lapply(IJK, strictIndexSequence)
        result <- array(data=RelevantNA(x$x), dim=lengths(fullindices))
        matches <- mapply(match, x=fullindices, table=strictindices)
        ok <- lapply(lapply(matches, is.na), "!")
        result[ok[[1]], ok[[2]], ok[[3]]] <-
          as.array(x)[matches[[1]][ok[[1]]],
                      matches[[2]][ok[[2]]],
                      matches[[3]][ok[[3]]]]
        if(drop) result <- result[,,,drop=TRUE]
        return(result)
      }
      IJK <- lapply(IJK, getElement, name="strict")
      I <- IJK[[1]]
      J <- IJK[[2]]
      K <- IJK[[3]]
      #' number of values to be returned along each margin
      newdims <- sapply(IJK, getElement, name="n")
      #' dimnames of return array
      newdn <- lapply(IJK, getElement, name="s")
      #' find all required data (not necessarily in required order)
      inI <- I$lo
      inJ <- J$lo
      inK <- K$lo
      df <- data.frame(i=x$i, j=x$j, k=x$k, x=x$x)
      use <- with(df, inI[i] & inJ[j] & inK[k])
      df <- df[use, ,drop=FALSE]
      #' contract sub-array to (1:n) * (1:m) * (1:l)
      df <- transform(df,
                      i = cumsum(inI)[i],
                      j = cumsum(inJ)[j],
                      k = cumsum(inK)[k])
      Imap <- I$map
      Jmap <- J$map
      Kmap <- K$map
      if(nrow(df) == 0 || (is.null(Imap) && is.null(Jmap) && is.null(Kmap))) {
        ## return values are already in correct position
        outdf <- df
      } else {
        #' invert map to determine output positions (reorder/repeat entries)
        snI <- seq_len(I$n)
        snJ <- seq_len(J$n)
        snK <- seq_len(K$n)
        imap <- Imap %orifnull% snI
        jmap <- Jmap %orifnull% snJ
        kmap <- Kmap %orifnull% snK
        whichi <- split(seq_along(imap), factor(imap, levels=snI))
        whichj <- split(seq_along(jmap), factor(jmap, levels=snJ))
        whichk <- split(seq_along(kmap), factor(kmap, levels=snK))
        dat.i <- whichi[df$i]
        dat.j <- whichj[df$j]
        dat.k <- whichk[df$k]
        stuff <- mapply(expandwithdata,
                        i=dat.i, j=dat.j, k=dat.k, x=df$x,
                        SIMPLIFY=FALSE)
        outdf <- rbindCompatibleDataFrames(stuff)
      }
      x <- sparse3Darray(i=outdf$i, j=outdf$j, k=outdf$k, x=outdf$x,
                         dims=newdims, dimnames=newdn)
      dimx <- newdims
      dn <- newdn
    }
    if(drop) {
      retain <- (dimx > 1)
      nretain <- sum(retain)
      if(nretain == 2) {
        #' result is a matrix
        retained <- which(retain)
        newi <- getElement(x, name=c("i","j","k")[ retained[1] ])
        newj <- getElement(x, name=c("i","j","k")[ retained[2] ])
        newdim <- dimx[retain]
        newdn <- dn[retain]
        return(sparseMatrix(i=newi, j=newj, x=x$x, dims=newdim, dimnames=newdn))
      } else if(nretain == 1) {
        #' sparse vector
        retained <- which(retain)
        newi <- getElement(x, name=c("i","j","k")[retained])
        #' ensure 'strict' 
        ord <- order(newi)
        newi <- newi[ord]
        newx <- x$x[ord]
        if(any(dup <- c(FALSE, diff(newi) == 0))) {
          retain <- !dup
          ii <- cumsum(retain)
          newi <- newi[retain]
          newx <- as(tapply(newx, ii, sum), typeof(newx))
        }
        x <- sparseVector(x=newx, i=newi, length=dimx[retained])
      } else if(nretain == 0) {
        #' single value
        x <- as.vector(as.array(x))
      }
    }
    return(x)
  }

  expandwithdata <- function(i, j, k, x) {
    z <- expand.grid(i=i, j=j, k=k)
    if(nrow(z) > 0)
      z$x <- x
    return(z)
  }

  Extract
})


rbindCompatibleDataFrames <- function(x) {
  #' faster version of Reduce(rbind, x) when entries are known to be compatible
  nama2 <- colnames(x[[1]])
  y <- vector(mode="list", length=length(nama2))
  names(y) <- nama2
  for(nam in nama2)
    y[[nam]] <- unlist(lapply(x, getElement, name=nam))
  return(as.data.frame(y))
}


"[<-.sparse3Darray" <- function(x, i, j, k, ..., value) {
  dimx <- dim(x)
  dn <- dimnames(x) %orifnull% rep(list(NULL), 3)
  #' interpret indices
  if(!missing(i) && length(dim(i)) == 2) {
    ## matrix index
    ijk <- as.matrix(i)
    if(!(missing(j) && missing(k)))
      stop("If i is a matrix, j and k should not be given", call.=FALSE)
    if(ncol(ijk) != 3)
      stop("If i is a matrix, it should have 3 columns", call.=FALSE)
    if(!all(inside3Darray(dimx, i)))
      stop("Some indices lie outside array limits", call.=FALSE)
    if(nrow(ijk) == 0)
      return(x) # no items to replace
    ## assemble data frame
    xdata <- data.frame(i=x$i, j=x$j, k=x$k, x=x$x)
    ## match xdata into ijk (not necessarily the first match in original order)
    m <- matchIntegerDataFrames(xdata[,1:3,drop=FALSE], ijk)
    ## ------- OLDER VERSION: --------
    ## convert triples of integers to character codes
    ## icode <- apply(ijk, 1, paste, collapse=",") << is too slow >>
    ## icode <- paste(ijk[,1], ijk[,2], ijk[,3], sep=",")
    ## xcode <- paste(x$i, x$j, x$k, sep=",")
    ##  m <- match(xcode, icode)
    ## -------------------------------
    ## remove any matches, retaining only data that do not match 'i'
    xdata <- xdata[is.na(m), , drop=FALSE]   # sic
    ## ensure replacement value is vector-like
    value <- as.vector(value)
    nv <- length(value)
    if(nv != nrow(i) && nv != 1)
      stop(paste("Number of items to replace", paren(nrow(i)),
                 "does not match number of items given", paren(nv)),
           call.=FALSE)
    vdata <- data.frame(i=ijk[,1], j=ijk[,2], k=ijk[,3], x=value)
    ## combine
    ydata <- rbind(xdata, vdata)
    y <- with(ydata, sparse3Darray(i=i,j=j,k=k,x=x,
                                   dims=dimx, dimnames=dn, strict=TRUE))
    return(y)
  }
  I <- grokIndexVector(if(missing(i)) NULL else i, dimx[1], dn[[1]])
  J <- grokIndexVector(if(missing(j)) NULL else j, dimx[2], dn[[2]])
  K <- grokIndexVector(if(missing(k)) NULL else k, dimx[3], dn[[3]])
  IJK <- list(I,J,K)
  if(!all(sapply(lapply(IJK, getElement, name="full"), is.null))) {
    warning("indices exceed array bounds; extending the array dimensions",
            call.=FALSE)
    fullindices   <- lapply(IJK, fullIndexSequence)
    ## strictindices <- lapply(IJK, strictIndexSequence)
    dnew <- pmax(dimx, sapply(fullindices, max))
    result <- array(data=RelevantZero(x$x), dim=dnew)
    result[cbind(x$i, x$j, x$k)] <- x$x
    result[fullindices[[1]], fullindices[[2]], fullindices[[3]]] <- value
    result <- as.sparse3Darray(result)
    return(result)
  }
  IJK <- lapply(IJK, getElement, name="strict")
  if(all(sapply(IJK, getElement, name="nind") == 0)) {
    # no elements are indexed
    return(x)
  }
  I <- IJK[[1]]
  J <- IJK[[2]]
  K <- IJK[[3]]
  #' extract current array entries
  xdata <- data.frame(i=x$i, j=x$j, k=x$k, x=x$x)
  #' identify data volume that will be overwritten
  inI <- I$lo
  inJ <- J$lo
  inK <- K$lo
  #' remove data that will be overwritten
  retain <- !with(xdata, inI[i] & inJ[j] & inK[k])
  xdata <- xdata[retain,,drop=FALSE]
  #' expected dimensions of 'value' implied by indices
  dimVshould <- sapply(IJK, getElement, name="nind")
  dimV <- dim(value)
  if(length(dimV) == 3) {
    #' both source and destination are 3D
    if(all(dimVshould == dimV)) {
      #' replace 3D block by 3D block of same dimensions
      value <- as.sparse3Darray(value)
      vdata <- data.frame(i=value$i, j=value$j, k=value$k, x=value$x)
      # determine positions of replacement data in original array
      vdata <- transform(vdata,
                         i=replacementIndex(i, I),
                         j=replacementIndex(j, J),
                         k=replacementIndex(k, K))
    } else
      stop(paste("Replacement value has wrong dimensions:",
                 paste(dimV, collapse="x"),
                 "instead of",
                 paste(dimVshould, collapse="x")),
           call.=FALSE)
  } else if(is.null(dimV)) {
    #' replacement value is a vector or sparseVector
    value <- as(value, "sparseVector")
    iv <- value@i
    xv <- value@x
    nv <- value@length
    collapsing <- (dimVshould == 1)
    realdim <- sum(!collapsing)
    if(nv == 1) {
      #' replacement value is a constant
      value <- as.vector(value[1])
      if(identical(value, RelevantZero(x$x))) {
        #' assignment causes relevant entries to be set to zero;
        #' these entries have already been deleted from 'xdata';
        #' nothing to add
        vdata <- data.frame(i=integer(0), j=integer(0), k=integer(0),
                            x=x$x[integer(0)])
      } else {
        #' replicate the constant
        vdata <- expand.grid(i=I$i, j=J$i, k=K$i, x=as.vector(value[1]))
      }
    } else if(realdim == 0) {
        stop(paste("Replacement value has too many entries:",
                   nv, "instead of 1"),
             call.=FALSE)
    } else if(realdim == 1) {
      theindex <- which(!collapsing)
      # target slice is one-dimensional
      if(nv != dimVshould[theindex]) 
        stop(paste("Replacement value has wrong number of entries:",
                   nv, "instead of", dimVshould[theindex]),
             call.=FALSE)
      newpos <- replacementIndex(iv, IJK[[theindex]])
      vdata <- switch(theindex,
                      data.frame(i=newpos, j=J$i,    k=K$i,     x=xv),
                      data.frame(i=I$i,    j=newpos, k=K$i,     x=xv),
                      data.frame(i=I$i,    j=J$i,    k=newpos,  x=xv))
    } else {
      # target slice is two-dimensional
      sdim <- dimVshould[!collapsing]
      sd1 <- sdim[1]
      sd2 <- sdim[2]
      if(nv != sd1)
        stop(paste("Length of replacement vector", paren(nv),
                   "does not match dimensions of array subset",
                   paren(paste(dimVshould, collapse="x"))),
             call.=FALSE)
      firstindex <- which(!collapsing)[1]
      secondindex <- which(!collapsing)[2]
      pos1 <- replacementIndex(iv, IJK[[firstindex]])
      pos2 <- replacementIndex(seq_len(sd2), IJK[[secondindex]])
      xv   <- rep(xv, sd2)
      pos2 <- rep(pos2, each=length(pos1))
      pos1 <- rep(pos1, sd2)
      pos3 <- if(length(pos1)) IJK[[which(collapsing)]]$i else integer(0)
      vdata <- data.frame(i=pos3, j=pos3, k=pos3, x=xv)
      vdata[,firstindex] <- pos1
      vdata[,secondindex] <- pos2
    }
  } else if(identical(dimVshould[dimVshould > 1],  dimV[dimV > 1])) {
    #' lower dimensional sets of the same dimension
    value <- value[drop=TRUE]
    dimV <- dim(value)
    dropping <- (dimVshould == 1)
    if(length(dimV) == 2) {
      value <- as(value, "TsparseMatrix")
      iv <- value@i + 1L
      jv <- value@j + 1L
      xv <- value@x
      firstindex <- which(!dropping)[1]
      secondindex <- which(!dropping)[2]
      pos1 <- replacementIndex(iv, IJK[[firstindex]])
      pos2 <- replacementIndex(jv, IJK[[secondindex]])
      pos3 <- if(length(pos1)) IJK[[which(dropping)]]$i else integer(0)
      vdata <- data.frame(i=pos3, j=pos3, k=pos3, x=xv)
      vdata[,firstindex] <- pos1
      vdata[,secondindex] <- pos2
    } else {
      value <- as(value, "sparseVector")
      iv <- value@i
      xv <- value@x
      vdata <- data.frame(i=if(dropping[1]) I$i else replacementIndex(iv, I),
                          j=if(dropping[2]) J$i else replacementIndex(iv, J),
                          k=if(dropping[3]) K$i else replacementIndex(iv, K),
                          x=xv)
    }
  } else
    stop(paste("Replacement value has wrong dimensions:",
               paste(dimV, collapse="x"),
               "instead of",
               paste(dimVshould, collapse="x")),
         call.=FALSE)
    
  ## combine
  if(nrow(vdata) > 0)
    xdata <- rbind(xdata, vdata)
  y <- with(xdata, sparse3Darray(i=i,j=j,k=k,x=x,
                                 dims=dimx, dimnames=dn, strict=TRUE))
  return(y)
}

bind.sparse3Darray <- function(A,B,along) {
  A <- as.sparse3Darray(A)
  B <- as.sparse3Darray(B)
  check.1.integer(along)
  stopifnot(along %in% 1:3)
  dimA <- dim(A)
  dimB <- dim(B)
  if(!all(dimA[-along] == dimB[-along]))
    stop("dimensions of A and B do not match")
  dimC <- dimA
  dimC[along] <- dimA[along] + dimB[along]
  # extract data
  Adf <- SparseEntries(A)
  Bdf <- SparseEntries(B)
  # realign 'B' coordinate
  Bdf[,along] <- Bdf[,along] + dimA[along]
  # combine
  C <- EntriesToSparse(rbind(Adf, Bdf), dimC)
  # add dimnames
  dnA <- dimnames(A)
  dnB <- dimnames(B)
  if(!is.null(dnA) || !is.null(dnB)) {
    if(length(dnA) != 3) dnA <- rep(list(NULL), 3)
    if(length(dnB) != 3) dnB <- rep(list(NULL), 3)
    dnC <- dnA
    dnC[[along]] <- c(dnA[[along]] %orifnull% rep("", dimA[along]),
                      dnB[[along]] %orifnull% rep("", dimB[along]))
    dimnames(C) <- dnC
  }
  return(C)
}


anyNA.sparse3Darray <- function(x, recursive=FALSE) {
  anyNA(x$x)
}

unionOfSparseIndices <- function(A, B) {
  #' A, B are data frames of indices i, j, k
  ijk <- unique(rbind(A, B))
  colnames(ijk) <- c("i", "j", "k")
  return(ijk)
}
  
Ops.sparse3Darray <- function(e1,e2=NULL){
  if(nargs() == 1L) {
    switch(.Generic,
           "!" = {
             result <- do.call(.Generic, list(as.array(e1)))
           },
           "-" = ,
           "+" = {
             result <- e1
             result$x <- do.call(.Generic, list(e1$x))
           },
           stop(paste("Unary", sQuote(.Generic),
                      "is undefined for sparse 3D arrays."), call.=FALSE))
    return(result)
  }
  # binary operation
  # Decide whether full or sparse
  elist <- list(e1, e2)
  isfull <- sapply(elist, inherits, what=c("matrix", "array"))
  if(any(isfull) &&
     any(sapply(lapply(elist[isfull], dim), prod) > 1)) {
    # full array
    n1 <- length(dim(e1))
    n2 <- length(dim(e2))
    e1 <- if(n1 == 3) as.array(e1) else
          if(n1 == 2) as.matrix(e1) else as.vector(as.matrix(as.array(e1)))
    e2 <- if(n2 == 3) as.array(e2) else
          if(n2 == 2) as.matrix(e2) else as.vector(as.matrix(as.array(e2)))
    result <- do.call(.Generic, list(e1, e2))
    return(result)
  }
  # sparse result (usually)
  e1 <- as.sparse3Darray(e1)
  e2 <- as.sparse3Darray(e2)
  dim1 <- dim(e1)
  dim2 <- dim(e2)
  mode1 <- typeof(e1$x)
  mode2 <- typeof(e2$x)
  zero1 <- vector(mode=mode1, length=1L)
  zero2 <- vector(mode=mode2, length=1L)
  
  if(prod(dim1) == 1) {
    ## e1 is constant
    e1 <- as.vector(as.array(e1))
    z12 <- do.call(.Generic, list(e1, zero2))
    if(!isRelevantZero(z12)) {
      # full matrix/array will be generated
      result <- do.call(.Generic, list(e1, as.array(e2)[drop=TRUE]))
    } else {
      # sparse 
      result <- e2
      result$x <- do.call(.Generic, list(e1, e2$x))
    }
    return(result)
  }

  if(prod(dim2) == 1) {
    ## e2 is constant
    e2 <- as.vector(as.array(e2))
    z12 <- do.call(.Generic, list(zero1, e2))
    if(!isRelevantZero(z12)) {
      # full matrix/array will be generated
      result <- do.call(.Generic, list(as.array(e1)[drop=TRUE], e2))
    } else {
      # sparse 
      result <- e1
      result$x <- do.call(.Generic, list(e1$x, e2))
    }
    return(result)
  }
  
  z12 <- do.call(.Generic, list(zero1, zero2))
  if(!isRelevantZero(z12)) {
    #' Result is an array
    e1 <- as.array(e1)
    e2 <- as.array(e2)
    result <- do.call(.Generic, list(e1, e2))
    return(result)
  }

  # Result is sparse
  if(identical(dim1, dim2)) {
    #' extents are identical
    ijk1 <- SparseIndices(e1)
    ijk2 <- SparseIndices(e2)
    if(identical(ijk1, ijk2)) {
      #' patterns of nonzero entries are identical
      ijk <- ijk1
      values <- do.call(.Generic, list(e1$x, e2$x))
    } else {			   
      #' different patterns of nonzero entries
      ijk <- unionOfSparseIndices(ijk1, ijk2)
      values <- as.vector(do.call(.Generic, list(e1[ijk], e2[ijk])))
    }			      
    dn <- dimnames(e1) %orifnull% dimnames(e2)
    result <- sparse3Darray(i=ijk$i, j=ijk$j, k=ijk$k, x=values,
                              dims=dim1, dimnames=dn, strict=TRUE)
    return(result)
  }

  drop1 <- (dim1 == 1)
  drop2 <- (dim2 == 1)
  if(!any(drop1 & !drop2) && identical(dim1[!drop2], dim2[!drop2])) {
    #' dim2 is a slice of dim1
    ijk1 <- data.frame(i=e1$i, j=e1$j, k=e1$k)
    ijk2 <- data.frame(i=e2$i, j=e2$j, k=e2$k)
    expanding <- which(drop2 & !drop1)
    if(length(expanding) == 1) {
      n <- dim1[expanding]
      m <- nrow(ijk2)
      ijk2 <- as.data.frame(lapply(ijk2, rep, times=n))
      ijk2[,expanding] <- rep(seq_len(n), each=m)
      ijk <- unionOfSparseIndices(ijk1, ijk2)
      ijkdrop <- ijk
      if(nrow(ijkdrop) > 0) ijkdrop[,expanding] <- 1
      xout <- do.call(.Generic, list(e1[ijk], e2[ijkdrop]))
      result <- sparse3Darray(i=ijk[,1L], j=ijk[,2L], k=ijk[,3L],
                              x=as.vector(xout),
                              dims=dim1, dimnames=dimnames(e1), strict=TRUE)
      return(result)
    }
  }

  if(!any(drop2 & !drop1) && identical(dim2[!drop1], dim1[!drop1])) {
    #' dim1 is a slice of dim2
    ijk1 <- data.frame(i=e1$i, j=e1$j, k=e1$k)
    ijk2 <- data.frame(i=e2$i, j=e2$j, k=e2$k)
    expanding <- which(drop1 & !drop2)
    if(length(expanding) == 1) {
      n <- dim2[expanding]
      m <- nrow(ijk1)
      ijk1 <- as.data.frame(lapply(ijk1, rep, times=n))
      ijk1[,expanding] <- rep(seq_len(n), each=m)
      ijk <- unionOfSparseIndices(ijk1, ijk2)
      ijkdrop <- ijk
      if(nrow(ijkdrop) > 0) ijkdrop[,expanding] <- 1L
      xout <- do.call(.Generic, list(e1[ijkdrop], e2[ijk]))
      result <- sparse3Darray(i=ijk[,1L], j=ijk[,2L], k=ijk[,3L],
                              x=as.vector(xout),
                              dims=dim2, dimnames=dimnames(e2), strict=TRUE)
      return(result)
    }
  }

  if(all(drop1[-1]) && dim1[1L] == dim2[1L]) {
    #' e1 is a (sparse) vector matching the first extent of e2
    if(.Generic %in% c("*", "&")) {
      # result is sparse
      ijk <- data.frame(i=e2$i, j=e2$j, k=e2$k)
      ones <- rep(1L, nrow(ijk))
      i11 <- data.frame(i=e2$i, j=ones, k=ones)
      xout <- do.call(.Generic, list(e1[i11], e2[ijk]))
      result <- sparse3Darray(i=ijk[,1L], j=ijk[,2L], k=ijk[,3L],
                              x=as.vector(xout),
                              dims=dim2, dimnames=dimnames(e2), strict=TRUE)
    } else {
      # result is full array
      e1 <- as.array(e1)[,,,drop=TRUE]
      e2 <- as.array(e2)
      result <- do.call(.Generic, list(e1, e2))
    }
    return(result)
  }
  
  stop(paste("Non-conformable arrays:",
             paste(dim1, collapse="x"), "and", paste(dim2, collapse="x")),
       call.=FALSE)
}

Math.sparse3Darray <- function(x, ...){
  z <- RelevantZero(x$x)
  fz <- do.call(.Generic, list(z))
  if(!isRelevantZero(fz)) {
    # result is a full array
    result <- do.call(.Generic, list(as.array(x), ...))
    return(result)
  }
  x$x <- do.call(.Generic, list(x$x))
  return(x)
}

Complex.sparse3Darray <- function(z) {
  oo <- RelevantZero(z$x)
  foo <- do.call(.Generic, list(z=oo))
  if(!isRelevantZero(foo)) {
    # result is a full array
    result <- do.call(.Generic, list(z=as.array(z)))
    return(result)
  }
  z$x <- do.call(.Generic, list(z=z$x))
  return(z)
}

Summary.sparse3Darray <- function(..., na.rm=FALSE) {
  argh <- list(...)
  is3D <- sapply(argh, inherits, what="sparse3Darray")
  if(any(is3D)) {
    xvalues   <- lapply(argh[is3D], getElement, name="x")
    fullsizes <- sapply(lapply(argh[is3D], dim), prod)
    argh[is3D] <- xvalues
    #' zero entry should be appended if and only if there are any empty cells
    zeroes <- lapply(xvalues, RelevantZero)
    zeroes <- zeroes[lengths(xvalues) < fullsizes]
    argh <- append(argh, zeroes)
  }
  rslt <- do.call(.Generic, append(argh, list(na.rm=na.rm)))
  return(rslt)
}


SparseIndices <- function(x) {
  #' extract indices of entries of sparse vector/matrix/array
  nd <- length(dim(x))
  if(nd > 3)
    stop("Arrays of more than 3 dimensions are not supported", call.=FALSE)
  if(nd == 0 || nd == 1) {
    x <- as(x, "sparseVector")
    df <- data.frame(i=x@i)
  } else if(nd == 2) {
    x <- as(x, "TsparseMatrix")
    df <- data.frame(i=x@i + 1L, j=x@j + 1L)
  } else if(nd == 3) {
    x <- as.sparse3Darray(x)
    df <- data.frame(i=x$i, j=x$j, k=x$k)
  }
  return(df)
}

SparseEntries <- function(x) {
  #' extract entries of sparse vector/matrix/array
  nd <- length(dim(x))
  if(nd > 3)
    stop("Arrays of more than 3 dimensions are not supported", call.=FALSE)
  if(nd == 0 || nd == 1) {
    x <- as(x, "sparseVector")
    df <- data.frame(i=x@i, x=x@x)
  } else if(nd == 2) {
    x <- as(x, "TsparseMatrix")
    df <- data.frame(i=x@i + 1L, j=x@j + 1L, x=x@x)
  } else if(nd == 3) {
    x <- as.sparse3Darray(x)
    df <- data.frame(i=x$i, j=x$j, k=x$k, x=x$x)
  }
  return(df)
}

EntriesToSparse <- function(df, dims) {
  #' convert data frame of indices and values
  #' to sparse vector/matrix/array
  nd <- length(dims)
  if(nd == 0)
    return(with(df, as(sum(x), typeof(x))))
  sn <- seq_len(nd)
  colnames(df)[sn] <- c("i","j","k")[sn]
  if(nd == 1) {
    #' sparse vector: duplicate entries not allowed
    df <- df[with(df, order(i)), , drop=FALSE]
    dup <- c(FALSE, with(df, diff(i) == 0))
    if(any(dup)) {
      #' accumulate values at the same array location
      first <- !dup
      newi <- cumsum(first)
      newx <- as(tapply(df$x, newi, sum), typeof(df$x))
      df <- data.frame(i=df$i[first], x=newx)
    }
    result <- with(df, sparseVector(i=i, x=x, length=dims))
  } else if(nd == 2) {
    result <- with(df, sparseMatrix(i=i, j=j, x=x, dims=dims))
  } else if(nd == 3) {
    result <- with(df, sparse3Darray(i=i, j=j, k=k, x=x, dims=dims))
  }
  return(result)
}

evalSparse3Dentrywise <- function(expr, envir) {
  ## DANGER: this assumes all sparse arrays in the expression
  ##         have the same pattern of nonzero elements!
  e <- as.expression(substitute(expr))
  ## get names of all variables in the expression
  varnames <- all.vars(e)
  allnames <- all.names(e, unique=TRUE)
  funnames <- allnames[!(allnames %in% varnames)]
  if(length(varnames) == 0)
    stop("No variables in this expression")
  ## get the values of the variables
  if(missing(envir)) {
    envir <- parent.frame() # WAS: sys.parent()
  } else if(is.list(envir)) {
    envir <- list2env(envir, parent=parent.frame())
  }
  vars <- mget(varnames, envir=envir, inherits=TRUE, ifnotfound=list(NULL))
  funs <- mget(funnames, envir=envir, inherits=TRUE, ifnotfound=list(NULL))
  ## find out which variables are sparse3Darray
  isSpud <- sapply(vars, inherits, what="sparse3Darray")
  if(!any(isSpud))
    stop("No sparse 3D arrays in this expression")
  spuds <- vars[isSpud]
  template <- spuds[[1L]]
  ## replace each array by its entries, and evaluate
  spudvalues <- lapply(spuds, getElement, name="x")
  ## minimal safety check
  if(length(unique(lengths(spudvalues))) > 1)
    stop("Different numbers of sparse entries", call.=FALSE)
  vars[isSpud] <- spudvalues
  v <- eval(e, append(vars, funs))
  ## reshape as 3D array
  result <- sparse3Darray(x=v,
  	                  i=template$i,
  	                  j=template$j,
  	                  k=template$k,
			  dims=dim(template),
			  dimnames=dimnames(template))
  return(result)
}