Merge branch 'master-dev'

Author: theogf

.gitignore

@@ -1,2 +1,3 @@
 *.json
 Manifest.toml
+coverage/

src/KernelFunctions.jl

@@ -3,7 +3,6 @@ module KernelFunctions
 export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa, kernelpdmat
 export get_params, set_params!
 
-
 export Kernel
 export ConstantKernel, WhiteKernel, ZeroKernel
 export SqExponentialKernel, ExponentialKernel, GammaExponentialKernel
@@ -13,6 +12,7 @@ export LinearKernel, PolynomialKernel
 export RationalQuadraticKernel, GammaRationalQuadraticKernel
 export KernelSum, KernelProduct
 
+export SelectTransform, ChainTransform, ScaleTransform, LowRankTransform, IdentityTransform, FunctionTransform
 
 
 using Distances, LinearAlgebra
@@ -25,16 +25,21 @@ const defaultobs = 2
 include("utils.jl")
 include("distances/dotproduct.jl")
 include("distances/delta.jl")
-include("transform/transform.jl")
 
 
+"""
+Abstract type defining a slice-wise transformation on an input matrix
+"""
+abstract type Transform end
 abstract type Kernel{T,Tr<:Transform} end
 
+include("transform/transform.jl")
 kernels = ["exponential","matern","polynomial","constant","rationalquad","exponentiated"]
 for k in kernels
     include(joinpath("kernels",k*".jl"))
 end
 include("matrix/kernelmatrix.jl")
+include("matrix/kernelpdmat.jl")
 include("kernels/kernelsum.jl")
 include("kernels/kernelproduct.jl")
 

src/generic.jl

