Merge pull request #52 from sharanry/gabor

Add Gabor Kernel

Author: theogf

src/KernelFunctions.jl

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

src/kernels/gabor.jl

@@ -0,0 +1,78 @@
+"""
+    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{K<:Kernel} <: BaseKernel
+    kernel::K
+    function GaborKernel(;ell=nothing, p=nothing)
+        k = _gabor(ell=ell, p=p)
+        new{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 == :kernel 
+        return getfield(k, v)
+    elseif v == :ell
+        kernel1 = k.kernel.kernels[1]
+    	if kernel1 isa TransformedKernel
+    	    return 1 ./ kernel1.transform.s[1]
+    	else
+    	    return 1.0
+    	end
+    elseif v == :p
+        kernel2 = k.kernel.kernels[2]
+        if kernel2 isa TransformedKernel
+            return 1 ./ kernel2.transform.s[1]
+        else
+            return 1.0
+        end
+    else
+        error("Invalid Property")
+    end
+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;
+    obsdim::Int=defaultobs)
+    kernelmatrix(κ.kernel, X; obsdim=obsdim)
+end
+
+function kernelmatrix(
+    κ::GaborKernel,
+    X::AbstractMatrix,
+    Y::AbstractMatrix;
+    obsdim::Int=defaultobs)
+    kernelmatrix(κ.kernel, X, Y; obsdim=obsdim)
+end
+
+function kerneldiagmatrix(
+    κ::GaborKernel,
+    X::AbstractMatrix;
+    obsdim::Int=defaultobs) #TODO Add test
+    kerneldiagmatrix(κ.kernel, X; obsdim=obsdim)
+end

src/trainable.jl

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

test/kernels/gabor.jl

@@ -0,0 +1,14 @@
+@testset "Gabor" begin
+    v1 = rand(3); v2 = rand(3)
+    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

test/runtests.jl

@@ -67,6 +67,7 @@ using KernelFunctions: metric
         include(joinpath("kernels", "exponential.jl"))
         include(joinpath("kernels", "exponentiated.jl"))
         include(joinpath("kernels", "fbm.jl"))
+        include(joinpath("kernels", "gabor.jl"))
         include(joinpath("kernels", "kernelproduct.jl"))
         include(joinpath("kernels", "kernelsum.jl"))
         include(joinpath("kernels", "matern.jl"))