covMan: Creator Function for 'covMan' Objects
In kergp: Gaussian Process Laboratory

Description Usage Arguments Details Note Examples

Creator function for covMan objects representing a covariance kernel entered manually.

   covMan(kernel, hasGrad = FALSE, acceptMatrix = FALSE,
          inputs = paste("x", 1:d, sep = ""), 
          d = length(inputs), parNames,
          par = NULL, parLower = NULL, parUpper = NULL,
          label = "covMan", ...)

`kernel`	A (semi-)positive definite function. This must be an object of class `"function"` with formal arguments named `"x1"`, `"x2"` and `"par"`. The first two formal arguments are locations vectors or matrices. The third formal is for the vector θ of all covariance parameters. An analytical gradient can be computed and returned as an attribute of the result with name `"gradient"`. See Details.
`hasGrad`	Logical indicating whether the `kernel` function returns the gradient w.r.t. the vector of parameters as a `"gradient"` attribute of the result. See Details
`acceptMatrix`	Logical indicating whether `kernel` admits matrices as arguments. Default is `FALSE`. See Examples below.
`inputs`	Character vector giving the names of the inputs used as arguments of `kernel`. Optional if `d` is given.
`d`	Integer specifying the spatial dimension (equal to the number of inputs). Optional if `inputs` is given.
`parNames`	Vector of character strings containing the parameter names.
`par, parLower, parUpper`	Optional numeric vectors containing the parameter values, lower bounds and upper bounds.
`label`	Optional character string describing the kernel.
`...`	Not used at this stage.

The formals and the returned value of the kernel function must be in accordance with the value of acceptMatrix.

When acceptMatrix is FALSE, the formal arguments x1 and x2 of kernel are numeric vectors with length d. The returned result is a numeric vector of length 1. The attribute named "gradient" of the returned value (if provided in accordance with the value of hasGrad) must then be a numeric vector with length equal to the number of covariance parameters. It must contain the derivative of the kernel value K(x1, x2, θ) with respect to the parameter vector θ.
When acceptMatrix is TRUE, the formals x1 and x2 are matrices with d columns and with n1 and n2 rows. The result is then a covariance matrix with n1 rows and n2 columns. The gradient attribute (if provided in accordance with the value of hasGrad) must be a list with length equal to the number of covariance parameters. The list element l must contain a numeric matrix with dimension c(n_1, n_2) which is the derivative of the covariance matrix w.r.t. the covariance parameter θ[l].

The kernel function must be symmetric with respect to its first two arguments, and it must be positive definite, which is not checked. If the function returns an object with a "gradient" attribute but hasGrad was set to FALSE, the gradient will not be used in optimization.

The name of the class was motivated by earlier stages in the development.

myCovMan <- 
      covMan(
         kernel = function(x1, x2, par) { 
         htilde <- (x1 - x2) / par[1]
         SS2 <- sum(htilde^2)
         d2 <- exp(-SS2)
         kern <- par[2] * d2
         d1 <- 2 * kern * SS2 / par[1]            
         attr(kern, "gradient") <- c(theta = d1,  sigma2 = d2)
         return(kern)
      },
      hasGrad = TRUE,    
      d = 1,
      label = "myGauss",
      parLower = c(theta = 0.0, sigma2 = 0.0),
      parUpper = c(theta = Inf, sigma2 = Inf),
      parNames = c("theta", "sigma2"),
      par = c(NA, NA)
      )
      
# Let us now code the same kernel in C
kernCode <- "
       SEXP kern, dkern;
       int nprotect = 0, d;
       double SS2 = 0.0, d2, z, *rkern, *rdkern;

       d = LENGTH(x1);
       PROTECT(kern = allocVector(REALSXP, 1)); nprotect++;
       PROTECT(dkern = allocVector(REALSXP, 2)); nprotect++;
       rkern = REAL(kern);
       rdkern = REAL(dkern);

       for (int i = 0; i < d; i++) {
         z = ( REAL(x1)[i] - REAL(x2)[i] ) / REAL(par)[0];
         SS2 += z * z; 
       }

       d2 = exp(-SS2);
       rkern[0] = REAL(par)[1] * d2;
       rdkern[1] =  d2; 
       rdkern[0] =  2 * rkern[0] * SS2 / REAL(par)[0];

       SET_ATTR(kern, install(\"gradient\"), dkern);
       UNPROTECT(nprotect);
       return kern;
     "
     
myCovMan  

## "inline" the C function into an R function: much more efficient! 

## Not run: 
require(inline)
kernC <- cfunction(sig = signature(x1 = "numeric", x2 = "numeric",
                                   par = "numeric"),
                    body = kernCode)
myCovMan <- covMan(kernel = kernC, hasGrad = TRUE, d = 1,
                   parNames = c("theta", "sigma2"))
myCovMan

## End(Not run)

## A kernel admitting matricial arguments
myCov <- covMan(
    
    kernel = function(x1, x2, par) { 
      # x1 : matrix of size n1xd
      # x2 : matrix of size n2xd
            
      d <- ncol(x1)
            
      SS2 <- 0
      for (j in 1:d){
        Aj <- outer(x1[, j], x2[, j], "-")
        Aj2 <- Aj^2
        SS2 <- SS2 + Aj2 / par[j]^2
      }
      D2 <- exp(-SS2)
      kern <- par[d + 1] * D2
    },
    acceptMatrix = TRUE,
    d = 2,
    label = "myGauss",
    parLower = c(theta1 = 0.0, theta2 = 0.0, sigma2 = 0.0),
    parUpper = c(theta1 = Inf, theta2 = Inf, sigma2 = Inf),
    parNames = c("theta1", "theta2", "sigma2"),
    par = c(NA, NA, NA)

)
      
coef(myCov) <- c(0.5, 1, 4)
show(myCov)

## computing the covariance kernel between two points
X <- matrix(c(0, 0), ncol = 2)
Xnew <- matrix(c(0.5, 1), ncol = 2)
colnames(X) <- colnames(Xnew) <- inputNames(myCov)
covMat(myCov, X)            ## same points
covMat(myCov, X, Xnew)      ## two different points

# computing covariances between sets of given locations
X <- matrix(c(0, 0.5, 0.7, 0, 0.5, 1), ncol = 2)
t <- seq(0, 1, length.out = 3)
Xnew <- as.matrix(expand.grid(t, t))
colnames(X) <- colnames(Xnew) <- inputNames(myCov)
covMat(myCov, X)         ## covariance matrix
covMat(myCov, X, Xnew)   ## covariances between design and new data