@@ -29,13 +29,6 @@ for kernel in [:ExponentialKernel,:SqExponentialKernel,:Matern32Kernel,:Matern52
     end
 end
 
-function set!(k::Kernel,x)
+function set_params!(k::Kernel,x)
     @error "Setting parameters to this kernel is either not possible or has not been implemented"
 end
-
-set_params!(k::Kernel{T,<:ScaleTransform{<:Base.RefValue{<:Tρ}}},ρ::AbstractVector{<:Tρ}) where {T,Tρ<:Real} = set!(k.transform,ρ[1])
-set_params!(k::Kernel{T,<:ScaleTransform{<:AbstractVector{<:Tρ}}},ρ::AbstractVector{<:Tρ}) where {T,Tρ<:Real} = set!(k.transform,ρ)
-set_params!(k::Kernel{T,<:LowRankTransform{<:AbstractMatrix{<:Tm}}},m::AbstractMatrix{<:Tm}) where {T,Tm<:Real} = set!(k.transform,m)
-
-get_params(k::Kernel{T,<:ScaleTransform{<:Base.RefValue{<:Tρ}}}) where {T,Tρ} = [k.transform.s[]]
-get_params(k::Kernel{T,<:ScaleTransform{<:AbstractVector{<:Tρ}}}) where {T,Tρ} = k.transform.s

src/kernels/kernelproduct.jl

@@ -24,7 +24,7 @@ Base.:*(k::Kernel,kp::KernelProduct) = KernelProduct(vcat(k,kp.kernels))
 Base.:*(kp::KernelProduct,k::Kernel) = KernelProduct(vcat(kp.kernels,k))
 
 Base.length(k::KernelProduct) = length(k.kernels)
-metric(k::KernelProduct) = getmetric.(k.kernels) #TODO Add test
+metric(k::KernelProduct) = metric.(k.kernels) #TODO Add test
 transform(k::KernelProduct) = transform.(k.kernels) #TODO Add test
 transform(k::KernelProduct,x::AbstractVecOrMat) = transform.(k.kernels,[x]) #TODO Add test
 transform(k::KernelProduct,x::AbstractVecOrMat,obsdim::Int) = transform.(k.kernels,[x],obsdim) #TODO Add test

src/matrix/kernelmatrix.jl

@@ -62,6 +62,13 @@ Calculate the kernel matrix of `X` (and `Y`) with respect to kernel `κ`.
 """
 kernelmatrix
 
+function kernelmatrix(
+        κ::Kernel,
+        X::AbstractVector{<:Real};
+        obsdim::Int=defaultobs
+        )
+        kernelmatrix(κ,reshape(X,1,:),obsdim=2)
+end
 
 function kernelmatrix(
         κ::Kernel,

src/matrix/kernelpdmat.jl

@@ -10,11 +10,13 @@ function kernelpdmat(
     K = kernelmatrix(κ,X,obsdim=obsdim)
     Kmax =maximum(K)
     α = eps(eltype(K))
-    while !isposdef(K+αI) && α < 0.01*Kmax
+    while !isposdef(K+α*I) && α < 0.01*Kmax
         α *= 2.0
     end
     if α >= 0.01*Kmax
-        @error "Adding noise on the diagonal was not sufficient to build a positive-definite matrix:\n - Check that your kernel parameters are not extreme\n - Check that your data is sufficiently sparse\n - Maybe use a different kernel"
+        throw(ErrorException("Adding noise on the diagonal was not sufficient to build a positive-definite matrix:\n\t- Check that your kernel parameters are not extreme\n\t- Check that your data is sufficiently sparse\n\t- Maybe use a different kernel"))
     end
-    return PDMat(K+αI)
+    return PDMat(K+α*I)
 end
+
+kernelpdmat(κ::Kernel,X::AbstractVector{<:Real};obsdim=defaultobs) = kernelpdmat(κ,reshape(X,1,:),obsdim=2)

src/transform/lowranktransform.jl

@@ -15,6 +15,10 @@ function set!(t::LowRankTransform{<:AbstractMatrix{T}},M::AbstractMatrix{T}) whe
     @assert size(t) == size(M) "Size of the given matrix $(size(M)) and the projection matrix $(size(t)) are not the same"
     t.proj .= M
 end
+set_params!(k::Kernel{T,<:LowRankTransform{<:AbstractMatrix{<:Tm}}},m::AbstractMatrix{<:Tm}) where {T,Tm<:Real} = set!(k.transform,m)
+
+get_params(k::Kernel{T,<:LowRankTransform}) where {T} = get_params(k.transform)
+get_params(t::LowRankTransform) = t.proj
 
 Base.size(tr::LowRankTransform,i::Int) = size(tr.proj,i)
 Base.size(tr::LowRankTransform) = size(tr.proj) #  TODO Add test
@@ -27,7 +31,7 @@ end
 
 function transform(t::LowRankTransform,x::AbstractVector{<:Real},obsdim::Int=defaultobs) #TODO Add test
     @assert size(t,2) == length(x) "Vector has wrong dimensions $(length(x)) compared to projection matrix"
-    t.proj*X
+    t.proj*x
 end
 
 _transform(t::LowRankTransform,X::AbstractVecOrMat{<:Real},obsdim::Int=defaultobs) = obsdim == 2 ? t.proj * X : X * t.proj'

src/transform/scaletransform.jl

@@ -31,17 +31,23 @@ end
 function set!(t::ScaleTransform{Base.RefValue{T}},ρ::T) where {T<:Real}
     t.s[] = ρ
 end
-
 function set!(t::ScaleTransform{AbstractVector{T}},ρ::AbstractVector{T}) where {T<:Real}
     @assert length(ρ) == dim(t) "Trying to set a vector of size $(length(ρ)) to ScaleTransform of dimension $(dim(t))"
     t.s .= ρ
 end
+set_params!(k::Kernel{T,<:ScaleTransform{<:Base.RefValue{<:Tρ}}},ρ::Tρ) where {T,Tρ<:Real} = set!(k.transform,ρ)
+set_params!(k::Kernel{T,<:ScaleTransform{<:Base.RefValue{<:Tρ}}},ρ::AbstractVector{<:Tρ}) where {T,Tρ<:Real} = set!(k.transform,ρ[1])
+set_params!(k::Kernel{T,<:ScaleTransform{<:AbstractVector{<:Tρ}}},ρ::AbstractVector{<:Tρ}) where {T,Tρ<:Real} = set!(k.transform,ρ)
+
+get_params(k::Kernel{T,<:ScaleTransform}) where {T} = get_params(k.transform)
+get_params(t::ScaleTransform{<:Base.RefValue}) = [t.s[]]
+get_params(t::ScaleTransform{<:AbstractVector}) = t.s
 
 dim(str::ScaleTransform{Base.RefValue{<:Real}}) = 1 #TODO Add test
 dim(str::ScaleTransform{<:AbstractVector{<:Real}}) = length(str.s)
 
 function transform(t::ScaleTransform{<:AbstractVector{<:Real}},X::AbstractMatrix{<:Real},obsdim::Int)
-    @boundscheck if dim(t) != size(X,!Bool(obsdim-1)+1)
+    @boundscheck if dim(t) != size(X,feature_dim(obsdim))
         throw(DimensionMismatch("Array has size $(size(X,!Bool(obsdim-1)+1)) on dimension $(!Bool(obsdim-1)+1)) which does not match the length of the scale transform length , $(dim(t)).")) #TODO Add test
     end
     _transform(t,X,obsdim)

src/transform/selecttransform.jl

@@ -10,29 +10,30 @@ Select the dimensions `dims` that the kernel is applied to.
 """
 struct SelectTransform{T<:AbstractVector{<:Int}} <: Transform
     select::T
-    dim_max::Int
 end
 
-function SelectTransform(dims::AbstractVector{T}) where {T<:Int}
+function SelectTransform(dims::V) where {V<:AbstractVector{T} where  {T<:Int}}
     @assert all(dims.>0) "Selective dimensions should all be positive integers"
-    SelectTransform{T}(dims,maximum(dims))
+    SelectTransform{V}(dims)
 end
 
-function set!(t::SelectTransform{<:AbstractVector{T}},s::AbstractVector{T}) where {T<:Real}
-    t.proj .= s
-end
+get_params(t::SelectTransform) = t.select
+get_params(k::Kernel{T,<:SelectTransform}) where {T} = get_params(k.transform)
+
+set!(t::SelectTransform{<:AbstractVector{T}},dims::AbstractVector{T}) where {T<:Int} = t.select .= dims
+set_params!(k::Kernel{T,<:SelectTransform{Td}},dims::AbstractVector{Td}) where {T,Td<:Int} = set!(k.transform,dims)
 
 Base.maximum(t::SelectTransform) = maximum(t.select)
 
 function transform(t::SelectTransform,X::AbstractMatrix{<:Real},obsdim::Int=defaultobs)
-    @boundscheck t.dim_max <= size(X,feature_dim(obsdim)) ?
-        throw(DimensionMismatch("The highest index $(t.dim_max) is higher then the feature dimension of X : $(size(X,feature_dim(obsdim)))")) : nothing
+    @boundscheck maximum(t) >= size(X,feature_dim(obsdim)) ?
+        throw(DimensionMismatch("The highest index $(maximum(t)) is higher then the feature dimension of X : $(size(X,feature_dim(obsdim)))")) : nothing
     @inbounds _transform(t,X,obsdim)
 end
 
 function transform(t::SelectTransform,x::AbstractVector{<:Real},obsdim::Int=defaultobs) #TODO Add test
-    @assert t.dim_max <= length(x) "The highest index $(t.dim_max) is higher then the vector length : $(length(x))"
-    return x[t.select]
+    @assert maximum(t) <= length(x) "The highest index $(maximum(t)) is higher then the vector length : $(length(x))"
+    return @inbounds x[t.select]
 end
 
 _transform(t::SelectTransform,X::AbstractMatrix{<:Real},obsdim::Int=defaultobs) = obsdim == 2 ? X[t.select,:] : X[:,t.select]

src/transform/transform.jl

@@ -1,12 +1,6 @@
 export Transform, IdentityTransform, ScaleTransform, LowRankTransform, FunctionTransform, ChainTransform
 export transform
 
-"""
-Abstract type defining a slice-wise transformation on an input matrix
-"""
-abstract type Transform end
-
-
 """
 ```julia
     transform(t::Transform, X::AbstractMatrix)

test/test_kernelmatrix.jl

@@ -1,6 +1,7 @@
 using Distances, LinearAlgebra
 using Test
 using KernelFunctions
+using PDMats
 
 dims = [10,5]
 
@@ -35,7 +36,6 @@ k = SqExponentialKernel()
     @testset "KernelSum" begin
         k1 = SqExponentialKernel()
         k2 = LinearKernel()
-        k3 =
         ks = k1 + k2
         w1 = 0.4; w2 = 1.2;
         ks2 = KernelSum([k1,k2],weights=[w1,w2])
@@ -49,12 +49,22 @@ k = SqExponentialKernel()
     @testset "KernelProduct" begin
         k1 = SqExponentialKernel()
         k2 = LinearKernel()
+        k3 = RationalQuadraticKernel()
         kp = k1 * k2
+        kp2 = k1 * k3
+        @test all(KernelFunctions.metric(kp).==[KernelFunctions.metric(k1),KernelFunctions.metric(k2)])
         @test all(kernelmatrix(kp,A) .≈ kernelmatrix(k1,A) .* kernelmatrix(k2,A))
         @test all(kernelmatrix(kp*k1,A) .≈ kernelmatrix(k1,A).^2 .* kernelmatrix(k2,A))
         @test all(kernelmatrix(k1*kp,A) .≈ kernelmatrix(k1,A).^2 .* kernelmatrix(k2,A))
         @test all(kernelmatrix(kp,A) .≈ kernelmatrix(k1,A) .* kernelmatrix(k2,A))
         @test all(kernelmatrix(kp,A,B) .≈ kernelmatrix(k1,A,B) .* kernelmatrix(k2,A,B))
         @test all(kernelmatrix(kp,A) .≈ kernelmatrix(k1,A) .* kernelmatrix(k2,A))
+        @test all(kerneldiagmatrix(kp,A) .== kerneldiagmatrix(k1,A) .* kerneldiagmatrix(k2,A))
+    end
+    @testset "PDMat" begin
+        for obsdim in [1,2]
+            @test all(kernelpdmat(k,A,obsdim=obsdim) .≈ PDMat(kernelmatrix(k,A,obsdim=obsdim)))
+            # @test_throws ErrorException kernelpdmat(k,ones(100,100),obsdim=obsdim)
+        end
     end
 end

test/test_transform.jl

@@ -6,27 +6,58 @@ seed!(42)
 
 dims = (10,5)
 X = rand(dims...)
-##
+x = rand(dims[1])
 s = 3.0
 v1 = vcat(3.0,4.0*ones(dims[2]-1))
 v2 = vcat(3.0,4.0*ones(dims[1]-1))
-t = ScaleTransform(s)
-vt1 = ScaleTransform(v1)
-vt2 = ScaleTransform(v2)
-@test all(KernelFunctions.transform(t,X).==s*X)
-@test all(KernelFunctions.transform(vt1,X,1).==v1'.*X)
-@test all(KernelFunctions.transform(vt2,X,2).==v2.*X)
-##
 P = rand(5,10)
-tp = LowRankTransform(P)
-@test all(KernelFunctions.transform(tp,X,2).==P*X)
-##
+sdims = [1,2,3]
 f(x) = sin.(x)
-tf = FunctionTransform(f)
-KernelFunctions.transform(tf,X,1)
-@test all(KernelFunctions.transform(tf,X,1).==f(X))
-##
-tchain = ChainTransform([t,tp,tf])
-@test all(KernelFunctions.transform(tchain,X,2).==f(P*(s*X)))
-@test all(KernelFunctions.transform(tchain,X,2).==
-            KernelFunctions.transform(tf∘tp∘t,X,2))
+
+@testset "Transform Test" begin
+    ## Test Scale Transform
+    @testset "ScaleTransform" begin
+        t = ScaleTransform(s)
+        vt1 = ScaleTransform(v1)
+        vt2 = ScaleTransform(v2)
+        @test all(KernelFunctions.transform(t,X).==s*X)
+        @test all(KernelFunctions.transform(vt1,X,1).==v1'.*X)
+        @test all(KernelFunctions.transform(vt2,X,2).==v2.*X)
+    end
+    ## Test LowRankTransform
+    @testset "LowRankTransform" begin
+        tp = LowRankTransform(P)
+        @test all(KernelFunctions.transform(tp,X,2).==P*X)
+        @test all(KernelFunctions.transform(tp,x).==P*x)
+        @test all(get_params(SqExponentialKernel(tp)).==P)
+        P2 = rand(5,10)
+        KernelFunctions.set!(tp,P2)
+        @test all(tp.proj.==P2)
+    end
+    ## Test FunctionTransform
+    @testset "FunctionTransform" begin
+        tf = FunctionTransform(f)
+        KernelFunctions.transform(tf,X,1)
+        @test all(KernelFunctions.transform(tf,X,1).==f(X))
+    end
+    ## Test SelectTransform
+    @testset "SelectTransform" begin
+        ts = SelectTransform(sdims)
+        @test all(KernelFunctions.transform(ts,X,2).==X[sdims,:])
+        @test all(KernelFunctions.transform(ts,x).==x[sdims])
+        @test all(get_params(SqExponentialKernel(ts)).==sdims)
+        sdims2 = [2,3,5]
+        KernelFunctions.set!(ts,sdims2)
+        @test all(ts.select.==sdims2)
+    end
+    ## Test ChainTransform
+    @testset "ChainTransform" begin
+        t = ScaleTransform(s)
+        tp = LowRankTransform(P)
+        tf = FunctionTransform(f)
+        tchain = ChainTransform([t,tp,tf])
+        @test all(KernelFunctions.transform(tchain,X,2).==f(P*(s*X)))
+        @test all(KernelFunctions.transform(tchain,X,2).==
+                    KernelFunctions.transform(tf∘tp∘t,X,2))
+    end
+end