R/sparse3Darray.R

Defines functions evalSparse3Dentrywise EntriesToSparse SparseEntries SparseIndices Summary.sparse3Darray Complex.sparse3Darray Math.sparse3Darray Ops.sparse3Darray unionOfSparseIndices anyNA.sparse3Darray bind.sparse3Darray rbindCompatibleDataFrames as.array.sparse3Darray aperm.sparse3Darray print.sparse3Darray dimnames.sparse3Darray dim.sparse3Darray as.sparse3Darray sparse3Darray

Documented in anyNA.sparse3Darray aperm.sparse3Darray as.array.sparse3Darray as.sparse3Darray bind.sparse3Darray Complex.sparse3Darray dimnames.sparse3Darray dim.sparse3Darray EntriesToSparse evalSparse3Dentrywise Math.sparse3Darray Ops.sparse3Darray print.sparse3Darray rbindCompatibleDataFrames sparse3Darray SparseEntries SparseIndices Summary.sparse3Darray unionOfSparseIndices

#'
#' 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.44 $  $Date: 2021/03/04 07:42:18 $
#'

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=newi[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)
}

Try the spatstat.sparse package in your browser

Any scripts or data that you put into this service are public.

spatstat.sparse documentation built on March 16, 2021, 9:08 a.m.