R/mvQuad.R

Defines functions quadrature rescale.NIGrid rescale plot.NIGrid print.NIGrid dim.NIGrid .print.NIGrid.size size.NIGrid size getWeights getNodes copyNIGrid createNIGrid .sgridnD .SpGrSeq .fgridnD .grid1D readRule

Documented in copyNIGrid createNIGrid dim.NIGrid getNodes getWeights plot.NIGrid print.NIGrid quadrature readRule rescale rescale.NIGrid size size.NIGrid

# !diagnostics style=[true]
# !diagnostics level=[all]



# GetPredefinedTypes <- function(print.types = FALSE){
#    rule.names <- names(QuadRules)
#    rules <- vector(mode = "list", length = length(rule.names))
#
#    for (i in seq_along(rule.names)){
#       tmp <- NULL
#       for (j in seq_len(length(QuadRules[[rule.names[i]]]))){
#          len <- length(QuadRules[[rule.names[i]]][[j]][[2]])
#          if (len==0) len <- NA
#          tmp <- c(tmp, len)
#       }
#       rules[[i]] <- list(levels = tmp)
#       if (print.types){
#          cat("-----", rule.names[i], "---------------------- \n",
#              "levels: \n ",
#              tmp, "\n\n" )
#       }
#    }
#    names(rules) <- rule.names
#    invisible(rules)
# }

.hardCodedTypes <- c("cNC1", "cNC2", "cNC3", "cNC4", "cNC5", "cNC6",
                    "oNC0", "oNC1", "oNC2", "oNC3")

.extGaussQuad <- c("GLe", "GLa", "GHe", "GHN")

.preDefinedTypes <- structure(list(nLe = structure(list(levels = c(1L, 3L, 3L, 7L, 7L, 7L, 15L,
                                                                   15L, 15L, 15L, 15L, 15L, 31L, 31L, 31L, 31L, 31L, 31L, 31L,
                                                                   31L, 31L, 31L, 31L, 31L, 63L)), .Names = "levels"),
                                   GKr = structure(list(levels = c(NA, NA, 3L, NA, 5L, NA, 7L, NA, 9L, NA, 11L,
                                                                   NA, 13L, NA, 15L, NA, 17L, NA, 19L, NA, 21L, NA, 23L,
                                                                   NA, 25L, NA, 27L, NA, 29L)), .Names = "levels"),
                                   nHe = structure(list(levels = c(1L, 3L, 3L, 7L, 9L, 9L, 9L, 9L, 17L, 19L,
                                                        19L, 19L, 19L, 19L, 19L, 31L, 33L, 35L, 35L, 35L, 35L,
                                                        35L, 35L, 35L, 35L)), .Names = "levels"),
                                   nHN = structure(list(levels = c(1L, 3L, 3L, 7L, 9L, 9L, 9L, 9L, 17L, 19L,
                                                                   19L, 19L, 19L, 19L, 19L, 31L, 33L, 35L, 35L, 35L, 35L,
                                                                   35L, 35L, 35L, 35L)), .Names = "levels"),
                                   Leja = structure(list(levels = 1:141), .Names = "levels")),
                              .Names = c("nLe", "GKr", "nHe", "nHN", "Leja"))


#' reads a quadrature-rule from a text file
#'
#' \code{readRule} reads a quadrature-rule from a text file
#'
#' @param file file name of the text file containing the quadrature rule
#'
#' @details The text file containing the quadrature rule has to be formatted in the following way:
#'
#' The first line have to declare the domain \code{initial.domain a b}, where a and b denotes the lower and upper-bound for the integration domain.
#' This can be either a number or '-Inf'/'Inf' (for example \code{initial.domain 0 1} or \code{initial.domain 0 Inf})
#'
#' Every following line contains one single node and weight belonging to one level of the rule (format: 'level' 'node' 'weight').
#' This example shows the use for the "midpoint-rule" (levels: 1 - 3).
#'
#' > \code{initial.domain 0 1}
#'
#' > \code{1 0.5 1}
#'
#' > \code{2 0.25 0.5}
#'
#' > \code{2 0.75 0.5}
#'
#' > \code{3 0.166666666666667 0.333333333333333}
#'
#' > \code{3 0.5 0.333333333333333}
#'
#' > \code{3 0.833333333333333 0.333333333333333}
#'
#'
#' @return Returns an object of class 'customRule', which can be used for creating a 'NIGrid' (\code{\link{createNIGrid}})
#'
#' @seealso
#' \code{\link{createNIGrid}}
#'
#' @examples
#' \dontrun{myRule <- readRule(file="midpoint_rule.txt")}
#' \dontrun{nw <- createNIGrid(d=1, type = myRule.txt, level = 2)}

