Merge branch 'master-dev'

Author: theogf

benchmark/kernelmatrix.jl

@@ -2,7 +2,7 @@ using KernelFunctions
 
 SUITE["KernelFunctions"] = BenchmarkGroup()
 
-kernelnames = ["SquaredExponentialKernel"]
+kernelnames = ["SqExponentialKernel"]
 kerneltypes = ["ARD","ISO"]
 kernels=Dict{String,Dict{String,KernelFunctions.Kernel}}()
 for k in kernelnames

dev/debugAD.jl

@@ -12,39 +12,39 @@ testfunction(k,A,B) = det(kernelmatrix(k,A,B))
 testfunction(k,A) = sum(kernelmatrix(k,A))
 k = MaternKernel(vl)
 KernelFunctions.kappa(k,3)
-testfunction(SquaredExponentialKernel(vl),A)
+testfunction(SqExponentialKernel(vl),A)
 testfunction(MaternKernel(vl),A)
 @which kernelmatrix(MaternKernel(vl),A,B)
 #For debugging
 @info "Running Zygote gradients"
 Zygote.refresh()
 ## Zygote
-Zygote.gradient(x->testfunction(SquaredExponentialKernel(x),A),vl)
+Zygote.gradient(x->testfunction(SqExponentialKernel(x),A),vl)
 Zygote.gradient(x->testfunction(MaternKernel(x),A),vl)
-Zygote.gradient(x->testfunction(SquaredExponentialKernel(x),A,B),vl)[1]
+Zygote.gradient(x->testfunction(SqExponentialKernel(x),A,B),vl)[1]
 Zygote.gradient(x->testfunction(MaternKernel(x),A,B),vl)[1]
-Zygote.gradient(x->testfunction(SquaredExponentialKernel(x),A,B),l)
+Zygote.gradient(x->testfunction(SqExponentialKernel(x),A,B),l)
 Zygote.gradient(x->testfunction(MaternKernel(x),A,B),l)
-Zygote.gradient(x->testfunction(SquaredExponentialKernel(x),A),l)
+Zygote.gradient(x->testfunction(SqExponentialKernel(x),A),l)
 Zygote.gradient(x->testfunction(MaternKernel(x),A),l)
 Zygote.gradient(x->testfunction(MaternKernel(x),A),l)
 Zygote.gradient(x->kernelmatrix(MaternKernel(x,1.0),A)[1],l)
 @info "Running Tracker gradients"
 ## Tracker
-# Tracker.gradient(x->testfunction(SquaredExponentialKernel(vl),x,B),A)
-# Tracker.gradient(x->testfunction(SquaredExponentialKernel(l),x[:,:]),A)
-# # Tracker.gradient(x->testfunction(SquaredExponentialKernel(x),A,B),vl)
-# Tracker.gradient(x->testfunction(SquaredExponentialKernel(x),A),vl)
-# Tracker.gradient(x->testfunction(SquaredExponentialKernel(x),A,B),l)
-# Tracker.gradient(x->testfunction(SquaredExponentialKernel(x),A),l)
+# Tracker.gradient(x->testfunction(SqExponentialKernel(vl),x,B),A)
+# Tracker.gradient(x->testfunction(SqExponentialKernel(l),x[:,:]),A)
+# # Tracker.gradient(x->testfunction(SqExponentialKernel(x),A,B),vl)
+# Tracker.gradient(x->testfunction(SqExponentialKernel(x),A),vl)
+# Tracker.gradient(x->testfunction(SqExponentialKernel(x),A,B),l)
+# Tracker.gradient(x->testfunction(SqExponentialKernel(x),A),l)
 
 @info "Running ForwardDiff gradients"
 ## ForwardDiff
-ForwardDiff.gradient(x->testfunction(SquaredExponentialKernel(x),A,B),vl) #✓
+ForwardDiff.gradient(x->testfunction(SqExponentialKernel(x),A,B),vl) #✓
 ForwardDiff.gradient(x->testfunction(MaternKernel(x),A,B),vl) #✓
