Remove _kernel, sort out kappa (JuliaGaussianProcesses#95)

* Removes _kernel

* Remove _scale

* Removes binary kappa

* Removes k(x, y) -> Matrix

* Fixes tests and adds eval

* Removes metric from Transformed

* Refactors eval

* Remove eval from MyKernel

* Improves style src/basekernels/fbm.jl

Co-Authored-By: Théo Galy-Fajou <[email protected]>

* Reverts removal of _scale

* Reverts change to transforms.

Co-Authored-By: Théo Galy-Fajou <[email protected]>

* Requires Julia 1.3

* Removes 1.0 CI

* Reverts rationalquad.jl changes

* Removes custom show

* Adds reference to JuliaGaussianProcesses#96

* Moves abstract types around

* Pins Zygote to 0.4.16

* Adds extra tests

Co-authored-by: Théo Galy-Fajou <[email protected]>

Author: willtebbutt

.travis.yml

@@ -7,7 +7,6 @@ os:
   - linux
   - osx
 julia:
-  - 1.0
   - 1.3
   - 1.4
   - nightly

Project.toml

@@ -21,7 +21,8 @@ SpecialFunctions = "0.8, 0.9, 0.10"
 StatsBase = "0.32, 0.33"
 StatsFuns = "0.8, 0.9"
 ZygoteRules = "0.2"
-julia = "1.0"
+Zygote = "= 0.4.16"
+julia = "1.3"
 
 [extras]
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"

docs/src/userguide.md

@@ -27,8 +27,9 @@ To premultiply the kernel by a variance, you can use `*` or create a `ScaledKern
 To compute the kernel function on two vectors you can call
 ```julia
   k = SqExponentialKernel()
-  x1 = rand(3); x2 = rand(3)
-  kappa(k,x1,x2) == k(x1,x2) # Syntactic sugar
+  x1 = rand(3)
+  x2 = rand(3)
+  k(x1,x2)
 ```
 
 ## Creating a kernel matrix

src/KernelFunctions.jl

@@ -4,7 +4,7 @@ KernelFunctions. [Github](https://github.com/JuliaGaussianProcesses/KernelFuncti
 """
 module KernelFunctions
 
-export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa
+export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!
 export transform
 export duplicate, set! # Helpers
 
@@ -42,6 +42,7 @@ const defaultobs = 2
 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

src/approximations/nystrom.jl

@@ -12,7 +12,7 @@ end
 function nystrom_sample(k::Kernel, X::AbstractMatrix, S::Vector{<:Integer}; obsdim::Integer=defaultobs)
     obsdim ∈ [1, 2] || throw(ArgumentError("`obsdim` should be 1 or 2 (see docs of kernelmatrix))"))
     Xₘ = obsdim == 1 ? X[S, :] : X[:, S]
-    C = k(Xₘ, X; obsdim=obsdim)
+    C = kernelmatrix(k, Xₘ, X; obsdim=obsdim)
     Cs = C[:, S]
     return (C, Cs)
 end

src/basekernels/cosine.jl

@@ -9,6 +9,7 @@ The cosine kernel is a stationary kernel for a sinusoidal given by
 struct CosineKernel <: SimpleKernel end
 
 kappa(κ::CosineKernel, d::Real) = cospi(d)
+
 metric(::CosineKernel) = Euclidean()
 
 Base.show(io::IO, ::CosineKernel) = print(io, "Cosine Kernel")

src/basekernels/exponential.jl

@@ -12,9 +12,11 @@ related form of the kernel or [`GammaExponentialKernel`](@ref) for a generalizat
 struct SqExponentialKernel <: SimpleKernel end
 
 kappa(κ::SqExponentialKernel, d²::Real) = exp(-d²)
-iskroncompatible(::SqExponentialKernel) = true
+
 metric(::SqExponentialKernel) = SqEuclidean()
 
+iskroncompatible(::SqExponentialKernel) = true
+
 Base.show(io::IO,::SqExponentialKernel) = print(io,"Squared Exponential Kernel")
 
 ## Aliases ##
@@ -33,9 +35,11 @@ The exponential kernel is a Mercer kernel given by the formula:
 struct ExponentialKernel <: SimpleKernel end
 
 kappa(κ::ExponentialKernel, d::Real) = exp(-d)
-iskroncompatible(::ExponentialKernel) = true
+
 metric(::ExponentialKernel) = Euclidean()
 
+iskroncompatible(::ExponentialKernel) = true
+
 Base.show(io::IO, ::ExponentialKernel) = print(io, "Exponential Kernel")
 
 ## Alias ##
@@ -60,7 +64,9 @@ struct GammaExponentialKernel{Tγ<:Real} <: SimpleKernel
 end
 
 kappa(κ::GammaExponentialKernel, d²::Real) = exp(-d²^first(κ.γ))
-iskroncompatible(::GammaExponentialKernel) = true
+
 metric(::GammaExponentialKernel) = SqEuclidean()
 
+iskroncompatible(::GammaExponentialKernel) = true
+
 Base.show(io::IO, κ::GammaExponentialKernel) = print(io, "Gamma Exponential Kernel (γ = ", first(κ.γ), ")")

src/basekernels/exponentiated.jl

@@ -9,7 +9,9 @@ The exponentiated kernel is a Mercer kernel given by:
 struct ExponentiatedKernel <: SimpleKernel end
 
 kappa(κ::ExponentiatedKernel, xᵀy::Real) = exp(xᵀy)
+
 metric(::ExponentiatedKernel) = DotProduct()
+
 iskroncompatible(::ExponentiatedKernel) = true
 
 Base.show(io::IO, ::ExponentiatedKernel) = print(io, "Exponentiated Kernel")

src/basekernels/fbm.jl

@@ -17,6 +17,16 @@ struct FBMKernel{T<:Real} <: BaseKernel
     end
 end
 
+function (κ::FBMKernel)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
+    modX = sum(abs2, x)
+    modY = sum(abs2, y)
+    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)^first(κ.h) + abs2(y)^first(κ.h) - abs2(x - y)^first(κ.h)) / 2
+
 Base.show(io::IO, κ::FBMKernel) = print(io, "Fractional Brownian Motion Kernel (h = ", first(κ.h), ")")
 
 const sqroundoff = 1e-15
@@ -65,13 +75,3 @@ function kernelmatrix!(
     K .= _fbm.(vec(modX), reshape(modY, 1, :), modXY, κ.h)
     return K
 end
-
-function kappa(κ::FBMKernel, x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
-    modX = sum(abs2, x)
-    modY = sum(abs2, y)
-    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)^first(κ.h) + abs2(y)^first(κ.h) - abs2(x-y)^first(κ.h))/2

src/basekernels/gabor.jl

@@ -15,6 +15,8 @@ struct GaborKernel{K<:Kernel} <: BaseKernel
     end
 end
 
+(κ::GaborKernel)(x, y) = κ.kernel(x ,y)
+
 function _gabor(; ell = nothing, p = nothing)
     if ell === nothing
         if p === nothing
@@ -53,8 +55,6 @@ end
 
 Base.show(io::IO, κ::GaborKernel) = print(io, "Gabor Kernel (ell = ", κ.ell, ", p = ", κ.p, ")")
 
-kappa(κ::GaborKernel, x, y) = kappa(κ.kernel, x ,y)
-
 function kernelmatrix(
     κ::GaborKernel,
     X::AbstractMatrix;

src/basekernels/maha.jl

@@ -17,6 +17,7 @@ struct MahalanobisKernel{T<:Real, A<:AbstractMatrix{T}} <: SimpleKernel
 end
 
 kappa(κ::MahalanobisKernel, d::T) where {T<:Real} = exp(-d)
+
 metric(κ::MahalanobisKernel) = SqMahalanobis(κ.P)
 
 Base.show(io::IO, κ::MahalanobisKernel) = print(io, "Mahalanobis Kernel (size(P) = ", size(κ.P), ")")

src/basekernels/matern.jl

@@ -40,6 +40,7 @@ The matern 3/2 kernel is a Mercer kernel given by the formula:
 struct Matern32Kernel <: SimpleKernel end
 
 kappa(κ::Matern32Kernel, d::Real) = (1 + sqrt(3) * d) * exp(-sqrt(3) * d)
+
 metric(::Matern32Kernel) = Euclidean()
 
 Base.show(io::IO, ::Matern32Kernel) = print(io, "Matern 3/2 Kernel")
@@ -55,6 +56,7 @@ The matern 5/2 kernel is a Mercer kernel given by the formula:
 struct Matern52Kernel <: SimpleKernel end
 
 kappa(κ::Matern52Kernel, d::Real) = (1 + sqrt(5) * d + 5 * d^2 / 3) * exp(-sqrt(5) * d)
+
 metric(::Matern52Kernel) = Euclidean()
 
 Base.show(io::IO, ::Matern52Kernel) = print(io, "Matern 5/2 Kernel")

src/basekernels/polynomial.jl

@@ -15,6 +15,7 @@ struct LinearKernel{Tc<:Real} <: SimpleKernel
 end
 
 kappa(κ::LinearKernel, xᵀy::Real) = xᵀy + first(κ.c)
+
 metric(::LinearKernel) = DotProduct()
 
 Base.show(io::IO, κ::LinearKernel) = print(io, "Linear Kernel (c = ", first(κ.c), ")")
@@ -37,7 +38,8 @@ struct PolynomialKernel{Td<:Real, Tc<:Real} <: SimpleKernel
     end
 end
 
-kappa(κ::PolynomialKernel, xᵀy::T) where {T<:Real} = (xᵀy + first(κ.c))^(first(κ.d))
+kappa(κ::PolynomialKernel, xᵀy::Real) = (xᵀy + first(κ.c))^(first(κ.d))
+
 metric(::PolynomialKernel) = DotProduct()
 
 Base.show(io::IO, κ::PolynomialKernel) = print(io, "Polynomial Kernel (c = ", first(κ.c), ", d = ", first(κ.d), ")")

src/generic.jl

@@ -3,34 +3,25 @@ Base.length(::Kernel) = 1
 Base.iterate(k::Kernel) = (k,nothing)
 Base.iterate(k::Kernel, ::Any) = nothing
 
-# default fallback for evaluating a kernel with two arguments (such as vectors etc)
-kappa(κ::Kernel, x, y) = kappa(κ, evaluate(metric(κ), x, y))
-kappa(κ::TransformedKernel, x, y) = kappa(kernel(κ), apply(κ.transform,x), apply(κ.transform,y))
-kappa(κ::TransformedKernel{<:BaseKernel,<:ScaleTransform}, x, y) = kappa(κ, _scale(κ.transform, metric(κ), x, y))
+printshifted(io::IO, o, shift::Int) = print(io, o)
+
+# See https://github.com/JuliaGaussianProcesses/KernelFunctions.jl/issues/96
 _scale(t::ScaleTransform, metric::Euclidean, x, y) =  first(t.s) * evaluate(metric, x, y)
 _scale(t::ScaleTransform, metric::Union{SqEuclidean,DotProduct}, x, y) =  first(t.s)^2 * evaluate(metric, x, y)
 _scale(t::ScaleTransform, metric, x, y) = evaluate(metric, apply(t, x), apply(t, y))
 
-printshifted(io::IO, o, shift::Int) = print(io, o)
-Base.show(io::IO, κ::Kernel) = print(io, nameof(typeof(κ)))
-
 ### Syntactic sugar for creating matrices and using kernel functions
 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)(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
-end
-
 for k in nameof.(subtypes(BaseKernel))
     @eval begin
         @deprecate($k(ρ::Real;args...),transform($k(args...),ρ))
         @deprecate($k(ρ::AbstractVector{<:Real};args...),transform($k(args...),ρ))
     end
 end
+
+# Fallback implementation of evaluate for `SimpleKernel`s.
+(k::SimpleKernel)(x, y) = kappa(k, evaluate(metric(k), x, y))

src/kernels/kernelproduct.jl

@@ -22,7 +22,7 @@ Base.:*(kp::KernelProduct,k::Kernel) = KernelProduct(vcat(kp.kernels,k))
 
 Base.length(k::KernelProduct) = length(k.kernels)
 
-kappa(κ::KernelProduct, x ,y) = prod(kappa(k, x, y) for k in κ.kernels)
+(κ::KernelProduct)(x, y) = prod(k(x, y) for k in κ.kernels)
 
 hadamard(x,y) = x.*y
 

src/kernels/kernelsum.jl

@@ -46,7 +46,7 @@ Base.:*(w::Real, k::KernelSum) = KernelSum(k.kernels, weights = w * k.weights) #
 
 Base.length(k::KernelSum) = length(k.kernels)
 
-kappa(κ::KernelSum, x, y) = sum(κ.weights[i] * kappa(κ.kernels[i], x, y) for i in 1:length(κ))
+(κ::KernelSum)(x, y) = sum(κ.weights[i] * κ.kernels[i](x, y) for i in 1:length(κ))
 
 function kernelmatrix(κ::KernelSum, X::AbstractMatrix; obsdim::Int = defaultobs)
     sum(κ.weights[i] * kernelmatrix(κ.kernels[i], X, obsdim = obsdim) for i in 1:length(κ))

src/kernels/scaledkernel.jl

@@ -15,7 +15,7 @@ end
 
 kappa(k::ScaledKernel, x) = first(k.σ²) * kappa(k.kernel, x)
 
-kappa(k::ScaledKernel, x, y) = first(k.σ²) * kappa(k.kernel, x, y)
+(k::ScaledKernel)(x, y) = first(k.σ²) * k.kernel(x, y)
 
 metric(k::ScaledKernel) = metric(k.kernel)
 

src/kernels/tensorproduct.jl

@@ -22,8 +22,8 @@ end
 
 Base.length(kernel::TensorProduct) = length(kernel.kernels)
 
-function kappa(kernel::TensorProduct, x, y)
-    return prod(kappa(k, xi, yi) for (k, xi, yi) in zip(kernel.kernels, x, y))
+function (kernel::TensorProduct)(x, y)
+    return prod(k(xi, yi) for (k, xi, yi) in zip(kernel.kernels, x, y))
 end
 
 # TODO: General implementation of `kernelmatrix` and `kerneldiagmatrix`

src/kernels/transformedkernel.jl

@@ -9,6 +9,12 @@ struct TransformedKernel{Tk<:Kernel,Tr<:Transform} <: Kernel
     transform::Tr
 end
 
+function (k::TransformedKernel)(x, y)
+    x′ = apply(k.transform, x)
+    y′ = apply(k.transform, y)
+    return k.kernel(x′, y′)
+end
+
 """
 ```julia
     transform(k::BaseKernel, t::Transform) (1)
@@ -29,10 +35,6 @@ transform(k::BaseKernel,ρ::AbstractVector) = TransformedKernel(k, ARDTransform(
 
 kernel(κ) = κ.kernel
 
-kappa(κ::TransformedKernel, x) = kappa(κ.kernel, x)
-
-metric(κ::TransformedKernel) = metric(κ.kernel)
-
 Base.show(io::IO, κ::TransformedKernel) = printshifted(io, κ, 0)
 
 function printshifted(io::IO, κ::TransformedKernel, shift::Int)

test/basekernels/exponential.jl

@@ -13,6 +13,7 @@
         @test GaussianKernel == SqExponentialKernel
         @test SEKernel == SqExponentialKernel
         @test repr(k) == "Squared Exponential Kernel"
+        @test KernelFunctions.iskroncompatible(k) == true
     end
     @testset "ExponentialKernel" begin
         k = ExponentialKernel()
@@ -22,6 +23,7 @@
         @test metric(ExponentialKernel()) == Euclidean()
         @test repr(k) == "Exponential Kernel"
         @test LaplacianKernel == ExponentialKernel
+        @test KernelFunctions.iskroncompatible(k) == true
     end
     @testset "GammaExponentialKernel" begin
         γ = 2.0
@@ -33,6 +35,7 @@
         @test metric(GammaExponentialKernel()) == SqEuclidean()
         @test metric(GammaExponentialKernel(γ=2.0)) == SqEuclidean()
         @test repr(k) == "Gamma Exponential Kernel (γ = $(γ))"
+        @test KernelFunctions.iskroncompatible(k) == true
 
         #Coherence :
         @test GammaExponentialKernel(γ=1.0)(v1,v2) ≈ SqExponentialKernel()(v1,v2)

test/basekernels/fbm.jl

@@ -8,7 +8,6 @@
     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()

test/basekernels/gabor.jl

@@ -5,8 +5,13 @@
     k = GaborKernel(ell=ell, p=p)
     @test k.ell ≈ ell atol=1e-5
     @test k.p ≈ p atol=1e-5
-    @test kappa(k,v1,v2) ≈ exp(-sqeuclidean(v1,v2) ./(k.ell.^2))*cospi(euclidean(v1,v2)./ k.p) atol=1e-5
-    @test kappa(k,v1,v2) ≈ kappa(transform(SqExponentialKernel(), 1/k.ell),v1,v2)*kappa(transform(CosineKernel(), 1/k.p), v1,v2) atol=1e-5
+
+    k_manual = exp(-sqeuclidean(v1, v2) / (k.ell^2)) * cospi(euclidean(v1, v2) / k.p)
+    @test k(v1,v2) ≈ k_manual atol=1e-5
+
+    lhs_manual = transform(SqExponentialKernel(), 1/k.ell)(v1,v2)
+    rhs_manual = transform(CosineKernel(), 1/k.p)(v1,v2)
+    @test k(v1,v2) ≈ lhs_manual * rhs_manual atol=1e-5
 
     k = GaborKernel()
     @test k.ell ≈ 1.0 atol=1e-5

test/basekernels/piecewisepolynomial.jl

@@ -11,7 +11,6 @@
 
     @test k2(v1, v2) ≈ k(v1, v2) atol=1e-5
 
-    @test k(v1, v2) ≈ kappa(k, v1, v2) atol=1e-5
     @test typeof(k(v1, v2)) <: Real
     @test size(kernelmatrix(k, m1, m2)) == (4, 4)
     @test size(kernelmatrix(k, m1)) == (4, 4)

test/kernels/custom.jl

@@ -4,19 +4,8 @@ struct MyKernel <: SimpleKernel end
 KernelFunctions.kappa(::MyKernel, d2::Real) = exp(-d2)
 KernelFunctions.metric(::MyKernel) = SqEuclidean()
 
-# some syntactic sugar
-(κ::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)
-
 @testset "custom" begin
     @test kappa(MyKernel(), 3) == kappa(SqExponentialKernel(), 3)
-    @test kappa(MyKernel(), 1, 3) == kappa(SqExponentialKernel(), 1, 3)
-    @test kappa(MyKernel(), [1, 2], [3, 4]) == kappa(SqExponentialKernel(), [1, 2], [3, 4])
     @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()([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])
 end

test/kernels/kernelproduct.jl

@@ -7,9 +7,9 @@
     k2 = SqExponentialKernel()
     k3 = RationalQuadraticKernel()
 
-    k = KernelProduct([k1,k2])
+    k = KernelProduct([k1, k2])
     @test length(k) == 2
-    @test kappa(k,v1,v2) == kappa(k1*k2,v1,v2)
-    @test kappa(k*k,v1,v2) ≈ kappa(k,v1,v2)^2
-    @test kappa(k*k3,v1,v2) ≈ kappa(k3*k,v1,v2)
+    @test k(v1, v2) == (k1 * k2)(v1, v2)
+    @test (k * k)(v1, v2) ≈ k(v1, v2)^2
+    @test (k * k3)(v1, v2) ≈ (k3 * k)(v1, v2)
 end