readRule <- function(file=NULL){
  if (is.null(file)) stop("filename necessary")

  tmp <- readLines(file)

  tmp.index <- which(grepl("initial.domain", tmp))
  if (length(tmp.index)!=1) stop("no or more than one 'initial.domain' defined in your file")
  initial.domain <- matrix(as.numeric(strsplit(tmp[tmp.index], split = " ", fixed = T)[[1]][-1]), ncol=2)
  tmp <- strsplit(tmp, split = " ", fixed = T)
  tmp <- do.call(rbind, tmp)
  colnames(tmp) <- c("l", "n", "w")

  suppressWarnings(levels <- unique(na.omit(as.numeric(tmp[,"l"]))))

  rule <- vector(mode = "list", length = max(levels))

  for (i in levels){
    tmp.l <- subset(tmp, tmp[,1]==i)
    rule[[i]] <- list(n = as.numeric(tmp.l[,2]),
                      w = as.numeric(tmp.l[,3]),
                      features = list(initial.domain=initial.domain))
  }

  class(rule) <- "costumRule"
  return(rule)
}


.grid1D <- function(type, level=1) {

   if (inherits(type, "costumRule")){
     if (!(level %in% c(1:length(type)))){
       stop(paste("degree (",level,") for user defined rule not supported \n") )
     } else {
       return(type[[level]])
     }
   }

  if (!inherits(type, "character")) stop("type of quadrature rule not appropriate defined")

   if (!((type %in% c(.hardCodedTypes, names(.preDefinedTypes), .extGaussQuad)))) {
      if (!existsFunction(type)) {
         stop("type of quadrature rule is not supported")
      } else {
         tmp.fun <- match.fun(type)
         return(tmp.fun(level))
      }
   }

   if (type %in% .hardCodedTypes) {
      quad.type <- substr(type,1,3)

      # closed-Newton-Cotes forumla
      if (quad.type=="cNC") {

         degree <- as.numeric(substr(type,4,4))
         if (!(degree %in% c(1:6)) ) stop("degree of closed-Newton-Cotes formula must between 1 and 6")

         if (degree == 1) weights <- c(1/2, 1/2)
         if (degree == 2) weights <- c(1/6, 4/6, 1/6)
         if (degree == 3) weights <- c(1/8, 3/8, 3/8, 1/8)
         if (degree == 4) weights <- c(7/90, 32/90, 12/90, 32/90, 7/90)
         if (degree == 5) weights <- c(19/288, 75/288, 50/288, 50/288, 75/288, 19/288)
         if (degree == 6) weights <- c(41/840, 216/840, 27/840, 272/840, 27/840, 216/840, 41/840)

         interv_length <- 1/level
         weights <- weights * interv_length
         steps <- interv_length/degree
         n <- seq(0,1,steps)
         w <- numeric(length(n))

         for (i in 1:level) {
            w[((i-1)*degree+1):(i*degree+1)] <- w[((i-1)*degree+1):(i*degree+1)] + weights
         }

         return(list(n = n,w = w, features=list(initial.domain=matrix(c(0,1), ncol = 2))))
      }

      # open-Newton-Cotes Forumla
      if (quad.type=="oNC") {
         degree <- as.numeric(substr(type,4,4))
         if (!(degree %in% c(0:3)) ) stop("degree of open-Newton-Cotes formula must between 0 and 3")

         if (degree == 0) {
            weights <- c(1)
            nodes <- c(1/2)
         }
         if (degree == 1) {
            weights <- c(1/2, 1/2)
            nodes <- c(1/4, 3/4)
         }
         if (degree == 2) {
            weights <- c(3/8, 2/8, 3/8)
            nodes <- c(1/6, 1/2, 5/6)
         }
         if (degree == 3) {
            weights <- c(13/48, 11/48, 11/48, 13/48)
            nodes <- c(1/8, 3/8, 5/8, 7/8)
         }
         interv_length <- 1/level
         weights <- weights * interv_length
         nodes <- nodes * interv_length

         n = numeric((degree+1)*level)
         w = rep(weights,level)

         for (i in 1:level){
            n[((i-1)*(degree+1)+1):(i*(degree+1))] <- (i-1)*interv_length + nodes
         }

         return(list(n = n, w = w, features=list(initial.domain=matrix(c(0,1), ncol = 2))))
      }
   }

  if (type %in% names(.preDefinedTypes)) {
    if (!(level %in% c(1:length(QuadRules[[type]]))) ){
      stop(paste("degree (",level,") for rule '", type, "' not supported \n") )
    } else {
      return(QuadRules[[type]][[level]])
    }
  }

  if (type %in% .extGaussQuad) {
    if (type=="GLe") {
      tmp.nw <- gauss.quad(level, kind="legendre")
      tmp.nw$nodes <- tmp.nw$nodes / 2 + 0.5
      tmp.nw$weights <- tmp.nw$weights / 2
      tmp.dom <- list(initial.domain=matrix(c(0,1), ncol = 2))
    }

    if (type=="GLa") {
      tmp.nw <- gauss.quad(level, kind="laguerre")
      tmp.dom <- list(initial.domain=matrix(c(0,Inf), ncol = 2))
    }

    if (type=="GHe") {
      tmp.nw <- gauss.quad(level, kind="hermite")
      tmp.nw$nodes <- tmp.nw$nodes * sqrt(2)
      tmp.nw$weights <- tmp.nw$weights * sqrt(2) / sqrt(exp(-(tmp.nw$nodes^2)))
      tmp.dom <- list(initial.domain=matrix(c(-Inf, Inf), ncol = 2))
    }
    if (type=="GHN") {
      tmp.nw <- gauss.quad(level, kind="hermite")
      tmp.nw$nodes <- tmp.nw$nodes * sqrt(2)
      tmp.nw$weights <- tmp.nw$weights * sqrt(2) / sqrt(exp(-(tmp.nw$nodes^2)))*dnorm(tmp.nw$nodes)
      tmp.dom <- list(initial.domain=matrix(c(-Inf, Inf), ncol = 2))
    }

    return(list(n = as.numeric(tmp.nw$nodes), w = as.numeric(tmp.nw$weights), features=tmp.dom))
  }

}

