R/lineqGP.R

In lineqGPR: Gaussian Process Regression Models with Linear Inequality Constraints

#' @title Hat Basis Functions for \code{"lineqGP"} Models
#' @description Evaluate the hat basis functions for \code{"lineqGP"} models.
#' @param x a vector (or matrix) with the input data.
#' @param u a vector (or matrix) with the locations of the knots.
#' @param d a number corresponding to the dimension of the input space.
#' @return A matrix with the hat basis functions. The basis functions are indexed by rows.
#'
#' @section Comments:
#' This function was tested mainly for 1D or 2D input spaces. It could change
#' in future versions for higher dimensions.
#' @author A. F. Lopez-Lopera.
#'
#' @references Lopez-Lopera, A. F., Bachoc, F., Durrande, N., and Roustant, O. (2017),
#' "Finite-dimensional Gaussian approximation with linear inequality constraints".
#' \emph{ArXiv e-prints}
#' \href{https://arxiv.org/abs/1710.07453}{[link]}
#'
#' Maatouk, H. and Bay, X. (2017),
#' "Gaussian process emulators for computer experiments with inequality constraints".
#' \emph{Mathematical Geosciences},
#' 49(5): 557-582.
#' \href{https://link.springer.com/article/10.1007/s11004-017-9673-2}{[link]}
#'
#' @examples
#' x <- seq(0, 1, 1e-3)
#' m <- 5
#' u <- seq(0, 1, 1/(m-1))
#' Phi <- basisCompute.lineqGP(x, u, d = 1)
#' matplot(Phi, type = "l", lty = 2, main = "Hat basis functions with m = 5")
#'
#' @export
basisCompute.lineqGP <- function(x, u, d = 1) {
  if (!is.matrix(x) || ncol(x) != d)
    x <- matrix(x, ncol = d)
  if (is.list(u) && length(u) == d) {
    m <- as.integer(lapply(u, length))
  } else if (is.double(u)) {
    u <- matrix(u, ncol = d)
    m <- nrow(u)
  } else {
    m <- nrow(u)
  }

  # precomputing some terms
  n <- nrow(x)
  delta <- 1/(m-1)

  # computing the hat basis functions
  if (d == 1){
    distAbs <- abs(outer(x[, 1]/delta, u[, 1]/delta, "-"))
    idx <- distAbs <= 1
    Phi <- matrix(0, n, m)
    Phi[idx] <- 1 - distAbs[idx]
  } else if (d >= 2){
    PhiList <- list()
    for (k in seq(d)) {
      PhiList[[k]] <- basisCompute.lineqGP(x[, k], u[[k]], d = 1)
      # distAbs <- abs(outer(x[, k]/delta[k], u[[k]]/delta[k], "-"))
      # idx <- distAbs <= 1
      # Phi_temp <- matrix(0, n, m[k])
      # Phi_temp[idx] <- 1 - distAbs[idx]
      # Phip[[k]] <- Phi_temp
    }
    Phi <- matrix(0, n, prod(m))
    for (i in seq(n)) {
      Phi_temp2 <- PhiList[[1]][i, ]
      for (k in 2:d) {
        Phi_temp2 <- Phi_temp2 %x% PhiList[[k]][i, ]
      }
      Phi[i, ] <- Phi_temp2
    }
  }
  return(Phi)
}