-ForwardDiff.gradient(x->testfunction(SquaredExponentialKernel(x),A),vl) #✓
+ForwardDiff.gradient(x->testfunction(SqExponentialKernel(x),A),vl) #✓
 ForwardDiff.gradient(x->testfunction(MaternKernel(x),A),vl) #✓
-ForwardDiff.gradient(x->testfunction(SquaredExponentialKernel(x[1]),A,B),[l])
+ForwardDiff.gradient(x->testfunction(SqExponentialKernel(x[1]),A,B),[l])
 ForwardDiff.gradient(x->testfunction(MaternKernel(x[1]),A,B),[l])
-ForwardDiff.gradient(x->testfunction(SquaredExponentialKernel(x[1]),A),[l])
+ForwardDiff.gradient(x->testfunction(SqExponentialKernel(x[1]),A),[l])
 ForwardDiff.gradient(x->testfunction(MaternKernel(x[1]),A),[l])

dev/matrixvsvectors.jl

@@ -14,7 +14,7 @@ timekf = similar(timestheno); memkf = similar(timestheno)
     B = randn(D,1001)
 
     # Standardised eq kernel with length-scale 0.1.
-    medkf = median(@benchmark KernelFunctions.kernelmatrix(SquaredExponentialKernel(0.01),$A,$B,obsdim=2))
+    medkf = median(@benchmark KernelFunctions.kernelmatrix(SqExponentialKernel(0.01),$A,$B,obsdim=2))
     timekf[i] = medkf.time/1e6; memkf[i] = medkf.memory/2^20
     medstheno = median(@benchmark pw(eq(; l=0.1), ColsAreObs($A), ColsAreObs($B)))
     timestheno[i] = medstheno.time/1e6; memstheno[i] = medstheno.memory/2^20

src/KernelFunctions.jl

@@ -1,7 +1,7 @@
 module KernelFunctions
 
 export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa
-export Kernel, SquaredExponentialKernel, MaternKernel, Matern32Kernel, Matern52Kernel
+export Kernel, SqExponentialKernel, MaternKernel, Matern32Kernel, Matern52Kernel
 
 export Transform, ScaleTransform
 
@@ -18,7 +18,7 @@ include("utils.jl")
 include("transform/transform.jl")
 include("kernelmatrix.jl")
 
-kernels = ["squaredexponential","matern"]
+kernels = ["sqexponential","matern"]
 for k in kernels
     include(joinpath("kernels",k*".jl"))
 end

src/generic.jl

@@ -1,11 +1,12 @@
 
 @inline metric(κ::Kernel) = κ.metric
 kernels =
