R/kernels.R

In kerndwd: Distance Weighted Discrimination (DWD) and Kernel Methods

## S4 object definitions and assigment/accessor functions for the slots.
##
## created  10.09.03 alexandros karatzoglou
## updated  23.08.05


setClass("kernel", representation("function", kpar = "list"))
setClass("kernelMatrix",
         representation("matrix"),
         prototype = structure(.Data = matrix()))

setClassUnion("listI", c("list", "numeric", "vector", "integer", "matrix"))
setClassUnion(
  "output",
  c(
    "matrix",
    "factor",
    "vector",
    "logical",
    "numeric",
    "list",
    "integer",
    "NULL"
  )
)
setClassUnion("input", c("matrix", "list"))
setClassUnion("kfunction", c("function", "character"))
setClassUnion("mpinput", c("matrix", "data.frame", "missing"))
setClassUnion("lpinput", c("list", "missing"))
setClassUnion("kpinput", c("kernelMatrix", "missing"))

if (!isGeneric("kpar")) {
  if (is.function("kpar"))
    fun <- kpar
  else
    fun <- function(object)
      standardGeneric("kpar")
  setGeneric("kpar", fun)
}
setGeneric("kpar<-", function(x, value)
  standardGeneric("kpar<-"))



setGeneric("as.kernelMatrix", function(x, center = FALSE)
  standardGeneric("as.kernelMatrix"))
setMethod("as.kernelMatrix", signature(x = "matrix"),
          function(x, center = FALSE)
          {
            if (center) {
              m <- dim(x)[1]
              x <- t(t(x - colSums(x) / m) -  rowSums(x) / m) + sum(x) / m ^ 2
            }
            
            return(new("kernelMatrix", .Data = x))
          })



## kernel functions
## Functions for computing a kernel value, matrix, matrix-vector
## product and quadratic form
##
## author : alexandros karatzoglou


## Define the kernel objects,
## functions with an additional slot for the kernel parameter list.
## kernel functions take two vector arguments and return a scalar (dot product)


rbfdot <- function(sigma = 1)
{
  rval <- function(x, y = NULL)
  {
    if (!is(x, "vector"))
      stop("x must be a vector")
    if (!is(y, "vector") && !is.null(y))
      stop("y must a vector")
    if (is(x, "vector") && is.null(y)) {
      return(1)
    }
    if (is(x, "vector") && is(y, "vector")) {
      if (!length(x) == length(y))
        stop("number of dimension must be the same on both data points")
      return(exp(sigma * (
        2 * crossprod(x, y) - crossprod(x) - crossprod(y)
      )))
      # sigma/2 or sigma ??
    }
  }
  return(new(
    "rbfkernel",
    .Data = rval,
    kpar = list(sigma = sigma)
  ))
}
setClass(
  "rbfkernel",
  prototype = structure(
    .Data = function() {
    },
    kpar = list()
  ),
  contains = c("kernel")
)

laplacedot <- function(sigma = 1)
{
  rval <- function(x, y = NULL)
  {
    if (!is(x, "vector"))
      stop("x must be a vector")
    if (!is(y, "vector") && !is.null(y))
      stop("y must a vector")
    if (is(x, "vector") && is.null(y)) {
      return(1)
    }
    if (is(x, "vector") && is(y, "vector")) {
      if (!length(x) == length(y))
        stop("number of dimension must be the same on both data points")
      return(exp(-sigma * sqrt(-(
        round(2 * crossprod(x, y) - crossprod(x) - crossprod(y), 9)
      ))))
    }
  }
  return(new(
    "laplacekernel",
    .Data = rval,
    kpar = list(sigma = sigma)
  ))
}

setClass(
  "laplacekernel",
  prototype = structure(
    .Data = function() {
    },
    kpar = list()
  ),
  contains = c("kernel")
)

besseldot <- function(sigma = 1,
                      order = 1,
                      degree = 1)
{
  rval <- function(x, y = NULL)
  {
    if (!is(x, "vector"))
      stop("x must be a vector")
    if (!is(y, "vector") && !is.null(y))
      stop("y must a vector")
    if (is(x, "vector") && is.null(y)) {
      return(1)
    }
    if (is(x, "vector") && is(y, "vector")) {
      if (!length(x) == length(y))
        stop("number of dimension must be the same on both data points")
      lim <- 1 / (gamma(order + 1) * 2 ^ (order))
      bkt <-
        sigma * sqrt(-(2 * crossprod(x, y) - crossprod(x) - crossprod(y)))
      if (bkt < 10e-5)
        res <- lim
      else
        res <- besselJ(bkt, order) * (bkt ^ (-order))
      return((res / lim) ^ degree)
    }
  }
  return(new(
    "besselkernel",
    .Data = rval,
    kpar = list(
      sigma = sigma ,
      order = order ,
      degree = degree
    )
  ))
}

setClass(
  "besselkernel",
  prototype = structure(
    .Data = function() {
    },
    kpar = list()
  ),
  contains = c("kernel")
)

anovadot <- function(sigma = 1, degree = 1)
{
  rval <- function(x, y = NULL)
  {
    if (!is(x, "vector"))
      stop("x must be a vector")
    if (!is(y, "vector") && !is.null(y))
      stop("y must a vector")
    if (is(x, "vector") && is.null(y)) {
      return(1)
    }
    if (is(x, "vector") && is(y, "vector")) {
      if (!length(x) == length(y))
        stop("number of dimension must be the same on both data points")
      
      res <- sum(exp(-sigma * (x - y) ^ 2))
      return((res) ^ degree)
    }
  }
  return(new(
    "anovakernel",
    .Data = rval,
    kpar = list(sigma = sigma , degree = degree)
  ))
}

setClass(
  "anovakernel",
  prototype = structure(
    .Data = function() {
    },
    kpar = list()
  ),
  contains = c("kernel")
)


splinedot <- function()
{
  rval <- function(x, y = NULL)
  {
    if (!is(x, "vector"))
      stop("x must be a vector")
    if (!is(y, "vector") && !is.null(y))
      stop("y must a vector")
    if (is(x, "vector") && is.null(y)) {
      return(1)
    }
    if (is(x, "vector") && is(y, "vector")) {
      if (!length(x) == length(y))
        stop("number of dimension must be the same on both data points")
      minv <- pmin(x, y)
      res <- 1 + x * y * (1 + minv) - ((x + y) / 2) * minv ^ 2 + (minv ^
                                                                    3) / 3
      fres <- prod(res)
      return(fres)
    }
  }
  return(new("splinekernel", .Data = rval, kpar = list()))
}

setClass(
  "splinekernel",
  prototype = structure(
    .Data = function() {
    },
    kpar = list()
  ),
  contains = c("kernel")
)



fourierdot <- function(sigma = 1)
{
  rval <- function(x, y = NULL)
  {
    if (!is(x, "vector"))
      stop("x must be a vector")
    if (!is(y, "vector") && !is.null(y))
      stop("y must a vector")
    if (is(x, "vector") && is.null(y)) {
      return(1)
    }
    if (is(x, "vector") && is(y, "vector")) {
      if (!length(x) == length(y))
        stop("number of dimension must be the same on both data points")
      res <- 	(1 - sigma ^ 2) / 2 * (1 - 2 * sigma * cos(x - y) + sigma ^ 2)
      fres <- prod(res)
      return(fres)
    }
  }
  return(new("fourierkernel", .Data = rval, kpar = list()))
}

setClass(
  "fourierkernel",
  prototype = structure(
    .Data = function() {
    },
    kpar = list(sigma = 1)
  ),
  contains = c("kernel")
)





tanhdot <- function(scale = 1, offset = 1)
{
  rval <- function(x, y = NULL)
  {
    if (!is(x, "vector"))
      stop("x must be a  vector")
    if (!is(y, "vector") && !is.null(y))
      stop("y must be a vector")
    if (is(x, "vector") && is.null(y)) {
      tanh(scale * crossprod(x) + offset)
    }
    if (is(x, "vector") && is(y, "vector")) {
      if (!length(x) == length(y))
        stop("number of dimension must be the same on both data points")
      tanh(scale * crossprod(x, y) + offset)
    }
  }
  return(new(
    "tanhkernel",
    .Data = rval,
    kpar = list(scale = scale, offset = offset)
  ))
}
setClass(
  "tanhkernel",
  prototype = structure(
    .Data = function() {
    },
    kpar = list()
  ),
  contains = c("kernel")
)