#' @title Linear Systems of Inequalities for \code{"lineqGP"} Models
#' @description Build the linear system of inequalities for \code{"lineqGP"} models.
#' @param m the number of linear inequality constraints.
#' @param constrType a character string corresponding to the type of the inequality constraint.
#' Options: "boundedness", "monotonicity", "convexity", "linear"
#' @param l the value (or vector) with the lower bound.
#' @param u the value (or vector) with the upper bound.
#' @param A a matrix containing the structure of the linear equations.
#' @param d the value with the input dimension.
#' @param lineqSysType a character string corresponding to the type of the
#' linear system. Options: \code{twosides}, \code{oneside} (see \code{\link{bounds2lineqSys}} for more details).
#' @param constrIdx for d > 1, a logical vector with the indices of active constrained dimensions.
#' @param rmInf If \code{TRUE}, inactive constraints are removed
#' (e.g. \eqn{-\infty \leq x \leq \infty}{-Inf \le x \le Inf}).
#' @return  A list with the linear system of inequalities: \code{list(A,l,u)} (\code{twosides}) or \code{list(M,g)} (\code{oneside}).
#'
#' @section Comments:
#' This function could change in future versions for more types of inequality
#' constraints in higher dimensions.
#' @seealso \code{\link{bounds2lineqSys}}
#' @author A. F. Lopez-Lopera.
#'
#' @references Lopez-Lopera, A. F., Bachoc, F., Durrande, N., and Roustant, O. (2017),
#' "Finite-dimensional Gaussian approximation with linear inequality constraints".
#' \emph{ArXiv e-prints}
#' \href{https://arxiv.org/abs/1710.07453}{[link]}
#'
#' @examples
#' linSys1 <- lineqGPSys(m = 5, constrType = "boundedness", l = 0, u = 1, lineqSysType = "twosides")
#' linSys1
#' linSys2 <- lineqGPSys(m = 5, constrType = "boundedness", l = 0, u = 1, lineqSysType = "oneside")
#' linSys2
#'
#' @export
lineqGPSys <- function(m = nrow(A),
                       constrType = c("boundedness", "monotonicity",
                                      "convexity", "linear"),
                       l = -Inf, u = Inf,  A = diag(m), d = length(m),
                       lineqSysType = "twosides", constrIdx = seq(length(m)),
                       rmInf = TRUE) {
  constrType <- match.arg(constrType)
  if (constrType != "linear") {
    if (d == 1) {
      switch(constrType,
             boundedness = {
               linSys <- bounds2lineqSys(m, l, u, lineqSysType = lineqSysType, rmInf = rmInf)
             }, monotonicity = {
                 if (length(l) == 1) l <- c(-Inf, rep(l, m-1))
                 if (length(u) == 1) u <- c(Inf, rep(u, m-1))
                 A <- diag(m)
                 if (m == 2) {
                   A[2, 1] <- -1
                 } else {
                   diag(A[-1, -ncol(A)]) <- -1
                 }
                 # diag(A[-1, -ncol(A)]) <- -1
                 linSys <- bounds2lineqSys(nrow(A), l, u, A, lineqSysType, rmInf)
             }, convexity = {
                 if (length(l) == 1) l <- c(-rep(Inf, 2), rep(l, m-2))
                 if (length(u) == 1) u <- c(rep(Inf, 2), rep(u, m-2))
                 A <- diag(m)
                 diag(A[-seq(2), -c(ncol(A)-1,ncol(A))]) <- 1
                 diag(A[-seq(2), -c(1,ncol(A))]) <- -2
                 linSys <- bounds2lineqSys(nrow(A), l, u, A, lineqSysType, rmInf)
             }
      )
    } else {
      # For the of case convexity for d >= 2, constrType == "convexity" is a
      # weak version of convexity
      temp <- sapply(m, lineqGPSys, constrType, l, u,
                     lineqSysType = "twosides", rmInf = FALSE)
      Abase <- temp[1, ]
      lbase <- temp[2, ]
      ubase <- temp[3, ]
      Idiag <- lapply(as.list(m), diag)
      ones <- lapply(Idiag, diag)
      
      Anames <- matrix(paste("Idiag[[",seq(d),"]]", sep = ""),
                       d, d, byrow = TRUE)
      diag(Anames) <- paste("Abase[[",seq(d),"]]", sep = "")
      lbnames <- matrix(paste("ones[[",seq(d),"]]", sep = ""),
                        d, d, byrow = TRUE)
      diag(lbnames) <- paste("lbase[[",seq(d),"]]", sep = "")
      ubnames <- matrix(paste("ones[[",seq(d),"]]", sep = ""),
                        d, d, byrow = TRUE)
      diag(ubnames) <- paste("ubase[[",seq(d),"]]", sep = "")
      A <- l <- u <- c()
      for (k in constrIdx) {
        Atemp <- eval(parse(text = paste(Anames[k, ], collapse = " %x% ")))
        A <- rbind(A, Atemp)
        lbtemp <- eval(parse(text = paste(lbnames[k, ], collapse = " %x% ")))
        l <- c(l, lbtemp)
        ubtemp <- eval(parse(text = paste(ubnames[k, ], collapse = " %x% ")))
        u <- c(u, ubtemp)
      }
    linSys <- bounds2lineqSys(nrow(A), l, u, A, lineqSysType, rmInf)
    }
  } else {
    linSys <- bounds2lineqSys(nrow(A), l, u, A, lineqSysType, rmInf)
  }
  return(linSys)
}