.fgridnD <- function(dim, type, level) {
  le <- numeric(dim)
  domain <- matrix(NA, nrow = dim, ncol = 2)
  storage <- vector(mode = "list", length = dim)
  for (i in 1:dim) {
    tmp <- .grid1D(type[[i]], level[i])
    le[i] <- length(tmp$w)
    domain[i, ] <- tmp$features$initial.domain
    storage[[i]] <- tmp
  }

  No.Points <- prod(le)
  w <- rep(1, No.Points)
  n <- matrix(NA, nrow = No.Points, ncol = dim)

  for (i in 1 :dim) {
    storage[[i]]
    n[ ,i] = rep(storage[[i]]$n, each = prod(le[1:i - 1]), times = (No.Points/prod(le[1:i])))
    w = w * rep(storage[[i]]$w, each = prod(le[1:i - 1]), times = (No.Points/prod(le[1:i])))
  }
  return(list(n = n, w = w, features = list(initial.domain = domain)))
}

.SpGrSeq <- function(dimension, norm, level.trans){
  seq.vec       <- rep(0, dimension)
  a             <- norm - dimension
  seq.vec[ 1 ]  <- a
  fs            <- matrix(NA, nrow=choose(norm - 1, dimension - 1), ncol=dimension )
  fs[1, ]       <- seq.vec
  cnt           <- 1
  cnt.seq       <- 1

  while ( seq.vec[ dimension ] < a ) {
    if ( cnt == dimension ) {
      for (i in seq(cnt-1, 1, -1) ) {
        cnt <- i
        if ( seq.vec[ i ] != 0 ) { break }
      }
    }
    seq.vec[ cnt ]  <- seq.vec[ cnt ] - 1
    cnt             <- cnt + 1
    seq.vec[ cnt ]  <- a - sum( seq.vec[1:(cnt-1)] )
    if ( cnt < dimension ) {
      seq.vec[ (cnt+1):dimension ] <- rep(0, dimension - cnt)
    }
    cnt.seq <- cnt.seq + 1
    fs[cnt.seq, ] = seq.vec
  }

  fs <- fs + 1
  fs <- level.trans(fs)

  return(fs)
}



.sgridnD <- function(dim, type, level, level.trans){
  minq <- max(0,level-dim)
  maxq <- level-1

  noSubGrids <- sum(choose(c(minq:maxq) + dim - 1, dim - 1))
  maxStorSize <- 100 # fine tuning?
  storage <- vector(mode="list", length = maxStorSize)

  nw.table <- data.table()
  nw.names <- c(paste("n", c(1:dim), sep="_"), "w.c")

  cnt <- 0
  w.c <- NULL
  for (q in (minq:maxq)) {
    bq  <- (-1)^(maxq-q)*choose(dim-1,dim+q-level)

    GridSeq <- .SpGrSeq(dim,dim+q,level.trans)
    for (i in 1:dim(GridSeq)[1]){
      res <- .fgridnD(dim,type,GridSeq[i,])

      cnt <- cnt + 1
      storage[[cnt]] <- data.table::data.table(res[[1]], bq*res[[2]])

      if (cnt == maxStorSize) {
        nw.table <- rbindlist(c(list(nw.table), storage))
        data.table::setnames(nw.table, nw.names)
        data.table::setkeyv(x = nw.table, cols = nw.names[-(dim+1)])
        nw.table[, w:=sum(w.c), by = eval(nw.names[-(dim+1)])]
        nw.table <- unique(nw.table, by = eval(nw.names[-(dim+1)]))
        nw.table[,w.c:=NULL]

        cnt <- 0
      }
    }
  }

  if (cnt > 0) {
    nw.table <- rbindlist(c(list(nw.table), storage[1:cnt]))
    data.table::setnames(nw.table, nw.names)
    data.table::setkeyv(x = nw.table, cols = nw.names[-(dim+1)])
    nw.table[, w:=sum(w.c), by = eval(nw.names[-(dim+1)])]
    nw.table <- unique(nw.table, by = eval(nw.names[-(dim+1)]))
    nw.table[,w.c:=NULL]
  }

  nw.table <- nw.table[w!=0]

  n <- nw.table[,nw.names[-(dim+1)], with=FALSE]
  w <- nw.table[, w]

  row.names(n) <- NULL
  colnames(n) <- NULL

  return(list(n = as.matrix(n), w = as.matrix(w), features=res[[3]]))
}


