[WIP] Rework on kernelmatrix to work with Vectors and more complex kernels

src/KernelFunctions.jl

@@ -10,6 +10,7 @@ export duplicate, set! # Helpers
 
 export Kernel
 export ConstantKernel, WhiteKernel, EyeKernel, ZeroKernel
+export CosineKernel
 export SqExponentialKernel, RBFKernel, GaussianKernel, SEKernel
 export LaplacianKernel, ExponentialKernel, GammaExponentialKernel
 export ExponentiatedKernel
@@ -43,6 +44,7 @@ Abstract type defining a slice-wise transformation on an input matrix
 abstract type Transform end
 abstract type Kernel end
 abstract type BaseKernel <: Kernel end
+abstract type SimpleKernel <: BaseKernel end
 
 include("utils.jl")
 include("distances/dotproduct.jl")

src/basekernels/constant.jl

@@ -7,14 +7,15 @@ Create a kernel that always returning zero
 ```
 The output type depends of `x` and `y`
 """
-struct ZeroKernel <: BaseKernel end
+struct ZeroKernel <: SimpleKernel end
 
 kappa(κ::ZeroKernel, d::T) where {T<:Real} = zero(T)
 
 metric(::ZeroKernel) = Delta()
 
 Base.show(io::IO, ::ZeroKernel) = print(io, "Zero Kernel")
 
+issimple(ZeroKernel) = true
 
 """
     WhiteKernel()
@@ -24,7 +25,7 @@ Base.show(io::IO, ::ZeroKernel) = print(io, "Zero Kernel")
 ```
 Kernel function working as an equivalent to add white noise. Can also be called via `EyeKernel()`
 """
-struct WhiteKernel <: BaseKernel end
+struct WhiteKernel <: SimpleKernel end
 
 """
     EyeKernel()
@@ -48,7 +49,7 @@ Kernel function always returning a constant value `c`
     κ(x,y) = c
 ```
 """
-struct ConstantKernel{Tc<:Real} <: BaseKernel
+struct ConstantKernel{Tc<:Real} <: SimpleKernel
     c::Vector{Tc}
     function ConstantKernel(;c::T=1.0) where {T<:Real}
         new{T}([c])

src/basekernels/cosine.jl

@@ -6,7 +6,7 @@ The cosine kernel is a stationary kernel for a sinusoidal given by
     κ(x,y) = cos( π * (x-y) )
 ```
 """
-struct CosineKernel <: BaseKernel end
+struct CosineKernel <: SimpleKernel end
 
 kappa(κ::CosineKernel, d::Real) = cospi(d)
 metric(::CosineKernel) = Euclidean()

src/basekernels/exponential.jl

@@ -9,7 +9,7 @@ Can also be called via `SEKernel`, `GaussianKernel` or `SEKernel`.
 See also [`ExponentialKernel`](@ref) for a
 related form of the kernel or [`GammaExponentialKernel`](@ref) for a generalization.
 """
-struct SqExponentialKernel <: BaseKernel end
+struct SqExponentialKernel <: SimpleKernel end
 
 kappa(κ::SqExponentialKernel, d²::Real) = exp(-d²)
 iskroncompatible(::SqExponentialKernel) = true
@@ -30,7 +30,7 @@ The exponential kernel is a Mercer kernel given by the formula:
     κ(x,y) = exp(-‖x-y‖)
 ```
 """
-struct ExponentialKernel <: BaseKernel end
+struct ExponentialKernel <: SimpleKernel end
 
 kappa(κ::ExponentialKernel, d::Real) = exp(-d)
 iskroncompatible(::ExponentialKernel) = true
@@ -51,7 +51,7 @@ The γ-exponential kernel is an isotropic Mercer kernel given by the formula:
 Where `γ > 0`, (the keyword `γ` can be replaced by `gamma`)
 For `γ = 1`, see `SqExponentialKernel` and `γ = 0.5`, see `ExponentialKernel`
 """
-struct GammaExponentialKernel{Tγ<:Real} <: BaseKernel
+struct GammaExponentialKernel{Tγ<:Real} <: SimpleKernel
     γ::Vector{Tγ}
     function GammaExponentialKernel(; gamma::T=2.0, γ::T=gamma) where {T<:Real}
         @check_args(GammaExponentialKernel, γ, γ >= zero(T), "γ > 0")

src/basekernels/exponentiated.jl

@@ -6,7 +6,7 @@ The exponentiated kernel is a Mercer kernel given by:
     κ(x,y) = exp(xᵀy)
 ```
 """
-struct ExponentiatedKernel <: BaseKernel end
+struct ExponentiatedKernel <: SimpleKernel end
 
 kappa(κ::ExponentiatedKernel, xᵀy::Real) = exp(xᵀy)
 metric(::ExponentiatedKernel) = DotProduct()

src/basekernels/fbm.jl

@@ -66,17 +66,6 @@ function kernelmatrix!(
     return K
 end
 
-## Apply kernel on two vectors ##
-function _kernel(
-        κ::FBMKernel,
-        x::AbstractVector,
-        y::AbstractVector;
-        obsdim::Int = defaultobs
-    )
-    @assert length(x) == length(y) "x and y don't have the same dimension!"
-    return kappa(κ, x, y)
-end
-
 function kappa(κ::FBMKernel, x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
     modX = sum(abs2, x)
     modY = sum(abs2, y)

src/basekernels/maha.jl

@@ -8,7 +8,7 @@ Mahalanobis distance-based kernel given by
 where the matrix P is the metric.
 
 """