#' @title Creation Method for the \code{"lineqGP"} S3 Class
#' @description Creation method for the \code{"lineqGP"} S3 class.
#' @param x a vector or matrix with the input data. The dimensions should be indexed by columns.
#' @param y a vector with the output data.
#' @param constrType a character string corresponding to the type of the inequality constraint.
#' Options: "boundedness", "monotonicity", "convexity", "linear";
#' Multiple constraints can be also defined, e.g. \code{constrType = c("boundedness", "monotonicity")}.
#' @return A list with the following elements.
#' \item{x,y,constrType}{see \bold{Arguments}.}
#' \item{d}{a number corresponding to the input dimension.}
#' \item{constrIdx}{for d > 1, a logical vector with the indices of active constrained dimensions.}
#' \item{localParam}{a list with specific parameters required for \code{"lineqGP"} models:
#' \code{m} (number of basis functions), \code{sampler}, and \code{samplingParams}.
#' See \code{\link{simulate.lineqGP}}.}
#' \item{kernParam}{a list with the kernel parameters: \code{par} (kernel parameters), \code{type}, \code{nugget}.
#' See \code{\link{kernCompute}}}
#' \item{bounds}{the limit values if \code{constrType = "boundedness"}.}
#' \item{(Lambda,lb,ub)}{the linear system of inequalities if \code{constrType = "linear"}.}
#'
#' @seealso \code{\link{augment.lineqGP}}, \code{\link{predict.lineqGP}}, \code{\link{simulate.lineqGP}}
#' @author A. F. Lopez-Lopera.
#' @references Lopez-Lopera, A. F., Bachoc, F., Durrande, N., and Roustant, O. (2017),
#' "Finite-dimensional Gaussian approximation with linear inequality constraints".
#' \emph{ArXiv e-prints}
#' \href{https://arxiv.org/abs/1710.07453}{[link]}
#'
#' @examples
#' # creating the model
#' sigfun <- function(x) return(1/(1+exp(-7*(x-0.5))))
#' x <- seq(0, 1, length = 5)
#' y <- sigfun(x)
#' model <- create(class = "lineqGP", x, y, constrType = "monotonicity")
#' model
#'
#' @method create lineqGP
#' @export
create.lineqGP <- function(x, y, constrType) {
  # changing the data as matrices
  if (!is.matrix(x))
    x <- as.matrix(x)
  if (!is.matrix(y) || ncol(y) != 1)
    y <- matrix(y)
  d <- ncol(x) # dimension of the input space

  # creating some lists for the model
  kernParam <- list(par = c(sigma2 = 1^2, theta = rep(0.2, d)),
                    type = "gaussian", nugget = 1e-7*sd(y))
  localParam <- list(m = rep(10*length(y), d), sampler = "HMC",
                     samplingParams = c(thinning = 1, burn.in = 1, scale = 0.1))
  # creating the full list for the model
  model <- list(x = x, y = y, constrType = constrType,
                d = d, constrIdx = seq(d), varnoise = 0,
                localParam = localParam, kernParam = kernParam)
  model$bounds <- c()
  if (any(constrType == "boundedness"))
    model$bounds <- rbind(model$bounds,
                          c(lower = min(y) - 0.05*abs(max(y) - min(y)),
                            upper = max(y) + 0.05*abs(max(y) - max(y))))
  if (any(constrType == "monotonicity"))
    model$bounds <- rbind(model$bounds, c(0, Inf))
  if (any(constrType == "convexity"))
    model$bounds <- rbind(model$bounds, c(0, Inf))

  if (any(constrType == "linear")) {
    # imposing positiveness constraints by default
    model$Lambda <- diag(prod(model$localParam$m))
    model$lb <- rep(0, nrow(model$Lambda))
    model$ub <- rep(Inf, nrow(model$Lambda))
  }
  return(model)
}

