R/kmBounded.R

In maatouk/constrKriging: Kriging with inequality constraints

#' @title Kriging model with boundedness constraints
#' @param design 1-column matrix of the design of experiments
#' @param response a vector containing the output values given by the real function at the design points
#' @param basis.size a value represents the number of the basis functions (descritization of 1D input set)
#' @param covtype an optimal character string specifying the covariance function to be used ("gauss" and "matern3_2" choice)
#' @param coef.cov a value corresponding to the length theta hyper-parameters of covariance function
#' @param coef.var a value specifying the variance parameter
#' @param nugget an optimal value used as nugget effect to solve the numerical inverse matrix problem
#' @param lower lower bound constraint
#' @param upper upper bound constraint

#' @examples 
#' model = kmBounded1D(design=c(0.1, 0.3, 0.5, 0.9), response=c(7, -8, 9, 15), lower=-10, upper = 18, coef.cov=1)
kmBounded1D <- function(design, response, 
                        basis.size = dim(design)[1]+2+10, 
                        covtype = "matern5_2",
                        coef.cov = 0.5*(max(design)-min(design)), # "LOO"
                        coef.var = var(response),
                        lower = min(response)-(max(response)-min(response))*.1,
                        upper = max(response)+(max(response)-min(response))*.1,
                        nugget = 1e-7*sd(response)) {
  
  if (!is.matrix(design)) design=matrix(design, ncol=1)
  
  
  #   if (coef.cov=="LOO") {
  #     object = kmBounded1D(design, response, 
  #                                      basis.size = dim(design)[1]+2+10, 
  #                                      covtype,
  #                                      coef.cov=0.5*(max(design)-min(design)),
  #                                      coef.var,
  #                                      lower,
  #                                      upper,
  #                                      nugget)
  #       
  #       theta=coef.cov_LOO(object)
  #     
  #     model = NULL
  #     while(is.null(model)) # auto raise nugget if needed
  #       try(model <- kmBounded1D(design, response, 
  #                                    basis.size, 
  #                                    covtype,
  #                                    coef.cov = theta,
  #                                    coef.var,
  #                                    lower,
  #                                    upper,
  #                                    nugget))
  #     
  #     return(model)
  #   }
  
  n <- nrow(design) # number of design points
  N <- basis.size   # discretization size
  p <- (N + 1) - n  # degree of freedom
  u <- seq(0, 1, by = 1/N) # discretization vector
  
  sig <- sqrt(coef.var)
  theta <- coef.cov
  
  
  if(covtype=='gauss'){
    # Gaussian covariance kernel 
    k <- function(x, xp, sig, theta){
      sig^2*exp(-(x-xp)^2/(2*theta^2))
    }
  }else if (covtype=='matern3_2'){
    # Matern 3/2 covariance kernel
    k <- function(x, xp, sig, theta){
      sig^2*(1+(sqrt(3)*abs(x-xp)/theta))*exp(-sqrt(3)*abs(x-xp)/theta)
    }
  }
  else if(covtype=="matern5_2"){
    # Matern 5/2 covariance kernel
    k <- function(x, xp, sig, theta){
      sig^2*(1+sqrt(5)*(abs(x-xp))/theta+(5*(x-xp)^2)/(3*theta^2))*exp(-sqrt(5)*(abs(x-xp))/theta)
    }
  }
  
  else stop('covtype', covtype, 'not supported')
  
  
  # Basis functions  
  phi <- function(x){
    ifelse(x >= -1 & x <= 1, 1-abs(x), 0)
  }
  phi0 <- function(x, N){
    phi(x*N)
  }
  phiN <- function(x, N){
    phi((x-u[N+1])*N)
  }
  phii <- function(x, i, N){
    phi((x-u[i+1])*N)
  }
  
  A <- matrix(data = NA, ncol = N+1, nrow = n)
  A[, 1] = phi0(design[, 1], N)
  A[, N+1] = phiN(design[, 1], N)
  for(j in 2 : N){
    A[, j] = phii(design[, 1], j-1, N)
  }
  
  fctGamma=function(.theta){
    Gamma <- matrix(data = NA, nrow = N+1, ncol = N+1)
    for(i in 1 : (N+1)){
      for(j in 1 : (N+1)){
        Gamma[i, j] = k(u[i], u[j], sig, .theta)
      }
    }
    Gamma <- Gamma + nugget * diag(N+1)
    return(Gamma)
  }
  Gamma=fctGamma(theta)
  
  invGamma <- chol2inv(chol(Gamma))
  
  Amat2 <- diag(N+1)
  Amat1 <- rbind(A, Amat2)
  Amat <- rbind(Amat1, -Amat2)
  
  zetoil <- solve.QP(invGamma, dvec=rep(0, N+1), Amat=t(Amat), bvec=c(response, rep(lower, N+1), rep(-upper, N+1)), meq=n)$solution
  
  return(structure(
    list(zetoil=zetoil, phi0=phi0, phiN=phiN, phii=phii, Amat=Amat, Gamma=Gamma, A=A, D=D, 
         fctGamma=fctGamma, invGamma=invGamma,
         call=list(design=design, response=response, basis.size=basis.size, covtype=covtype,
                   coef.cov=coef.cov, coef.var=coef.var, nugget=nugget, lower=lower,
                   upper=upper)), class="kmBounded1D"))
}




#' @title Apply basis functions to given data in design space
#' @param model km* model
#' @param newdata data in design space
#' @method Phi1D kmBounded
Phi1D.kmBounded1D <- function(model, newdata){
  N <- model$call$basis.size  
  x <- newdata
  
  v <- matrix(NA, nrow = length(x), ncol = N+1)
  v[, 1] <- model$phi0(x, N = N)
  v[, N+1] <- model$phiN(x, N = N)
  for(j in 2 : N){
    v[, j] = model$phii(x, j-1, N = N)
  }
  
  return(v)
}