Added tests for the new methods. And corrected bugs

Author: theogf

src/matrix/kernelmatrix.jl

@@ -47,10 +47,10 @@ function kernelmatrix!(
     κ::Kernel,
     X::AbstractVector
     )
-    if !check_dims(K, X, X, feature_dim(obsdim), obsdim)
+    if (size(K, 1) != size(K, 2)) || (length(X) != size(K, 1))
         throw(DimensionMismatch("Dimensions of the target array K $(size(K)) are not consistent with X $(size(X))"))
     end
-    map!(κ, K, X, X')
+    K .= κ.(X, X')
 end
 
 ## Wrapper for vector of reals
@@ -98,7 +98,10 @@ function kernelmatrix!(
     X::AbstractVector,
     Y::AbstractVector
     )
-    map!(K, κ, X, Y')
+    if (size(K, 1) != length(X)) || (size(K, 2) != length(Y))
+        throw(DimensionMismatch("Dimensions of the target array K $(size(K)) are not consistent with X $(size(X)) and Y $(size(Y))"))
+    end
+    K .= κ.(X, Y')
 end
 
 """
@@ -162,6 +165,15 @@ function kernelmatrix(
     _kernelmatrix(κ, X, Y, obsdim)
 end
 
+function kernelmatrix(κ::Kernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim::Int = defaultobs)
+    @assert obsdim ∈ [1, 2] "obsdim should be 1 or 2 (see docs of `kernelmatrix`))"
+    if obsdim == 1
+        kernelmatrix(κ, ColVecs(X'), ColVecs(Y'))
+    else
+        kernelmatrix(κ, ColVecs(X), ColVecs(Y))
+    end
+end
+
 @inline _kernelmatrix(κ::SimpleKernel, X, Y, obsdim) =
     map(x -> kappa(κ, x), pairwise(metric(κ), X, Y, dims = obsdim))
 

src/utils.jl

@@ -25,6 +25,7 @@ end
 
 Base.size(D::ColVecs) = (size(D.X, 2),)
 Base.getindex(D::ColVecs, i::Int) = view(D.X, :, i)
+Base.getindex(D::ColVecs, i::CartesianIndex{1}) = view(D.X, :, i)
 Base.getindex(D::ColVecs, i) = ColVecs(view(D.X, :, i))
 
 # Take highest Float among possibilities

test/matrix/kernelmatrix.jl

@@ -2,28 +2,50 @@
 
     rng = MersenneTwister(123456)
     dims = [10,5]
-
-    A = rand(rng, dims...)
-    B = rand(rng, dims...)
+    vA = [rand(rng, dims[1]) for _ in 1:dims[2]]
+    A = hcat(vA...)
+    vB = [rand(rng, dims[1]) for _ in 1:dims[2]]
+    B = hcat(vB...)
+    x = rand(rng, dims[1])
+    X = collect(reshape(x, 1, :))
+    y = rand(rng, dims[2])
+    Y = collect(reshape(y, 1 , :))
+    KX = zeros(dims[1], dims[1])
+    KXY = zeros(dims[1], dims[2])
     C = rand(rng, 8, 9)
     K = [zeros(dims[1],dims[1]),zeros(dims[2],dims[2])]
     Kdiag = [zeros(dims[1]),zeros(dims[2])]
     s = rand(rng)
     k = SqExponentialKernel()
+    struct baseSE <: KernelFunctions.BaseKernel end
+    (k::baseSE)(x, y) = exp(-evaluate(SqEuclidean(), x, y))
+    newk = baseSE()
     kt = transform(SqExponentialKernel(),s)
 
     @testset "Kernel Matrix Operations" begin
         @testset "Inplace Kernel Matrix" begin
+            @test kernelmatrix!(KX, k, x) ≈ kernelmatrix!(KX, k, X)
+            @test kernelmatrix!(KXY, k, x, y) ≈ kernelmatrix!(KXY, k, X, Y)
+            @test kernelmatrix!(K[2], k, vA) ≈ kernelmatrix(k, A) atol = 1e-5
+            @test kernelmatrix!(K[2], k, vA, vB) ≈ kernelmatrix(k, A, B) atol = 1e-5
             for obsdim in [1,2]
+                @show obsdim
                 @test kernelmatrix!(K[obsdim], k, A, B, obsdim = obsdim) == kernelmatrix(k, A, B, obsdim = obsdim)
                 @test kernelmatrix!(K[obsdim], k, A, obsdim = obsdim) == kernelmatrix(k, A, obsdim = obsdim)
                 @test kerneldiagmatrix!(Kdiag[obsdim], k, A, obsdim = obsdim) == kerneldiagmatrix(k, A, obsdim = obsdim)
                 @test_throws DimensionMismatch kernelmatrix!(K[obsdim], k, A, C, obsdim=obsdim)
                 @test_throws DimensionMismatch kernelmatrix!(K[obsdim], k, C, obsdim=obsdim)
                 @test_throws DimensionMismatch kerneldiagmatrix!(Kdiag[obsdim], k, C, obsdim=obsdim)
+                @test kernelmatrix!(K[obsdim], newk, A, B, obsdim = obsdim) ≈ kernelmatrix(k, A, B, obsdim = obsdim)
+                @test kernelmatrix!(K[obsdim], newk, A, obsdim = obsdim) ≈ kernelmatrix(k, A, obsdim = obsdim)
+                @test kerneldiagmatrix!(Kdiag[obsdim], newk, A, obsdim = obsdim) ≈ kerneldiagmatrix(k, A, obsdim = obsdim)
             end
         end
         @testset "Kernel matrix" begin
+            @test kernelmatrix(k, x) ≈ kernelmatrix(k, X)
+            @test kernelmatrix(k, x, y) ≈ kernelmatrix(k, X, Y)
+            @test kernelmatrix(k, vA) ≈ kernelmatrix(k, A) atol = 1e-5
+            @test kernelmatrix(k, vA, vB) ≈ kernelmatrix(k, A, B) atol = 1e-5
             for obsdim in [1,2]
                 @test kernelmatrix(k,A,B,obsdim=obsdim) == kappa.(k,pairwise(KernelFunctions.metric(k),A,B,dims=obsdim))
                 @test kernelmatrix(k,A,obsdim=obsdim) == kappa.(k,pairwise(KernelFunctions.metric(k),A,dims=obsdim))
@@ -32,6 +54,9 @@
                 @test k(A,obsdim=obsdim) == kernelmatrix(k,A,obsdim=obsdim)
                 # @test KernelFunctions._kernel(k,1.0,2.0) == KernelFunctions._kernel(k,[1.0],[2.0])
                 @test_throws DimensionMismatch kernelmatrix(k,A,C,obsdim=obsdim)
+                @test kernelmatrix!(K[obsdim], newk, A, B, obsdim = obsdim) ≈ kernelmatrix(k, A, B, obsdim = obsdim)
+                @test kernelmatrix!(K[obsdim], newk, A, obsdim = obsdim) ≈ kernelmatrix(k, A, obsdim = obsdim)
+                @test kerneldiagmatrix!(Kdiag[obsdim], newk, A, obsdim = obsdim) ≈ kerneldiagmatrix(k, A, obsdim = obsdim)
             end
         end
         @testset "Transformed Kernel Matrix Operations" begin