-struct MahalanobisKernel{T<:Real, A<:AbstractMatrix{T}} <: BaseKernel
+struct MahalanobisKernel{T<:Real, A<:AbstractMatrix{T}} <: SimpleKernel
     P::A
     function MahalanobisKernel(P::AbstractMatrix{T}) where {T<:Real}
         LinearAlgebra.checksquare(P)

src/basekernels/matern.jl

@@ -7,7 +7,7 @@ The matern kernel is a Mercer kernel given by the formula:
 ```
 For `ν=n+1/2, n=0,1,2,...` it can be simplified and you should instead use [`ExponentialKernel`](@ref) for `n=0`, [`Matern32Kernel`](@ref), for `n=1`, [`Matern52Kernel`](@ref) for `n=2` and [`SqExponentialKernel`](@ref) for `n=∞`.
 """
-struct MaternKernel{Tν<:Real} <: BaseKernel
+struct MaternKernel{Tν<:Real} <: SimpleKernel
     ν::Vector{Tν}
     function MaternKernel(;nu::T=1.5, ν::T=nu) where {T<:Real}
         @check_args(MaternKernel, ν, ν > zero(T), "ν > 0")
@@ -37,7 +37,7 @@ The matern 3/2 kernel is a Mercer kernel given by the formula:
     κ(x,y) = (1+√(3)‖x-y‖)exp(-√(3)‖x-y‖)
 ```
 """
-struct Matern32Kernel <: BaseKernel end
+struct Matern32Kernel <: SimpleKernel end
 
 kappa(κ::Matern32Kernel, d::Real) = (1 + sqrt(3) * d) * exp(-sqrt(3) * d)
 metric(::Matern32Kernel) = Euclidean()
@@ -52,7 +52,7 @@ The matern 5/2 kernel is a Mercer kernel given by the formula:
     κ(x,y) = (1+√(5)‖x-y‖ + 5/3‖x-y‖^2)exp(-√(5)‖x-y‖)
 ```
 """
-struct Matern52Kernel <: BaseKernel end
+struct Matern52Kernel <: SimpleKernel end
 
 kappa(κ::Matern52Kernel, d::Real) = (1 + sqrt(5) * d + 5 * d^2 / 3) * exp(-sqrt(5) * d)
 metric(::Matern52Kernel) = Euclidean()

src/basekernels/periodic.jl

@@ -8,7 +8,7 @@ Periodic Kernel as described in http://www.inference.org.uk/mackay/gpB.pdf eq. 4
     κ(x,y) = exp( - 0.5 sum_i(sin (π(x_i - y_i))/r_i))
 ```
 """
-struct PeriodicKernel{T} <: BaseKernel
+struct PeriodicKernel{T} <: SimpleKernel
     r::Vector{T}
     function PeriodicKernel(; r::AbstractVector{T} = ones(Float64, 1)) where {T<:Real}
         @assert all(r .> 0)

src/basekernels/piecewisepolynomial.jl

@@ -43,16 +43,6 @@ function kappa(
     return _piecewisepolynomial(κ, r, j)
 end
 
-function _kernel(
-    κ::PiecewisePolynomialKernel,
-    x::AbstractVector,
-    y::AbstractVector;
-    obsdim::Int = defaultobs,
-)
-    @assert length(x) == length(y) "x and y don't have the same dimension!"
-    return kappa(κ,x,y)
-end
-
 function kernelmatrix(
     κ::PiecewisePolynomialKernel{V},
     X::AbstractMatrix;
@@ -62,6 +52,17 @@ function kernelmatrix(
     return map(r->_piecewisepolynomial(κ, r, j), pairwise(metric(κ), X; dims=obsdim))
 end
 
+function kernelmatrix(
+    κ::PiecewisePolynomialKernel{V},
+    X::AbstractMatrix,
+    Y::AbstractMatrix;
+    obsdim::Int = defaultobs
+) where {V}
+    j = div(size(X, feature_dim(obsdim)), 2) + V + 1
+    return map(r->_piecewisepolynomial(κ, r, j), pairwise(metric(κ), X, Y; dims=obsdim))
+end
+
+
 function _kernelmatrix(κ::PiecewisePolynomialKernel{V}, X, Y, obsdim) where {V}
     j = div(size(X, feature_dim(obsdim)), 2) + V + 1
     return map(r->_piecewisepolynomial(κ, r, j), pairwise(metric(κ), X, Y; dims=obsdim))

src/basekernels/polynomial.jl

@@ -7,7 +7,7 @@ The linear kernel is a Mercer kernel given by
 ```
 Where `c` is a real number
 """
-struct LinearKernel{Tc<:Real} <: BaseKernel
+struct LinearKernel{Tc<:Real} <: SimpleKernel
     c::Vector{Tc}
     function LinearKernel(;c::T=0.0) where {T}
         new{T}([c])
@@ -28,7 +28,7 @@ The polynomial kernel is a Mercer kernel given by
 ```
 Where `c` is a real number, and `d` is a shape parameter bigger than 1. For `d = 1` see [`LinearKernel`](@ref)
 """
