R/GridData.R

In IRSN/smint: Smooth Multivariate Interpolation for Gridded and Scattered Data

##*****************************************************************************
#' @title S4 class for Grid
#' 
#' @description This class is intended handle Grid object,
#'     with support for useful methods like levels, ...
#'
#' @slot dim number of dimensions
#' 
#' @slot dimNames name of dimensions
#' 
#' @slot levels cardinal of discrete values
#' 
#' @slot index
#' 
#' @author Yves Deville
#'
#' @name Grid-class
#' @rdname Grid-class
#' 
#' @export
setClass("Grid",    	
         representation(
            dim = "integer",   
            dimNames = "character",
            levels = "list",
            index = "array"),
         validity = function(object) {
            if (length(object@dimNames) != object@dim) {
                return("incorrect length for 'dimNames'")
            }
            if (length(object@levels) != object@dim) {
                return("incorrect length for 'levels'")
            }
            if (!all(dim(object@index) == lapply(object@levels, length))) {
                return("'index' dimension conflicts with 'levels'")
            }
            if (any(sapply(object@levels, function(x) any(duplicated(x)) ))) {
                return("levels must be distinct'")
            }
       
            TRUE
         }
)

##*****************************************************************************
#' Create a Grid object.
#' 
#' @title Create a new \code{Grid} object
#'
#' @param levels A list with the levels of the variables. The length
#' of the list gives the spatial dimension. The \eqn{i}-th element
#' must be a numeric vector containing the distinct values of the
#' corresponding variable.
#'
#' @param nlevels Integer vector giving the number of levels by
#' dimension. This formal argument is used only when \code{levels} is
#' missing. In this case, the nodes for dimension \eqn{i} will be
#' regularly spaced between \code{0.0} and \code{1.0}.
#'
#' @param dimNames If the list \code{levels} does not have names,
#' then the character vector \code{dimNames} will be used.
#'
#' @param index An array with integer values and dimension
#' corresponding to the spatial dimensions. Each value represent the
#' number of the corresponding node in a "flat" representation
#' (data.frame or matrix), and will be used in tasks such as
#' interpolation.
#' 
#' @export
#'
#' @return An object with S4 class \code{"Grid"}. If needed, this
#' object can be coerced to a data frame or to a matrix by using the
#' methods \code{as.data.frame} or \code{as.matrix}.
#'
#' @author Yves Deville
#'
#' @note If \code{index} is not provided, the vectors in the
#' \code{levels} list can be unsorted, and they will be
#' sorted. However, when a non-default \code{index} is provided, the
#' levels must be sorted. This rule is intended to reduce the risk of
#' error.
#'
#' @examples
#' myGD1 <- Grid(levels = list("X" = c(0, 1), "Y" = c(0.0, 0.5, 1.0)))
#' ## the same with less effort
#' myGD2 <- Grid(nlevels = c("X" = 2, "Y" = 3))
#' nlevels(myGD2)
#' levels(myGD2)
#' 
Grid <- function(levels,
                     nlevels = NULL,
                     dimNames = NULL,
                     index = NULL) {
  
  if (missing(levels)) {
    if ( (length(nlevels) == 1) && (!is.null(dimNames)) ) {
      nlevels <- rep(nlevels, length.out = length(dimnames))
      names(nlevels) <- dimNames
    }
    dimNames <- names(nlevels)
    levels <- list()
    for (i in 1L:length(nlevels)) {
      levels[[i]] <- seq(from = 0.0, to = 1.0, length.out = nlevels[i])
    }
    names(levels) <- names(nlevels)
  } else {
    if (!is.null(nlevels)){
      warning("'levels' is given so 'nlevels' is ignored")
    }
    if (!is.list(levels)) stop("'levels' must be a list")
    levels <- lapply(levels, as.numeric)
    if (!is.null(index) && any(sapply(levels, is.unsorted))) {
      stop("'index' can be provided only when levels are sorted")
    }
    levels <- lapply(levels, sort)
  }
  
  dim <- length(levels)
  nLevels <- sapply(levels, length)
  levChar <- lapply(levels, round_levels)
  levChar <- lapply(levChar, as.character)
      
  if (is.null(dimNames)) {
    if (!is.null(nm <- names(nLevels))) dimNames <- nm
    else {
      dimNames <- paste("X", 1L:dim, sep = "") 
      names(levels) <- dimNames
    }
  }
  names(levChar) <- dimNames
  
  n <- prod(nLevels)
  names(nLevels) <- NULL ## to avoid having names in the 'dim' attribute
  if (is.null(index)) index <- array(1L:n, dim = nLevels, dimnames = levChar)

  myGD <- new("Grid",
              dim = as.integer(dim),
              dimNames = as.character(dimNames),
              levels = levels,
              index = index)
  myGD
  
}