#' creates a grid for numerical integration.
#'
#' \code{createNIGrid} Creates a grid for multivariate numerical integration.
#' The Grid can be based on different quadrature- and construction-rules.
#'
#' @param dim number of dimensions
#' @param type quadrature rule (see Details)
#' @param level accuracy level (typically number of grid points for the underlying 1D quadrature rule)
#' @param ndConstruction character vector which denotes the construction rule
#' for multidimensional grids.
#' \describe{
#' \item{\code{product}}{for product rule, returns a "full grid" (default)}
#' \item{\code{sparse}}{for combination technique, leads to a regular "sparse grid".}}
#'
#' @param level.trans logical variable denotes either to take the levels as number
#' of grid points (FALSE = default) or to transform in that manner that number of
#' grid points = 2^(levels-1) (TRUE). Alternatively \code{level.trans} can be a function, which takes (n x d)-matrix and returns
#' a matrix with the same dimensions (see the example; this feature is particularly useful for the 'sparse' construction rule,
#' to account for different importance of the dimensions).
#'
#' @details The following quadrature rules are supported (build-in).
#' \describe{
#'  \item{\code{cNC1, cNC2, ..., cNC6}}{closed Newton-Cotes Formula of degree 1-6 (1=trapezoidal-rule; 2=Simpson's-rule; ...),
#'  initial interval of integration: [0, 1]}
#'  \item{\code{oNC0, oNC1, ..., oNC3}}{open Newton-Cote Formula of degree 0-3 (0=midpoint-rule; ...),
#'  initial interval of integration: [0, 1]}
#'  \item{\code{GLe, GKr}}{Gauss-Legendre and Gauss-Kronrod rule for an initial interval of integration: [0, 1]}
#'  \item{\code{nLe}}{nested Gauss-Legendre rule for an initial interval of integration: [0, 1] (Knut Petras (2003). Smolyak cubature of given polynomial degree with few nodes for increasing dimension. Numerische Mathematik 93, 729-753)}
#'  \item{\code{GLa}}{Gauss-Laguerre rule for an initial interval of integration: [0, INF)}
#'  \item{\code{GHe}}{Gauss-Hermite rule for an initial interval of integration: (-INF, INF)}
#'  \item{\code{nHe}}{nested Gauss-Hermite rule for an initial interval of integration: (-INF, INF) (A. Genz and B. D. Keister (1996). Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight." Journal of Computational and Applied Mathematics 71, 299-309)}
#'  \item{\code{GHN}, \code{nHN}}{(nested) Gauss-Hermite rule as before but weights are multiplied by the standard normal density (\eqn{\hat(w)_i = w_i * \phi(x_i)}).}
#'  \item{\code{Leja}}{Leja-Points for an initial interval of integration: [0, 1]}}
#' The argument \code{type} and \code{level} can also be vector-value, different for each dimension (the later only for "product rule"; see examples)
#'
#' @return Returns an object of class 'NIGrid'. This object is basically an environment
#' containing nodes and weights and a list of features for this special grid. This
#' grid can be used for numerical integration (via \code{\link{quadrature}})
#' @references
#' Philip J. Davis, Philip Rabinowitz (1984): Methods of Numerical Integration
#'
#' F. Heiss, V. Winschel (2008): Likelihood approximation by numerical integration on sparse grids, Journal of Econometrics
#'
#' H.-J. Bungartz, M. Griebel (2004): Sparse grids, Acta Numerica
#' @seealso
#' \code{\link[=rescale.NIGrid]{rescale}}, \code{\link{quadrature}}, \code{\link[=print.NIGrid]{print}}, \code{\link[=plot.NIGrid]{plot}} and \code{\link[=size.NIGrid]{size}}
#' @examples
#' ## 1D-Grid --> closed Newton-Cotes Formula of degree 1 (trapeziodal-rule)
#' myGrid <- createNIGrid(dim=1, type="cNC1", level=10)
#' print(myGrid)
#' ## 2D-Grid --> nested Gauss-Legendre rule
#' myGrid <- createNIGrid(dim=2, type=c("GLe","nLe"), level=c(4, 7))
#' rescale(myGrid, domain = rbind(c(-1,1),c(-1,1)))
#' plot(myGrid)
#' print(myGrid)
#' myFun <- function(x){
#'    1-x[,1]^2*x[,2]^2
#' }
#' quadrature(f = myFun, grid = myGrid)
#' ## level transformation
#' levelTrans <- function(x){
#'   tmp <- as.matrix(x)
#'   tmp[, 2] <- 2*tmp[ ,2]
#'   return(tmp)
#' }
#' nw <- createNIGrid(dim=2, type="cNC1", level = 3,
#'    level.trans = levelTrans, ndConstruction = "sparse")
#' plot(nw)
createNIGrid <- function(dim=NULL, type=NULL, level=NULL,
                         ndConstruction="product", level.trans=NULL){
## to-do:
#     - adaptive grid (CW, 2015-01-05)

   if (is.null(dim) | (dim%%1)!=0 | dim < 1) {
      stop("'dim' is not appropriate defined")
   }

   if (!inherits(type, "list")){
     if (inherits(type, "character")) type <- as.list(type)
     if (inherits(type, "costumRule")) type <- list(type)
   }

   if (is.null(type) | (length(type) > 1 & length(type) < dim)) {
      stop("'type' is not appropriate defined")
   }

   if (length(type) < dim){
      type <- replicate(type[[1]], n=dim, FALSE)
   }

   if (is.null(level) | any(level%%1!=0)){
      stop("'level' is not appropriate defined")
   }

   if (length(level) < dim){
      level <- rep(level, dim)
   }

   level <- matrix(level, ncol = dim)

   if (is.logical(level.trans)) {
      if (level.trans == TRUE){
         level.trans <- function(x){2^(x-1)}
      }
   }

   if (!is.function(level.trans)){
      level.trans <- function(x){x}
   }

   if (dim > 1){
      if (ndConstruction=="product"){
         nw <- .fgridnD(dim = dim, type = type, level = level.trans(level))
      }
      if (ndConstruction=="sparse"){
         nw <- .sgridnD(dim = dim, type = type, level = level[1], level.trans = level.trans)
      }
   } else {
      nw <- .grid1D(type = type[[1]], level = level.trans(level))
   }

   object=new.env(parent=globalenv())

   object$dim <- dim

   type.text <- lapply(type,
                       function(tmp){
                         if (inherits(tmp, "costumRule")) return("costum")
                         if (inherits(tmp, "character"))  return(tmp)
                       })
   object$type <- unlist(type.text)
   object$level <- level
   object$level.trans <- level.trans

   object$ndConstruction <- ndConstruction

   if (exists("features", where=nw)){
     object$features <- c(type="static", move=0L, nw[["features"]])
   } else {
     object$features <- list(type="static", move=0L, initial.domain = matrix(c(NA, NA), ncol=2))
   }

   if (!exists("n", where=nw) | !exists("n", where=nw)) stop("this quadrature rule doesn't provide nodes (or/and) weights")

   object$nodes <- as.matrix(nw[["n"]])
   object$weights <- as.matrix(nw[["w"]])
   class(object) <- "NIGrid"
   return(object)
}