setClass(
  "polykernel",
  prototype = structure(
    .Data = function() {
    },
    kpar = list()
  ),
  contains = c("kernel")
)

polydot <- function(degree = 1,
                    scale = 1,
                    offset = 1)
{
  rval <- function(x, y = NULL)
  {
    if (!is(x, "vector"))
      stop("x must be a vector")
    if (!is(y, "vector") && !is.null(y))
      stop("y must be a vector")
    if (is(x, "vector") && is.null(y)) {
      (scale * crossprod(x) + offset) ^ degree
    }
    
    if (is(x, "vector") && is(y, "vector")) {
      if (!length(x) == length(y))
        stop("number of dimension must be the same on both data points")
      (scale * crossprod(x, y) + offset) ^ degree
    }
    
  }
  return(new(
    "polykernel",
    .Data = rval,
    kpar = list(
      degree = degree,
      scale = scale,
      offset = offset
    )
  ))
}

setClass(
  "vanillakernel",
  prototype = structure(
    .Data = function() {
    },
    kpar = list()
  ),
  contains = c("kernel")
)

vanilladot <- function()
{
  rval <- function(x, y = NULL)
  {
    if (!is(x, "vector"))
      stop("x must be a vector")
    if (!is(y, "vector") && !is.null(y))
      stop("y must be a vector")
    if (is(x, "vector") && is.null(y)) {
      crossprod(x)
    }
    
    if (is(x, "vector") && is(y, "vector")) {
      if (!length(x) == length(y))
        stop("number of dimension must be the same on both data points")
      crossprod(x, y)
    }
    
  }
  return(new("vanillakernel", .Data = rval, kpar = list()))
}


setMethod("show", signature(object = "kernel"),
          function(object)
          {
            switch(
              class(object),
              "rbfkernel" = cat(
                paste(
                  "Gaussian Radial Basis kernel function.",
                  "\n",
                  "Hyperparameter :" ,
                  "sigma = ",
                  kpar(object)$sigma,
                  "\n"
                )
              ),
              "laplacekernel" = cat(
                paste(
                  "Laplace kernel function.",
                  "\n",
                  "Hyperparameter :" ,
                  "sigma = ",
                  kpar(object)$sigma,
                  "\n"
                )
              ),
              "besselkernel" = cat(
                paste(
                  "Bessel kernel function.",
                  "\n",
                  "Hyperparameter :" ,
                  "sigma = ",
                  kpar(object)$sigma,
                  "order = ",
                  kpar(object)$order,
                  "degree = ",
                  kpar(object)$degree,
                  "\n"
                )
              ),
              "anovakernel" = cat(
                paste(
                  "Anova RBF kernel function.",
                  "\n",
                  "Hyperparameter :" ,
                  "sigma = ",
                  kpar(object)$sigma,
                  "degree = ",
                  kpar(object)$degree,
                  "\n"
                )
              ),
              "tanhkernel" = cat(
                paste(
                  "Hyperbolic Tangent kernel function.",
                  "\n",
                  "Hyperparameters :",
                  "scale = ",
                  kpar(object)$scale,
                  " offset = ",
                  kpar(object)$offset,
                  "\n"
                )
              ),
              "polykernel" = cat(
                paste(
                  "Polynomial kernel function.",
                  "\n",
                  "Hyperparameters :",
                  "degree = ",
                  kpar(object)$degree,
                  " scale = ",
                  kpar(object)$scale,
                  " offset = ",
                  kpar(object)$offset,
                  "\n"
                )
              ),
              "vanillakernel" = cat(paste(
                "Linear (vanilla) kernel function.", "\n"
              )),
              "splinekernel" = cat(paste("Spline kernel function.", "\n")),
            )
          })

## create accesor function as in "S4 Classses in 15 pages more or less", well..

if (!isGeneric("kpar")) {
  if (is.function(kpar))
    fun <- kpar
  else
    fun <- function(object)
      standardGeneric("kpar")
  setGeneric("kpar", fun)
}

setMethod("kpar", "kernel", function(object)
  object@kpar)




## Functions that return usefull kernel calculations (kernel matrix etc.)

## kernelMatrix function takes two or three arguments

kernelMatrix <- function(kernel, x, y = NULL)
{
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  if (!is(x, "matrix"))
    stop("x must be a matrix")
  if (!is(y, "matrix") && !is.null(y))
    stop("y must be a matrix")
  n <- nrow(x)
  res1 <- matrix(rep(0, n * n), ncol = n)
  if (is.null(y)) {
    for (i in 1:n) {
      for (j in i:n) {
        res1[i, j]  <- kernel(x[i, ], x[j, ])
      }
    }
    res1 <- res1 + t(res1)
    diag(res1) <- diag(res1) / 2
    
    
  }
  if (is(y, "matrix")) {
    m <- dim(y)[1]
    res1 <- matrix(0, dim(x)[1], dim(y)[1])
    for (i in 1:n) {
      for (j in 1:m) {
        res1[i, j] <- kernel(x[i, ], y[j, ])
      }
    }
  }
  
  return(as.kernelMatrix(res1))
}

setGeneric("kernelMatrix", function(kernel, x, y = NULL)
  standardGeneric("kernelMatrix"))



kernelMatrix.rbfkernel <- function(kernel, x, y = NULL)
{
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  if (!is(y, "matrix") &&
      !is.null(y))
    stop("y must be a matrix or a vector")
  sigma = kpar(kernel)$sigma
  n <- dim(x)[1]
  dota <- rowSums(x * x) / 2
  if (is(x, "matrix") && is.null(y)) {
    res <- crossprod(t(x))
    for (i in 1:n)
      res[i, ] <- exp(2 * sigma * (res[i, ] - dota - rep(dota[i], n)))
    return(as.kernelMatrix(res))
  }
  if (is(x, "matrix") && is(y, "matrix")) {
    if (!(dim(x)[2] == dim(y)[2]))
      stop("matrixes must have the same number of columns")
    m <- dim(y)[1]
    dotb <- rowSums(y * y) / 2
    res <- x %*% t(y)
    for (i in 1:m)
      res[, i] <- exp(2 * sigma * (res[, i] - dota - rep(dotb[i], n)))
    return(as.kernelMatrix(res))
  }
}
setMethod("kernelMatrix",
          signature(kernel = "rbfkernel"),
          kernelMatrix.rbfkernel)

kernelMatrix.laplacekernel <- function(kernel, x, y = NULL)
{
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  if (!is(y, "matrix") &&
      !is.null(y))
    stop("y must be a matrix or a vector")
  sigma = kpar(kernel)$sigma
  n <- dim(x)[1]
  dota <- rowSums(x * x) / 2
  if (is(x, "matrix") && is.null(y)) {
    res <- crossprod(t(x))
    for (i in 1:n)
      res[i, ] <-
        exp(-sigma * sqrt(round(-2 * (
          res[i, ] - dota - rep(dota[i], n)
        ), 9)))
    return(as.kernelMatrix(res))
  }
  if (is(x, "matrix") && is(y, "matrix")) {
    if (!(dim(x)[2] == dim(y)[2]))
      stop("matrixes must have the same number of columns")
    m <- dim(y)[1]
    dotb <- rowSums(y * y) / 2
    res <- x %*% t(y)
    for (i in 1:m)
      res[, i] <-
      exp(-sigma * sqrt(round(-2 * (
        res[, i] - dota - rep(dotb[i], n)
      ), 9)))
    return(as.kernelMatrix(res))
  }
}
setMethod("kernelMatrix",
          signature(kernel = "laplacekernel"),
          kernelMatrix.laplacekernel)