##*****************************************************************************
#' Random drawing of a \code{Grid} object
#'
#' The nodes are drawn by first choosing their number for each
#' spatial dimension using a Poisson distribution, and then by
#' drawing a sample of the uniform distribution with the suitable
#' size.
#' @title Random drawing of a \code{Grid} object
#'
#' @param dim The spatial dimension.
#'
#' @param nxMin Minimum number of nodes. Integer vector of length 1
#' or \code{dim}.example(Grid).
#'
#' @param nxE Expected number of nodes. Numeric vector of length 1
#' or \code{dim}.
#'
#' @param nlevels Integer vector giving the number of levels for each
#' spatial dimension.
#'
#' @param dimNames Names of the spatial dimensions.
#'
#' @export
#' 
#' @importFrom stats runif
#' @importFrom stats rpois
#' 
#' @return An object with S4 class \code{"Grid"}.
#'
#' @author Yves Deville
#'
#' @note This function is \emph{not} not for drawing random design
#' points from within the object.
#'
#' @seealso \code{\link{Grid-class}}.
#'
#' @examples
#' set.seed(1234)
#' ## specify dim: number of levels are random!
#' GD <- randGrid(dim = 4)
#' plot(GD)
#' ## specify number of levels (hence dim)
#' GD2 <- randGrid(nlevels = c(2, 3, 4))
#' plot(GD2)
#' 
randGrid <- function(dim = 1L + rpois(1L, lambda = 3),
                      nxMin = 1L, nxE = 4L,
                      nlevels, dimNames) {
   
   if (missing(nlevels)) {
     if (missing(dim))  {
       dim <- 1L + rpois(1L, lambda = 3)
     } 
     nlevels <- nxMin + rpois(dim, lambda = nxE - nxMin)
   } else {
     dim <- length(nlevels)
   }
   
   if (missing(dimNames)) {
      if (!is.null(nm <- names(nlevels))) dimNames <- nm
      else dimNames <- paste("X", 1L:dim, sep = "") 
   }
   
   levels <- list()
   for (i in 1L:dim) {
      levels[[i]] <- sort(runif(nlevels[i]))
   }
   names(levels) <- dimNames
   n <- prod(nlevels)
   
   levChar <- lapply(levels, round_levels)
   levChar <- lapply(levChar, as.character)
   
   index <- array(1L:n, dim = nlevels, dimnames = levChar)
   
   myGD <- new("Grid",
               dim = as.integer(dim),
               dimNames = as.character(dimNames),
               levels = levels,
               index = index)
   myGD
   
}

##*****************************************************************************
## Basic extraction methods
##*****************************************************************************

#' @title Get dimension of given Grid
#' @author Yves Deville
#' @param x Grid object
#' @return Number of dimensions/columns of the Grid object.
#' @method dim Grid
#' @exportMethod dim
#' @aliases dim,Grid-method
setMethod("dim", 
          signature = signature(x = "Grid"),
          definition = function(x){ x@dim })

#' @title Get dimension names of given Grid
#' @author Yves Deville
#' @param x Grid object
#' @return Names of dimensions/columns of the Grid object.
#' @method dimnames Grid
#' @exportMethod dimnames
#' @aliases dimnames,Grid-method
setMethod("dimnames", 
          signature = signature(x = "Grid"),
          definition = function(x){ x@dimNames })