#' copies an NIGrid-object
#'
#' \code{copyNIGrid} copies an NIGrid-object
#'
#' @param object1 original NIGrid-object
#' @param object2 destination; if NULL \code{copyNIGrid} returns a NIGrid-object
#' otherwise the \code{object2} will be overwritten.
#'
#' @return Returns a NIGrid-object or NULL
#' @examples
#' myGrid <- createNIGrid(dim=2, type="GHe", level=5)
#' myGrid.copy <- copyNIGrid(myGrid)
copyNIGrid <- function(object1, object2=NULL){
   if (is.null(object2)) {
      object2 = new.env(parent = globalenv())
      class(object2) = class(object1)
      nullFlag = TRUE
   }
   elements = ls(envir = object1)
   for (index in 1:length(elements)) {
      assign(elements[[index]], get(elements[[index]], envir = object1,
                                  inherits = FALSE), envir = object2)
   }
   if (nullFlag) {
      return(object2)
   } else {
      return(NULL)
   }
}

#' @name getNodes and getWeights
#' @rdname getNodesWeights
#'
#' @title get nodes and weights from an NIGrid-object
#' @description \code{getNodes} and \code{getWeights} extract the (potentially rescaled) nodes and weights
#'  out of an NIGrid-Object
#'
#' @param grid object of class \code{NIGrid}
#' @return Returns the nodes or weights of the given grid
#' @seealso
#' \code{\link{createNIGrid}}
#' @examples
#' myGrid <- createNIGrid(dim=2, type="cNC1", level=3)
#' getNodes(myGrid)
getNodes <- function(grid){
   if (!is(grid, "NIGrid")) { stop(" 'object' argument must be of class 'NIGrid' .") }
   if (grid$features$move==0) return(grid$nodes)

   if (grid$features$move==1) {
      n <- grid$nodes
      for (d in 1:grid$dim){
         n[,d] <- (grid$nodes[,d]-grid$features$initial.domain[d,1])*
            diff(grid$features$domain[d,])/diff(grid$features$initial.domain[d,]) +
            grid$features$domain[d,1]
      }
   }

   if (grid$features$move==2) {
      C <- as.matrix(grid$features$C)
      m <- grid$features$m


      if (grid$dim > 1) {
        if (grid$features$dec.type==0) {
          A <- diag(sqrt(diag(C)))
        }
        if (grid$features$dec.type==1) {
          res.dec <- eigen(C)
          A <- res.dec[[2]] %*% diag(sqrt(res.dec[[1]]))
        }
        if (grid$features$dec.type==2) {
          A <- t(chol(C))
        }
      } else {
        A <- as.matrix(sqrt(C))
      }

      n <- t(A %*% t(grid$nodes))
      for (d in 1:grid$dim){
         n[, d] <- n[, d] + m[d]
      }



   }
   return(n)
}


