First type and functions

Author: theogf

src/KernelFunctions.jl

@@ -0,0 +1,16 @@
+module KernelFunctions
+
+using Distances, LinearAlgebra
+
+const defaultobs = 2
+abstract type Kernel{T} where {T<:Real} end
+
+include("kernelmatrix.jl")
+include("kernels/common.jl")
+
+kernels = ("squaredexponential")
+for k in kernels
+    include(joinpath("kernels",k*".jl"))
+end
+
+end

src/common.jl

@@ -0,0 +1,8 @@
+
+"""Get method for the kernel metric"""
+@inline metric(κ::Kernel) = κ.metric
+"""Apply functions of a kernel on a distance"""
+@inline (κ::Kernel)(d::Real) = kappa(κ,d)
+
+
+@inline (κ::Kernel)(x::AbstractVector{<:Real},y::AbstractVector{<:Real}) = kappa(κ,evaluate(κ.(metric),x,y))

src/kernelmatrix.jl

@@ -0,0 +1,85 @@
+"""
+```
+    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.
+"""
+function kernelmatrix!(
+        K::Matrix{T},
+        κ::Kernel{T},
+        X::AbstractMatrix{T},
+        Y::AbstractMatrix{T};
+        obsdim::Integer = defaultobs
+    ) where {T<:Real}
+    basematrix!(σ, K, basefunction(κ), κ.α, X, Y)
+    kappamatrix!(κ, K)
+end
+
+function kernelmatrix(
+        κ::Kernel{T},
+        X::AbstractMatrix{T};
+        obsdim::Int = defaultobs,
+        symmetrize::Bool = true
+    ) where {T<:Real}
+    return symmetric_kappamatrix!(κ,pairwise(basefunction(κ),X,dims=obsdim),symmetrize)
+end
+
+function kernelmatrix(
+        κ::Kernel{T},
+        X::AbstractMatrix{T},
+        Y::AbstractMatrix{T};
+        obsdim::Int = defaultobs
+    ) where {T<:Real}
+    kappamatrix!(κ, pairwise(basefunction(κ), X, Y, dims=obsdim))
+end
+
+
+# Convenience Methods ======================================================================
+
+"""
+    kernel(κ::Kernel, x, y)
+
+Apply the kernel `κ` to ``x`` and ``y`` where ``x`` and ``y`` are vectors or scalars of
+some subtype of ``Real``.
+"""
+function kernel(κ::Kernel{T}, x::Real, y::Real) where {T}
+    kernel(κ, T(x), T(y))
+end
+
+function kernel(
+        κ::Kernel{T},
+        x::AbstractArray{T1},
+        y::AbstractArray{T2};
+        obsdim::Int = defaultobs
+    ) where {T,T1<:Real,T2<:Real}
+    kappamatrix!(κ, pairwise(metric(κ),X,Y,dims=obsdim))
+end
+
+"""
+```
+    kernelmatrix(κ::Kernel, X::Matrix ; obsdim::Int=2, symmetrize::Bool)
+```
+Calculate the kernel matrix of `X` with respect to kernel `κ`.
+"""
+function kernelmatrix(
+        κ::Kernel{T},
+        X::AbstractMatrix{T1};
+        obsdim::Int = defaultobs,
+        symmetrize::Bool = true
+    ) where {T,T1}
+    return symmetric_kappamatrix!(κ,pairwise(basefunction(κ),X,dims=obsdim),symmetrize)
+end
+
+"""
+    kernelmatrix(κ::Kernel, X::Matrix, Y::Matrix; obsdim::Int=2)
+
+Calculate the base matrix of `X` and `Y` with respect to kernel `κ`.
+"""
+function kernelmatrix(
+        κ::Kernel{T},
+        X::AbstractMatrix{T1},
+        Y::AbstractMatrix{T2};
+        obsdim=defaultobs
+    ) where {T,T1,T2}
+    kappamatrix!(κ, pairwise(basefunction(κ), X, Y, dims=dim(σ)))
+end

src/kernels/squaredexponential.jl

@@ -0,0 +1,48 @@
+@doc raw"""
+    SquaredExponentialKernel([α=1])
+
+The squared exponential kernel is an isotropic Mercer kernel given by the formula:
+
+```
+    κ(x,y) = exp(α‖x-y‖²)   α > 0
+```
+
+where `α` is a positive scaling parameter. See also [`ExponentialKernel`](@ref) for a
+related form of the kernel or [`GammaExponentialKernel`](@ref) for a generalization.
+
+# Examples
+
+```jldoctest; setup = :(using MLKernels)
+julia> SquaredExponentialKernel()
+SquaredExponentialKernel{Float64}(1.0)
+
+julia> SquaredExponentialKernel(2.0f0)
+SquaredExponentialKernel{Float32}(2.0)
+```
+"""
+struct SquaredExponentialKernel{T<:Real,A} <: Kernel{T}
+    α::A
+    metric::SemiMetric
+    function SquaredExponentialKernel{T}(α::A=T(1)) where {A<:Union{Real,AbstractVector{<:Real}},T<:Real}
+        @check_args(SquaredExponentialKernel, α, all(α .> zero(T)), "α > 0")
+        if A <: Real
+            return new{eltype(A),A}(α,SqEuclidean())
+        else
+            return new{eltype(A),A}(α,WeightedSqEuclidean(α))
+        end
+    end
+end
+
+function SquaredExponentialKernel(α::Union{T,AbstractVector{T}}=1.0) where {T<:Real}
+    SquaredExponentialKernel{promote_float(T)}(α)
+end
+
+@inline kappa(κ::SquaredExponentialKernel{T,<:Real}, d²::T) where {T} = exp(-κ.α*d²)
+@inline kappa(κ::SquaredExponentialKernel{T}, d²::T) where {T} = exp(-d²)
+
+function convert(
+        ::Type{K},
+        κ::SquaredExponentialKernel
+    ) where {K>:SquaredExponentialKernel{T,A} where {T,A}}
+    return SquaredExponentialKernel{T}(T.(κ.α))
+end

test/runtests.jl

@@ -0,0 +1,4 @@
+using Test
+using KernelFunctions
+
+@test 1+1 == 2