kernelMatrix.besselkernel <- function(kernel, x, y = NULL)
{
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  if (!is(y, "matrix") &&
      !is.null(y))
    stop("y must be a matrix or a vector")
  sigma = kpar(kernel)$sigma
  nu = kpar(kernel)$order
  ni = kpar(kernel)$degree
  n <- dim(x)[1]
  lim <- 1 / (gamma(nu + 1) * 2 ^ (nu))
  dota <- rowSums(x * x) / 2
  if (is(x, "matrix") && is.null(y)) {
    res <- crossprod(t(x))
    for (i in 1:n) {
      xx <- sigma * sqrt(round(-2 * (res[i, ] - dota - rep(dota[i], n)), 9))
      res[i, ] <- besselJ(xx, nu) * (xx ^ (-nu))
      res[i, which(xx < 10e-5)] <- lim
    }
    return(as.kernelMatrix((res / lim) ^ ni))
  }
  if (is(x, "matrix") && is(y, "matrix")) {
    if (!(dim(x)[2] == dim(y)[2]))
      stop("matrixes must have the same number of columns")
    m <- dim(y)[1]
    dotb <- rowSums(y * y) / 2
    res <- x %*% t(y)
    for (i in 1:m) {
      xx <- sigma * sqrt(round(-2 * (res[, i] - dota - rep(dotb[i], n)), 9))
      res[, i] <- besselJ(xx, nu) * (xx ^ (-nu))
      res[which(xx < 10e-5), i] <- lim
    }
    return(as.kernelMatrix((res / lim) ^ ni))
  }
}
setMethod("kernelMatrix",
          signature(kernel = "besselkernel"),
          kernelMatrix.besselkernel)


kernelMatrix.anovakernel <- function(kernel, x, y = NULL)
{
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  if (!is(y, "matrix") &&
      !is.null(y))
    stop("y must be a matrix or a vector")
  sigma = kpar(kernel)$sigma
  degree = kpar(kernel)$degree
  n <- dim(x)[1]
  if (is(x, "matrix") && is.null(y)) {
    a <- matrix(0,  dim(x)[2], n)
    res <- matrix(0, n , n)
    for (i in 1:n)
    {
      a[rep(TRUE, dim(x)[2]), rep(TRUE, n)] <- x[i, ]
      res[i, ] <- colSums(exp(-sigma * (a - t(x)) ^ 2)) ^ degree
    }
    return(as.kernelMatrix(res))
  }
  if (is(x, "matrix") && is(y, "matrix")) {
    if (!(dim(x)[2] == dim(y)[2]))
      stop("matrixes must have the same number of columns")
    
    m <- dim(y)[1]
    b <- matrix(0, dim(x)[2], m)
    res <- matrix(0, dim(x)[1], m)
    for (i in 1:n)
    {
      b[rep(TRUE, dim(x)[2]), rep(TRUE, m)] <- x[i, ]
      res[i, ] <- colSums(exp(-sigma * (b - t(y)) ^ 2)) ^ degree
    }
    return(as.kernelMatrix(res))
  }
}
setMethod("kernelMatrix",
          signature(kernel = "anovakernel"),
          kernelMatrix.anovakernel)


kernelMatrix.polykernel <- function(kernel, x, y = NULL)
{
  if (!is(y, "matrix") && !is.null(y))
    stop("y must be a matrix")
  scale = kpar(kernel)$scale
  offset = kpar(kernel)$offset
  degree = kpar(kernel)$degree
  if (is(x, "matrix") && is.null(y))
  {
    res <- (scale * crossprod(t(x)) + offset) ^ degree
    return(as.kernelMatrix(res))
  }
  if (is(x, "matrix") && is(y, "matrix")) {
    if (!(dim(x)[2] == dim(y)[2]))
      stop("matrixes must have the same number of columns")
    res <- (scale * crossprod(t(x), t(y)) + offset) ^ degree
    return(as.kernelMatrix(res))
  }
}
setMethod("kernelMatrix",
          signature(kernel = "polykernel"),
          kernelMatrix.polykernel)

kernelMatrix.vanilla <- function(kernel, x, y = NULL)
{
  if (!is(y, "matrix") && !is.null(y))
    stop("y must be a matrix")
  if (is(x, "matrix") && is.null(y)) {
    res <- crossprod(t(x))
    return(as.kernelMatrix(res))
  }
  if (is(x, "matrix") && is(y, "matrix")) {
    if (!(dim(x)[2] == dim(y)[2]))
      stop("matrixes must have the same number of columns")
    res <- crossprod(t(x), t(y))
    return(as.kernelMatrix(res))
  }
}
setMethod("kernelMatrix",
          signature(kernel = "vanillakernel"),
          kernelMatrix.vanilla)

kernelMatrix.tanhkernel <- function(kernel, x, y = NULL)
{
  if (!is(y, "matrix") && !is.null(y))
    stop("y must be a matrix")
  if (is(x, "matrix") && is.null(y)) {
    scale = kpar(kernel)$scale
    offset = kpar(kernel)$offset
    res <- tanh(scale * crossprod(t(x)) + offset)
    return(as.kernelMatrix(res))
  }
  if (is(x, "matrix") && is(y, "matrix")) {
    if (!(dim(x)[2] == dim(y)[2]))
      stop("matrixes must have the same number of columns")
    res <- tanh(scale * crossprod(t(x), t(y)) + offset)
    return(as.kernelMatrix(res))
  }
}
setMethod("kernelMatrix",
          signature(kernel = "tanhkernel"),
          kernelMatrix.tanhkernel)


kernelMatrix.splinekernel <- function(kernel, x, y = NULL)
{
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  if (!is(y, "matrix") &&
      !is.null(y))
    stop("y must be a matrix or a vector")
  sigma = kpar(kernel)$sigma
  degree = kpar(kernel)$degree
  n <- dim(x)[1]
  if (is(x, "matrix") && is.null(y)) {
    a <- matrix(0,  dim(x)[2], n)
    res <- matrix(0, n , n)
    x <- t(x)
    for (i in 1:n)
    {
      dr <- 	x + x[, i]
      dp <-   x * x[, i]
      dm <-   pmin(x, x[, i])
      res[i, ] <-
        apply((1 + dp + dp * dm - (dr / 2) * dm ^ 2 + (dm ^ 3) / 3), 2, prod)
    }
    return(as.kernelMatrix(res))
  }
  if (is(x, "matrix") && is(y, "matrix")) {
    if (!(dim(x)[2] == dim(y)[2]))
      stop("matrixes must have the same number of columns")
    
    m <- dim(y)[1]
    b <- matrix(0, dim(x)[2], m)
    res <- matrix(0, dim(x)[1], m)
    x <- t(x)
    y <- t(y)
    for (i in 1:n)
    {
      dr <- 	y + x[, i]
      dp <-   y * x[, i]
      dm <-   pmin(y, x[, i])
      res[i, ] <-
        apply((1 + dp + dp * dm - (dr / 2) * dm ^ 2 + (dm ^ 3) / 3), 2, prod)
    }
    return(as.kernelMatrix(res))
  }
}
setMethod("kernelMatrix",
          signature(kernel = "splinekernel"),
          kernelMatrix.splinekernel)


## kernelMult computes kernel matrix  - vector product
## function computing <x,x'> * z (<x,x'> %*% z)


kernelMult <- function(kernel,
                       x,
                       y = NULL,
                       z,
                       blocksize = 128)
{
  #  if(is.function(kernel)) ker <- deparse(substitute(kernel))
  #      kernel <- do.call(kernel, kpar)
  
  if (!is(x, "matrix"))
    stop("x must be a matrix")
  if (!is(y, "matrix") && !is.null(y))
    stop("y must be a matrix")
  if (!is(z, "matrix") &&
      !is(z, "vector"))
    stop("z must ba a matrix or a vector")
  n <- nrow(x)
  
  if (is.null(y))
  {
    ## check if z,x match
    z <- as.matrix(z)
    if (is.null(y) && !dim(z)[1] == n)
      stop("z columns/length do not match x columns")
    
    res1 <- matrix(rep(0, n * n), ncol = n)
    
    for (i in 1:n)
    {
      for (j in i:n)
      {
        res1[j, i] <- kernel(x[i, ], x[j, ])
      }
    }
    res1 <- res1 + t(res1)
    diag(res1) <- diag(res1) / 2
  }
  if (is(y, "matrix"))
  {
    m <- dim(y)[1]
    z <- as.matrix(z)
    
    if (!dim(z)[1] == m)
      stop("z has wrong dimension")
    res1 <- matrix(rep.int(0, m * n), ncol = m)
    for (i in 1:n)
    {
      for (j in 1:m)
      {
        res1[i, j] <- kernel(x[i, ], y[j, ])
      }
    }
  }
  return(res1 %*% z)
}

setGeneric("kernelMult", function(kernel,
                                  x,
                                  y = NULL,
                                  z,
                                  blocksize = 256)
  standardGeneric("kernelMult"))

kernelMult.character <-
  function(kernel,
           x,
           y = NULL,
           z,
           blocksize = 256)
  {
    return(x %*% z)
  }