#' @rdname getNodesWeights
#' @examples
#' getWeights(myGrid)
getWeights <- function(grid){
   if (!is(grid, "NIGrid")) { stop(" 'object' argument must be of class 'NIGrid' .") }
   if (grid$features$move==0) return(grid$weights)

   if (grid$features$move==1) {
      w <- grid$weights
      for (d in 1:grid$dim){
         w <- w * diff(grid$features$domain[d,])/diff(grid$features$initial.domain[d,])
      }
   }

   if (grid$features$move==2) {
      C <- as.matrix(grid$features$C)

      if (grid$features$dec.type==0) {
         w <- grid$weights * prod(sqrt(diag(C)))
      }
      if (grid$features$dec.type==1) {
         res.dec <- eigen(C)
         w <- grid$weights * prod(sqrt(res.dec[[1]]))
      }
      if (grid$features$dec.type==2) {
         A <- t(chol(C))
         w <- grid$weights * prod(diag(A))
      }
   }

   return(w)
}

#' @name size (size.NIGrid)
#' @rdname size.NIGrid
#'
#' @title returns the size of an NIGrid-object
#' @description Returns the size of an NIGrid-object
#' @param ... other arguments passed to the specific method
#'
#' @return Returns the grid size in terms of dimensions, number of grid points and used memory
#' @examples
#' myGrid <- createNIGrid(dim=2, type="GHe", level=5)
#' size(myGrid)
#' @export
size <- function(object, ...) UseMethod("size")

#' @rdname size.NIGrid
#' @param object a grid of type \code{NIGrid}
#' @export
size.NIGrid <- function(object, ...){
   size.dim <- object$dim
   size.no  <- length(object$weights)
   size.mem <- object.size(object$nodes)+object.size(object$weights)
   size.object <- list(dim=size.dim, gridpoints=size.no, memory=size.mem)
   class(size.object) <- "NIGrid.size"
   return(size.object)
}

.print.NIGrid.size <- function(x, ...){
   cat("grid size \n",
       " # dimensions: ", x[[1]],
       "\t # gridpoints: ", x[[2]],
       "\t used memory:", format(x[[3]], units="auto"), "\n", sep="")
}

#' @rdname size.NIGrid
#' @param x object of type \code{NIGrid}
#' @examples
#' dim(myGrid)
#' @export
dim.NIGrid <- function(x){
   return(x$dim)
}


#' @rdname print.NIGrid
#' @name print (print.NIGrid)
#' @title prints characteristic information for an NIGrid-object
#' @description Prints characteristic information for an NIGrid-object
#' @param x a grid of type \code{NIGrid}
#' @param ... further arguments passed to or from other methods
#'
#' @return Prints the information for an NIGrid-object (i.a. grid size (dimensions, grid points, memory usage),
#' type and support)
#' @examples
#' myGrid <- createNIGrid(dim=2, type="GHe", level=5)
#' print(myGrid)
#' @export
print.NIGrid <- function(x, ...){
   tmp <- size(x)
   cat("Grid for Numerical Integration  \n",
       " # dimensions: ", tmp[[1]],
       "\t # gridpoints: ", tmp[[2]],
       "\t used memory:", format(tmp[[3]], units="auto"), "\n\n", sep="")

   if (x$dim > 1) cat("nD-Construction principle:  '", x$ndConstruction, "'\n", sep="")
   tmp.level <- as.vector(x$level)
   if (length(unique(x$type))==1 & length(unique(tmp.level))==1) {
      cat(" same quadrature rule for each dimension \n")
      cat(" \t type: ", unique(x$type), "\n")
      cat(" \t level: ", unique(tmp.level), "\n")
      cat(" \t initial domain: ", x$features$initial.domain[1,])
   } else {
      cat(" individual quadrature rule for each dimension \n")
      for (d in 1:x$dim){
         cat(" dim =", d, "--> \t type:", x$type[d], "\t level: ",
             tmp.level[d], "\t initial domain:", x$features$initial.domain[d,], "\n")
      }
   }

   if (x$features$move!=0){
      if (x$features$move == 1){
         cat("\n Grid rescaled. New Domain: \n")
         tmp <- as.matrix(x$features$domain)
         colnames(tmp) <- c("a", "b")
         rownames(tmp) <- paste("dim ", 1:x$dim, ": ", sep="")
         print(tmp)
      }

      if (x$features$move == 2){
         cat("\n Grid rescaled.")
         cat("\n mean-vector: ", x$features$m)
         cat("\n covariance matrix: \n")
         print(x$features$C)
      }
   }
   cat("\n")
}

