WIP : Refactoring with Transform

Author: theogf

src/KernelFunctions.jl

@@ -3,13 +3,16 @@ module KernelFunctions
 export kernelmatrix, kernelmatrix!, kappa
 export Kernel, SquaredExponentialKernel
 
+export Transform, ScaleTransform
+
 using Distances, LinearAlgebra
 
 const defaultobs = 2
-abstract type Kernel{T<:Real} end
+abstract type Kernel{T,Tr} end
 
 include("utils.jl")
 include("common.jl")
+include("transform/transform.jl")
 include("kernelmatrix.jl")
 
 kernels = ["squaredexponential"]

src/common.jl

@@ -4,5 +4,8 @@
 """Apply functions of a kernel on a distance"""
 # @inline (κ::K)(d::Real) where {K<:Kernel} = kappa(κ,d)
 
+@inline transform(κ::Kernel) = κ.transform
+@inline transform(κ::Kernel,x::AbstractVecOrMat) = transform(κ.transform,x)
+@inline transform(κ::Kernel,x::AbstractVecOrMat,obsdim::Int) = transform(κ.transform,x,obsdim)
 
 # @inline (κ::Kernel)(x::AbstractVector{<:Real},y::AbstractVector{<:Real}) = kappa(κ,evaluate(κ.(metric),x,y))

src/kernelmatrix.jl

@@ -1,16 +1,16 @@
 
-function _kappamatrix!(κ::Kernel{T}, P::AbstractMatrix{T₁}) where {T<:Real,T₁<:Real}
+function _kappamatrix!(κ::Kernel, 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},
+        κ::Kernel,
         P::AbstractMatrix{T₁},
         symmetrize::Bool
-    ) where {T<:Real,T₁<:Real}
+    ) where {T₁<:Real}
     if !((n = size(P,1)) == size(P,2))
         throw(DimensionMismatch("Pairwise matrix must be square."))
     end