#' @title Augmenting Method for the \code{"lineqGP"} S3 Class
#' @description Augmenting method for the \code{"lineqGP"} S3 class.
#' @param x an object with class \code{lineqGP}.
#' @param ... further arguments passed to or from other methods.
#' @return An expanded \code{"lineqGP"} object with the following additional elements.
#' \item{Phi}{a matrix corresponding to the hat basis functions.
#' The basis functions are indexed by rows.}
#' \item{Gamma}{the covariance matrix of the Gassian vector \eqn{\boldsymbol{\xi}}{\xi}.}
#' \item{(Lambda,lb,ub)}{the linear system of inequalities.}
#' \item{...}{further parameters passed to or from other methods.}
#'
#' @details Some paramaters of the finite-dimensional GP with linear inequality
#' constraints are computed. Here, \eqn{\boldsymbol{\xi}}{\xi} is a centred Gaussian
#' vector with covariance \eqn{\boldsymbol{\Gamma}}{\Gamma}, s.t.
#' \eqn{\boldsymbol{\Phi} \boldsymbol{\xi} = \boldsymbol{y}}{\Phi \xi = y}
#' (interpolation constraints) and
#' \eqn{\boldsymbol{l} \leq \boldsymbol{\Lambda} \boldsymbol{\xi} \leq \boldsymbol{u}}{lb \le \Lambda \xi \le ub}
#' (inequality constraints).
#' @seealso \code{\link{create.lineqGP}}, \code{\link{predict.lineqGP}},
#'          \code{\link{simulate.lineqGP}}
#' @author A. F. Lopez-Lopera.
#'
#' @references Lopez-Lopera, A. F., Bachoc, F., Durrande, N., and Roustant, O. (2017),
#' "Finite-dimensional Gaussian approximation with linear inequality constraints".
#' \emph{ArXiv e-prints}
#' \href{https://arxiv.org/abs/1710.07453}{[link]}
#'
#' @examples
#' # creating the model
#' sigfun <- function(x) return(1/(1+exp(-7*(x-0.5))))
#' x <- seq(0, 1, length = 5)
#' y <- sigfun(x)
#' model <- create(class = "lineqGP", x, y, constrType = "monotonicity")
#'
#' # updating and expanding the model
#' model$localParam$m <- 30
#' model$kernParam$par <- c(1, 0.2)
#' model2 <- augment(model)
#' image(model2$Gamma, main = "covariance matrix")
#'
#' @importFrom broom augment
#' @export
augment.lineqGP <- function(x, ...) {
  model <- x
  if (!("nugget" %in% names(model$kernParam)))
    model$kernParam$nugget <- 0
  if ("bounds" %in% names(model)) {
    bounds <- model$bounds
  } else {
    bounds <- c(0, Inf)
  }
  # passing some terms from the model
  x <- model$x
  m <- model$localParam$m

  # computing the kernel matrix for the prior
  if (model$d == 1) {
    u <- matrix(seq(0, 1, by = 1/(m-1)), ncol = 1) # discretization vector
    Gamma <- kernCompute(u, u, model$kernParam$type, model$kernParam$par)
  } else if (model$d >= 2) {
    if (length(m) == 1){
      model$localParam$m <- rep(m, model$d)
      names(model$localParam$m) <- paste("m", seq(model$d), sep = "")
      mflag <- TRUE
    }
    m <- model$localParam$m
    u <- vector("list", model$d) # list with d sets of knots
    for (j in 1:model$d)
      u[[j]] <- matrix(seq(0, 1, by = 1/(m[j]-1)))

    AllGammas <- vector("list", model$d) # list with the d covariance matrices
    if (length(model$kernParam$par) == 2 && length(unique(m)) == 1) {
      uBase <- u[[1]]
      GammaBase <- kernCompute(uBase, uBase,
                               model$kernParam$type,
                               model$kernParam$par)
      for (j in 1:model$d) {
        AllGammas[[j]] <- GammaBase
      }
    } else {
      if (length(model$kernParam$par[-1]) != model$d) {
        model$kernParam$par <- c(model$kernParam$par[1],
                                 rep(model$kernParam$par[2], model$d))
        names(model$kernParam$par) <- c(names(model$kernParam$par)[1],
                                        paste(names(model$kernParam$par)[2],
                                              seq(model$d), sep = ""))
      }
      for (j in 1:model$d) {
        AllGammas[[j]] <- kernCompute(u[[j]], u[[j]],
                                      model$kernParam$type,
                                      model$kernParam$par[c(1, j+1)])
      }
    }
    expr <- paste("AllGammas[[", seq(model$d), "]]", collapse = " %x% ")
    Gamma <- eval(parse(text = expr))
  }
  model$Gamma <- Gamma# + model$kernParam$nugget*diag(nrow(Gamma))
  model$Phi <- basisCompute.lineqGP(x, u, ncol(x))

  # precomputing the linear system for the QP solver and MCMC samplers
  mt <- nrow(model$Gamma)
  M <- g <- numeric()
  Lambda <- lb <- ub <- numeric()
  mvec <- rep(0, length(model$constrType))
  for (i in seq(length(model$constrType))) {
    if (model$constrType[i] == "linear") {
      if (!("Lambda" %in% names(model)))
        stop('matrix Lambda is not defined')
      lsys <- lineqGPSys(nrow(model$Lambda), model$constrType[i], model$lb, model$ub,
                         model$Lambda, lineqSysType = "oneside")
      lsys2 <- lineqGPSys(nrow(model$Lambda), model$constrType[i], model$lb, model$ub,
                          model$Lambda, rmInf = FALSE)
    } else {
      bounds <- matrix(bounds, ncol = 2)
      lsys <- lineqGPSys(m, model$constrType[i], bounds[i,1], bounds[i,2],
                         d = model$d, constrIdx = model$constrIdx,
                         lineqSysType = "oneside")
      lsys2 <- lineqGPSys(m, model$constrType[i], bounds[i,1], bounds[i,2],
                          d = model$d, constrIdx = model$constrIdx,
                          rmInf = FALSE)
    }
    # oneside linear structure for QP.solver: M = [Lambda,-Lambda] and g = [-lb,ub]
    M <- rbind(M, lsys$M)
    g <- rbind(g, -matrix(lsys$g))
    # twosides linear structure (Lambda, lb, ub) for MCMC samplers
    Lambda <- rbind(Lambda, lsys2$A)
    lb <- c(lb, lsys2$l)
    ub <- c(ub, lsys2$u)
    # extra term required for HMC sampler
    mvec[i] <- nrow(lsys2$A)
  }
  # adding the parameters to the model structure
  model$Lambda <- Lambda
  model$lb <- lb
  model$ub <- ub
  model$lineqSys$M <- M # for QP solve
  model$lineqSys$g <- g # for QP solve
  model$localParam$mvec <- mvec # for HMC sampler
  model$u <- u
  return(model)
}