#' @rdname plot.NIGrid
#' @name plot (plot.NIGrid)
#' @title plots an NIGrid-object
#' @description
#' Plots the grid points of an NIGrid-object
#'
#' @param x a grid of type \code{NIGrid}
#' @param plot.dimension vector of length 1, 2 or 3. with the dimensions to be plotted
#' (see examples)
#' @param ... arguments passed to the default plot command
#'
#' @examples
#' myGrid <- createNIGrid(dim=4, type=c("GHe", "cNC1", "GLe", "oNC1"),
#'                        level=c(3,4,5,6))
#' plot(myGrid) ## dimension 1-min(3,dim(myGrid)) are plotted
#' ## Free arranged plots
#' plot(myGrid, plot.dimension=c(4,2,1))
#' plot(myGrid, plot.dimension=c(1,2))
#' plot(myGrid, plot.dimension=c(3))
#' @export
plot.NIGrid <- function(x, plot.dimension=NULL, ...){
   p.dim <- 1:min(x$dim, 3)

   arg.list <- list(...)

   if (!is.null(plot.dimension)){
      if (length(plot.dimension) >= 4) {
         cat("only first, second and third stated dimension is used")
         plot.dimension <- as.numeric(plot.dimension[1:3])
      }

      if (any(!(plot.dimension %in% (1:x$dim)))){
         stop("wrong dimensions stated to be plotted")
      }
      p.dim <- plot.dimension
   }

   if (length(p.dim)==1){
      x.tmp <- as.numeric(getNodes(x)[, p.dim])
      y.tmp <- rep(1,length(x.tmp))
      if (length(which("xlab" == names(arg.list))) == 0) arg.list = c(arg.list, xlab=paste("dim", p.dim))
      if (length(which("axes" == names(arg.list))) == 0) {
         arg.list = c(arg.list, axes=FALSE)
         createXaxis = TRUE
      }
      do.call(plot, c(list(x=x.tmp, y=y.tmp, ylab=""), arg.list))
      if (createXaxis) axis(side = 1)
   }
   if (length(p.dim)==2){
      tmp <- getNodes(x)
      x.tmp <- as.numeric(tmp[, p.dim[1]])
      y.tmp <- as.numeric(tmp[, p.dim[2]])

      if (length(which("xlab" == names(arg.list))) == 0) arg.list = c(arg.list, xlab=paste("dim", p.dim[1]))
      if (length(which("ylab" == names(arg.list))) == 0) arg.list = c(arg.list, ylab=paste("dim", p.dim[2]))
      do.call(plot, c(list(x=x.tmp, y=y.tmp), arg.list))
   }
   if (length(p.dim)==3){
     if (requireNamespace("rgl", quietly = TRUE)) {
       tmp <- getNodes(x)
       x.tmp <- as.numeric(tmp[, p.dim[1]])
       y.tmp <- as.numeric(tmp[, p.dim[2]])
       z.tmp <- as.numeric(tmp[, p.dim[3]])
       if (length(which("xlab" == names(arg.list))) == 0) arg.list = c(arg.list, xlab=paste("dim", p.dim[1]))
       if (length(which("ylab" == names(arg.list))) == 0) arg.list = c(arg.list, ylab=paste("dim", p.dim[2]))
       if (length(which("zlab" == names(arg.list))) == 0) arg.list = c(arg.list, zlab=paste("dim", p.dim[3]))
       do.call(rgl::plot3d, c(list(x=x.tmp, y=y.tmp, z=z.tmp), arg.list))
     } else {
       stop("wrong dimensions stated to be plotted (no 'rgl'-package for 3d-plot)")
     }
   }
}