setMethod(
  "kernelMult",
  signature(kernel = "character", x = "kernelMatrix"),
  kernelMult.character
)



kernelMult.rbfkernel <-
  function(kernel,
           x,
           y = NULL,
           z,
           blocksize = 256)
  {
    if (!is(y, "matrix") &&
        !is.null(y) && !is(y, "vector"))
      stop("y must be a matrix or a vector")
    if (!is(z, "matrix") &&
        !is(z, "vector"))
      stop("z must be a matrix or a vector")
    if (is(x, "vector"))
      x <- as.matrix(x)
    if (is(y, "vector"))
      y <- as.matrix(y)
    sigma <- kpar(kernel)$sigma
    n <- dim(x)[1]
    m <- dim(x)[2]
    nblocks <- floor(n / blocksize)
    lowerl <- 1
    upperl <- 0
    dota <- as.matrix(rowSums(x ^ 2))
    
    if (is.null(y))
    {
      z <- as.matrix(z)
      if (!dim(z)[1] == n)
        stop("z rows must equal x rows")
      
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      if (nblocks > 0)
      {
        dotab <- rep(1, blocksize) %*% t(dota)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          res[lowerl:upperl, ] <-
            exp(sigma * (2 * x[lowerl:upperl, ] %*% t(x) - dotab - dota[lowerl:upperl] %*%
                           t(rep.int(1, n)))) %*% z
          lowerl <- upperl + 1
        }
      }
      if (lowerl <= n)
        res[lowerl:n, ] <-
        exp(sigma * (
          2 * x[lowerl:n, ] %*% t(x) - rep.int(1, n + 1 - lowerl) %*% t(dota) - dota[lowerl:n] %*%
            t(rep.int(1, n))
        )) %*% z
      
    }
    if (is(y, "matrix"))
    {
      n2 <- dim(y)[1]
      z <- as.matrix(z)
      
      if (!dim(z)[1] == n2)
        stop("z length must equal y rows")
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      dotb <- as.matrix(rowSums(y * y))
      
      if (nblocks > 0)
      {
        dotbb <-  rep(1, blocksize) %*% t(dotb)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          res[lowerl:upperl, ] <-
            exp(sigma * (2 * x[lowerl:upperl, ] %*% t(y) - dotbb - dota[lowerl:upperl] %*%
                           t(rep.int(1, n2)))) %*% z
          lowerl <- upperl + 1
        }
      }
      if (lowerl <= n)
        res[lowerl:n, ] <-
        exp(sigma * (
          2 * x[lowerl:n, ] %*% t(y) - rep.int(1, n + 1 - lowerl) %*% t(dotb) - dota[lowerl:n] %*%
            t(rep.int(1, n2))
        )) %*% z
    }
    return(res)
  }
setMethod("kernelMult",
          signature(kernel = "rbfkernel"),
          kernelMult.rbfkernel)


kernelMult.laplacekernel <-
  function(kernel,
           x,
           y = NULL,
           z,
           blocksize = 256)
  {
    if (!is(y, "matrix") &&
        !is.null(y) && !is(y, "vector"))
      stop("y must be a matrix or a vector")
    if (!is(z, "matrix") &&
        !is(z, "vector"))
      stop("z must be a matrix or a vector")
    if (is(x, "vector"))
      x <- as.matrix(x)
    if (is(y, "vector"))
      y <- as.matrix(y)
    sigma <- kpar(kernel)$sigma
    n <- dim(x)[1]
    m <- dim(x)[2]
    nblocks <- floor(n / blocksize)
    lowerl <- 1
    upperl <- 0
    dota <- as.matrix(rowSums(x ^ 2))
    
    if (is.null(y))
    {
      z <- as.matrix(z)
      if (!dim(z)[1] == n)
        stop("z rows must equal x rows")
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      if (nblocks > 0)
      {
        dotab <- rep(1, blocksize) %*% t(dota)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          res[lowerl:upperl, ] <-
            exp(-sigma * sqrt(-round(
              2 * x[lowerl:upperl, ] %*% t(x) - dotab - dota[lowerl:upperl] %*% t(rep.int(1, n)), 9
            ))) %*% z
          lowerl <- upperl + 1
          
        }
      }
      if (lowerl <= n)
        res[lowerl:n, ] <-
        exp(-sigma * sqrt(-round(
          2 * x[lowerl:n, ] %*% t(x) - rep.int(1, n + 1 - lowerl) %*% t(dota) - dota[lowerl:n] %*%
            t(rep.int(1, n)),
          9
        ))) %*% z
      
    }
    if (is(y, "matrix"))
    {
      n2 <- dim(y)[1]
      z <- as.matrix(z)
      
      if (!dim(z)[1] == n2)
        stop("z length must equal y rows")
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      dotb <- as.matrix(rowSums(y * y))
      
      if (nblocks > 0)
      {
        dotbb <-  rep(1, blocksize) %*% t(dotb)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          res[lowerl:upperl, ] <-
            exp(-sigma * sqrt(-round(
              2 * x[lowerl:upperl, ] %*% t(y) - dotbb - dota[lowerl:upperl] %*% t(rep.int(1, n2)), 9
            ))) %*% z
          lowerl <- upperl + 1
        }
      }
      if (lowerl <= n)
        res[lowerl:n, ] <-
        exp(-sigma * sqrt(-round(
          2 * x[lowerl:n, ] %*% t(y) - rep.int(1, n + 1 - lowerl) %*% t(dotb) - dota[lowerl:n] %*%
            t(rep.int(1, n2)),
          9
        ))) %*% z
    }
    return(res)
  }
setMethod("kernelMult",
          signature(kernel = "laplacekernel"),
          kernelMult.laplacekernel)



kernelMult.besselkernel <-
  function(kernel,
           x,
           y = NULL,
           z,
           blocksize = 256)
  {
    if (!is(y, "matrix") &&
        !is.null(y) && !is(y, "vector"))
      stop("y must be a matrix")
    if (!is(z, "matrix") &&
        !is(z, "vector"))
      stop("z must be a matrix or a vector")
    if (is(x, "vector"))
      x <- as.matrix(x)
    if (is(y, "vector"))
      y <- as.matrix(y)
    sigma <- kpar(kernel)$sigma
    nu <- kpar(kernel)$order
    ni <- kpar(kernel)$degree
    n <- dim(x)[1]
    m <- dim(x)[2]
    nblocks <- floor(n / blocksize)
    lowerl <- 1
    upperl <- 0
    lim <- 1 / (gamma(nu + 1) * 2 ^ (nu))
    dota <- as.matrix(rowSums(x ^ 2))
    
    if (is.null(y))
    {
      z <- as.matrix(z)
      if (!dim(z)[1] == n)
        stop("z rows must equal x rows")
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      if (nblocks > 0)
      {
        dotab <- rep(1, blocksize) %*% t(dota)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          
          xx <-
            sigma * sqrt(-round(2 * x[lowerl:upperl, ] %*% t(x) - dotab - dota[lowerl:upperl] %*%
                                  t(rep.int(1, n)), 9))
          res1 <- besselJ(xx, nu) * (xx ^ (-nu))
          res1[which(xx < 10e-5)] <- lim
          
          res[lowerl:upperl, ] <- ((res1 / lim) ^ ni) %*% z
          lowerl <- upperl + 1
        }
      }
      if (lowerl <= n)
      {
        xx <-
          sigma * sqrt(-round(
            2 * x[lowerl:n, ] %*% t(x) - rep.int(1, n + 1 - lowerl) %*% t(dota) - dota[lowerl:n] %*%
              t(rep.int(1, n)),
            9
          ))
        res1 <- besselJ(xx, nu) * (xx ^ (-nu))
        res1[which(xx < 10e-5)] <- lim
        res[lowerl:n, ] <- ((res1 / lim) ^ ni) %*% z
      }
    }
    if (is(y, "matrix"))
    {
      n2 <- dim(y)[1]
      z <- as.matrix(z)
      
      if (!dim(z)[1] == n2)
        stop("z length must equal y rows")
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      dotb <- as.matrix(rowSums(y * y))
      
      if (nblocks > 0)
      {
        dotbb <-  rep(1, blocksize) %*% t(dotb)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          xx <-
            sigma * sqrt(-round(2 * x[lowerl:upperl, ] %*% t(y) - dotbb - dota[lowerl:upperl] %*%
                                  t(rep.int(1, n2)), 9))
          res1 <- besselJ(xx, nu) * (xx ^ (-nu))
          res1[which(xx < 10e-5)] <- lim
          
          res[lowerl:upperl, ] <- ((res1 / lim) ^ ni) %*% z
          lowerl <- upperl + 1
        }
      }
      if (lowerl <= n)
      {
        xx <-
          sigma * sqrt(-round(
            2 * x[lowerl:n, ] %*% t(y) - rep.int(1, n + 1 - lowerl) %*% t(dotb) - dota[lowerl:n] %*%
              t(rep.int(1, n2)),
            9
          ))
        res1 <- besselJ(xx, nu) * (xx ^ (-nu))
        res1[which(xx < 10e-5)] <- lim
        res[lowerl:n, ] <- ((res1 / lim) ^ ni) %*% z
      }
    }
    return(res)
  }