#' @title Prediction Method for the \code{"lineqGP"} S3 Class
#' @description Prediction method for the \code{"lineqGP"} S3 class.
#' @param object an object with class \code{"lineqGP"}.
#' @param xtest a vector (or matrix) with the test input design.
#' @param ... further arguments passed to or from other methods.
#' @return A \code{"lineqGP"} object with the following elements.
#' \item{Lambda}{a matrix corresponding to the linear set of inequality constraints.}
#' \item{lb}{the lower bound vector of the inequalities constraints.}
#' \item{ub}{the upper bound vector of the inequalities constraints.}
#' \item{Phi.test}{a matrix corresponding to the hat basis functions evaluated
#' at \code{xtest}. The basis functions are indexed by rows.}
#' \item{mu}{the unconstrained GP mean predictor.}
#' \item{Sigma}{the unconstrained GP prediction conditional covariance matrix.}
#' \item{xi.map}{the GP maximum a posteriori (MAP) predictor given the inequality constraints.}
#'
#' @details The posterior paramaters of the finite-dimensional GP with linear inequality
#' constraints are computed. Here, \eqn{\boldsymbol{\xi}}{\xi} is a centred Gaussian
#' vector with covariance \eqn{\boldsymbol{\Gamma}}{\Gamma}, s.t.
#' \eqn{\boldsymbol{\Phi} \boldsymbol{\xi} = \boldsymbol{y}}{\Phi \xi = y}
#' (interpolation constraints) and
#' \eqn{\boldsymbol{l} \leq \boldsymbol{\Lambda} \boldsymbol{\xi} \leq \boldsymbol{u}}{lb \le \Lambda \xi \le ub}
#' (inequality constraints).
#' @seealso \code{\link{create.lineqGP}}, \code{\link{augment.lineqGP}},
#'          \code{\link{simulate.lineqGP}}
#' @author A. F. Lopez-Lopera.
#'
#' @references Lopez-Lopera, A. F., Bachoc, F., Durrande, N., and Roustant, O. (2017),
#' "Finite-dimensional Gaussian approximation with linear inequality constraints".
#' \emph{ArXiv e-prints}
#' \href{https://arxiv.org/abs/1710.07453}{[link]}
#'
#' @examples
#' # creating the model
#' sigfun <- function(x) return(1/(1+exp(-7*(x-0.5))))
#' x <- seq(0, 1, length = 5)
#' y <- sigfun(x)
#' model <- create(class = "lineqGP", x, y, constrType = "monotonicity")
#'
#' # updating and expanding the model
#' model$localParam$m <- 30
#' model$kernParam$par <- c(1, 0.2)
#' model <- augment(model)
#'
#' # predictions from the model
#' xtest <- seq(0, 1, length = 100)
#' ytest <- sigfun(xtest)
#' pred <- predict(model, xtest)
#' plot(xtest, ytest, type = "l", lty = 2, main = "Kriging predictions")
#' lines(xtest, pred$Phi.test %*% pred$mu, type = "l", col = "blue")
#' lines(xtest, pred$Phi.test %*% pred$xi.map, type = "l", col = "red")
#' legend("right", c("ytest", "mean", "mode"), lty = c(2,1,1),
#'        col = c("black", "blue", "red"))
#'
#' @importFrom quadprog solve.QP
#' @export
predict.lineqGP <- function(object, xtest, ...) {
  model <- augment(object)
  if (!is.matrix(xtest))
    xtest <- matrix(xtest, ncol = model$d)

  # passing some terms from the model
  pred <- list()
  class(pred) <- class(model)
  pred$Lambda <- model$Lambda
  pred$lb <- model$lb
  pred$ub <- model$ub

  # precomputing some terms
  pred$Phi.test <- basisCompute.lineqGP(xtest, model$u, d = model$d)

  # # computing the linear system for the QP solver
  # A <- rbind(model$Phi, model$lineqSys$M)
  # b <- rbind(matrix(model$y, ncol = 1), model$lineqSys$g)

  # computing the conditional mean vector and conditional covariance matrix
  # given the interpolation points
  GammaPhit <- model$Gamma %*% t(model$Phi)
  PhiGammaPhit <- model$Phi %*% GammaPhit 
  if (model$varnoise != 0)
    PhiGammaPhit <- PhiGammaPhit + model$varnoise*diag(nrow(model$Phi))
  invPhiGammaPhit <- chol2inv(chol(PhiGammaPhit))
  pred$mu <- GammaPhit %*% invPhiGammaPhit %*% model$y
  pred$Sigma <- model$Gamma - GammaPhit %*% invPhiGammaPhit %*% t(GammaPhit)

  # computing the conditional mode vector given the interpolation points
  # and the inequality constraints
  # invGamma <- chol2inv(chol(model$Gamma))
  # d <- matrix(0, nrow = nrow(model$Gamma))
  # pred$xi.map <- solve.QP(invGamma, d, t(A), b, nrow(model$x))$solution
  if (min(eigen(pred$Sigma, symmetric = TRUE)$values) <= 0) # numerical stability
    pred$Sigma <- pred$Sigma + model$kernParam$nugget*diag(nrow(pred$Sigma))
  invSigma <- chol2inv(chol(pred$Sigma))
  pred$xi.map <- solve.QP(invSigma, t(pred$mu) %*% invSigma,
                          t(model$lineqSys$M), model$lineqSys$g)$solution
  return(pred)
}