#' @rdname rescale.NIGrid
#' @name rescale (rescale.NIGrid)
#' @title moves, rescales and/or rotates a multidimensional grid.
#' @description
#' \code{rescale.NIGrid} manipulates a grid for more efficient numerical integration with
#' respect to a given domain (bounded integral) or vector
#' of means and covariance matrix (unbounded integral).
#'
#' @param object an initial grid of type \code{NIGrid}
#' @param ... further arguments passed to or from other methods
#' @return This function modifies the "support-attribute" of the grid. The
#' recalculation of the nodes and weights is done when the \code{\link{getNodes}} or \code{\link{getWeights}}
#' are used.
#' @references Peter Jaeckel (2005): A note on multivariate Gauss-Hermite quadrature
#' @seealso
#' \code{\link{quadrature}}, \code{\link{createNIGrid}}
#' @export
rescale <- function(object, ...) UseMethod("rescale")

#' @rdname rescale.NIGrid
#' @param domain a (d x 2)-matrix with the boundaries for each dimension
#' @param m vector of means
#' @param C covariance matrix
#' @param dec.type type of covariance decomposition (\cite{Peter Jaeckel (2005)})
#' @examples
#' C = matrix(c(2,0.9,0.9,2),2)
#' m = c(-.5, .3)
#' par(mfrow=c(3,1))
#'
#' myGrid <- createNIGrid(dim=2, type="GHe", level=5)
#'
#' rescale(myGrid, m=m, C=C, dec.type=0)
#' plot(myGrid, col="red")
#'
#' rescale(myGrid, m=m, C=C, dec.type=1)
#' plot(myGrid, col="green")
#'
#' rescale(myGrid, m=m, C=C, dec.type=2)
#' plot(myGrid, col="blue")
#' @export
rescale.NIGrid <- function(object, domain=NULL, m=NULL, C=NULL, dec.type=0, ...){
   # check plausibility --> bounded vs. unbounded integral (dimension wise)

   if (is.null(domain) & is.null(m) & is.null(C)) {
      object$features[["move"]] <- 0L
      object$features[["domain"]] <- NULL
      object$features[["m"]] <- NULL
      object$features[["C"]] <- NULL
      object$features[["dec.type"]] <- NULL
      return(invisible(NULL))
   }

   if (!is.null(domain) & (!is.null(m) | !is.null(C))) stop("rescaling not appropriate defined")
   if (is.null(domain) & (is.null(m) | is.null(C))) stop("rescaling not appropriate defined")

   if (!is.null(domain)){
      if (any(is.infinite(object$features[["initial.domain"]]))) stop("rescaling not appropriate defined (bounded domain for unbounded rule)")
      object$features[["move"]] <- 1L
      object$features[["domain"]] <- matrix(domain, ncol=2)
      return(invisible(NULL))
   }

   if (is.null(domain)) {
      C <- as.matrix(C)
      C.dim <- dim(C)
      if (C.dim[1]==C.dim[2]) {
         C.dim <- C.dim[1]
      } else stop("(in rescale) covariance matrix C is not quadratic")

      if (!(all(C == t(C)) & all(eigen(C)$values > 0))) {
         # check for symmetry and positive semi definitness
         stop("(in rescale) no appropriate covariance matrix C")
      }

      if (!length(m)==C.dim) stop("(in rescale) vector of mean not appropriate")

      if (!object$dim==C.dim) stop("(in rescale) initial grid not appropriate")
      if (C.dim==1) dec.type<-0


      object$features[["move"]] <- 2L
      object$features[["m"]] <- m
      object$features[["C"]] <- C
      object$features[["dec.type"]] <- dec.type
      return(invisible(NULL))
   }
}


#' @title computes the approximated Integral
#' @description
#'  \code{quadrature} computes the integral for a given function based on an NIGrid-object
#'
#' @param f a function which takes the x-values as a (n x d) matrix as a first argument
#' @param grid a grid of type \code{NIGrid}
#' @param ... further arguments for the function \code{f}
#' @return The approximated value of the integral
#' @seealso \code{\link{createNIGrid}}, \code{\link[=rescale.NIGrid]{rescale}}
#' @examples
#' myGrid <- createNIGrid(dim=2, type="GLe", level=5)
#' rescale(myGrid, domain=rbind(c(-1,1),c(-1,1)))
#' plot(myGrid, col="blue")
#' myFun <- function(x){
#'    1 - x[,1]^2 * x[,2]^2
#' }
#' quadrature(myFun, myGrid)
quadrature <- function(f, grid=NULL, ...){
## to-do:
#    - add adaptivity (CW; 2015-01-05)
#    - add error estimate (CW; 2015-01-05)
   f <- match.fun(f)
   ff <- function(x) f(x, ...)

   if (is.null(grid)) stop("quadrature not appropriate specified")

   if (!is.null(grid) & !is(grid, "NIGrid")) stop("'Grid' isn't of class NIGrid")

   # place to implement adaptivity
   y <- ff(getNodes(grid))
   w <- getWeights(grid)
   if (length(y) != length(w)) stop("function 'f' fits not to the supported Grid")
   return(sum(y*w))
}

Try the mvQuad package in your browser

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

mvQuad documentation built on Sept. 19, 2023, 5:06 p.m.