setMethod("kernelMult",
          signature(kernel = "besselkernel"),
          kernelMult.besselkernel)

kernelMult.anovakernel <-
  function(kernel,
           x,
           y = NULL,
           z,
           blocksize = 256)
  {
    if (!is(y, "matrix") &&
        !is.null(y) && !is(y, "vector"))
      stop("y must be a matrix or a vector")
    if (!is(z, "matrix") &&
        !is(z, "vector"))
      stop("z must be a matrix or a vector")
    if (is(x, "vector"))
      x <- as.matrix(x)
    if (is(y, "vector"))
      y <- as.matrix(y)
    sigma <- kpar(kernel)$sigma
    degree <- kpar(kernel)$degree
    n <- dim(x)[1]
    m <- dim(x)[2]
    nblocks <- floor(n / blocksize)
    lowerl <- 1
    upperl <- 0
    
    
    if (is.null(y))
    {
      z <- as.matrix(z)
      
      if (!dim(z)[1] == n)
        stop("z rows must equal x rows")
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      
      if (nblocks > 0)
      {
        a <- matrix(0, m, blocksize)
        re <- matrix(0, n, blocksize)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          for (j in 1:n)
          {
            a[rep(TRUE, m), rep(TRUE, blocksize)] <- x[j, ]
            re[j, ] <-
              colSums(exp(-sigma * (a - t(x[lowerl:upperl, ])) ^ 2)) ^ degree
          }
          res[lowerl:upperl, ] <- t(re) %*% z
          lowerl <- upperl + 1
          
        }
      }
      if (lowerl <= n) {
        a <- matrix(0, m, n - lowerl + 1)
        re <- matrix(0, n, n - lowerl + 1)
        for (j in 1:n)
        {
          a[rep(TRUE, m), rep(TRUE, n - lowerl + 1)] <- x[j, ]
          re[j, ] <-
            colSums(exp(-sigma * (a - t(x[lowerl:n, , drop = FALSE])) ^ 2)) ^ degree
        }
        res[lowerl:n, ] <- t(re) %*% z
      }
    }
    if (is(y, "matrix"))
    {
      n2 <- dim(y)[1]
      nblocks <- floor(n2 / blocksize)
      
      z <- as.matrix(z)
      
      if (!dim(z)[1] == n2)
        stop("z length must equal y rows")
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      
      if (nblocks > 0)
      {
        b <- matrix(0, m, blocksize)
        re <- matrix(0, n, blocksize)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          for (j in 1:n)
          {
            b[rep(TRUE, dim(x)[2]), rep(TRUE, blocksize)] <- x[j, ]
            re[j, ] <-
              colSums(exp(-sigma * (b - t(y[lowerl:upperl, ])) ^ 2) ^ degree)
          }
          res[, 1] <- res[, 1] + re %*% z[lowerl:upperl, ]
          lowerl <- upperl + 1
        }
      }
      if (lowerl <= n)
      {
        b <- matrix(0, dim(x)[2], n2 - lowerl + 1)
        re <- matrix(0, n, n2 - lowerl + 1)
        for (i in 1:n)
        {
          b[rep(TRUE, dim(x)[2]), rep(TRUE, n2 - lowerl + 1)] <- x[i, ]
          re[i, ] <-
            colSums(exp(-sigma * (b - t(y[lowerl:n2, , drop = FALSE])) ^ 2) ^ degree)
        }
        
        res[, 1] <- res[, 1] + re %*% z[lowerl:n2]
      }
    }
    return(res)
  }
setMethod("kernelMult",
          signature(kernel = "anovakernel"),
          kernelMult.anovakernel)



kernelMult.splinekernel <-
  function(kernel,
           x,
           y = NULL,
           z,
           blocksize = 256)
  {
    if (!is(y, "matrix") &&
        !is.null(y) && !is(y, "vector"))
      stop("y must be a matrix or a vector")
    if (!is(z, "matrix") &&
        !is(z, "vector"))
      stop("z must be a matrix or a vector")
    n <- dim(x)[1]
    m <- dim(x)[2]
    if (is(x, "vector"))
      x <- as.matrix(x)
    if (is(y, "vector"))
      y <- as.matrix(y)
    nblocks <- floor(n / blocksize)
    lowerl <- 1
    upperl <- 0
    
    if (is.null(y))
    {
      z <- as.matrix(z)
      
      if (!dim(z)[1] == n)
        stop("z rows must equal x rows")
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      x <- t(x)
      if (nblocks > 0)
      {
        re <- matrix(0, dim(z)[1], blocksize)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          
          for (j in lowerl:upperl)
          {
            dr <-  x + x[, j]
            dp <-  x * x[, j]
            dm <-  pmin(x, x[, j])
            re[, j - (i - 1) * blocksize] <-
              apply((1 + dp + dp * dm - (dr / 2) * dm ^ 2 + (dm ^ 3) / 3), 2, prod)
          }
          res[lowerl:upperl, ] <- crossprod(re, z)
          lowerl <- upperl + 1
        }
      }
      if (lowerl <= n) {
        a <- matrix(0, m, n - lowerl + 1)
        re <- matrix(0, dim(z)[1], n - lowerl + 1)
        for (j in lowerl:(n - lowerl + 1))
        {
          dr <- x + x[, j]
          dp <- x * x[, j]
          dm <-  pmin(x, x[, j])
          re[, j - nblocks * blocksize] <-
            apply((1 + dp + dp * dm - (dr / 2) * dm ^ 2 + (dm ^ 3) / 3), 2, prod)
        }
        res[lowerl:n, ] <- crossprod(re, z)
      }
    }
    if (is(y, "matrix"))
    {
      n2 <- dim(y)[1]
      nblocks <- floor(n2 / blocksize)
      z <- as.matrix(z)
      
      if (!dim(z)[1] == n2)
        stop("z length must equal y rows")
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      x <- t(x)
      y <- t(y)
      if (nblocks > 0)
      {
        re <- matrix(0, dim(z)[1], blocksize)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          
          for (j in lowerl:upperl)
          {
            dr <- y + x[, j]
            dp <- y * x[, j]
            dm <- pmin(y, x[, j])
            re[, j - (i - 1) * blocksize] <-
              apply((1 + dp + dp * dm - (dr / 2) * dm ^ 2 + (dm ^ 3) / 3), 2, prod)
          }
          res[lowerl:upperl] <- crossprod(re, z)
          lowerl <- upperl + 1
        }
      }
      if (lowerl <= n)
      {
        b <- matrix(0, dim(x)[2], n - lowerl + 1)
        re <- matrix(0, dim(z)[1], n - lowerl + 1)
        for (j in lowerl:(n - lowerl + 1))
        {
          dr <- y + x[, j]
          dp <- y * x[, j]
          dm <-  pmin(y, x[, j])
          re[, j - nblocks * blocksize] <-
            apply((1 + dp + dp * dm - (dr / 2) * dm ^ 2 + (dm ^ 3) / 3), 2, prod)
        }
        res[lowerl:n] <-  crossprod(re, z)
      }
    }
    return(res)
  }
setMethod("kernelMult",
          signature(kernel = "splinekernel"),
          kernelMult.splinekernel)


