Add Gabor Kernel

src/KernelFunctions.jl

@@ -14,7 +14,7 @@ export ExponentiatedKernel
 export MaternKernel, Matern32Kernel, Matern52Kernel
 export LinearKernel, PolynomialKernel
 export RationalQuadraticKernel, GammaRationalQuadraticKernel
-export MahalanobisKernel
+export MahalanobisKernel, GaborKernel
 export KernelSum, KernelProduct
 export TransformedKernel, ScaledKernel
 
@@ -45,7 +45,7 @@ include("distances/dotproduct.jl")
 include("distances/delta.jl")
 include("transform/transform.jl")
 
-for k in ["exponential","matern","polynomial","constant","rationalquad","exponentiated","cosine","maha"]
+for k in ["exponential","matern","polynomial","constant","rationalquad","exponentiated","cosine","maha","gabor"]
     include(joinpath("kernels",k*".jl"))
 end
 include("kernels/transformedkernel.jl")

src/kernels/gabor.jl

@@ -0,0 +1,74 @@
+import Base.getproperty
+
+"""
+    GaborKernel(; ell::Real=1.0, p::Real=1.0)
+
+Gabor kernel with length scale ell and period p. Given by
+```math
+    κ(x,y) =  h(x-z), h(t) = exp(-sum(t.^2./(ell.^2)))*cos(pi*sum(t./p))
+```
+
+"""
+struct GaborKernel{T<:Real, K<:Kernel} <: BaseKernel
+    κ::K
+    function GaborKernel(;ell=nothing, p=nothing)
+        k = _gabor(ell=ell, p=p)
+        if ell == nothing ell=1.0 end
+        if p == nothing p=1.0 end
+        new{Union{typeof(ell),typeof(p)}, typeof(k)}(k)
+    end
+end
+
+function _gabor(; ell = nothing, p = nothing)
+    if ell === nothing
+        if p === nothing
+            return SqExponentialKernel() * CosineKernel()
+        else
+            return SqExponentialKernel() * transform(CosineKernel(), 1 ./ p)
+        end
+    elseif p === nothing
+        return transform(SqExponentialKernel(), 1 ./ ell) * CosineKernel()
+    else
+        return transform(SqExponentialKernel(), 1 ./ ell) * transform(CosineKernel(), 1 ./ p)
+    end
+end
+
+function Base.getproperty(k::GaborKernel, v::Symbol)
+    if v == :κ 
+        return getfield(k, v)
+    elseif v == :ell && typeof(k.κ.kernels[1]) <: SqExponentialKernel
+        return 1.0
+    elseif v == :ell && typeof(k.κ.kernels[1]) <: TransformedKernel
+        return 1 ./ k.κ.kernels[1].transform.s[1]
+    elseif v == :p && typeof(k.κ.kernels[2]) <: CosineKernel
+        return 1.0
+    elseif v == :p && typeof(k.κ.kernels[2]) <: TransformedKernel
+        return 1 ./ k.κ.kernels[2].transform.s[1]
+    else
+        error("Invalid Property")
+    end
+end
+
+kappa(κ::GaborKernel, x, y) where {T<:Real} = kappa(κ.κ, x ,y)
+
+function kernelmatrix(
+    κ::GaborKernel,
+    X::AbstractMatrix;
+    obsdim::Int=defaultobs)
+    kernelmatrix(κ.κ, X; obsdim=obsdim)
+end
+
+function kernelmatrix(
+    κ::GaborKernel,
+    X::AbstractMatrix,
+    Y::AbstractMatrix;
+    obsdim::Int=defaultobs)
+    kernelmatrix(κ.κ, X, Y; obsdim=obsdim)
+end
+
+function kerneldiagmatrix(
+    κ::GaborKernel,
+    X::AbstractMatrix;
+    obsdim::Int=defaultobs) #TODO Add test
+    kerneldiagmatrix(κ.κ, X; obsdim=obsdim)
+end

src/trainable.jl

@@ -18,6 +18,8 @@ trainable(k::RationalQuadraticKernel) = (k.α,)
 
 trainable(k::MahalanobisKernel) = (k.P,)
 
+trainable(k::GaborKernel) = (k.κ,)
+
 #### Composite kernels
 
 trainable(κ::KernelProduct) = κ.kernels

test/test_kernels.jl

@@ -68,6 +68,15 @@ x = rand()*2; v1 = rand(3); v2 = rand(3); id = IdentityTransform()
             @test k(v1,v2) ≈ exp(dot(v1,v2))
         end
     end
+    @testset "Gabor" begin
+        ell = abs(rand())
+        p = abs(rand())
+        k = GaborKernel(ell=ell, p=p)
+        @test k.ell == ell
+        @test k.p == p
+        @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
+    end
     @testset "Matern" begin
         @testset "MaternKernel" begin
             ν = 2.0