@@ -38,7 +38,12 @@ function kernelmatrix!(
         _kappamatrix!(κ, pairwise!(K, metric(κ), X, Y, dims=obsdim))
 end
 
-
+"""
+```
+    kernelmatrix!(K::Matrix, κ::Kernel, X::Matrix; obsdim::Integer=2, symmetrize::Bool=true)
+```
+In-place version of `kernelmatrix` where pre-allocated matrix `K` will be overwritten with the kernel matrix.
+"""
 function kernelmatrix!(
         K::Matrix{T₁},
         κ::Kernel{T},
@@ -70,24 +75,33 @@ function kernel(
         obsdim::Int = defaultobs
     ) where {T,T₁<:Real,T₂<:Real}
     # TODO Verify dimensions
-    kappa(κ, evaluate(metric(κ),x,y))
+    kappa(κ, evaluate(metric(κ),transform(κ,x),transform(κ,y)))
 end
 
 """
 ```
-    kernelmatrix(κ::Kernel, X::Matrix ; obsdim::Int=2, symmetrize::Bool)
+    kernelmatrix(κ::Kernel, X::Matrix ; obsdim::Int=2, symmetrize::Bool=true)
 ```
 Calculate the kernel matrix of `X` with respect to kernel `κ`.
 """
 function kernelmatrix(
-        κ::Kernel{T},
-        X::AbstractMatrix{T₁};
+        κ::Kernel{T,<:Transform{A}},
+        X::AbstractMatrix;
         obsdim::Int = defaultobs,
         symmetrize::Bool = true
-    ) where {T,T₁<:Real}
-    Tₛ = typeof(zero(eltype(X))*zero(T))
-    m = size(X,obsdim)
-    return kernelmatrix!(Matrix{promote_float(T,T₁)}(undef,m,m),κ,X,obsdim=obsdim,symmetrize=symmetrize)
+    ) where {T,A}
+    # Tₖ = typeof(zero(eltype(X))*zero(T))
+    # m = size(X,obsdim)
+    K = map(x->kappa(κ,x),pairwise(metric(κ),transform(κ,X,obsdim),dims=obsdim))
+    # K = Matrix{Tₖ}(undef,m,m)
+    # for i in 1:m
+    #     tx = transform(κ,@view X[i,:])
+    #     for j in 1:i
+    #         K[i,j] = kappa(κ,kernel(κ,tx,transform(@view X[j,:])))
+    #     end
+    # end
+    return K
+    # return kernelmatrix!(Matrix{Tₖ}(undef,m,m),κ,X,obsdim=obsdim,symmetrize=symmetrize)
 end
 
 """
@@ -102,10 +116,10 @@ function kernelmatrix(
         Y::AbstractMatrix{T₂};
         obsdim=defaultobs
     ) where {T,T₁<:Real,T₂<:Real}
-    Tₛ = typeof(zero(eltype(X))*zero(eltype(Y))*zero(T))
+    Tₖ = typeof(zero(eltype(X))*zero(eltype(Y))*zero(T))
     m = size(X,obsdim)
     n = size(Y,obsdim)
-    kernelmatrix!(Matrix{Tₛ}(undef,m,n),κ,X,Y,obsdim=obsdim)
+    kernelmatrix!(Matrix{Tₖ}(undef,m,n),κ,X,Y,obsdim=obsdim)
 end
 
 
@@ -117,8 +131,30 @@ Calculate the diagonal matrix of `X` with respect to kernel `κ`
 """
 function kerneldiagmatrix(
         κ::Kernel{T},
-        X::AbstractMatrix{T₁}
+        X::AbstractMatrix{T₁};
+        obsdim::Int = 2
+        ) where {T,T₁}
+        n = size(X,obsdim)
+        Tₖ = typeof(zero(T)*zero(eltype(X)))
+        K = Vector{Tₖ}(undef,n)
+        kerneldiagmatrix!(K,κ,X,obsdim=obsdim)
+        return K
+end
+
+function kerneldiagmatrix!(
+        K::AbstractVector{T₁},
+        κ::Kernel{T},
+        X::AbstractMatrix{T₂};
+        obsdim::Int = 2
         ) where {T,T₁,T₂}
-        @error "Not implemented yet"
-        #TODO
+        if obsdim == 1
+            for i in eachindex(K)
+                @inbounds @views K[i] = kernel(κ, X[i,:],X[i,:])
+            end
+        else
+            for i in eachindex(K)
+                @inbounds @views K[i] = kernel(κ,X[:,i],X[:,i])
+            end
+        end
+        return K
 end

src/kernels/squaredexponential.jl

@@ -1,10 +1,10 @@
-@doc raw"""
+"""
     SquaredExponentialKernel([α=1])
 
 The squared exponential kernel is an isotropic Mercer kernel given by the formula:
 
 ```
-    κ(x,y) = exp(α‖x-y‖²)   α > 0
+    κ(x,y) = exp(-‖x-y‖²)
 ```
 
 where `α` is a positive scaling parameter. See also [`ExponentialKernel`](@ref) for a
@@ -23,29 +23,31 @@ julia> SquaredExponentialKernel([2.0,3.0])
 SquaredExponentialKernel{Float64,Array{Float64}}(1.0)
 ```
 """
-struct SquaredExponentialKernel{T<:Real,A} <: Kernel{T}
-    α::A
+struct SquaredExponentialKernel{T,Tr<:Transform} <: Kernel{T,Tr}
+    transform::Tr
     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
+    function SquaredExponentialKernel{T,Tr}(transform::Tr) where {T,Tr<:Transform}
+        return new{T,Tr}(transform,SqEuclidean())
     end
 end
 
-function SquaredExponentialKernel(α::Union{T,AbstractVector{T}}=1.0) where {T<:Real}
-    SquaredExponentialKernel{promote_float(T)}(α)
+function SquaredExponentialKernel(α::T=1.0) where {T<:Real}
+    SquaredExponentialKernel{T,ScaleTransform{T}}(ScaleTransform(α))
 end
 
-@inline kappa(κ::SquaredExponentialKernel{T,<:Real}, d²::Real) where {T} = exp(-κ.α*d²)
-@inline kappa(κ::SquaredExponentialKernel{T}, d²::Real) where {T} = exp(-d²)
+function SquaredExponentialKernel(α::A) where {A<:AbstractVector{<:Real}}
+    SquaredExponentialKernel{eltype(A),ScaleTransform{A}}(ScaleTransform(α))
+end
 
-function convert(
-        ::Type{K},
-        κ::SquaredExponentialKernel
-    ) where {K>:SquaredExponentialKernel{T,A} where {T,A}}
-    return SquaredExponentialKernel{T}(T.(κ.α))
+function SquaredExponentialKernel(t::T) where {T<:Transform}
+    SquaredExponentialKernel{eltype(t),T}(t)
 end
+
+@inline kappa(κ::SquaredExponentialKernel, d²::Real) where {T} = exp(-d²)
+
+# function convert(
+#         ::Type{K},
+#         κ::SquaredExponentialKernel
+#     ) where {K>:SquaredExponentialKernel{T,A} where {T,A}}
+#     return SquaredExponentialKernel{T}(T.(κ.α))
+# end

src/transform/transform.jl

@@ -0,0 +1,35 @@
+abstract type Transform{T} end
+
+struct TransformChain{T} <: Transform{T}
+end
+
+
+
+struct InputTransform{T} <: Transform{T}
+
+end
+
+struct ScaleTransform{T<:Union{Real,AbstractVector{<:Real}}} <: Transform{T}
+    s::T
+end
+
+
+function ScaleTransform(s::T=1.0) where {T<:Real}
+    @check_args(ScaleTransform, s, s > zero(T), "s > 0")
+    ScaleTransform{T}(s)
+end
+
+function ScaleTransform(s::T,dims::Integer) where {T<:Real}
+    @check_args(ScaleTransform, s, s > zero(T), "s > 0")
+    ScaleTransform{Vector{T}}(fill(s,dims))
+end
+
+function ScaleTransform(s::A) where {A<:AbstractVector{<:Real}}
+    @check_args(ScaleTransform, s, all(s.>zero(eltype(A))), "s > 0")
+    ScaleTransform{A}(s)
+end
+
+transform(t::ScaleTransform{<:AbstractVector{<:Real}},x::AbstractVector{<:Real}) = t.s.*x
+transform(t::ScaleTransform{<:AbstractVector{<:Real}},X::AbstractMatrix{<:Real},obsdim::Int) = obsdim == 1 ? t.s'.*X : t.s.*X
+
+transform(t::ScaleTransform{<:Real},x::AbstractVecOrMat) = t.s*x

test/constructors.jl

@@ -7,7 +7,8 @@ vl = [l,l]
 ## SquaredExponentialKernel
 @testset "SquaredExponentialKernel" begin
     @test KernelFunctions.metric(SquaredExponentialKernel(l)) == SqEuclidean()
-    @test KernelFunctions.metric(SquaredExponentialKernel(vl)) == WeightedSqEuclidean(vl)
+    @test KernelFunctions.transform(SquaredExponentialKernel(l)) == ScaleTransform(l)
+    @test KernelFunctions.transform(SquaredExponentialKernel(vl)) == ScaleTransform(vl)
 end
 
 SquaredExponentialKernel(l)