Fix merge conflicts

Author: sharanry

Project.toml

@@ -1,6 +1,6 @@
 name = "KernelFunctions"
 uuid = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"

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 = $(κ.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/approximations/nystrom.jl

@@ -0,0 +1,9 @@
+@testset "nystrom" begin
+    dims = [10,5]
+    X = rand(dims...)
+    k = SqExponentialKernel()
+    for obsdim in [1, 2]
+        @test kernelmatrix(k, X; obsdim=obsdim) ≈ kernelmatrix(nystrom(k, X, 1.0; obsdim=obsdim))
+        @test kernelmatrix(k, X; obsdim=obsdim) ≈ kernelmatrix(nystrom(k, X, collect(1:dims[obsdim]); obsdim=obsdim))
+    end
+end

test/distances/delta.jl

@@ -0,0 +1,10 @@
+@testset "delta" begin
+    A = rand(10,5)
+    B = rand(20,5)
+    d = KernelFunctions.Delta()
+    @test pairwise(d,A,dims=1) == Matrix(I,size(A,1),size(A,1))
+    @test pairwise(d,A,B,dims=1) == zeros(size(A,1),size(B,1))
+    @test d(1,2) == 0
+    @test d(1,1) == 1
+    @test_throws DimensionMismatch d(rand(3),rand(4))
+end

test/distances/dotproduct.jl

@@ -0,0 +1,8 @@
+@testset "dotproduct" begin
+    A = rand(10,5)
+    B = rand(20,5)
+    d = KernelFunctions.DotProduct()
+    @test diag(pairwise(d,A,dims=2)) == [dot(A[:,i],A[:,i]) for i in 1:size(A,2)]
+    @test_throws DimensionMismatch d(rand(3),rand(4))
+    @test d(3.0,2.0) == 6.0
+end

test/test_generic.jl renamed to test/generic.jl

@@ -1,8 +1,5 @@
-using KernelFunctions
-
-k = SqExponentialKernel()
-
-@testset "Generic functions to test" begin
+@testset "generic" begin
+    k = SqExponentialKernel()
     @test length(k) == 1
     @test iterate(k) == (k,nothing)
     @test iterate(k,1) == nothing

test/kernels/constant.jl

@@ -0,0 +1,25 @@
+@testset "constant" begin
+    @testset "ZeroKernel" begin
+        k = ZeroKernel()
+        @test eltype(k) == Any
+        @test kappa(k,2.0) == 0.0
+        @test KernelFunctions.metric(ZeroKernel()) == KernelFunctions.Delta()
+    end
+    @testset "WhiteKernel" begin
+        k = WhiteKernel()
+        @test eltype(k) == Any
+        @test kappa(k,1.0) == 1.0
+        @test kappa(k,0.0) == 0.0
+        @test EyeKernel == WhiteKernel
+        @test metric(WhiteKernel()) == KernelFunctions.Delta()
+    end
+    @testset "ConstantKernel" begin
+        c = 2.0
+        k = ConstantKernel(c=c)
+        @test eltype(k) == Any
+        @test kappa(k,1.0) == c
+        @test kappa(k,0.5) == c
+        @test metric(ConstantKernel()) == KernelFunctions.Delta()
+        @test metric(ConstantKernel(c=2.0)) == KernelFunctions.Delta()
+    end
+end

test/kernels/cosine.jl

@@ -0,0 +1,14 @@
+@testset "cosine" begin
+    rng = MersenneTwister(123456)
+    x = rand(rng)*2
+    v1 = rand(rng, 3)
+    v2 = rand(rng, 3)
+
+    k = CosineKernel()
+    @test eltype(k) == Any
+    @test kappa(k, 1.0) ≈ -1.0 atol=1e-5
+    @test kappa(k, 2.0) ≈ 1.0 atol=1e-5
+    @test kappa(k, 1.5) ≈ 0.0 atol=1e-5
+    @test kappa(k,x) ≈ cospi(x) atol=1e-5
+    @test k(v1, v2) ≈ cospi(sqrt(sum(abs2.(v1-v2)))) atol=1e-5
+end

test/kernels/custom.jl

@@ -0,0 +1,22 @@
+# minimal definition of a custom kernel
+struct MyKernel <: Kernel 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/exponential.jl

@@ -0,0 +1,34 @@
+@testset "exponential" begin
+    rng = MersenneTwister(123456)
+    x = rand(rng)*2
+    v1 = rand(rng, 3)
+    v2 = rand(rng, 3)
+    @testset "SqExponentialKernel" begin
+        k = SqExponentialKernel()
+        @test kappa(k,x) ≈ exp(-x)
+        @test k(v1,v2) ≈ exp(-norm(v1-v2)^2)
+        @test kappa(SqExponentialKernel(),x) == kappa(k,x)
+        @test metric(SqExponentialKernel()) == SqEuclidean()
+    end
+    @testset "ExponentialKernel" begin
+        k = ExponentialKernel()
+        @test kappa(k,x) ≈ exp(-x)
+        @test k(v1,v2) ≈ exp(-norm(v1-v2))
+        @test kappa(ExponentialKernel(),x) == kappa(k,x)
+        @test metric(ExponentialKernel()) == Euclidean()
+    end
+    @testset "GammaExponentialKernel" begin
+        γ = 2.0
+        k = GammaExponentialKernel(γ=γ)
+        @test kappa(k,x) ≈ exp(-(x)^(γ))
+        @test k(v1,v2) ≈ exp(-norm(v1-v2)^(2γ))
+        @test kappa(GammaExponentialKernel(),x) == kappa(k,x)
+        @test GammaExponentialKernel(gamma=γ).γ == [γ]
+        @test metric(GammaExponentialKernel()) == SqEuclidean()
+        @test metric(GammaExponentialKernel(γ=2.0)) == SqEuclidean()
+
+        #Coherence :
+        @test KernelFunctions._kernel(GammaExponentialKernel(γ=1.0),v1,v2) ≈ KernelFunctions._kernel(SqExponentialKernel(),v1,v2)
+        @test KernelFunctions._kernel(GammaExponentialKernel(γ=0.5),v1,v2) ≈ KernelFunctions._kernel(ExponentialKernel(),v1,v2)
+    end
+end

test/kernels/exponentiated.jl

@@ -0,0 +1,12 @@
+@testset "exponentiated" begin
+    rng = MersenneTwister(123456)
+    x = rand(rng)*2
+    v1 = rand(rng, 3)
+    v2 = rand(rng, 3)
+
+    k = ExponentiatedKernel()
+    @test kappa(k,x) ≈ exp(x)
+    @test kappa(k,-x) ≈ exp(-x)
+    @test k(v1,v2) ≈ exp(dot(v1,v2))
+    @test metric(ExponentiatedKernel()) == KernelFunctions.DotProduct()
+end

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/kernels/gabor.jl

@@ -0,0 +1,13 @@
+@testset "Gabor" begin
+    ell = abs(rand())
+    p = abs(rand())
+    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 = GaborKernel()
+    @test k.ell ≈ 1.0 atol=1e-5
+    @test k.p ≈ 1.0 atol=1e-5
+end