Merge pull request #57 from theogf/syntacticsugarallkernels

Created concrete types to call syntactic sugar on all kernels

Author: theogf

src/generic.jl

@@ -15,10 +15,13 @@ printshifted(io::IO,κ::Kernel,shift::Int) = print(io,"$κ")
 Base.show(io::IO,κ::Kernel) = print(io,nameof(typeof(κ)))
 
 ### Syntactic sugar for creating matrices and using kernel functions
-for k in subtypes(BaseKernel)
-    if  k ∈ [FBMKernel] continue end #for kernels without `metric` or `kappa`
+function concretetypes(k, ktypes::Vector)
+    isempty(subtypes(k)) ? push!(ktypes, k) : concretetypes.(subtypes(k), Ref(ktypes))
+    return ktypes
+end
+
+for k in concretetypes(Kernel, [])
     @eval begin
-        @inline (κ::$k)(d::Real) = kappa(κ,d) #TODO Add test
         @inline (κ::$k)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = 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)

src/kernels/fbm.jl

@@ -10,28 +10,30 @@ For `h=1/2`, this is the Wiener Kernel, for `h>1/2`, the increments are
 positively correlated and for `h<1/2` the increments are negatively correlated.
 """
 struct FBMKernel{T<:Real} <: BaseKernel
-    h::T
+    h::Vector{T}
     function FBMKernel(; h::T=0.5) where {T<:Real}
-        @assert h<=1.0 && h>=0.0 "FBMKernel: Given Hurst index h is invalid."
-        return new{T}(h)
+        @assert 0.0 <= h <= 1.0 "FBMKernel: Given Hurst index h is invalid."
+        return new{T}([h])
     end
 end
 
-Base.show(io::IO, κ::FBMKernel) = print(io, "Fractional Brownian Motion Kernel (h = $(k.h))")
+Base.show(io::IO, κ::FBMKernel) = print(io, "Fractional Brownian Motion Kernel (h = $(first(k.h)))")
+
+const sqroundoff = 1e-15
 
 _fbm(modX, modY, modXY, h) = (modX^h + modY^h - modXY^h)/2
 
 function kernelmatrix(κ::FBMKernel, X::AbstractMatrix; obsdim::Int = defaultobs)
     @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
-    modX = sum(abs2, X; dims = 3 - obsdim)
-    modXX = pairwise(SqEuclidean(), X, dims = obsdim)
+    modX = sum(abs2, X; dims = feature_dim(obsdim))
+    modXX = pairwise(SqEuclidean(sqroundoff), X, dims = obsdim)
     return _fbm.(vec(modX), reshape(modX, 1, :), modXX, κ.h)
 end
 
 function kernelmatrix!(K::AbstractMatrix, κ::FBMKernel, X::AbstractMatrix; obsdim::Int = defaultobs)
     @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
-    modX = sum(abs2, X; dims = 3 - obsdim)
-    modXX = pairwise(SqEuclidean(), X, dims = obsdim)
+    modX = sum(abs2, X; dims = feature_dim(obsdim))
+    modXX = pairwise(SqEuclidean(sqroundoff), X, dims = obsdim)
     K .= _fbm.(vec(modX), reshape(modX, 1, :), modXX, κ.h)
     return K
 end
@@ -43,9 +45,9 @@ function kernelmatrix(
     obsdim::Int = defaultobs,
 )
     @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
-    modX = sum(abs2, X, dims=3-obsdim)
-    modY = sum(abs2, Y, dims=3-obsdim)
-    modXY = pairwise(SqEuclidean(), X, Y,dims=obsdim)
+    modX = sum(abs2, X, dims = feature_dim(obsdim))
+    modY = sum(abs2, Y, dims = feature_dim(obsdim))
+    modXY = pairwise(SqEuclidean(sqroundoff), X, Y,dims = obsdim)
     return _fbm.(vec(modX), reshape(modY, 1, :), modXY, κ.h)
 end
 
@@ -57,9 +59,9 @@ function kernelmatrix!(
     obsdim::Int = defaultobs,
 )
     @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
-    modX = sum(abs2, X, dims=3-obsdim)
-    modY = sum(abs2, Y, dims=3-obsdim)
-    modXY = pairwise(SqEuclidean(), X, Y,dims=obsdim)
+    modX = sum(abs2, X, dims = feature_dim(obsdim))
+    modY = sum(abs2, Y, dims = feature_dim(obsdim))
+    modXY = pairwise(SqEuclidean(sqroundoff), X, Y,dims = obsdim)
     K .= _fbm.(vec(modX), reshape(modY, 1, :), modXY, κ.h)
     return K
 end
@@ -72,23 +74,15 @@ function _kernel(
         obsdim::Int = defaultobs
     )
     @assert length(x) == length(y) "x and y don't have the same dimension!"
-    return κ(x,y)
+    return kappa(κ, x, y)
 end
 
-#Syntactic Sugar
-function (κ::FBMKernel)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
+function kappa(κ::FBMKernel, x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
     modX = sum(abs2, x)
     modY = sum(abs2, y)
-    modXY = sqeuclidean(x, y)
-    return (modX^κ.h + modY^κ.h - modXY^κ.h)/2
+    modXY = evaluate(SqEuclidean(sqroundoff), x, y)
+    h = first(κ.h)
+    return (modX^h + modY^h - modXY^h)/2
 end
 
-(κ::FBMKernel)(x::Real, y::Real) = (abs2(x)^κ.h + abs2(y)^κ.h - abs2(x-y)^κ.h)/2
-
-function (κ::FBMKernel)(X::AbstractMatrix{<:Real}, Y::AbstractMatrix{<:Real}; obsdim::Integer=defaultobs)
-    return kernelmatrix(κ, X, Y, obsdim=obsdim)
-end
-
-function (κ::FBMKernel)(X::AbstractMatrix{<:Real}; obsdim::Integer=defaultobs)
-    return kernelmatrix(κ, X, obsdim=obsdim)
-end
+(κ::FBMKernel)(x::Real, y::Real) = (abs2(x)^first(κ.h) + abs2(y)^first(κ.h) - abs2(x-y)^first(κ.h))/2

src/matrix/kernelmatrix.jl

@@ -1,9 +1,8 @@
 """
