R/kernel_Cubic.R

In CollinErickson/GauPro: Gaussian Process Fitting

#' Cubic Kernel R6 class
#'
#' @docType class
#' @importFrom R6 R6Class
#' @export
#' @useDynLib GauPro, .registration = TRUE
#' @importFrom Rcpp evalCpp
#' @importFrom stats optim
# @keywords data, kriging, Gaussian process, regression
#' @return Object of \code{\link{R6Class}} with methods for fitting GP model.
#' @format \code{\link{R6Class}} object.
#' @examples
#' k1 <- Cubic$new(beta=runif(6)-.5)
#' plot(k1)
#'
#' n <- 12
#' x <- matrix(seq(0,1,length.out = n), ncol=1)
#' y <- sin(2*pi*x) + rnorm(n,0,1e-1)
#' gp <- GauPro_kernel_model$new(X=x, Z=y, kernel=Cubic$new(1),
#'                               parallel=FALSE, restarts=0)
#' gp$predict(.454)
Cubic <- R6::R6Class(
  classname = "GauPro_kernel_Cubic",
  inherit = GauPro_kernel_beta,
  public = list(
    #' @description Calculate covariance between two points
    #' @param x vector.
    #' @param y vector, optional. If excluded, find correlation
    #' of x with itself.
    #' @param beta Correlation parameters.
    #' @param s2 Variance parameter.
    #' @param params parameters to use instead of beta and s2.
    k = function(x, y=NULL, beta=self$beta, s2=self$s2, params=NULL) {
      if (!is.null(params)) {
        # lenpar <- length(params)
        # beta <- params[1:(lenpar-1)]
        # logs2 <- params[lenpar]

        lenparams <- length(params)
        if (self$beta_est) {
          beta <- params[1:self$beta_length]
        } else {
          beta <- self$beta
        }
        if (self$s2_est) {
          logs2 <- params[lenparams]
        } else {
          logs2 <- self$logs2
        }

        s2 <- 10^logs2
      } else {
        if (is.null(beta)) {beta <- self$beta}
        if (is.null(s2)) {s2 <- self$s2}
      }
      theta <- 10^beta
      if (is.null(y)) {
        if (is.matrix(x)) {
          if (self$useC && Sys.info()[['sysname']] == "Windows") {
            val <- s2 * corr_cubic_matrix_symC(x, theta)
          } else {
            val <- outer(1:nrow(x), 1:nrow(x),
                         Vectorize(function(i,j){self$kone(x[i,],x[j,],
                                                           theta=theta, s2=s2)}))
          }
          return(val)
        } else {
          return(s2 * 1)
        }
      }
      if (is.matrix(x) & is.matrix(y)) {
        if (self$useC && Sys.info()[['sysname']] == "Windows") {
          s2 * corr_cubic_matrixC(x, y, theta)
        } else {
          outer(1:nrow(x), 1:nrow(y),
                Vectorize(function(i,j){self$kone(x[i,],y[j,],
                                                  theta=theta, s2=s2)}))
        }
      } else if (is.matrix(x) & !is.matrix(y)) {
        if (self$useC && Sys.info()[['sysname']] == "Windows") {
          s2 * corr_cubic_matrixvecC(x, y, theta)
        } else {
          apply(x, 1, function(xx) {self$kone(xx, y, theta=theta, s2=s2)})
        }
      } else if (is.matrix(y)) {
        if (self$useC && Sys.info()[['sysname']] == "Windows") {
          s2 * corr_cubic_matrixvecC(y, x, theta)
        } else {
          apply(y, 1, function(yy) {self$kone(yy, x, theta=theta, s2=s2)})
        }
      } else {
        self$kone(x, y, theta=theta, s2=s2)
      }
    },
    #' @description Find covariance of two points
    #' @param x vector
    #' @param y vector
    #' @param beta correlation parameters on log scale
    #' @param theta correlation parameters on regular scale
    #' @param s2 Variance parameter
    kone = function(x, y, beta, theta, s2) {
      if (missing(theta)) {theta <- 10^beta}
      h <- x-y
      d <- h/theta
      r <- ifelse(abs(d) <= 0.5,
                  1-6*d^2+6*abs(d)^3,
                  ifelse(abs(d) <= 1,
                         2*(1-abs(d))^3,
                         0))
      prod(r) * s2
    },
    #' @description Derivative of covariance with respect to parameters
    #' @param params Kernel parameters
    #' @param X matrix of points in rows
    #' @param C_nonug Covariance without nugget added to diagonal
    #' @param C Covariance with nugget
    #' @param nug Value of nugget
    dC_dparams = function(params=NULL, X, C_nonug, C, nug) {
      n <- nrow(X)
      # stop("cubic dC_dparams not implemented")

      lenparams <- length(params)
      if (lenparams > 0) {
        if (self$beta_est) {
          beta <- params[1:self$beta_length]
        } else {
          beta <- self$beta
        }
        if (self$s2_est) {
          logs2 <- params[lenparams]
        } else {
          logs2 <- self$logs2
        }
      } else {
        beta <- self$beta
        logs2 <- self$logs2
      }

      # lenparams <- length(params)
      # beta <- params[1:(lenparams - 1)]
      theta <- 10^beta
      log10 <- log(10)
      # logs2 <- params[lenparams]
      s2 <- 10 ^ logs2

      # if (is.null(params)) {params <- c(self$beta, self$logs2)}
      if (missing(C_nonug)) { # Assume C missing too, must have nug
        C_nonug <- self$k(x=X, params=params)
        C <- C_nonug + diag(nug*s2, nrow(C_nonug))
      }

      lenparams_D <- self$beta_length*self$beta_est + self$s2_est
      if (self$useC && Sys.info()[['sysname']] == "Windows") {
        dC_dparams <- kernel_cubic_dC(X, theta, C_nonug, self$s2_est,
                                      self$beta_est, lenparams_D, s2*nug, s2)
      } else {
        dC_dparams <- array(dim=c(lenparams_D, n, n), data = 0)

        # Deriv for logs2
        if (self$s2_est) {
          dC_dparams[lenparams_D,,] <- C * log10 # Deriv for logs2
        }

        # Deriv for beta
        if (self$beta_est) {
          for (i in seq(1, n-1, 1)) {
            for (j in seq(i+1, n, 1)) {
              h <- X[i,] - X[j,] #x-y
              d <- h/theta
              dabs <- abs(d)
              r <- ifelse(abs(d) <= 0.5,
                          1-6*d^2+6*abs(d)^3,
                          ifelse(abs(d) <= 1,
                                 2*(1-abs(d))^3,
                                 0))
              # prod(r)
              for (k in 1:length(beta)) {
                drk_ddk <- sign(d[k]) * ifelse(dabs[k] <= 0.5,
                                               -12*dabs[k]+18*dabs[k]^2,
                                               ifelse(dabs[k] <= 1,
                                                      -6*(1-dabs[k])^2, 0))
                dC_dparams[k,i,j] <- (s2 * log10 * (-h[k]) / theta[k] *
                                        drk_ddk * prod(r[-k]))
                dC_dparams[k,j,i] <- dC_dparams[k,i,j]
              }
            }
          }
          for (i in seq(1, n, 1)) { # Get diagonal set to zero
            for (k in 1:length(beta)) {
              dC_dparams[k,i,i] <- 0
            }
          }
        }
      }
      return(dC_dparams)
    },
    #' @description Derivative of covariance with respect to X
    #' @param XX matrix of points
    #' @param X matrix of points to take derivative with respect to
    #' @param theta Correlation parameters
    #' @param beta log of theta
    #' @param s2 Variance parameter
    dC_dx = function(XX, X, theta, beta=self$beta, s2=self$s2) {
      # stop("cubic dC_dparams not implemented")
      if (missing(theta)) {theta <- 10^beta}
      if (!is.matrix(XX)) {stop()}
      D <- ncol(XX)
      if (ncol(X) != D) {stop()}
      n <- nrow(X)
      nn <- nrow(XX)
      dC_dx <- array(NA, dim=c(nn, D, n))
      for (i in 1:nn) {
        for (k in 1:n) {
          h <- XX[i,] - X[k,]
          d <- h/theta
          r <- ifelse(abs(d) <= 0.5,
                      1-6*d^2+6*abs(d)^3,
                      ifelse(abs(d) <= 1,
                             2*(1-abs(d))^3,
                             0))
          # k <- prod(r) * s2

          dr_dd <- ifelse(abs(d) <= 0.5,
                          # 1-6*d^2+6*abs(d)^3,
                          -12*d + 18*d^2*sign(d),
                          ifelse(abs(d) <= 1,
                                 # 2*(1-abs(d))^3,
                                 -6*(1-abs(d))^2*sign(d),
                                 0))
          for (j in 1:D) {
            prodrj <- if (D==1) {1} else {prod(r[-j])}
            dC_dx[i, j, k] <- s2 * prodrj * dr_dd[j] / theta[j]
          }
        }
      }
      dC_dx
    },
    #' @description Print this object
    print = function() {
      cat('GauPro kernel: Cubic\n')
      cat('\tD    =', self$D, '\n')
      cat('\tbeta =', signif(self$beta, 3), '\n')
      cat('\ts2   =', self$s2, '\n')
    }
  )
)

#' @rdname Cubic
#' @export
#' @param beta Initial beta value
#' @param s2 Initial variance
#' @param D Number of input dimensions of data
#' @param beta_lower Lower bound for beta
#' @param beta_upper Upper bound for beta
#' @param beta_est Should beta be estimated?
#' @param s2_lower Lower bound for s2
#' @param s2_upper Upper bound for s2
#' @param s2_est Should s2 be estimated?
#' @param useC Should C code used? Much faster.
k_Cubic <- function(beta, s2=1, D,
                    beta_lower=-8, beta_upper=6, beta_est=TRUE,
                    s2_lower=1e-8, s2_upper=1e8, s2_est=TRUE,
                    useC=TRUE) {
  Cubic$new(
    beta=beta,
    s2=s2,
    D=D,
    beta_lower=beta_lower,
    beta_upper=beta_upper,
    beta_est=beta_est,
    s2_lower=s2_lower,
    s2_upper=s2_upper,
    s2_est=s2_est,
    useC=useC
  )
}