#' @title Set dimension names of Grid
#' @author Yves Deville
#' @param x Grid object
#' @param value New dimension names of Grid
#' @return Modified Grid object
#' @method dimnames Grid
#' @exportMethod dimnames<-
#' @aliases dimnames<-,Grid-method
setMethod("dimnames<-", 
          signature = signature(x = "Grid"),
          definition = function(x, value){
              if (length(value) != x@dim) {
                  stop(sprintf("'value' has length %d but must have length %d",
                               length(value), x@dim))
              }
              x@dimNames[] <- names(x@levels) <- value
              L <- dimnames(x@index)
              names(L) <- value
              dimnames(x@index) <- L
              x
          }
          )

#' @title Get levels of Grid
#' @author Yves Deville
#' @param x Grid object
#' @return Levels of the Grid object.
#' @method levels Grid
#' @exportMethod levels
#' @aliases levels,Grid-method
setMethod("levels", 
          signature = signature(x = "Grid"),
          definition = function(x){ x@levels })

#' @title Get number of levels of Grid
#' @author Yves Deville
#' @param x Grid object
#' @return Number of levels of the Grid object.
#' @method nlevels Grid
#' @exportMethod nlevels
#' @aliases nlevels,Grid-method
setMethod("nlevels", 
          signature = signature(x = "Grid"),
          definition = function(x){
            sapply(x@levels, length)
          })

## XXX problem : define the 'index' function
## setMethod("index", 
##           signature = signature(object = "Grid"),
##           definition = function(object){ object@index })

##*****************************************************************************
## Coercion S4 methods
##*****************************************************************************

#' @title Coerce Grid in data.frame
#' @author Yves Deville
#' @param x Grid object
#' @param jitter if TRUE, add a small white noise to data. Defaults to FALSE.
#' @param ... Not used.
#' @return Coerced data.frame object
#' @method as.data.frame Grid
#' @exportMethod as.data.frame
#' @importFrom stats rnorm
#' @aliases as.data.frame,Grid-method
setMethod("as.data.frame", 
          signature = signature(x = "Grid"),
          definition = function(x, jitter = FALSE, ...){
             n <- prod(dim(x@index))
             lev <- x@levels
             names(lev) <- x@dimNames
             df <- expand.grid(lev)
             dni <- Matrix::invPerm(x@index)
             df <- df[dni, ]
             rownames(df) <- as.character(1L:n)
             if (jitter) {
                 dx <- lapply(lev, function(x) min(diff(x)) / 30)
                 for (j in 1:ncol(df)) {
                     if (is.finite(dx[[j]])) {
                         df[ , j] <-  df[ , j] + rnorm(n, mean = 0, sd = dx[[j]])
                     }
                 }
             }
             df
         }
)

#' @title Coerce Grid in matrix
#' @author Yves Deville
#' @param x Grid object
#' @param jitter if TRUE, add a small white noise to data. Defaults to FALSE.
#' @param ... Not used.
#' @return Coerced matrix object
#' @method as.matrix Grid
#' @exportMethod as.matrix
#' @aliases as.matrix,Grid-method
setMethod("as.matrix", 
          signature = signature(x = "Grid"),
          definition = function(x, jitter = FALSE, ...){
              as.matrix(as.data.frame(x, jitter = jitter))
          })

##*****************************************************************************
## Classical methods
##*****************************************************************************          
#' @title Get length of Grid
#' @author Yves Deville
#' @param x Grid object
#' @return Length of the Grid object.
#' @method length Grid
#' @exportMethod length
#' @aliases length,Grid-method
setMethod("length", 
          signature = signature(x = "Grid"),
          definition = function(x){
              length(x@index)
          })

#' @title Plot Grid as pairs
#' @author Yves Deville
#' @param x Grid object
#' @param y optional data to add in plot alongside x (should have same length)
#' @param jitter if TRUE, add a small white noise to data. Defaults to FALSE.
#' @param ... Optional arguments passed to 'pairs()'
#' @return 'pairs' object
#' @method plot Grid
#' @importFrom graphics pairs
#' @exportMethod plot
#' @aliases plot,Grid-method
setMethod("plot", 
          signature = signature(x = "Grid"),
          definition = function(x, y, jitter = FALSE, ...){
              df <- as.data.frame(x, jitter = jitter)
              if (!missing(y)) {
                  df <- cbind(df, y)
              }
              pairs(df, ...)
          })