-```
-    kernelmatrix!(K::Matrix, κ::Kernel, X::Matrix; obsdim::Integer=2)
-    kernelmatrix!(K::Matrix, κ::Kernel, X::Matrix, Y::Matrix; obsdim::Integer=2)
-```
-In-place version of `kernelmatrix` where pre-allocated matrix `K` will be overwritten with the kernel matrix.
+    kernelmatrix!(K::Matrix, κ::Kernel, X::Matrix; obsdim::Integer = 2)
+    kernelmatrix!(K::Matrix, κ::Kernel, X::Matrix, Y::Matrix; obsdim::Integer = 2)
+
+In-place version of [`kernelmatrix`](@ref) where pre-allocated matrix `K` will be overwritten with the kernel matrix.
 """
 kernelmatrix!
 
@@ -21,7 +20,7 @@ function kernelmatrix!(
         map!(x->kappa(κ,x),K,pairwise(metric(κ),X,dims=obsdim))
 end
 
-kernelmatrix!(K::Matrix, κ::TransformedKernel, X::AbstractMatrix; obsdim::Int = defaultobs) =
+kernelmatrix!(K::AbstractMatrix, κ::TransformedKernel, X::AbstractMatrix; obsdim::Int = defaultobs) =
         kernelmatrix!(K, kernel(κ), apply(κ.transform, X, obsdim = obsdim), obsdim = obsdim)
 
 function kernelmatrix!(
@@ -61,13 +60,12 @@ _kernel(κ::TransformedKernel, x::AbstractVector, y::AbstractVector; obsdim::Int
         _kernel(kernel(κ), apply(κ.transform, x), apply(κ.transform, y), obsdim = obsdim)
 
 """
