Merge pull request #50 from JuliaGaussianProcesses/periodickernel

Add Periodic Kernel

Author: theogf

Project.toml

@@ -15,7 +15,7 @@ ZygoteRules = "700de1a5-db45-46bc-99cf-38207098b444"
 
 [compat]
 Compat = "2.2, 3"
-Distances = "0.8"
+Distances = "0.8.2"
 Requires = "1.0.1"
 SpecialFunctions = "0.8, 0.9, 0.10"
 StatsBase = "0.32, 0.33"

docs/src/api.md

@@ -28,6 +28,7 @@ LinearKernel
 PolynomialKernel
 RationalQuadraticKernel
 GammaRationalQuadraticKernel
+PeriodicKernel
 ZeroKernel
 ConstantKernel
 WhiteKernel

docs/src/kernels.md

@@ -91,6 +91,14 @@ The [Polynomial Kernel](@ref KernelFunctions.PolynomialKernel) is defined as
   k(x,x';c,d) = \left(\langle x,x'\rangle + c\right)^d
 ```
 
+## Periodic Kernels
+
+### PeriodicKernel
+
+```math
+  k(x,x';r) = \exp\left(-0.5 \sum_i (sin (π(x_i - x'_i))/r_i)^2\right)
+```
+
 ## Constant Kernels
 
 ### ConstantKernel

src/KernelFunctions.jl

@@ -16,6 +16,7 @@ export MaternKernel, Matern32Kernel, Matern52Kernel
 export LinearKernel, PolynomialKernel
 export RationalQuadraticKernel, GammaRationalQuadraticKernel
 export MahalanobisKernel, GaborKernel, PiecewisePolynomialKernel
+export PeriodicKernel
 export KernelSum, KernelProduct
 export TransformedKernel, ScaledKernel
 
@@ -44,9 +45,10 @@ abstract type BaseKernel <: Kernel end
 include("utils.jl")
 include("distances/dotproduct.jl")
 include("distances/delta.jl")
+include("distances/sinus.jl")
 include("transform/transform.jl")
 
-for k in ["exponential","matern","polynomial","constant","rationalquad","exponentiated","cosine","maha","fbm","gabor","piecewisepolynomial"]
+for k in ["exponential","matern","polynomial","constant","rationalquad","exponentiated","cosine","maha","fbm","gabor","periodic","piecewisepolynomial"]
     include(joinpath("kernels",k*".jl"))
 end
 include("kernels/transformedkernel.jl")

src/distances/delta.jl

@@ -5,8 +5,8 @@ end
     @boundscheck if length(a) != length(b)
         throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
     end
-    return a==b
+    return a == b
 end
 
-@inline (dist::Delta)(a::AbstractArray,b::AbstractArray) = Distances._evaluate(dist,a,b)
-@inline (dist::Delta)(a::Number,b::Number) = a==b
+@inline (dist::Delta)(a::AbstractArray, b::AbstractArray) = Distances._evaluate(dist, a, b)
+@inline (dist::Delta)(a::Number,b::Number) = a == b

src/distances/dotproduct.jl

@@ -1,12 +1,13 @@
-struct DotProduct <: Distances.PreMetric
-end
+struct DotProduct <: Distances.PreMetric end
+# struct DotProduct <: Distances.UnionSemiMetric end
 
-@inline function Distances._evaluate(::DotProduct,a::AbstractVector{T},b::AbstractVector{T}) where {T}
+@inline function Distances._evaluate(::DotProduct, a::AbstractVector{T}, b::AbstractVector{T}) where {T}
     @boundscheck if length(a) != length(b)
         throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
     end
     return dot(a,b)
 end
 
-@inline (dist::DotProduct)(a::AbstractArray,b::AbstractArray) = Distances._evaluate(dist,a,b)
-@inline (dist::DotProduct)(a::Number,b::Number) = a*b
+@inline Distances.eval_op(::DotProduct, a::Real, b::Real) = a * b
+@inline (dist::DotProduct)(a::AbstractArray,b::AbstractArray) = Distances._evaluate(dist, a, b)
+@inline (dist::DotProduct)(a::Number,b::Number) = a * b

src/distances/sinus.jl

@@ -0,0 +1,16 @@
+struct Sinus{T} <: Distances.SemiMetric
+# struct Sinus{T} <: Distances.UnionSemiMetric
+    r::Vector{T}
+end
+
+Distances.parameters(d::Sinus) = d.r
+@inline Distances.eval_op(::Sinus, a::Real, b::Real, p::Real) = abs2(sinpi(a - b) / p)
+@inline (dist::Sinus)(a::AbstractArray, b::AbstractArray) = Distances._evaluate(dist, a, b)
+@inline (dist::Sinus)(a::Number, b::Number) = abs2(sinpi(a - b) / first(dist.r))
+
+@inline function Distances._evaluate(d::Sinus, a::AbstractVector{T}, b::AbstractVector{T}) where {T}
+    @boundscheck if (length(a) != length(b)) || length(a) != length(d.r)
+        throw(DimensionMismatch("Dimensions of the inputs are not matching : a = $(length(a)), b = $(length(b)), r = $(length(d.r))"))
+    end
+    return sum(abs2, sinpi.(a - b) ./ d.r)
+end

src/kernels/periodic.jl

@@ -0,0 +1,27 @@
+"""
+    PeriodicKernel(r::AbstractVector)
+    PeriodicKernel(dims::Int)
+    PeriodicKernel(T::DataType, dims::Int)
+
+Periodic Kernel as described in http://www.inference.org.uk/mackay/gpB.pdf eq. 47.
+```
+    κ(x,y) = exp( - 0.5 sum_i(sin (π(x_i - y_i))/r_i))
+```
+"""
+struct PeriodicKernel{T} <: BaseKernel
+    r::Vector{T}
+    function PeriodicKernel(; r::AbstractVector{T} = ones(Float64, 1)) where {T<:Real}
+        @assert all(r .> 0)
+        new{T}(r)
+    end
+end
+
+PeriodicKernel(dims::Int) = PeriodicKernel(Float64, dims)
+
+PeriodicKernel(T::DataType, dims::Int = 1) = PeriodicKernel(r = ones(T, dims))
+
+metric(κ::PeriodicKernel) = Sinus(κ.r)
+
+kappa(κ::PeriodicKernel, d::Real) = exp(- 0.5d)
+
+Base.show(io::IO, κ::PeriodicKernel) = print(io, "Periodic Kernel, length(r) = $(length(κ.r))")

src/trainable.jl

@@ -14,6 +14,8 @@ trainable(k::MaternKernel) = (k.ν,)
 
 trainable(k::LinearKernel) = (k.c,)
 
+trainable(k::PeriodicKernel) = (k.r,)
+
 trainable(k::PolynomialKernel) = (k.d, k.c)
 
 trainable(k::RationalQuadraticKernel) = (k.α,)

src/zygote_adjoints.jl

@@ -1,5 +1,13 @@
 @adjoint function evaluate(s::DotProduct, x::AbstractVector, y::AbstractVector)
-  dot(x,y), Δ -> begin
-    (nothing, Δ.*y, Δ.*x)
+  dot(x, y), Δ -> begin
+    (nothing, Δ .* y, Δ .* x)
   end
 end
+
+# @adjoint function evaluate(s::Sinus, x::AbstractVector, y::AbstractVector)
+#   d = evaluate(s, x, y)
+#   s = sum(sin.(π*(x-y)))
+#   d, Δ -> begin
+#     (Sinus(Δ ./ s.r), 2Δ .* cos.(x - y) * d, -2Δ .* cos.(x - y) * d)
+#   end
+# end

test/distances/sinus.jl

@@ -0,0 +1,9 @@
+@testset "sinus" begin
+    A = rand(10)
+    B = rand(10)
+    p = rand(10)
+    d = KernelFunctions.Sinus(p)
+    @test Distances.parameters(d) == p
+    @test evaluate(d, A, B) == sum(abs2.(sinpi.(A - B) ./ p))
+    @test d(3.0, 2.0) == abs2(sinpi(3.0 - 2.0) / first(p))
+end

test/kernels/periodic.jl

@@ -0,0 +1,10 @@
+@testset "Periodic Kernel" begin
+    x = rand()*2; v1 = rand(3); v2 = rand(3);
+    r = rand(3)
+    k = PeriodicKernel(r = r)
+    @test kappa(k, x) ≈ exp(-0.5x)
+    @test k(v1, v2) ≈ exp(-0.5 * sum(abs2, sinpi.(v1 - v2) ./ r))
+    @test k(v1, v2) == k(v2, v1)
+    @test PeriodicKernel(3)(v1, v2) == PeriodicKernel(r = ones(3))(v1, v2)
+    @test_nowarn println(k)
+end

test/runtests.jl

@@ -49,6 +49,7 @@ using KernelFunctions: metric
     @testset "distances" begin
         include(joinpath("distances", "dotproduct.jl"))
         include(joinpath("distances", "delta.jl"))
+        include(joinpath("distances", "sinus.jl"))
     end
 
     @testset "transform" begin
@@ -71,6 +72,7 @@ using KernelFunctions: metric
         include(joinpath("kernels", "kernelproduct.jl"))
         include(joinpath("kernels", "kernelsum.jl"))
         include(joinpath("kernels", "matern.jl"))
+        include(joinpath("kernels", "periodic.jl"))
         include(joinpath("kernels", "polynomial.jl"))
         include(joinpath("kernels", "piecewisepolynomial.jl"))
         include(joinpath("kernels", "rationalquad.jl"))

test/trainable.jl

@@ -1,5 +1,5 @@
 @testset "trainable" begin
-    ν = 2.0; c = 3.0; d = 2.0; γ = 2.0; α = 2.5; h = 0.5
+    ν = 2.0; c = 3.0; d = 2.0; γ = 2.0; α = 2.5; h = 0.5; r = rand(3)
 
     kc = ConstantKernel(c=c)
     @test all(params(kc) .== params([c]))
@@ -22,6 +22,9 @@
     kp = PolynomialKernel(c=c, d=d)
     @test all(params(kp) .== params([d], [c]))
 
+    kpe = PeriodicKernel(r = r)
+    @test all(params(kpe) .== params(r))
+
     kr = RationalQuadraticKernel(α=α)
     @test all(params(kr) .== params([α]))