#' @title Print Grid object informations
#' @author Yves Deville
#' @param object Grid object
#' @return Nothing.
#' @method show Grid
#' @exportMethod show
#' @aliases show,Grid-method
setMethod("show", 
          signature = signature(object = "Grid"),
          definition = function(object){
             cat("Grid Data object\n")
             cat(sprintf("  o dimension : %d\n", dim(object)))
             nx <- sapply(levels(object), length)
             nn <- sprintf("(%s)", dimnames(object))
             nn <- paste(nx, nn)
             nn <- paste(nn, collapse = ", ")
             cat(sprintf("  o dim names : %s\n",
                         paste(dimnames(object), collapse = ", ")))
             cat(sprintf("  o number of nodes : %s\n", nn))
             cat(sprintf("  o total number of nodes : %d\n", prod(nx)))
          })

show <- function(object, ...) {
   UseMethod("show")
}

##*****************************************************************************
#' @title Permutation method
#' @description The generalised transpose changes the order of the factors 
#'    keeping the order of the observations. The transformation to a data frame is
#'    not efficient in term of computations, but is easily understood.
#' @author Yves Deville
#' @param a Grid object
#' @param perm new order indexes
#' @param ... Not used.
#' @return Permuted Grid.
#' @method aperm Grid
#' @exportMethod aperm
#' @aliases aperm,Grid-method
##*****************************************************************************          
## setMethod("aperm",
##           signature = signature(a = "Grid"),
##           definition = function(a, perm, ...){
##              indexNew <- aperm(a@index, perm = perm, ...)
##              dimNamesNew <- a@dimNames[perm]
##              levelsNew <- list()
##              for (j in 1L:dim(a)) {
##                  nm <- dimNamesNew[j]
##                  levelsNew[[nm]] <- a@levels[[perm[j]]] 
##              }
##              new("Grid",
##                  dim = dim(a),
##                  dimNames = dimNamesNew,
##                  levels = levelsNew,
##                  index = indexNew)      
##          })
setMethod("aperm",
          signature = signature(a = "Grid"),
          definition = function(a, perm, ...){
              if (!setequal(sort(as.integer(perm)), 1L:dim(a)))
                  stop("'perm' must be a valid permutation of ",
                       dim(a), " items")
              des <- as.data.frame(a)
              as.Grid(des[ , perm])
          })

##*****************************************************************************
## Coercion S3 methods
##*****************************************************************************
#' Coercion to \code{Grid}.
#'
#' The dimensions of the \code{Grid} object are matched to the \eqn{d}
#' columns of the object given in \code{object}, in the same
#' order. Each column is coerced into a factor. Each combination
#' of the \eqn{d} different levels must be found exactly one in the
#' data frame or matrix, and the the row number of that combination is
#' stored in the \code{index} slot of the object.
#'
#' @title Coercion to \code{Grid}
#'
#' @param object An object to be coerced into a \code{Grid} object.
#'
#' @param ... Not used yet.
#'
#' @export
#' 
#' @return Object with S4 class \code{"Grid"}.
#' 
#' @note A numeric object is coerced into a \code{Grid} object
#' with dimension \eqn{d =1}.
#'
#' @examples
#' set.seed(1234)
#' GDnum <- as.Grid(runif(8))
#' GDlist <- as.Grid(list("X" = runif(8), "Y" = runif(4)))
#' df <- expand.grid(X = runif(6), Y = runif(4))
#' GDdf <- as.Grid(df)
#' GDmat <- as.Grid(as.matrix(df))
as.Grid <- function(object, ...) {
   UseMethod("as.Grid")
}

#' @title Coercion to \code{Grid}
#' @method as.Grid default
#' @param object An object to be coerced into a \code{Grid} object.
#' @param ... Not used yet.
#' @rdname as.Grid
#' @export
as.Grid.default <- function(object, ...) {
  as.Grid(as.data.frame(object), ...)
}

#' @title Coercion to \code{Grid}
#' @method as.Grid matrix
#' @param object An object to be coerced into a \code{Grid} object.
#' @param ... Not used yet.
#' @rdname as.Grid
#' @export
as.Grid.matrix <- function(object, ...) {
  as.Grid(as.data.frame(object), ...)
}