-```
-    kernelmatrix(κ::Kernel, X::Matrix ; obsdim::Int=2)
-    kernelmatrix(κ::Kernel, X::Matrix, Y::Matrix; obsdim::Int=2)
-```
+    kernelmatrix(κ::Kernel, X::Matrix; obsdim::Int = 2)
+    kernelmatrix(κ::Kernel, X::Matrix, Y::Matrix; obsdim::Int = 2)
+
 Calculate the kernel matrix of `X` (and `Y`) with respect to kernel `κ`.
-`obsdim=1` means the matrix `X` (and `Y`) has size #samples x #dimension
-`obsdim=2` means the matrix `X` (and `Y`) has size #dimension x #samples
+`obsdim = 1` means the matrix `X` (and `Y`) has size #samples x #dimension
+`obsdim = 2` means the matrix `X` (and `Y`) has size #dimension x #samples
 """
 kernelmatrix
 
@@ -109,12 +107,11 @@ kernelmatrix(κ::TransformedKernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim
         kernelmatrix(kernel(κ), apply(κ.transform, X, obsdim = obsdim), apply(κ.transform, Y, obsdim = obsdim), obsdim = obsdim)
 
 """
-```
-    kerneldiagmatrix(κ::Kernel, X::Matrix; obsdim::Int=2)
-```
+    kerneldiagmatrix(κ::Kernel, X::Matrix; obsdim::Int = 2)
+
 Calculate the diagonal matrix of `X` with respect to kernel `κ`
-`obsdim=1` means the matrix `X` has size #samples x #dimension
-`obsdim=2` means the matrix `X` has size #dimension x #samples
+`obsdim = 1` means the matrix `X` has size #samples x #dimension
+`obsdim = 2` means the matrix `X` has size #dimension x #samples
 """
 function kerneldiagmatrix(
         κ::Kernel,
@@ -130,10 +127,9 @@ function kerneldiagmatrix(
 end
 
 """
