test utils revamp (#159)

* Fix style

* Fix convention

* First pass over test set implementation

* Add standardised tests to BaseKernels

* Test composite kernels

* Fix some tests

* Fix maha

* Fix sm

* Fix up maha

* Remove redundant file

* Move existing test utils over to module

* Add Gamma Exponential kernel reference

* Update src/matrix/kernelpdmat.jl

Co-authored-by: Théo Galy-Fajou <[email protected]>

* Remove repeated code

* Warn about breaking change

* Update src/test_utils.jl

Co-authored-by: David Widmann <[email protected]>

* Bump patch

* Fix up tests

* Remove dead space

* Fix rational quadratic parameter test

* Fix some style issues

* Add extra parameter check

* Update src/basekernels/rationalquad.jl

Co-authored-by: David Widmann <[email protected]>

* Tweak check

* Fix RQ convention to match EQ

* Refactor tests

* Fix nn issues

* Fix weird printing issue

* Update src/test_utils.jl

Co-authored-by: David Widmann <[email protected]>

* Update test/kernels/kernelsum.jl

Co-authored-by: Théo Galy-Fajou <[email protected]>

* Update src/test_utils.jl

Co-authored-by: Théo Galy-Fajou <[email protected]>

* Test FBM kernel

* Fix up Gabor

* Loosen dof bound

* Perturb test

Co-authored-by: Will Tebbutt <[email protected]>
Co-authored-by: Théo Galy-Fajou <[email protected]>
Co-authored-by: David Widmann <[email protected]>

Author: 4 people

Project.toml

@@ -8,10 +8,12 @@ Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 Functors = "d9f16b24-f501-4c13-a1f2-28368ffc5196"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 ZygoteRules = "700de1a5-db45-46bc-99cf-38207098b444"
 
 [compat]

src/KernelFunctions.jl

@@ -44,6 +44,8 @@ export NystromFact, nystrom
 
 export spectral_mixture_kernel, spectral_mixture_product_kernel
 
+export ColVecs, RowVecs
+
 export MOInput
 export IndependentMOKernel, LatentFactorMOKernel
 
@@ -108,6 +110,8 @@ include(joinpath("mokernels", "slfm.jl"))
 
 include("zygote_adjoints.jl")
 
+include("test_utils.jl")
+
 function __init__()
     @require Kronecker="2c470bb0-bcc8-11e8-3dad-c9649493f05e" begin
         include(joinpath("matrix", "kernelkroneckermat.jl"))

src/basekernels/gabor.jl

@@ -57,17 +57,10 @@ end
 
 Base.show(io::IO, κ::GaborKernel) = print(io, "Gabor Kernel (ell = ", κ.ell, ", p = ", κ.p, ")")
 
-function kernelmatrix(κ::GaborKernel, X::AbstractMatrix; obsdim::Int=defaultobs)
-    return kernelmatrix(κ.kernel, X; obsdim=obsdim)
-end
+kernelmatrix(κ::GaborKernel, x::AbstractVector) = kernelmatrix(κ.kernel, x)
 
-function kernelmatrix(
-    κ::GaborKernel, X::AbstractMatrix, Y::AbstractMatrix;
-    obsdim::Int=defaultobs,
-)
-    return kernelmatrix(κ.kernel, X, Y; obsdim=obsdim)
+function kernelmatrix(κ::GaborKernel, x::AbstractVector, y::AbstractVector)
+    return kernelmatrix(κ.kernel, x, y)
 end
 
-function kerneldiagmatrix(κ::GaborKernel, X::AbstractMatrix; obsdim::Int=defaultobs) #TODO Add test
-    return kerneldiagmatrix(κ.kernel, X; obsdim=obsdim)
-end
+kerneldiagmatrix(κ::GaborKernel, x::AbstractVector) = kerneldiagmatrix(κ.kernel, x)

src/basekernels/nn.jl

@@ -42,13 +42,13 @@ function kernelmatrix(::NeuralNetworkKernel, x::RowVecs, y::RowVecs)
     X_2 = sum(x.X .* x.X; dims=2)
     Y_2 = sum(y.X .* y.X; dims=2)
     XY = x.X * y.X'
-    return asin.(XY ./ sqrt.((X_2 .+ 1)' * (Y_2 .+ 1)))
+    return asin.(XY ./ sqrt.((X_2 .+ 1) * (Y_2 .+ 1)'))
 end
 
 function kernelmatrix(::NeuralNetworkKernel, x::RowVecs)
     X_2_1 = sum(x.X .* x.X; dims=2) .+ 1
     XX = x.X * x.X'
-    return asin.(XX ./ sqrt.(X_2_1' * X_2_1))
+    return asin.(XX ./ sqrt.(X_2_1 * X_2_1'))
 end
 
 Base.show(io::IO, κ::NeuralNetworkKernel) = print(io, "Neural Network Kernel")

src/basekernels/periodic.jl

@@ -26,4 +26,6 @@ metric(κ::PeriodicKernel) = Sinus(κ.r)
 
 kappa(κ::PeriodicKernel, d::Real) = exp(- 0.5d)
 
-Base.show(io::IO, κ::PeriodicKernel) = print(io, "Periodic Kernel (length(r) = ", length(κ.r), ")")
+function Base.show(io::IO, κ::PeriodicKernel)
+    print(io, "Periodic Kernel, length(r) = $(length(κ.r))")
+end

src/basekernels/rationalquad.jl

@@ -1,48 +1,62 @@
 """
-    RationalQuadraticKernel(; α = 2.0)
+    RationalQuadraticKernel(; α=2.0)
 
 The rational-quadratic kernel is a Mercer kernel given by the formula:
 ```
-    κ(x,y)=(1+||x−y||²/α)^(-α)
+    κ(x, y) = (1 + ||x − y||² / (2α))^(-α)
 ```
-where `α` is a shape parameter of the Euclidean distance. Check [`GammaRationalQuadraticKernel`](@ref) for a generalization.
+where `α` is a shape parameter of the Euclidean distance. Check
+[`GammaRationalQuadraticKernel`](@ref) for a generalization.
 """
 struct RationalQuadraticKernel{Tα<:Real} <: SimpleKernel
     α::Vector{Tα}
     function RationalQuadraticKernel(;alpha::T=2.0, α::T=alpha) where {T}
-        @check_args(RationalQuadraticKernel, α, α > zero(T), "α > 1")
+        @check_args(RationalQuadraticKernel, α, α > zero(T), "α > 0")
         return new{T}([α])
     end
 end
 
 @functor RationalQuadraticKernel
 
-kappa(κ::RationalQuadraticKernel, d²::T) where {T<:Real} = (one(T)+d²/first(κ.α))^(-first(κ.α))
+function kappa(κ::RationalQuadraticKernel, d²::T) where {T<:Real}
+    return (one(T) + d² / (2 * first(κ.α)))^(-first(κ.α))
+end
+
 metric(::RationalQuadraticKernel) = SqEuclidean()
 
-Base.show(io::IO, κ::RationalQuadraticKernel) = print(io, "Rational Quadratic Kernel (α = ", first(κ.α), ")")
+function Base.show(io::IO, κ::RationalQuadraticKernel)
+    print(io, "Rational Quadratic Kernel (α = $(first(κ.α)))")
+end
 
 """
-`GammaRationalQuadraticKernel([ρ=1.0[,α=2.0[,γ=2.0]]])`
+`GammaRationalQuadraticKernel([α=2.0 [, γ=2.0]])`
+
 The Gamma-rational-quadratic kernel is an isotropic Mercer kernel given by the formula:
 ```
-    κ(x,y)=(1+ρ^(2γ)||x−y||^(2γ)/α)^(-α)
+    κ(x, y) = (1 + ||x−y||^γ / α)^(-α)
 ```
 where `α` is a shape parameter of the Euclidean distance and `γ` is another shape parameter.
 """
 struct GammaRationalQuadraticKernel{Tα<:Real, Tγ<:Real} <: SimpleKernel
     α::Vector{Tα}
     γ::Vector{Tγ}
-    function GammaRationalQuadraticKernel(;alpha::Tα=2.0, gamma::Tγ=2.0, α::Tα=alpha, γ::Tγ=gamma) where {Tα<:Real, Tγ<:Real}
-        @check_args(GammaRationalQuadraticKernel, α, α > one(Tα), "α > 1")
-        @check_args(GammaRationalQuadraticKernel, γ, γ >= one(Tγ), "γ >= 1")
+    function GammaRationalQuadraticKernel(
+        ;alpha::Tα=2.0, gamma::Tγ=2.0, α::Tα=alpha, γ::Tγ=gamma,
+    ) where {Tα<:Real, Tγ<:Real}
+        @check_args(GammaRationalQuadraticKernel, α, α > zero(Tα), "α > 0")
+        @check_args(GammaRationalQuadraticKernel, γ, zero(γ) < γ <= 2, "0 < γ <= 2")
         return new{Tα, Tγ}([α], [γ])
     end
 end
 
 @functor GammaRationalQuadraticKernel
 
-kappa(κ::GammaRationalQuadraticKernel, d²::T) where {T<:Real} = (one(T)+d²^first(κ.γ)/first(κ.α))^(-first(κ.α))
+function kappa(κ::GammaRationalQuadraticKernel, d²::Real)
+    return (one(d²) + d²^(first(κ.γ) / 2) / first(κ.α))^(-first(κ.α))
+end
+
 metric(::GammaRationalQuadraticKernel) = SqEuclidean()
 
-Base.show(io::IO, κ::GammaRationalQuadraticKernel) = print(io, "Gamma Rational Quadratic Kernel (α = ", first(κ.α), ", γ = ", first(κ.γ), ")")
+function Base.show(io::IO, κ::GammaRationalQuadraticKernel)
+    print(io, "Gamma Rational Quadratic Kernel (α = $(first(κ.α)), γ = $(first(κ.γ)))")
+end

src/basekernels/sm.jl

@@ -54,7 +54,7 @@ function spectral_mixture_kernel(
     γs::AbstractMatrix{<:Real},
     ωs::AbstractMatrix{<:Real}
 )
-    spectral_mixture_kernel(SqExponentialKernel(), αs, γs, ωs)
+    return spectral_mixture_kernel(SqExponentialKernel(), αs, γs, ωs)
 end
 
 """
@@ -95,14 +95,14 @@ function spectral_mixture_product_kernel(
         throw(DimensionMismatch("The dimensions of αs, γs, ans ωs do not match"))
     end
     return TensorProduct(spectral_mixture_kernel(h, α, reshape(γ, 1, :), reshape(ω, 1, :))
-               for (α, γ, ω) in zip(eachrow(αs), eachrow(γs), eachrow(ωs)))
+        for (α, γ, ω) in zip(eachrow(αs), eachrow(γs), eachrow(ωs)))
 end
 
 function spectral_mixture_product_kernel(
     αs::AbstractMatrix{<:Real},
     γs::AbstractMatrix{<:Real},
     ωs::AbstractMatrix{<:Real}
 )
-    spectral_mixture_product_kernel(SqExponentialKernel(), αs, γs, ωs)
+    return spectral_mixture_product_kernel(SqExponentialKernel(), αs, γs, ωs)
 end
 

src/test_utils.jl

@@ -0,0 +1,135 @@
+module TestUtils
+
+const __ATOL = 1e-9
+
+using LinearAlgebra
+using KernelFunctions
+using Random
+using Test
+
+"""
+    test_interface(
+        k::Kernel,
+        x0::AbstractVector,
+        x1::AbstractVector,
+        x2::AbstractVector;
+        atol=__ATOL,
+    )
+
+Run various consistency checks on `k` at the inputs `x0`, `x1`, and `x2`.
+`x0` and `x1` should be of the same length with different values, while `x0` and `x2` should
+be of different lengths.
+
+    test_interface([rng::AbstractRNG], k::Kernel, T::Type{<:AbstractVector}; atol=__ATOL)
+
+`test_interface` offers certain types of test data generation to make running these tests
+require less code for common input types. For example, `Vector{<:Real}`, `ColVecs{<:Real}`,
+and `RowVecs{<:Real}` are supported. For other input vector types, please provide the data
+manually. 
+"""
+function test_interface(
+    k::Kernel,
+    x0::AbstractVector,
+    x1::AbstractVector,
+    x2::AbstractVector;
+    atol=__ATOL,
+)
+    # TODO: uncomment the tests of ternary kerneldiagmatrix.
+
+    # Ensure that we have the required inputs.
+    @assert length(x0) == length(x1)
+    @assert length(x0) ≠ length(x2)
+
+    # Check that kerneldiagmatrix basically works.
+    # @test kerneldiagmatrix(k, x0, x1) isa AbstractVector
+    # @test length(kerneldiagmatrix(k, x0, x1)) == length(x0)
+
+    # Check that pairwise basically works.
+    @test kernelmatrix(k, x0, x2) isa AbstractMatrix
+    @test size(kernelmatrix(k, x0, x2)) == (length(x0), length(x2))
+
+    # Check that elementwise is consistent with pairwise.
+    # @test kerneldiagmatrix(k, x0, x1) ≈ diag(kernelmatrix(k, x0, x1)) atol=atol
+
+    # Check additional binary elementwise properties for kernels.
+    # @test kerneldiagmatrix(k, x0, x1) ≈ kerneldiagmatrix(k, x1, x0)
+    @test kernelmatrix(k, x0, x2) ≈ kernelmatrix(k, x2, x0)' atol=atol
+
+    # Check that unary elementwise basically works.
+    @test kerneldiagmatrix(k, x0) isa AbstractVector
+    @test length(kerneldiagmatrix(k, x0)) == length(x0)
+
+    # Check that unary pairwise basically works.
+    @test kernelmatrix(k, x0) isa AbstractMatrix
+    @test size(kernelmatrix(k, x0)) == (length(x0), length(x0))
+    @test kernelmatrix(k, x0) ≈ kernelmatrix(k, x0)' atol=atol
+
+    # Check that unary elementwise is consistent with unary pairwise.
+    @test kerneldiagmatrix(k, x0) ≈ diag(kernelmatrix(k, x0)) atol=atol
+
+    # Check that unary pairwise produces a positive definite matrix (approximately).
+    @test eigmin(Matrix(kernelmatrix(k, x0))) > -atol
+
+    # Check that unary elementwise / pairwise are consistent with the binary versions.
+    # @test kerneldiagmatrix(k, x0) ≈ kerneldiagmatrix(k, x0, x0) atol=atol
+    @test kernelmatrix(k, x0) ≈ kernelmatrix(k, x0, x0) atol=atol
+
+    # Check that basic kernel evaluation succeeds and is consistent with `kernelmatrix`.
+    @test k(first(x0), first(x1)) isa Real
+    @test kernelmatrix(k, x0, x2) ≈ [k(xl, xr) for xl in x0, xr in x2]
+
+    tmp = Matrix{Float64}(undef, length(x0), length(x2))
+    @test kernelmatrix!(tmp, k, x0, x2) ≈ kernelmatrix(k, x0, x2)
+
+    tmp_square = Matrix{Float64}(undef, length(x0), length(x0))
+    @test kernelmatrix!(tmp_square, k, x0) ≈ kernelmatrix(k, x0)
+
+    tmp_diag = Vector{Float64}(undef, length(x0))
+    @test kerneldiagmatrix!(tmp_diag, k, x0) ≈ kerneldiagmatrix(k, x0)
+end
+
+function test_interface(
+    rng::AbstractRNG, k::Kernel, ::Type{Vector{T}}; kwargs...
+) where {T<:Real}
+    test_interface(k, randn(rng, T, 3), randn(rng, T, 3), randn(rng, T, 2); kwargs...)
+end
+
+function test_interface(
+    rng::AbstractRNG, k::Kernel, ::Type{<:ColVecs{T}}; dim_in=2, kwargs...,
+) where {T<:Real}
+    test_interface(
+        k,
+        ColVecs(randn(rng, T, dim_in, 3)),
+        ColVecs(randn(rng, T, dim_in, 3)),
+        ColVecs(randn(rng, T, dim_in, 2));
+        kwargs...,
+    )
+end
+
+function test_interface(
+    rng::AbstractRNG, k::Kernel, ::Type{<:RowVecs{T}}; dim_in=2, kwargs...,
+) where {T<:Real}
+    test_interface(
+        k,
+        RowVecs(randn(rng, T, 3, dim_in)),
+        RowVecs(randn(rng, T, 3, dim_in)),
+        RowVecs(randn(rng, T, 2, dim_in));
+        kwargs...,
+    )
+end
+
+function test_interface(k::Kernel, T::Type{<:AbstractVector}; kwargs...)
+    test_interface(Random.GLOBAL_RNG, k, T; kwargs...)
+end
+
+function test_interface(rng::AbstractRNG, k::Kernel, T::Type{<:Real}; kwargs...)
+    test_interface(rng, k, Vector{T}; kwargs...)
+    test_interface(rng, k, ColVecs{T}; kwargs...)
+    test_interface(rng, k, RowVecs{T}; kwargs...)
+end
+
+function test_interface(k::Kernel, T::Type{<:Real}=Float64; kwargs...)
+    test_interface(Random.GLOBAL_RNG, k, T; kwargs...)
+end
+
+end # module

src/transform/lineartransform.jl

@@ -29,7 +29,7 @@ end
 (t::LinearTransform)(x::Real) = vec(t.A * x)
 (t::LinearTransform)(x::AbstractVector{<:Real}) = t.A * x
 
-_map(t::LinearTransform, x::AbstractVector{<:Real}) = ColVecs(t.A * x')
+_map(t::LinearTransform, x::AbstractVector{<:Real}) = ColVecs(t.A * collect(x'))
 _map(t::LinearTransform, x::ColVecs) = ColVecs(t.A * x.X)
 _map(t::LinearTransform, x::RowVecs) = RowVecs(x.X * t.A')
 

src/utils.jl

@@ -70,6 +70,8 @@ struct RowVecs{T, TX<:AbstractMatrix{T}, S} <: AbstractVector{S}
     end
 end
 
+RowVecs(x::AbstractVector) = RowVecs(reshape(x, :, 1))
+
 Base.size(D::RowVecs) = (size(D.X, 1),)
 Base.getindex(D::RowVecs, i::Int) = view(D.X, i, :)
 Base.getindex(D::RowVecs, i::CartesianIndex{1}) = view(D.X, i, :)