#' @title Coercion to \code{Grid}
#' @method as.Grid numeric
#' @param object An object to be coerced into a \code{Grid} object.
#' @param ... Not used yet.
#' @rdname as.Grid
#' @export
as.Grid.numeric <- function(object, ...) {
  as.Grid(as.data.frame(object), ...)
}

#' @title Coercion to \code{Grid}
#' @method as.Grid data.frame
#' @param object An object to be coerced into a \code{Grid} object.
#' @param ... Not used yet.
#' @rdname as.Grid
#' @export
as.Grid.data.frame <- function(object, ...) {
   
   X <- lapply(object, as.factor)
   levels <- lapply(X, levels)
   levels <- lapply(levels, as.numeric)
   Table <- tapply(object[ , 1L], X, length)
   if (any(is.na(Table)) || any(Table != 1)) {
       stop("each combination of the levels must appear",
            " exactly once in the data frame 'object'")
   }
   n <- prod(sapply(levels, length))
   index <- tapply(1L:n, X, function(x) x)
   ## L <- dimnames(index)
   ## names(L) <- names(object)
   ## dimnames(index) <- L
   myGD <- new("Grid",
               dim = ncol(object),
               dimNames = names(object),
               levels = levels,
               index = index)
   myGD
   
}

#' @title Coercion to \code{Grid}
#' @method as.Grid list
#' @param object An object to be coerced into a \code{Grid} object.
#' @param ... Not used yet.
#' @rdname as.Grid
#' @export
as.Grid.list <- function(object, ...) {
  dim <- length(object)
  if (is.null(names(object))) {
    names(object) <- paste("X", 1L:dim, sep = "")
  }
  object <- lapply(object, as.numeric)
  object <- lapply(object, sort)
  nx <- sapply(object, length)
  index <- array(1L:prod(nx), dim = nx)
  myGD <- new("Grid",
               dim = dim,
               dimNames = names(object),
               levels = object,
               index = index)
   myGD
}
  
##*****************************************************************************
#' Apply a function to a \code{Grid} object.
#'
#' The function is applied to each vector combination of the levels by first
#' using the \code{as.matrix} coercion.
#'
#' @title Apply a function to a \code{Grid} object.
#' 
#' @param object An object with S4 class \code{"Grid"}
#'
#' @param fun The function to be applied.
#'
#' @param ... Further arguments to be passed to \code{fun}.
#' 
#' @export
#'
#' @return A vector of values.
#'
#' @seealso \code{\link{array_Grid}} to reshape the result as an array.
#' 
#' @note The result values are given in the order specified by the
#' \code{Grid} object.
#' 
#' @examples
#' myGD <- Grid(levels = list("X" = c(0, 1),
#'                            "Y" = c(0.0, 0.5, 1.0),
#'                            "Z" = c(0.0, 0.2, 0.4, 1.0)))
#' ## gaussian function
#' apply_Grid(myGD, fun = function(x) exp(-x[1]^2 - 3 * x[2]^2 - 2 * x[3]^2))
#' 
#' 
apply_Grid <- function(object, fun, ...) {
    
  X <- as.matrix(object)
  fTest <- fun(X[1, ], ...)
  
  if (length(fTest) > 1L) {
    stop("'fun' returning a vector value is  not implemented yet")
  }

  Y <- apply(X = X, MARGIN = 1, FUN = fun, ...)
  
  ## is as.array is TRUE,
  if (FALSE) {
    dim(Y) <- nlevels(object)
    ## nms <- lapply(levels(object), as.character)
    dimnames(Y) <- dimnames(object@index)
  } 
  
  Y
  
}