-struct PolynomialKernel{Td<:Real, Tc<:Real} <: BaseKernel
+struct PolynomialKernel{Td<:Real, Tc<:Real} <: SimpleKernel
     d::Vector{Td}
     c::Vector{Tc}
     function PolynomialKernel(; d::Td=2.0, c::Tc=0.0) where {Td<:Real, Tc<:Real}

src/basekernels/rationalquad.jl

@@ -7,7 +7,7 @@ The rational-quadratic kernel is a Mercer kernel given by the formula:
 ```
 where `α` is a shape parameter of the Euclidean distance. Check [`GammaRationalQuadraticKernel`](@ref) for a generalization.
 """
-struct RationalQuadraticKernel{Tα<:Real} <: BaseKernel
+struct RationalQuadraticKernel{Tα<:Real} <: SimpleKernel
     α::Vector{Tα}
     function RationalQuadraticKernel(;alpha::T=2.0, α::T=alpha) where {T}
         @check_args(RationalQuadraticKernel, α, α > zero(T), "α > 1")
@@ -28,7 +28,7 @@ The Gamma-rational-quadratic kernel is an isotropic Mercer kernel given by the f
 ```
 where `α` is a shape parameter of the Euclidean distance and `γ` is another shape parameter.
 """
-struct GammaRationalQuadraticKernel{Tα<:Real, Tγ<:Real} <: BaseKernel
+struct GammaRationalQuadraticKernel{Tα<:Real, Tγ<:Real} <: SimpleKernel
     α::Vector{Tα}
     γ::Vector{Tγ}
     function GammaRationalQuadraticKernel(;alpha::Tα=2.0, gamma::Tγ=2.0, α::Tα=alpha, γ::Tγ=gamma) where {Tα<:Real, Tγ<:Real}

src/generic.jl

@@ -22,7 +22,7 @@ end
 
 for k in concretetypes(Kernel, [])
     @eval begin
-        @inline (κ::$k)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = kappa(κ, x, y)
+        @inline (κ::$k)(x, y) = kappa(κ, x, y)
         @inline (κ::$k)(X::AbstractMatrix{T}, Y::AbstractMatrix{T}; obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ, X, Y, obsdim=obsdim)
         @inline (κ::$k)(X::AbstractMatrix{T}; obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ, X, obsdim=obsdim)
     end

src/kernels/scaledkernel.jl

@@ -15,6 +15,8 @@ end
 
 kappa(k::ScaledKernel, x) = first(k.σ²) * kappa(k.kernel, x)
 
+kappa(k::ScaledKernel, x, y) = first(k.σ²) * kappa(k.kernel, x, y)
+
 metric(k::ScaledKernel) = metric(k.kernel)
 
 Base.:*(w::Real, k::Kernel) = ScaledKernel(k, w)

src/kernels/tensorproduct.jl

@@ -22,7 +22,6 @@ end
 
 Base.length(kernel::TensorProduct) = length(kernel.kernels)
 
-(kernel::TensorProduct)(x, y) = kappa(kernel, x, y)
 function kappa(kernel::TensorProduct, x, y)
     return prod(kappa(k, xi, yi) for (k, xi, yi) in zip(kernel.kernels, x, y))
 end

src/kernels/transformedkernel.jl

@@ -39,3 +39,17 @@ function printshifted(io::IO, κ::TransformedKernel, shift::Int)
     printshifted(io, κ.kernel, shift)
     print(io,"\n" * ("\t" ^ (shift + 1)) * "- $(κ.transform)")
 end
+
+# Kernel matrix operations
+
+kernelmatrix!(K::AbstractMatrix, κ::TransformedKernel, X::AbstractMatrix; obsdim::Int = defaultobs) =
+    kernelmatrix!(K, kernel(κ), apply(κ.transform, X, obsdim = obsdim), obsdim = obsdim)
+
+kernelmatrix!(K::AbstractMatrix, κ::TransformedKernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim::Int = defaultobs) =
+    kernelmatrix!(K, kernel(κ), apply(κ.transform, X, obsdim = obsdim), apply(κ.transform, Y, obsdim = obsdim), obsdim = obsdim)
+
+kernelmatrix(κ::TransformedKernel, X::AbstractMatrix; obsdim::Int = defaultobs) =
+    kernelmatrix(kernel(κ), apply(κ.transform, X, obsdim = obsdim), obsdim = obsdim)
+
+kernelmatrix(κ::TransformedKernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim::Int = defaultobs) =
+    kernelmatrix(kernel(κ), apply(κ.transform, X, obsdim = obsdim), apply(κ.transform, Y, obsdim = obsdim), obsdim = obsdim)