-for k in [:SquaredExponentialKernel,:MaternKernel,:Matern32Kernel,:Matern52Kernel]
+for k in [:SqExponentialKernel,:MaternKernel,:Matern32Kernel,:Matern52Kernel]
     eval(quote
         @inline (κ::$k)(d::Real) = kappa(κ,d)
         @inline (κ::$k)(x::AbstractVector{T},y::AbstractVector{T}) where {T} = kernel(κ,evaluate(κ.(metric),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
 ### Transform generics

src/kernels/matern.jl

@@ -1,26 +1,26 @@
 """
-    MaternKernel([[ρ=1],ν=3/2])
+    MaternKernel([ρ=1.0,[ν=1.0]])
 
 The matern kernel is an isotropic Mercer kernel given by the formula:
 
 ```
-    κ(x,y) = 2^{1-ν}/Γ(ν)*(√(2ν)‖x-y‖/ρ)^ν K_ν(√(2ν)‖x-y‖/ρ)
+    κ(x,y) = 2^{1-ν}/Γ(ν)*(√(2ν)‖x-y‖)^ν K_ν(√(2ν)‖x-y‖)
 ```
 
-For `ν=n+1/2, n=0,1,2,...` it can be simplified (it will be converted automatically).
-`ρ` is a lengthscale parameter.
+For `ν=n+1/2, n=0,1,2,...` it can be simplified and you should instead use `ExponentialKernel` for `n=0`, `Matern32Kernel`, for `n=1`, Matern52Kernel for `n=2` and `SqExponentialKernel` for `n=∞`.
+`ρ` is the lengthscale parameter(s) or the transform object.
 
 # Examples
 
 ```jldoctest; setup = :(using KernelFunctions)
 julia> MaternKernel()
-Matern3_2Kernel{Float64,Float64}(1.0)
+MaternKernel{Float64,Float64}(1.0,1.0)
 
 julia> MaternKernel(2.0f0,3.0)
 MaternKernel{Float32,Float32}(2.0,3.0)
 
-julia> MaternKernel([2.0,3.0],5/2)
-Matern5_2Kernel{Float64,Array{Float64}}([2.0,3.0])
+julia> MaternKernel([2.0,3.0],2.5)
+MaternKernel{Float64,Array{Float64}}([2.0,3.0],2.5)
 ```
 """
 struct MaternKernel{T,Tr<:Transform} <: Kernel{T,Tr}
@@ -34,47 +34,17 @@ end
 
 function MaternKernel(ρ::T₁=1.0,ν::T₂=1.5) where {T₁<:Real,T₂<:Real}
     @check_args(MaternKernel, ν, ν > zero(T₂), "ν > 0")
-    if ν == 0.5
-        ExponentialKernel{T₁,ScaleTransform{T₁}}(ScaleTransform(ρ))
-    elseif ν == 1.5
-        Matern32Kernel{T₁,ScaleTransform{T₁}}(ScaleTransform(ρ))
-    elseif ν == 2.5
-        Matern52Kernel{T₁,ScaleTransform{T₁}}(ScaleTransform(ρ))
-    elseif ν == Inf
-        SquaredExponentialKernel{T₁,ScaleTransform{T₁}}(ScaleTransform(ρ))
-    else
-        MaternKernel{T₁,ScaleTransform{T₁}}(ScaleTransform(ρ),ν)
-    end
+    MaternKernel{T₁,ScaleTransform{T₁}}(ScaleTransform(ρ),ν)
 end
 
 function MaternKernel(ρ::A,ν::T=1.5) where {A<:AbstractVector{<:Real},T<:Real}
     @check_args(MaternKernel, ν, ν > zero(T), "ν > 0")
-    if ν == 0.5
-        ExponentialKernel{eltype(A),ScaleTransform{A}}(ScaleTransform(ρ))
-    elseif ν == 1.5
-        Matern32Kernel{eltype(A),ScaleTransform{A}}(ScaleTransform(ρ))
-    elseif ν == 2.5
-        Matern52Kernel{eltype(A),ScaleTransform{A}}(ScaleTransform(ρ))
-    elseif ν == Inf
-        SquaredExponentialKernel{eltype(A),ScaleTransform{A}}(ScaleTransform(ρ))
-    else
-        MaternKernel{eltype(A),ScaleTransform{A}}(ScaleTransform(ρ),ν)
-    end
+    MaternKernel{eltype(A),ScaleTransform{A}}(ScaleTransform(ρ),ν)
 end
 
-function MaternKernel(t::T₁,ν::T₂=1.5) where {T₁<:Transform,T₂<:Real}
-    @check_args(MaternKernel, ν, ν > zero(T₂), "ν > 0")
-    if ν == 0.5
-        ExponentialKernel{eltype(t),T₁}(t)
-    elseif ν == 1.5
-        Matern32Kernel{eltype(t),T₁}(t)
-    elseif ν == 2.5
-        Matern52Kernel{eltype(t),T₁}(t)
-    elseif ν == Inf
-        SquaredExponentialKernel{eltype(t),T₁}(t)
-    else
-        MaternKernel{eltype(t),T₁}(t,ν)
-    end
+function MaternKernel(t::Tr,ν::T=1.5) where {Tr<:Transform,T<:Real}
+    @check_args(MaternKernel, ν, ν > zero(T), "ν > 0")
+    MaternKernel{eltype(t),Tr}(t,ν)
 end
 
 @inline kappa(κ::MaternKernel, d::Real) where {T} = exp((1.0-κ.ν)*logtwo - lgamma(κ.ν) - κ.ν*log(sqrt(2κ.ν)*d))*besselk(κ.ν,sqrt(2κ.ν)*d)
@@ -96,8 +66,8 @@ function Matern32Kernel(ρ::A) where {A<:AbstractVector{<:Real}}
     Matern32Kernel{eltype(A),ScaleTransform{A}}(ScaleTransform(ρ))
 end
 
-function Matern32Kernel(t::Transform)
-    Matern52Kernel{eltype(A),ScaleTransform{A}}(t)
+function Matern32Kernel(t::Tr) where {Tr<:Transform}
+    Matern52Kernel{eltype(Tr),Tr}(t)
 end
 
 @inline kappa(κ::Matern32Kernel, d::T) where {T<:Real} = (1+sqrt(3)*d)*exp(-sqrt(3)*d)
@@ -118,8 +88,8 @@ function Matern52Kernel(ρ::A) where {A<:AbstractVector{<:Real}}
     Matern52Kernel{eltype(A),ScaleTransform{A}}(ScaleTransform(ρ))
 end
 
-function Matern52Kernel(t::Transform)
-    Matern52Kernel{eltype(A),ScaleTransform{A}}(t)
+function Matern52Kernel(t::Tr) where {Tr<:Transform}
+    Matern52Kernel{eltype(Tr),Tr}(t)
 end
 
 @inline kappa(κ::Matern52Kernel, d::Real) where {T} = (1+sqrt(5)*d+5*d^2/3)*exp(-sqrt(5)*d)

src/kernels/sqexponential.jl

@@ -0,0 +1,53 @@
+"""
+    SqExponentialKernel([α=1])
+
+The squared exponential kernel is an isotropic Mercer kernel given by the formula:
+
+```
+    κ(x,y) = exp(-‖x-y‖²)
+```
+
+See also [`ExponentialKernel`](@ref) for a
+related form of the kernel or [`GammaExponentialKernel`](@ref) for a generalization.
+
+# Examples
+
+```jldoctest; setup = :(using KernelFunctions)
+julia> SqExponentialKernel()
+SqExponentialKernel{Float64,Float64}(1.0)
+
+julia> SqExponentialKernel(2.0f0)
+SqExponentialKernel{Float32,Float32}(2.0)
+
+julia> SqExponentialKernel([2.0,3.0])
+SqExponentialKernel{Float64,Array{Float64}}(1.0)
+```
+"""
+struct SqExponentialKernel{T,Tr<:Transform} <: Kernel{T,Tr}
+    transform::Tr
+    metric::SemiMetric
+    function SqExponentialKernel{T,Tr}(transform::Tr) where {T,Tr<:Transform}
+        return new{T,Tr}(transform,SqEuclidean())
+    end
+end
+
+function SqExponentialKernel(α::T=1.0) where {T<:Real}
+    SqExponentialKernel{T,ScaleTransform{T}}(ScaleTransform(α))
+end
+
+function SqExponentialKernel(α::A) where {A<:AbstractVector{<:Real}}
+    SqExponentialKernel{eltype(A),ScaleTransform{A}}(ScaleTransform(α))
+end
+
+function SqExponentialKernel(t::T) where {T<:Transform}
+    SqExponentialKernel{eltype(t),T}(t)
+end
+
+@inline kappa(κ::SqExponentialKernel, d²::Real) where {T} = exp(-d²)
+
+# function convert(
+#         ::Type{K},
+#         κ::SqExponentialKernel
+#     ) where {K>:SqExponentialKernel{T,A} where {T,A}}
+#     return SqExponentialKernel{T}(T.(κ.α))
+# end

src/kernels/squaredexponential.jl

@@ -3,23 +3,30 @@ using KernelFunctions, Test, Distances
 # test type conversion
 l = 2.0
 vl = [l,l]
+s = ScaleTransform(3.0)
 
-## SquaredExponentialKernel
-@testset "SquaredExponentialKernel" begin
-    @test KernelFunctions.metric(SquaredExponentialKernel(l)) == SqEuclidean()
-    @test KernelFunctions.transform(SquaredExponentialKernel(l)) == ScaleTransform(l)
-    @test KernelFunctions.transform(SquaredExponentialKernel(vl)) == ScaleTransform(vl)
+## SqExponentialKernel
+@testset "SqExponentialKernel" begin
+    @test KernelFunctions.metric(SqExponentialKernel(l)) == SqEuclidean()
+    @test KernelFunctions.transform(SqExponentialKernel(l)) == ScaleTransform(l)
+    @test KernelFunctions.transform(SqExponentialKernel(vl)) == ScaleTransform(vl)
+    @test KernelFunctions.transform(SqExponentialKernel(s)) == s
 end
 
+## MaternKernel
+
 @testset "MaternKernel" begin
     @test KernelFunctions.metric(MaternKernel(l)) == Euclidean()
-    @test KernelFunctions.metric(MaternKernel(l,1.5)) == Euclidean()
-    @test KernelFunctions.metric(MaternKernel(l,2.5)) == Euclidean()
+    @test KernelFunctions.metric(Matern32Kernel(l)) == Euclidean()
+    @test KernelFunctions.metric(Matern52Kernel(l)) == Euclidean()
     @test KernelFunctions.transform(MaternKernel(l)) == ScaleTransform(l)
+    @test KernelFunctions.transform(Matern32Kernel(l)) == ScaleTransform(l)
+    @test KernelFunctions.transform(Matern52Kernel(l)) == ScaleTransform(l)
     @test KernelFunctions.transform(MaternKernel(vl)) == ScaleTransform(vl)
-    @test isa(MaternKernel(),Matern32Kernel)
-    @test isa(MaternKernel(1.0,1.0),MaternKernel)
-    @test isa(MaternKernel(1.0,1.5),Matern32Kernel)
-    @test isa(MaternKernel(1.0,2.5),Matern52Kernel)
-    @test isa(MaternKernel(1.0,Inf),SquaredExponentialKernel)
+    @test KernelFunctions.transform(Matern32Kernel(vl)) == ScaleTransform(vl)
+    @test KernelFunctions.transform(Matern52Kernel(vl)) == ScaleTransform(vl)
+    @test KernelFunctions.transform(MaternKernel(s)) == s
+    @test KernelFunctions.transform(Matern32Kernel(s)) == s
+    @test KernelFunctions.transform(Matern52Kernel(s)) == s
+
 end

test/constructors.jl

@@ -7,7 +7,7 @@ A = rand(dims...)
 B = rand(dims...)
 K = [zeros(dims[1],dims[1]),zeros(dims[2],dims[2])]
 Kdiag = [zeros(dims[1]),zeros(dims[2])]
-kernels = [SquaredExponentialKernel(),MaternKernel(1.0,1.0),Matern32Kernel(),Matern52Kernel()]
+kernels = [SqExponentialKernel(),MaternKernel(),Matern32Kernel(),Matern52Kernel()]
 @testset "Inplace Kernel Matrix" begin
     for k in kernels
         @testset "$k" begin

test/kernelmatrix.jl

@@ -7,7 +7,7 @@ dims = [10,5]
 A = rand(dims...)
 B = rand(dims...)
 K = [zeros(dims[1],dims[1]),zeros(dims[2],dims[2])]
-kernels = [SquaredExponentialKernel,MaternKernel]
+kernels = [SqExponentialKernel,MaternKernel,Matern32Kernel,Matern52Kernel]
 l = 2.0
 vl = l*ones(dims[1])
 testfunction(k,A,B) = det(kernelmatrix(k,A,B))