##*****************************************************************************
#' Sample from/in a design object.
#'
#' From the object, \code{size} 'new' design points are sampled or
#' drawn. When \code{atSample} is \code{FALSE}, the range of each
#' column is computed or extracted from the object \code{x}, and then
#' independent drawings from the corresponding uniform distributions
#' are used.
#' 
#' @title Sample from/in a design object 
#'
#' @usage sampleIn(X, size = 1L, atSample = FALSE, ...)
#' 
#' @param X The object, typically a design matrix with class
#' \code{matrix} or an object with class \code{"Grid"}.
#'
#' @param size Number of sample points.
#'
#' @param atSample Logical. If \code{TRUE} the new sample points are
#' obtained using \code{sample}
#'
#' @param ... Other arguments to be passed to the \code{sample} function
#' when \code{atSample} is \code{TRUE}, typically \code{replace = TRUE}
#' is needed when the object \code{x} contains less than the target
#' number given in \code{size}.
#' 
#' @export
#'
#' @return A matrix with \code{size} rows. When \code{atSample} is
#' \code{TRUE}, the matrix has an attribute named \code{index} which
#' gives the position of the sampled items in the original grid.
#'
#' @note When \code{x} is of class \code{"Grid"}, it may be the case
#' that the drawings with \code{atSample = FALSE} have no meaning since some
#' of the variables may be discrete variables and not continuous ones. 
#' 
#' @seealso The \code{\link{sample}} function.
#'
#' @examples
#' g <- Grid(levels = list("Temp" = c(400, 450, 500),
#'                         "Bore" = c(0, 600, 1200, 1800, 2400)))
#' X1 <- sampleIn(g, size = 4)
#' X2 <- sampleIn(g, size = 4, atSample = TRUE)
#' 
#' ## this must be zero
#' sum(abs(as.matrix(g)[attr(X2, "index"), ] - X2))
#'
#'
#'
##*****************************************************************************
## Sample
##*****************************************************************************
sampleIn <-  function(X, size = 1L, atSample = FALSE, ...) standardGeneric("sampleIn")
if (!isGeneric("sampleIn")) {
    setGeneric("sampleIn", sampleIn)
}

#' @title Sample in matrix
#' @param X design matrix
#' @param size Number of sample points.
#' @param atSample should new sample points obtained using \code{sample}
#' @param ... optional args passed to 'sample()'
#' @method sampleIn matrix
#' @exportMethod sampleIn
#' @rdname sampleIn
#' @importFrom stats runif
setMethod("sampleIn", 
          signature = signature(X = "matrix"),
          definition = function(X, size = 1L, atSample = FALSE, ...) {
              if (atSample) {
                  ind <- sample(x = 1L:nrow(X), size = size, ...)
                  XNew <- X[ind, , drop = FALSE]
                  attr(XNew, which = "index") <- ind
              } else {
                  XRange <- apply(X = X, MARGIN = 2, FUN = range)
                  XNew <- matrix(NA, nrow = size, ncol = ncol(X))
                  colnames(XNew) <- colnames(X)
                  for (j in 1:ncol(X)) {
                      XNew[ , j] <- runif(size, min = XRange[1, j],  max = XRange[2, j])
                  }
              }
              XNew   
          }
      )

#' @title Sample in data.frame
#' @param X design data.frame
#' @param size Number of sample points.
#' @param atSample should new sample points obtained using \code{sample}
#' @param ... optional args passed to 'sample()'
#' @method sampleIn data.frame
#' @rdname sampleIn
#' @exportMethod sampleIn
#' @importFrom stats runif
setMethod("sampleIn", 
          signature = signature(X = "data.frame"),
          definition = function(X, size = 1L, atSample = FALSE, ...) {
              sampleIn(X = as.matrix(X), size = size, atSample = atSample, ...)
          }
          )

#' @title Sample in Grid
#' @param X design Grid
#' @param size Number of sample points.
#' @param atSample should new sample points obtained using \code{sample}
#' @param ... optional args passed to 'sample()'
#' @method sampleIn Grid
#' @rdname sampleIn
#' @exportMethod sampleIn
#' @importFrom stats runif
setMethod("sampleIn", 
          signature = signature(X = "Grid"),
          definition = function(X, size = 1L, atSample = FALSE, ...) {
              if (atSample) {
                  xMat <- as.matrix(X)
                  ind <- sample(x = 1L:nrow(xMat), size = size, ...)
                  XNew <- xMat[ind, , drop = FALSE]
                  attr(XNew, which = "index") <- ind
              } else {
                  XRange <- as.data.frame(lapply(levels(X), range))
                  XNew <- matrix(NA, nrow = size, ncol = dim(X))
                  colnames(XNew) <- dimnames(X)
                  for (j in 1L:dim(X)) {
                      XNew[ , j] <- runif(size, min = XRange[1, j],  max = XRange[2, j])
                  }
              }
              XNew   
          }
      )