kernelMult.polykernel <-
  function(kernel,
           x,
           y = NULL,
           z,
           blocksize = 256)
  {
    if (!is(y, "matrix") &&
        !is.null(y) && !is(y, "vector"))
      stop("y must be a matrix")
    if (!is(z, "matrix") &&
        !is(z, "vector"))
      stop("z must be a matrix or a vector")
    if (is(x, "vector"))
      x <- as.matrix(x)
    if (is(y, "vector"))
      y <- as.matrix(y)
    degree <- kpar(kernel)$degree
    scale <- kpar(kernel)$scale
    offset <- kpar(kernel)$offset
    n <- dim(x)[1]
    m <- dim(x)[2]
    nblocks <- floor(n / blocksize)
    lowerl <- 1
    upperl <- 0
    if (is.null(y))
    {
      z <- as.matrix(z)
      
      if (!dim(z)[1] == n)
        stop("z rows must equal x rows")
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      if (nblocks > 0)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          res[lowerl:upperl, ] <-
            ((scale * x[lowerl:upperl, ] %*% t(x) + offset) ^ degree) %*% z
          lowerl <- upperl + 1
        }
      if (lowerl <= n)
        res[lowerl:n, ] <-
        ((scale * x[lowerl:n, ] %*% t(x) + offset) ^ degree) %*% z
    }
    if (is(y, "matrix"))
    {
      n2 <- dim(y)[1]
      z <- as.matrix(z)
      
      if (!dim(z)[1] == n2)
        stop("z length must equal y rows")
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      
      if (nblocks > 0)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          res[lowerl:upperl, ] <-
            ((scale * x[lowerl:upperl, ] %*% t(y) + offset) ^ degree) %*% z
          lowerl <- upperl + 1
        }
      if (lowerl <= n)
        res[lowerl:n, ] <-
        ((scale * x[lowerl:n, ] %*% t(y) + offset) ^ degree) %*% z
    }
    return(res)
  }
setMethod("kernelMult",
          signature(kernel = "polykernel"),
          kernelMult.polykernel)


kernelMult.tanhkernel <-
  function(kernel,
           x,
           y = NULL,
           z,
           blocksize = 256)
  {
    if (!is(y, "matrix") &&
        !is.null(y) && !is(y, "vector"))
      stop("y must be a matrix or a vector")
    if (!is(z, "matrix") &&
        !is(z, "vector"))
      stop("z must be a matrix or a vector")
    if (is(x, "vector"))
      x <- as.matrix(x)
    if (is(y, "vector"))
      y <- as.matrix(y)
    scale <- kpar(kernel)$scale
    offset <- kpar(kernel)$offset
    n <- dim(x)[1]
    m <- dim(x)[2]
    nblocks <- floor(n / blocksize)
    lowerl <- 1
    upperl <- 0
    
    
    if (is.null(y))
    {
      z <- as.matrix(z)
      
      if (!dim(z)[1] == n)
        stop("z rows must equal x rows")
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      if (nblocks > 0)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          res[lowerl:upperl, ] <-
            tanh(scale * x[lowerl:upperl, ] %*% t(x) + offset) %*% z
          lowerl <- upperl + 1
        }
      if (lowerl <= n)
        res[lowerl:n, ] <-
        tanh(scale * x[lowerl:n, ] %*% t(x) + offset) %*% z
    }
    if (is(y, "matrix"))
    {
      n2 <- dim(y)[1]
      z <- as.matrix(z)
      
      if (!dim(z)[1] == n2)
        stop("z length must equal y rows")
      res <- matrix(rep(0, dim(z)[2] * n), ncol = dim(z)[2])
      
      if (nblocks > 0)
        for (i in 1:nblocks)
        {
          upperl = upperl + blocksize
          res[lowerl:upperl, ] <-
            tanh(scale * x[lowerl:upperl, ] %*% t(y) + offset) %*% z
          lowerl <- upperl + 1
        }
      if (lowerl <= n)
        res[lowerl:n, ] <-
        tanh(scale * x[lowerl:n, ] %*% t(y) + offset) %*% z
    }
    return(res)
  }
setMethod("kernelMult",
          signature(kernel = "tanhkernel"),
          kernelMult.tanhkernel)


kernelMult.vanillakernel <-
  function(kernel,
           x,
           y = NULL,
           z,
           blocksize = 256)
  {
    if (!is(y, "matrix") &&
        !is.null(y) && !is(y, "vector"))
      stop("y must be a matrix or vector")
    if (!is(z, "matrix") &&
        !is(z, "vector"))
      stop("z must be a matrix or a vector")
    n <- dim(x)[1]
    m <- dim(x)[2]
    if (is(x, "vector"))
      x <- as.matrix(x)
    if (is(y, "vector"))
      y <- as.matrix(y)
    if (is.null(y))
    {
      z <- as.matrix(z)
      if (!dim(z)[1] == n)
        stop("z rows must equal x rows")
      res <- t(crossprod(crossprod(x, z), t(x)))
    }
    if (is(y, "matrix"))
    {
      n2 <- dim(y)[1]
      z <- as.matrix(z)
      
      if (!dim(z)[1] == n2)
        stop("z length must equal y rows")
      res <- t(crossprod(crossprod(y, z), t(x)))
    }
    return(res)
  }
setMethod("kernelMult",
          signature(kernel = "vanillakernel"),
          kernelMult.vanillakernel)


## kernelPol return the quadratic form of a kernel matrix
## kernelPol returns the scalar product of x y componentwise with polarities
## of z and k

kernelPol <- function(kernel, x, y = NULL, z, k = NULL)
{
  if (!is(x, "matrix"))
    stop("x must be a matrix")
  if (!is(y, "matrix") && !is.null(y))
    stop("y must be a matrix")
  if (!is(z, "matrix") &&
      !is(z, "vector"))
    stop("z must ba a matrix or a vector")
  n <- nrow(x)
  z <- as.matrix(z)
  
  
  if (!dim(z)[1] == n)
    stop("z must have the length equal to x colums")
  res1 <- matrix(rep(0, n * n), ncol = n)
  if (is.null(y))
  {
    for (i in 1:n)
    {
      for (j in i:n)
      {
        res1[i, j] <- kernel(x[i, ], x[j, ]) * z[j] * z[i]
      }
    }
    res1 <- res1 + t(res1)
    diag(res1) <- diag(res1) / 2
  }
  if (is(x, "matrix") && is(y, "matrix")) {
    m <- dim(y)[1]
    if (is.null(k))
      stop("k not specified!")
    k <- as.matrix(k)
    if (!dim(x)[2] == dim(y)[2])
      stop("matrixes must have the same number of columns")
    if (!dim(z)[2] == dim(k)[2])
      stop("z and k vectors must have the same number of columns")
    if (!dim(x)[1] == dim(z)[1])
      stop("z and x must have the same number of rows")
    if (!dim(y)[1] == dim(k)[1])
      stop("y and k must have the same number of rows")
    res1 <- matrix(0, dim(x)[1], dim(y)[1])
    for (i in 1:n)
    {
      for (j in 1:m)
      {
        res1[i, j] <- kernel(x[i, ], y[j, ]) * z[i] * k[j]
      }
    }
  }
  return(res1)
}

setGeneric("kernelPol", function(kernel, x, y = NULL, z, k = NULL)
  standardGeneric("kernelPol"))


kernelPol.rbfkernel <- function(kernel, x, y = NULL, z, k = NULL)
{
  if (!is(y, "matrix") &&
      !is.null(y) &&
      !is(y, "vector"))
    stop("y must be a matrix a vector or NULL")
  if (!is(z, "matrix") &&
      !is(z, "vector"))
    stop("z must be a matrix or a vector")
  if (!is(k, "matrix") &&
      !is(k, "vector") &&
      !is.null(k))
    stop("k must be a matrix or a vector")
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  sigma <- kpar(kernel)$sigma
  n <- dim(x)[1]
  dota <- rowSums(x * x) / 2
  z <- as.matrix(z)
  if (!dim(z)[1] == n)
    stop("z must have the length equal to x colums")
  if (is.null(y))
  {
    if (is(z, "matrix") && !dim(z)[1] == n)
      stop("z must have size equal to x colums")
    res <- crossprod(t(x))
    for (i in 1:n)
      res[i, ] <-
      z[i, ] * (exp(2 * sigma * (res[i, ] - dota - rep(dota[i], n))) * z)
    return(res)
  }
  if (is(y, "matrix"))
  {
    if (is.null(k))
      stop("k not specified!")
    m <- dim(y)[1]
    k <- as.matrix(k)
    if (!dim(k)[1] == m)
      stop("k must have equal rows to y")
    if (!dim(x)[2] == dim(y)[2])
      stop("matrixes must have the same number of columns")
    dotb <- rowSums(y * y) / 2
    res <- x %*% t(y)
    for (i in 1:m)
      #2*sigma or sigma
      res[, i] <-
      k[i, ] * (exp(2 * sigma * (res[, i] - dota - rep(dotb[i], n))) * z)
    return(res)
  }
}
setMethod("kernelPol",
          signature(kernel = "rbfkernel"),
          kernelPol.rbfkernel)