-```
-    kerneldiagmatrix!(K::AbstractVector,κ::Kernel, X::Matrix; obsdim::Int=2)
-```
-In place version of `kerneldiagmatrix`
+    kerneldiagmatrix!(K::AbstractVector,κ::Kernel, X::Matrix; obsdim::Int = 2)
+
+In place version of [`kerneldiagmatrix`](@ref)
 """
 function kerneldiagmatrix!(
         K::AbstractVector,

src/trainable.jl

@@ -4,6 +4,8 @@ import .Flux.trainable
 
 trainable(k::ConstantKernel) = (k.c,)
 
+trainable(k::FBMKernel) = (k.h,)
+
 trainable(k::GammaExponentialKernel) = (k.γ,)
 
 trainable(k::GammaRationalQuadraticKernel) = (k.α, k.γ)

test/kernels/custom.jl

@@ -5,7 +5,6 @@ KernelFunctions.kappa(::MyKernel, d2::Real) = exp(-d2)
 KernelFunctions.metric(::MyKernel) = SqEuclidean()
 
 # some syntactic sugar
-(κ::MyKernel)(d::Real) = kappa(κ, d)
 (κ::MyKernel)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = kappa(κ, x, y)
 (κ::MyKernel)(X::AbstractMatrix{<:Real}, Y::AbstractMatrix{<:Real}; obsdim = 2) = kernelmatrix(κ, X, Y; obsdim = obsdim)
 (κ::MyKernel)(X::AbstractMatrix{<:Real}; obsdim = 2) = kernelmatrix(κ, X; obsdim = obsdim)
@@ -17,7 +16,6 @@ KernelFunctions.metric(::MyKernel) = SqEuclidean()
     @test kernelmatrix(MyKernel(), [1 2; 3 4], [5 6; 7 8]) == kernelmatrix(SqExponentialKernel(), [1 2; 3 4], [5 6; 7 8])
     @test kernelmatrix(MyKernel(), [1 2; 3 4]) == kernelmatrix(SqExponentialKernel(), [1 2; 3 4])
 
-    @test MyKernel()(3) == SqExponentialKernel()(3)
     @test MyKernel()([1, 2], [3, 4]) == SqExponentialKernel()([1, 2], [3, 4])
     @test MyKernel()([1 2; 3 4], [5 6; 7 8]) == SqExponentialKernel()([1 2; 3 4], [5 6; 7 8])
     @test MyKernel()([1 2; 3 4]) == SqExponentialKernel()([1 2; 3 4])

test/kernels/fbm.jl

@@ -0,0 +1,17 @@
+@testset "FBM" begin
+    h = 0.3
+    k = FBMKernel(h = h)
+    v1 = rand(3); v2 = rand(3)
+    @test k(v1,v2) ≈ (sqeuclidean(v1, zero(v1))^h + sqeuclidean(v2, zero(v2))^h - sqeuclidean(v1-v2, zero(v1-v2))^h)/2 atol=1e-5
+
+    # kernelmatrix tests
+    m1 = rand(3,3)
+    m2 = rand(3,3)
+    @test kernelmatrix(k, m1, m1) ≈ kernelmatrix(k, m1) atol=1e-5
+    @test kernelmatrix(k, m1, m2) ≈ k(m1, m2) atol=1e-5
+
+
+    x1 = rand()
+    x2 = rand()
+    @test kernelmatrix(k, x1*ones(1,1), x2*ones(1,1))[1] ≈ k(x1, x2) atol=1e-5
+end

test/runtests.jl

@@ -62,17 +62,18 @@ using KernelFunctions: metric
     end
 
     @testset "kernels" begin
+        include(joinpath("kernels", "constant.jl"))
+        include(joinpath("kernels", "cosine.jl"))
         include(joinpath("kernels", "exponential.jl"))
+        include(joinpath("kernels", "exponentiated.jl"))
+        include(joinpath("kernels", "fbm.jl"))
+        include(joinpath("kernels", "kernelproduct.jl"))
+        include(joinpath("kernels", "kernelsum.jl"))
         include(joinpath("kernels", "matern.jl"))
         include(joinpath("kernels", "polynomial.jl"))
-        include(joinpath("kernels", "constant.jl"))
         include(joinpath("kernels", "rationalquad.jl"))
-        include(joinpath("kernels", "exponentiated.jl"))
-        include(joinpath("kernels", "cosine.jl"))
-        include(joinpath("kernels", "transformedkernel.jl"))
         include(joinpath("kernels", "scaledkernel.jl"))
-        include(joinpath("kernels", "kernelsum.jl"))
-        include(joinpath("kernels", "kernelproduct.jl"))
+        include(joinpath("kernels", "transformedkernel.jl"))
 
         # Legacy tests that don't correspond to anything meaningful in src. Unclear how
         # helpful these are.

test/trainable.jl

@@ -1,17 +1,27 @@
 @testset "trainable" begin
-    ν = 2.0; c = 3.0; d = 2.0; γ = 2.0; α = 2.5
+    ν = 2.0; c = 3.0; d = 2.0; γ = 2.0; α = 2.5; h = 0.5
+
     kc = ConstantKernel(c=c)
     @test all(params(kc) .== params([c]))
-    km = MaternKernel(ν=ν)
-    @test all(params(km) .== params([ν]))
-    kl = LinearKernel(c=c)
-    @test all(params(kl) .== params([c]))
+
+    kfbm = FBMKernel(h = h)
+    @test all(params(kfbm) .== params([h]))
+
     kge = GammaExponentialKernel(γ=γ)
     @test all(params(kge) .== params([γ]))
+
     kgr = GammaRationalQuadraticKernel(γ=γ, α=α)
     @test all(params(kgr) .== params([α], [γ]))
+
+    kl = LinearKernel(c=c)
+    @test all(params(kl) .== params([c]))
+
+    km = MaternKernel(ν=ν)
+    @test all(params(km) .== params([ν]))
+
     kp = PolynomialKernel(c=c, d=d)
     @test all(params(kp) .== params([d], [c]))
+
     kr = RationalQuadraticKernel(α=α)
     @test all(params(kr) .== params([α]))