##=======================================================================
#' Drop the dimensions with a unique level in a \code{Grid} object.
#'
#' @title Drop the dimensions with a unique level in a Grid
#' object
#' 
#' @param x An object with class \code{Grid} 
#'
#' @param ... Not used yet.
#' 
#' @export
#' 
#' @importFrom methods new
#'
#' @return An object of class \code{Grid} but with the fixed
#' dimensions removed.
#'
#' @seealso The usual \code{\link{drop}} function for arrays.
#' 
#' @examples
#' myGD <- Grid(nlevels = c(3, 1, 4, 1, 5))
#' myGD
#' plot(myGD)
#' myGD1 <- drop_Grid(myGD)
##  plot(myGD1)
drop_Grid <- function(x, ...) {
    nl <- nlevels(x) 
    fixed <- (nlevels(x) == 1L)
    if (any(fixed)) {
        levNew <- list()
        for (i in 1:dim(x)) {
            if (!fixed[i]) {
                nm <- x@dimNames[[i]]
                levNew[[nm]] <- x@levels[[i]]
            }
        }
        ## print(levNew)
        ind <- order(x@index)
        dim(ind) <- nl[!fixed]
        xNew <- new("Grid",
                    dim = x@dim - sum(fixed),   
                    dimNames = dimnames(x)[!fixed],
                    levels = levNew,
                    index =  ind)
        
        return(xNew)
        
    } else {
        return(x)
    }
    
}

##======================================================================
#' Find boundary points in a \code{Grid} object.
#'
#' When \code{type} is \code{"data.frame"}, the returned object is a
#' data frame containing only the boundary points. When \code{type}
#' is \code{"array"} the result is a logical array with its elements
#' in correspondence with the \code{index} slot of the
#' object. Finally, when \code{type} is \code{"index"}, the result is
#' an integer vector giving the indices of the boundary points in the
#' order of the nodes defined by the object. This order is the one
#' used by the \code{as.data.frame} coercion.
#' 
#' @title Find boundary points in a \code{Grid} object
#' 
#' @param x An object with class \code{"Grid"}.
#'
#' @param type The wanted type of result.
#' 
#' @export
#'
#' @return An object describing the boundary points. See
#' \bold{Details}.
#'
#' @note Remind that when using the \code{plot} method for a
#' \code{Grid} object, some nodes are generally hidden because
#' several points have the same projection when they are shown in a
#' two-dimensional scatterplot.
#'
#' When one or more of the levels has length \eqn{\leq 2}{<= 2}, all
#' points of the grid are boundary points!
#'
#' @examples
#' ## define a Grid object
#' myGD <- Grid(nlevels = c(3, 4))
#' bd <- boundary_Grid(myGD, type = "index")
#'
#' ## use a different color for boundary points
#' cols <- rep("black", length(myGD))
#' cols[bd] <- "red"
#' plot(myGD, col = cols, pch = 16, cex = 2, main = "Boundary points")
#'
#' ## repeat this after a generalised transposition
#' myGD2 <- aperm(myGD, perm = c(2, 1))
#' bd2 <- boundary_Grid(myGD2, type = "index")
#' cols2 <- rep("black", length(myGD2))
#' cols2[bd2] <- "red"
#' plot(myGD2, col = cols2, pch = 16, cex = 2, main = "Boundary points")
#'
#' ## 3-dimensional
#' myGD3 <- Grid(nlevels = c("x" = 3, "y"= 4, "z" = 6))
#' bd3 <- boundary_Grid(myGD3, type = "index")
#' cols3 <- rep("black", length(myGD3))
#' cols3[bd3] <- "red"
#' plot(myGD3, jitter = TRUE, col = cols3, pch = 16, cex = 2, main = "Boundary points")
boundary_Grid <- function(x, type = c("data.frame", "array", "index")) {
    
    result <- match.arg(type)
    nl <- nlevels(x)
    n <- prod(nl)
    logA <- array(FALSE, dim = nlevels(x), dimnames = levels(x))
    
    for (i in 1:dim(x)) {
        si <- slice.index(x@index, i)
        logA[si == 1] <- TRUE
        logA[si == nl[i]] <- TRUE
    }
    
    if (result == "array") {
        return(logA)
    } else {
        dni <- as.vector(x@index)
        ## ind is the vector of number in the R-canonical order
        ind <- dni[logA]
        if (result == "index") return(ind)
        else if (result == "data.frame") {
            df <- as.data.frame(x)[ind, ]
            return(df)
        }
    }
    
}