kernelPol.laplacekernel <- function(kernel, x, y = NULL, z, k = NULL)
{
  if (!is(y, "matrix") &&
      !is.null(y) &&
      !is(y, "vector"))
    stop("y must be a matrix, vector or NULL")
  if (!is(z, "matrix") &&
      !is(z, "vector"))
    stop("z must be a matrix or a vector")
  if (!is(k, "matrix") &&
      !is(k, "vector") &&
      !is.null(k))
    stop("k must be a matrix or a vector")
  sigma <- kpar(kernel)$sigma
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  n <- dim(x)[1]
  dota <- rowSums(x * x) / 2
  z <- as.matrix(z)
  if (!dim(z)[1] == n)
    stop("z must have the length equal to x colums")
  if (is.null(y))
  {
    if (is(z, "matrix") && !dim(z)[1] == n)
      stop("z must have size equal to x colums")
    res <- crossprod(t(x))
    for (i in 1:n)
      res[i, ] <-
      z[i, ] * (exp(-sigma * sqrt(-round(
        2 * (res[i, ] - dota - rep(dota[i], n)), 9
      ))) * z)
    return(res)
  }
  if (is(y, "matrix"))
  {
    if (is.null(k))
      stop("k not specified!")
    m <- dim(y)[1]
    k <- as.matrix(k)
    if (!dim(k)[1] == m)
      stop("k must have equal rows to y")
    if (!dim(x)[2] == dim(y)[2])
      stop("matrixes must have the same number of columns")
    dotb <- rowSums(y * y) / 2
    res <- x %*% t(y)
    for (i in 1:m)
      #2*sigma or sigma
      res[, i] <-
      k[i, ] * (exp(-sigma * sqrt(-round(
        2 * (res[, i] - dota - rep(dotb[i], n)), 9
      ))) * z)
    return(res)
  }
}
setMethod("kernelPol",
          signature(kernel = "laplacekernel"),
          kernelPol.laplacekernel)


kernelPol.besselkernel <- function(kernel, x, y = NULL, z, k = NULL)
{
  if (!is(y, "matrix") &&
      !is.null(y) && !is(y, "vector"))
    stop("y must be a matrix or NULL")
  if (!is(z, "matrix") &&
      !is(z, "vector"))
    stop("z must be a matrix or a vector")
  if (!is(k, "matrix") &&
      !is(k, "vector") &&
      !is.null(k))
    stop("k must be a matrix or a vector")
  sigma <- kpar(kernel)$sigma
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  nu <- kpar(kernel)$order
  ni <- kpar(kernel)$degree
  n <- dim(x)[1]
  lim <- 1 / (gamma(nu + 1) * 2 ^ nu)
  dota <- rowSums(x * x) / 2
  z <- as.matrix(z)
  
  if (!dim(z)[1] == n)
    stop("z must have the length equal to x colums")
  if (is.null(y))
  {
    if (is(z, "matrix") && !dim(z)[1] == n)
      stop("z must have size equal to x colums")
    res <- crossprod(t(x))
    for (i in 1:n)
    {
      xx <- sigma * sqrt(-round(2 * (res[i, ] - dota - rep(dota[i], n)), 9))
      res[i, ] <- besselJ(xx, nu) * (xx ^ (-nu))
      res[i, which(xx < 10e-5)] <- lim
      res[i, ] <- z[i, ] * (((res[i, ] / lim) ^ ni) * z)
    }
    return(res)
  }
  if (is(y, "matrix"))
  {
    if (is.null(k))
      stop("k not specified!")
    m <- dim(y)[1]
    if (!dim(k)[1] == m)
      stop("k must have equal rows to y")
    k <- as.matrix(k)
    if (!dim(x)[2] == dim(y)[2])
      stop("matrixes must have the same number of columns")
    dotb <- rowSums(y * y) / 2
    res <- x %*% t(y)
    for (i in 1:m) {
      #2*sigma or sigma
      xx <-
        sigma * sqrt(-round(2 * (res[, i] - dota - rep(dotb[i], n)), 9))
      res[, i] <- besselJ(xx, nu) * (xx ^ (-nu))
      res[which(xx < 10e-5), i] <- lim
      res[, i] <- k[i, ] * (((res[, i] / lim) ^ ni) * z)
    }
    return(res)
  }
}
setMethod("kernelPol",
          signature(kernel = "besselkernel"),
          kernelPol.besselkernel)


kernelPol.anovakernel <- function(kernel, x, y = NULL, z, k = NULL)
{
  if (!is(y, "matrix") &&
      !is.null(y) && !is(y, "vector"))
    stop("y must be a matrix or NULL")
  if (!is(z, "matrix") &&
      !is(z, "vector"))
    stop("z must be a matrix or a vector")
  if (!is(k, "matrix") &&
      !is(k, "vector") &&
      !is.null(k))
    stop("k must be a matrix or a vector")
  sigma <- kpar(kernel)$sigma
  degree <- kpar(kernel)$degree
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  n <- dim(x)[1]
  z <- as.matrix(z)
  if (!dim(z)[1] == n)
    stop("z must have the length equal to x colums")
  if (is.null(y))
  {
    if (is(z, "matrix") && !dim(z)[1] == n)
      stop("z must have size equal to x colums")
    a <- matrix(0, dim(x)[2], n)
    res <- matrix(0, n, n)
    for (i in 1:n)
    {
      a[rep(TRUE, dim(x)[2]), rep(TRUE, n)] <- x[i, ]
      res[i, ] <-
        z[i, ] * ((colSums(exp(
          -sigma * (a - t(x)) ^ 2
        )) ^ degree) * z)
    }
    return(res)
  }
  if (is(y, "matrix"))
  {
    if (is.null(k))
      stop("k not specified!")
    m <- dim(y)[1]
    k <- as.matrix(k)
    if (!dim(k)[1] == m)
      stop("k must have equal rows to y")
    if (!dim(x)[2] == dim(y)[2])
      stop("matrixes must have the same number of columns")
    
    b <- matrix(0, dim(x)[2], m)
    res <- matrix(0, dim(x)[1], m)
    for (i in 1:n)
    {
      b[rep(TRUE, dim(x)[2]), rep(TRUE, m)] <- x[i, ]
      res[i, ] <-
        z[i, ] * ((colSums(exp(
          -sigma * (b - t(y)) ^ 2
        )) ^ degree) * k)
    }
    return(res)
  }
}
setMethod("kernelPol",
          signature(kernel = "anovakernel"),
          kernelPol.anovakernel)


kernelPol.splinekernel <- function(kernel, x, y = NULL, z, k = NULL)
{
  if (!is(y, "matrix") &&
      !is.null(y) && !is(y, "vector"))
    stop("y must be a matrix or NULL")
  if (!is(z, "matrix") &&
      !is(z, "vector"))
    stop("z must be a matrix or a vector")
  if (!is(k, "matrix") &&
      !is(k, "vector") &&
      !is.null(k))
    stop("k must be a matrix or a vector")
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  sigma <- kpar(kernel)$sigma
  degree <- kpar(kernel)$degree
  n <- dim(x)[1]
  z <- as.vector(z)
  if (!(length(z) == n))
    stop("z must have the length equal to x colums")
  if (is.null(y))
  {
    res <- kernelMatrix(kernel, x)
    return(unclass(z * t(res * z)))
  }
  if (is(y, "matrix"))
  {
    if (is.null(k))
      stop("k not specified!")
    m <- dim(y)[1]
    k <- as.vector(k)
    if (!(length(k) == m))
      stop("k must have length equal to rows of y")
    
    res <- kernelMatrix(kernel, x, y)
    return(unclass(k * t(res * z)))
  }
}
setMethod("kernelPol",
          signature(kernel = "splinekernel"),
          kernelPol.splinekernel)


