All needed components for kernel matrices

Author: theogf

src/KernelFunctions.jl

@@ -1,14 +1,18 @@
 module KernelFunctions
 
+export kernelmatrix, kernelmatrix!, kappa
+export Kernel, SquaredExponentialKernel
+
 using Distances, LinearAlgebra
 
 const defaultobs = 2
-abstract type Kernel{T} where {T<:Real} end
+abstract type Kernel{T<:Real} end
 
+include("utils.jl")
+include("common.jl")
 include("kernelmatrix.jl")
-include("kernels/common.jl")
 
-kernels = ("squaredexponential")
+kernels = ["squaredexponential"]
 for k in kernels
     include(joinpath("kernels",k*".jl"))
 end

src/common.jl

@@ -2,7 +2,7 @@
 """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 (κ::K)(d::Real) where {K<:Kernel} = kappa(κ,d)
 
 
-@inline (κ::Kernel)(x::AbstractVector{<:Real},y::AbstractVector{<:Real}) = kappa(κ,evaluate(κ.(metric),x,y))
+# @inline (κ::Kernel)(x::AbstractVector{<:Real},y::AbstractVector{<:Real}) = kappa(κ,evaluate(κ.(metric),x,y))

src/kernelmatrix.jl

@@ -1,44 +1,61 @@
+
+function _kappamatrix!(κ::Kernel{T}, P::AbstractMatrix{T}) where {T<:Real}
+    for i in eachindex(P)
+        @inbounds P[i] = kappa(κ, P[i])
+    end
+    P
+end
+
+function _symmetric_kappamatrix!(
+        κ::Kernel{T},
+        P::AbstractMatrix{T},
+        symmetrize::Bool
+    ) where {T<:Real}
+    if !((n = size(P,1)) == size(P,2))
+        throw(DimensionMismatch("Pairwise matrix must be square."))
+    end
+    for j = 1:n, i = (1:j)
+        @inbounds P[i,j] = kappa(κ, P[i,j])
+    end
+    symmetrize ? LinearAlgebra.copytri!(P, 'U') : P
+end
+
+
 """
 ```
     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},
+        K::Matrix{T₁},
         κ::Kernel{T},
-        X::AbstractMatrix{T},
-        Y::AbstractMatrix{T};
-        obsdim::Integer = defaultobs
-    ) where {T<:Real}
-    basematrix!(σ, K, basefunction(κ), κ.α, X, Y)
-    kappamatrix!(κ, K)
+        X::AbstractMatrix{T₂},
+        Y::AbstractMatrix{T₃};
+        obsdim::Int = defaultobs
+        ) where {T,T₁,T₂,T₃}
+        #TODO Check dimension consistency
+        _kappamatrix!(κ, pairwise!(K,metric(κ), X, Y, dims=obsdim))
 end
 
-function kernelmatrix(
+
+function kernelmatrix!(
+        K::Matrix{T₁},
         κ::Kernel{T},
-        X::AbstractMatrix{T};
+        X::AbstractMatrix{T₂};
         obsdim::Int = defaultobs,
         symmetrize::Bool = true
-    ) where {T<:Real}
-    return symmetric_kappamatrix!(κ,pairwise(basefunction(κ),X,dims=obsdim),symmetrize)
+        ) where {T,T₁<:Real,T₂<:Real}
+        #TODO Check dimension consistency
+        _symmetric_kappamatrix!(κ,pairwise!(metric(κ),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)
-
+```
+    kernel(κ::Kernel, x, y; obsdim=2)
+```
 Apply the kernel `κ` to ``x`` and ``y`` where ``x`` and ``y`` are vectors or scalars of
 some subtype of ``Real``.
 """
@@ -48,11 +65,12 @@ end
 
 function kernel(
         κ::Kernel{T},
-        x::AbstractArray{T1},
-        y::AbstractArray{T2};
+        x::AbstractArray{T₁},
+        y::AbstractArray{T₂};
         obsdim::Int = defaultobs
-    ) where {T,T1<:Real,T2<:Real}
-    kappamatrix!(κ, pairwise(metric(κ),X,Y,dims=obsdim))
+    ) where {T,T₁<:Real,T₂<:Real}
+    # TODO Verify dimensions
+    _kappamatrix!(κ, pairwise(metric(κ),X,Y,dims=obsdim))
 end
 
 """
@@ -63,23 +81,39 @@ Calculate the kernel matrix of `X` with respect to kernel `κ`.
 """
 function kernelmatrix(
         κ::Kernel{T},
-        X::AbstractMatrix{T1};
+        X::AbstractMatrix{T₁};
         obsdim::Int = defaultobs,
         symmetrize::Bool = true
-    ) where {T,T1}
-    return symmetric_kappamatrix!(κ,pairwise(basefunction(κ),X,dims=obsdim),symmetrize)
+    ) where {T,T₁<:Real}
+    return _symmetric_kappamatrix!(κ,pairwise(metric(κ),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};
+        X::AbstractMatrix{T₁},
+        Y::AbstractMatrix{T₂};
         obsdim=defaultobs
-    ) where {T,T1,T2}
-    kappamatrix!(κ, pairwise(basefunction(κ), X, Y, dims=dim(σ)))
+    ) where {T,T₁<:Real,T₂<:Real}
+    _kappamatrix!(κ, pairwise(metric(κ), X, Y, dims=obsdim))
+end
+
+
+"""
+```
+    kerneldiagmatrix(κ::Kernel, X::Matrix; obsdim::Int=2)
+```
+Calculate the diagonal matrix of `X` with respect to kernel `κ`
+"""
+function kerneldiagmatrix(
+        κ::Kernel{T},
+        X::AbstractMatrix{T₁}
+        ) where {T,T₁,T₂}
+        @error "Not implemented yet"
+        #TODO
 end

src/utils.jl

@@ -0,0 +1,18 @@
+## Macro for checking
+macro check_args(K, param, cond, desc=string(cond))
+    quote
+        if !($(esc(cond)))
+            throw(ArgumentError(string(
+                $(string(K)), ": ", $(string(param)), " = ", $(esc(param)), " does not ",
+                "satisfy the constraint ", $(string(desc)), ".")))
+        end
+    end
+end
+
+function promote_float(Tₖ::DataType...)
+    if length(Tₖ) == 0
+        return Float64
+    end
+    T = promote_type(Tₖ...)
+    return T <: Real ? T : Float64
+end