##======================================================================
#' Subgrid by selection of levels in one dimension.
#'
#' @title Subgrid by selection of levels in one dimension.
#'
#' @param X An object with class \code{"Grid"}.
#'
#' @param subset An expression concerning \emph{one} of the dimension of \code{X}.
#' This can be compared to an expression in the \code{subset} method.
#' 
#' @param type The wanted type of result. When \code{type} is
#' \code{"index"}, the returned result is an integer vector
#' containing the indices of the subgrid nodes in the original
#' grid. When \code{type} is \code{"Grid"} the result is the subgrid
#' as an object of this class. Finally when \code{type} is
#' \code{"both"} a list is returned containing the two previous types
#' of result under the names \code{index} and \code{Grid}.
#' 
#' @param drop Logical. If \code{TRUE} and if only one level is selected
#' with \code{subset} then the corresponding dimension of the (flat)
#' \code{Grid} result will be dropped.
#' 
#' @export
#' 
#' @importFrom methods is
#'
#' @return A vector of integer indices, a \code{Grid} object or a list
#' embedding these two elements, depending on the value of \code{type}.
#'
#' @seealso The \code{\link{subset}} method.
#'
#' @note The new \code{Grid} returned (if any) uses consecutive
#' numbers for the nodes between \eqn{1} and the number of new nodes. So it
#' no longer keeps trace of the nodes position in \code{X}.  However
#' the surviving nodes are kept in the order that they had in
#' \code{X}.
#'
#' When a response vector or matrix say \code{Y} is related to the
#' original grid \code{X}, the index vector returned by
#' \code{subgrid} can be used to select the corresponding responses.
#' 
#' @examples
#' myGD <- Grid(levels = list("Temp" = c(250, 260, 280, 290, 300),
#'                            "Press" = c(3, 4, 6, 10)))
#' myGD2 <- subset_Grid(myGD, subset = Temp > 280) 
#' 
subset_Grid <- function(X, subset,
                    type = c("index", "Grid", "both"), drop = TRUE) {
  
    type <- match.arg(type)
    e <- substitute(subset)
    if (!is(X, "Grid")) {
        stop("'X' must be of class \"Grid\"")
    }
    av <- all.vars(e)
    if (length(av) > 1L) {
        stop("'subset' must concern only one variable")
    }
    j <- match(av, names(levels(X)))
    if (is.na(j)) {
        stop("variable '", av, "'does not match one of the dimension")
    }
    df <- as.data.frame(levels(X)[[j]])
    colnames(df) <- av
    ind <- eval(e, df)

    if (sum(ind) == 0) stop("empty subgrid")
    
    ind <- (1L:nlevels(X)[[j]])[ind]
    ## slice.index(X@index, j)
    ind2 <- X@index[slice.index(X@index, j) %in% ind]
    
    if (type == "index") return(ind2)
    else {
        Levels <- levels(X)
        Index <- Matrix::invPerm(order(ind2))
        
        if ((length(ind) == 1) && drop) {
            Levels[[j]] <- NULL
            dim(Index) <- sapply(Levels, length)
            X2 <- new("Grid",
                      dim = dim(X) - 1L,   
                      dimNames = dimnames(X)[-j],
                      levels = Levels,
                      index = Index)
                      
        } else {
            Levels[[j]] <- Levels[[j]][ind]
            dim(Index) <- sapply(Levels, length)
            X2 <- new("Grid",
                      dim = dim(X),
                      dimNames = dimnames(X),
                      levels = Levels,
                      index = Index)
        }
        if (type == "Grid") return(X2)
        else {
            return(list(index = ind2, Grid = X2))
        }
    }

}