kernelPol.polykernel <- function(kernel, x, y = NULL, z, k = NULL)
{
  if (!is(y, "matrix") &&
      !is.null(y) && !is(y, "vector"))
    stop("y must be a matrix or NULL")
  if (!is(z, "matrix") &&
      !is(z, "vector"))
    stop("z must be a matrix or a vector")
  if (!is(k, "matrix") &&
      !is(k, "vector") &&
      !is.null(k))
    stop("k must be a matrix or a vector")
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  degree <- kpar(kernel)$degree
  scale <- kpar(kernel)$scale
  offset <- kpar(kernel)$offset
  n <- dim(x)[1]
  
  if (is(z, "matrix"))
  {
    z <- as.vector(z)
  }
  m <- length(z)
  
  if (!(m == n))
    stop("z must have the length equal to x colums")
  if (is.null(y))
  {
    res <- z * t(((scale * crossprod(t(
      x
    )) + offset) ^ degree) * z)
    return(res)
  }
  if (is(y, "matrix"))
  {
    if (is.null(k))
      stop("k not specified!")
    m <- dim(y)[1]
    k <- as.vector(k)
    if (!(length(k) == m))
      stop("k must have length equal to rows of y")
    if (!dim(x)[2] == dim(y)[2])
      stop("matrixes must have the same number of columns")
    res <- k * t(((scale * x %*% t(y) + offset) ^ degree) * z)
    return(res)
  }
}
setMethod("kernelPol",
          signature(kernel = "polykernel"),
          kernelPol.polykernel)


kernelPol.tanhkernel <- function(kernel, x, y = NULL, z, k = NULL)
{
  if (!is(y, "matrix") &&
      !is.null(y) &&
      !is(y, "vector"))
    stop("y must be a matrix, vector or NULL")
  if (!is(z, "matrix") &&
      !is(z, "vector"))
    stop("z must be a matrix or a vector")
  if (!is(k, "matrix") &&
      !is(k, "vector") &&
      !is.null(k))
    stop("k must be a matrix or a vector")
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  scale <- kpar(kernel)$scale
  offset <- kpar(kernel)$offset
  n <- dim(x)[1]
  if (is(z, "matrix"))
  {
    z <- as.vector(z)
  }
  m <- length(z)
  
  if (!(m == n))
    stop("z must have the length equal to x colums")
  if (is.null(y))
  {
    res <- z * t(tanh(scale * crossprod(t(x)) + offset) * z)
    return(res)
  }
  if (is(y, "matrix"))
  {
    if (is.null(k))
      stop("k not specified!")
    m <- dim(y)[1]
    k <- as.vector(k)
    if (!(length(k) == m))
      stop("k must have length equal rows to y")
    if (!dim(x)[2] == dim(y)[2])
      stop("matrixes x, y must have the same number of columns")
    res <- k * t(tanh(scale * x %*% t(y) + offset) * z)
    return(res)
  }
}
setMethod("kernelPol",
          signature(kernel = "tanhkernel"),
          kernelPol.tanhkernel)


kernelPol.vanillakernel <- function(kernel, x, y = NULL, z, k = NULL)
{
  if (!is(y, "matrix") &&
      !is.null(y) &&
      !is(y, "vector"))
    stop("y must be a matrix, vector or NULL")
  if (!is(z, "matrix") &&
      !is(z, "vector"))
    stop("z must be a matrix or a vector")
  if (!is(k, "matrix") &&
      !is(k, "vector") &&
      !is.null(k))
    stop("k must be a matrix or a vector")
  n <- dim(x)[1]
  if (is(z, "matrix"))
  {
    z <- as.vector(z)
  }
  m <- length(z)
  
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  if (!(m == n))
    stop("z must have the length equal to x colums")
  if (is.null(y))
  {
    res <- z * t(crossprod(t(x)) * z)
    return(res)
  }
  if (is(y, "matrix"))
  {
    if (is.null(k))
      stop("k not specified!")
    m <- dim(y)[1]
    k <- as.vector(k)
    if (!length(k) == m)
      stop("k must have length equal rows to y")
    if (!dim(x)[2] == dim(y)[2])
      stop("matrixes x, y must have the same number of columns")
    for (i in 1:m)
      res <- k * t(x %*% t(y) * z)
    return(res)
  }
}
setMethod("kernelPol",
          signature(kernel = "vanillakernel"),
          kernelPol.vanillakernel)

## kernelFast returns the kernel matrix, its usefull in algorithms
## which require iterative kernel matrix computations

kernelFast <- function(kernel, x, y, a)
{
  return(kernelMatrix(kernel, x, y))
}
setGeneric("kernelFast", function(kernel, x, y, a)
  standardGeneric("kernelFast"))



kernelFast.rbfkernel <- function(kernel, x, y, a)
{
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  if (!is(y, "matrix"))
    stop("y must be a matrix or a vector")
  sigma = kpar(kernel)$sigma
  n <- dim(x)[1]
  dota <- a / 2
  if (is(x, "matrix") && is(y, "matrix")) {
    if (!(dim(x)[2] == dim(y)[2]))
      stop("matrixes must have the same number of columns")
    m <- dim(y)[1]
    dotb <- rowSums(y * y) / 2
    res <- x %*% t(y)
    for (i in 1:m)
      res[, i] <- exp(2 * sigma * (res[, i] - dota - rep(dotb[i], n)))
    return(res)
  }
}
setMethod("kernelFast",
          signature(kernel = "rbfkernel"),
          kernelFast.rbfkernel)

kernelFast.laplacekernel <- function(kernel, x, y, a)
{
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  if (!is(y, "matrix"))
    stop("y must be a matrix or a vector")
  sigma = kpar(kernel)$sigma
  n <- dim(x)[1]
  dota <- a / 2
  if (is(x, "matrix") && is(y, "matrix")) {
    if (!(dim(x)[2] == dim(y)[2]))
      stop("matrixes must have the same number of columns")
    m <- dim(y)[1]
    dotb <- rowSums(y * y) / 2
    res <- x %*% t(y)
    for (i in 1:m)
      res[, i] <-
      exp(-sigma * sqrt(round(-2 * (
        res[, i] - dota - rep(dotb[i], n)
      ), 9)))
    return(res)
  }
}
setMethod("kernelFast",
          signature(kernel = "laplacekernel"),
          kernelFast.laplacekernel)

kernelFast.besselkernel <- function(kernel, x, y, a)
{
  if (is(x, "vector"))
    x <- as.matrix(x)
  if (is(y, "vector"))
    y <- as.matrix(y)
  if (!is(y, "matrix"))
    stop("y must be a matrix or a vector")
  sigma = kpar(kernel)$sigma
  nu = kpar(kernel)$order
  ni = kpar(kernel)$degree
  n <- dim(x)[1]
  lim <- 1 / (gamma(nu + 1) * 2 ^ (nu))
  dota <- a / 2
  if (is(x, "matrix") && is(y, "matrix")) {
    if (!(dim(x)[2] == dim(y)[2]))
      stop("matrixes must have the same number of columns")
    m <- dim(y)[1]
    dotb <- rowSums(y * y) / 2
    res <- x %*% t(y)
    for (i in 1:m) {
      xx <- sigma * sqrt(round(-2 * (res[, i] - dota - rep(dotb[i], n)), 9))
      res[, i] <- besselJ(xx, nu) * (xx ^ (-nu))
      res[which(xx < 10e-5), i] <- lim
    }
    return((res / lim) ^ ni)
  }
}
setMethod("kernelFast",
          signature(kernel = "besselkernel"),
          kernelFast.besselkernel)


kernelFast.anovakernel <- function(kernel, x, y, a)
{
  return(kernelMatrix(kernel, x, y))
}
setMethod("kernelFast",
          signature(kernel = "anovakernel"),
          kernelFast.anovakernel)


kernelFast.polykernel <- function(kernel, x, y, a)
{
  return(kernelMatrix(kernel, x, y))
}
setMethod("kernelFast",
          signature(kernel = "polykernel"),
          kernelFast.polykernel)

kernelFast.vanilla <- function(kernel, x, y, a)
{
  return(kernelMatrix(kernel, x, y))
}
setMethod("kernelFast",
          signature(kernel = "vanillakernel"),
          kernelFast.vanilla)

kernelFast.tanhkernel <- function(kernel, x, y, a)
{
  return(kernelMatrix(kernel, x, y))
}
setMethod("kernelFast",
          signature(kernel = "tanhkernel"),
          kernelFast.tanhkernel)

kernelFast.splinekernel <- function(kernel, x, y, a)
{
  return(kernelMatrix(kernel, x, y))
}
setMethod("kernelFast",
          signature(kernel = "splinekernel"),
          kernelFast.splinekernel)