#' @title Simulation Method for the \code{"lineqGP"} S3 Class
#' @description Simulation method for the \code{"lineqGP"} S3 class.
#' @param object an object with class \code{"lineqGP"}.
#' @param nsim	the number of simulations.
#' @param seed see \code{\link{simulate}}.
#' @param xtest a vector (or matrix) with the test input design.
#' @param ... further arguments passed to or from other methods.
#' @return A \code{"lineqGP"} object with the following elements.
#' \item{x}{a vector (or matrix) with the training input design.}
#' \item{y}{the training output vector at \code{x}.}
#' \item{xtest}{a vector (or matrix) with the test input design.}
#' \item{Phi.test}{a matrix corresponding to the hat basis functions evaluated
#' at \code{xtest}. The basis functions are indexed by rows.}
#' \item{xi.sim}{the posterior sample-path of the finite-dimensional Gaussian vector.}
#' \item{ysim}{the posterior sample-path of the observed GP.
#' Note: \code{ysim = Phi.test \%*\% xi.sim}.}
#'
#' @details The posterior sample-path of the finite-dimensional GP with linear inequality
#' constraints are computed. Here, \eqn{\boldsymbol{\xi}}{\xi} is a centred Gaussian
#' vector with covariance \eqn{\boldsymbol{\Gamma}}{\Gamma}, s.t.
#' \eqn{\boldsymbol{\Phi} \boldsymbol{\xi} = \boldsymbol{y}}{\Phi \xi = y}
#' (interpolation constraints) and
#' \eqn{\boldsymbol{l} \leq \boldsymbol{\Lambda} \boldsymbol{\xi} \leq \boldsymbol{u}}{lb \le \Lambda \xi \le ub}
#' (inequality constraints).
#' @seealso \code{\link{create.lineqGP}}, \code{\link{augment.lineqGP}},
#'          \code{\link{predict.lineqGP}}
#' @author A. F. Lopez-Lopera.
#'
#' @references Lopez-Lopera, A. F., Bachoc, F., Durrande, N., and Roustant, O. (2017),
#' "Finite-dimensional Gaussian approximation with linear inequality constraints".
#' \emph{ArXiv e-prints}
#' \href{https://arxiv.org/abs/1710.07453}{[link]}
#'
#' @examples
#' # creating the model
#' sigfun <- function(x) return(1/(1+exp(-7*(x-0.5))))
#' x <- seq(0, 1, length = 5)
#' y <- sigfun(x)
#' model <- create(class = "lineqGP", x, y, constrType = "monotonicity")
#'
#' # updating and expanding the model
#' model$localParam$m <- 30
#' model$kernParam$par <- c(1, 0.2)
#' model <- augment(model)
#'
#' # sampling from the model
#' xtest <- seq(0, 1, length = 100)
#' ytest <- sigfun(xtest)
#' sim.model <- simulate(model, nsim = 50, seed = 1, xtest = xtest)
#' mu <- apply(sim.model$ysim, 1, mean)
#' qtls <- apply(sim.model$ysim, 1, quantile, probs = c(0.05, 0.95))
#' matplot(xtest, t(qtls), type = "l", lty = 1, col = "gray90",
#'         main = "Constrained Kriging model")
#' polygon(c(xtest, rev(xtest)), c(qtls[2,], rev(qtls[1,])), col = "gray90", border = NA)
#' lines(xtest, ytest, lty = 2)
#' lines(xtest, mu, type = "l", col = "darkgreen")
#' points(x, y, pch = 20)
#' legend("right", c("ytrain", "ytest", "mean", "confidence"), lty = c(NaN,2,1,NaN),
#'        pch = c(20,NaN,NaN,15), col = c("black", "black", "darkgreen", "gray80"))
#'
#' @importFrom stats simulate
#' @export
simulate.lineqGP <- function(object, nsim = 1, seed = NULL, xtest, ...) {
  model <- augment(object)
  pred <- predict(model, xtest)

  # computing the transformed conditional mode and covariance matrix given
  # the interpolation points and the inequality constraints
  eta.mu <- as.vector(pred$Lambda %*% pred$mu)
  eta.map <- as.vector(pred$Lambda %*% pred$xi.map)
  Sigma.eta <- pred$Lambda %*% pred$Sigma %*% t(pred$Lambda)
  if (min(eigen(Sigma.eta, symmetric = TRUE)$values) <= 0)
    Sigma.eta <- Sigma.eta + model$kernParam$nugget*diag(nrow(Sigma.eta))

  # listing control terms
  control <- as.list(unlist(model$localParam$samplingParam))
  control$mvec <- model$localParam$mvec # for HMC
  control$constrType <- model$constrType # for HMC

  # sampling from the truncated multinormal
  tmvPar <- list(mu = eta.mu, map = eta.map, Sigma = Sigma.eta, lb = pred$lb, ub = pred$ub)
  class(tmvPar) <- model$localParam$sampler
  set.seed(seed)
  eta <- tmvrnorm(tmvPar, nsim, control)

  # passing some terms to the simulated model
  simModel <- list()
  simModel$x <- model$x
  simModel$y <- model$y
  simModel$xtest <- xtest
  simModel$Phi.test <- pred$Phi.test
  simModel$xi.map <- pred$xi.map
  simModel$xi.sim <- qr.solve(pred$Lambda, eta)
  simModel$ysim <- pred$Phi.test %*% simModel$xi.sim
  class(simModel) <- class(model)
  return(simModel)
}