test utils revamp

Project.toml

@@ -1,17 +1,19 @@
 name = "KernelFunctions"
 uuid = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
-version = "0.7.2"
+version = "0.8.0"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 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

@@ -4,6 +4,10 @@ KernelFunctions. [Github](https://github.com/JuliaGaussianProcesses/KernelFuncti
 """
 module KernelFunctions
 
+@warn "SqExponentialKernel changed convention in version 0.8.0. This kernel now divides the
+ squared distance by 2 to better align itself with standard practice. This warning will
+ be removed in 0.9.0."
+
 export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!
 export transform
 export duplicate, set! # Helpers
@@ -31,6 +35,8 @@ export NystromFact, nystrom
 
 export spectral_mixture_kernel, spectral_mixture_product_kernel
 
+export ColVecs, RowVecs
+
 export MOInput
 export IndependentMOKernel, LatentFactorMOKernel
 
@@ -53,35 +59,53 @@ abstract type Kernel end
 abstract type SimpleKernel <: Kernel end
 
 include("utils.jl")
-include("distances/pairwise.jl")
-include("distances/dotproduct.jl")
-include("distances/delta.jl")
-include("distances/sinus.jl")
-include("transform/transform.jl")
-
-for f in readdir(joinpath(@__DIR__, "basekernels"))
- endswith(f, ".jl") && include(joinpath("basekernels", f))
-end
-
-include("kernels/transformedkernel.jl")
-include("kernels/scaledkernel.jl")
-include("matrix/kernelmatrix.jl")
-include("kernels/kernelsum.jl")
-include("kernels/kernelproduct.jl")
-include("kernels/tensorproduct.jl")
-include("approximations/nystrom.jl")
+include(joinpath("distances", "pairwise.jl"))
+include(joinpath("distances", "dotproduct.jl"))
+include(joinpath("distances", "delta.jl"))
+include(joinpath("distances", "sinus.jl"))
+include(joinpath("transform", "transform.jl"))
+
+include(joinpath("basekernels", "constant.jl"))
+include(joinpath("basekernels", "cosine.jl"))
+include(joinpath("basekernels", "exponential.jl"))
+include(joinpath("basekernels", "exponentiated.jl"))
+include(joinpath("basekernels", "fbm.jl"))
+include(joinpath("basekernels", "gabor.jl"))
+include(joinpath("basekernels", "maha.jl"))
+include(joinpath("basekernels", "matern.jl"))
+include(joinpath("basekernels", "nn.jl"))
+include(joinpath("basekernels", "periodic.jl"))
+include(joinpath("basekernels", "piecewisepolynomial.jl"))
+include(joinpath("basekernels", "polynomial.jl"))
+include(joinpath("basekernels", "rationalquad.jl"))
+include(joinpath("basekernels", "sm.jl"))
+include(joinpath("basekernels", "wiener.jl"))
+
+include(joinpath("kernels", "transformedkernel.jl"))
+include(joinpath("kernels", "scaledkernel.jl"))
+include(joinpath("matrix", "kernelmatrix.jl"))
+include(joinpath("kernels", "kernelsum.jl"))
+include(joinpath("kernels", "kernelproduct.jl"))
+include(joinpath("kernels", "tensorproduct.jl"))
+include(joinpath("approximations", "nystrom.jl"))
 include("generic.jl")
 
-include("mokernels/mokernel.jl")
-include("mokernels/moinput.jl")
-include("mokernels/independent.jl")
-include("mokernels/slfm.jl")
+include(joinpath("mokernels", "mokernel.jl"))
+include(joinpath("mokernels", "moinput.jl"))
+include(joinpath("mokernels", "independent.jl"))
+include(joinpath("mokernels", "slfm.jl"))
 
 include("zygote_adjoints.jl")
 
+include("test_utils.jl")
+
 function __init__()
- @require Kronecker="2c470bb0-bcc8-11e8-3dad-c9649493f05e" include("matrix/kernelkroneckermat.jl")
- @require PDMats="90014a1f-27ba-587c-ab20-58faa44d9150" include("matrix/kernelpdmat.jl")
+ @require Kronecker="2c470bb0-bcc8-11e8-3dad-c9649493f05e" begin
+ include(joinpath("matrix", "kernelkroneckermat.jl"))
+ end
+ @require PDMats="90014a1f-27ba-587c-ab20-58faa44d9150" begin
+ include(joinpath("matrix", "kernelpdmat.jl"))
+ end
 end
 
 end

src/basekernels/exponential.jl

@@ -11,7 +11,7 @@ related form of the kernel or [`GammaExponentialKernel`](@ref) for a generalizat
 """
 struct SqExponentialKernel <: SimpleKernel end
 
-kappa(κ::SqExponentialKernel, d²::Real) = exp(-d²)
+kappa(κ::SqExponentialKernel, d²::Real) = exp(-d² / 2)
 
 metric(::SqExponentialKernel) = SqEuclidean()
 
@@ -48,12 +48,15 @@ const LaplacianKernel = ExponentialKernel
 """
  GammaExponentialKernel(; γ = 2.0)
 
-The γ-exponential kernel is an isotropic Mercer kernel given by the formula:
+The γ-exponential kernel [1] is an isotropic Mercer kernel given by the formula:
 ```
- κ(x,y) = exp(-‖x-y‖^(2γ))
+ κ(x,y) = exp(-‖x-y‖^γ)
 ```
 Where `γ > 0`, (the keyword `γ` can be replaced by `gamma`)
-For `γ = 1`, see `SqExponentialKernel` and `γ = 0.5`, see `ExponentialKernel`
+For `γ = 2`, see `SqExponentialKernel` and `γ = 1`, see `ExponentialKernel`.
+
+[1] - Gaussian Processes for Machine Learning, Carl Edward Rasmussen and Christopher K. I.
+ Williams, MIT Press, 2006.
 """
 struct GammaExponentialKernel{Tγ<:Real} <: SimpleKernel
  γ::Vector{Tγ}
@@ -65,10 +68,12 @@ end
 
 @functor GammaExponentialKernel
 
-kappa(κ::GammaExponentialKernel, d²::Real) = exp(-d²^first(κ.γ))
+kappa(κ::GammaExponentialKernel, d::Real) = exp(-d^first(κ.γ))
 
-metric(::GammaExponentialKernel) = SqEuclidean()
+metric(::GammaExponentialKernel) = Euclidean()
 
 iskroncompatible(::GammaExponentialKernel) = true
 
-Base.show(io::IO, κ::GammaExponentialKernel) = print(io, "Gamma Exponential Kernel (γ = ", first(κ.γ), ")")
+function Base.show(io::IO, κ::GammaExponentialKernel)
+ print(io, "Gamma Exponential Kernel (γ = ", first(κ.γ), ")")
+end

src/basekernels/gabor.jl

@@ -11,7 +11,7 @@ struct GaborKernel{K<:Kernel} <: Kernel
  kernel::K
  function GaborKernel(;ell=nothing, p=nothing)
  k = _gabor(ell=ell, p=p)
- new{typeof(k)}(k)
+ return new{typeof(k)}(k)
  end
 end
 
@@ -57,24 +57,17 @@ end
 
 Base.show(io::IO, κ::GaborKernel) = print(io, "Gabor Kernel (ell = ", κ.ell, ", p = ", κ.p, ")")
 
-function kernelmatrix(
- κ::GaborKernel,
- X::AbstractMatrix;
- obsdim::Int=defaultobs)
- kernelmatrix(κ.kernel, X; obsdim=obsdim)
+function kernelmatrix(κ::GaborKernel, X::AbstractMatrix; obsdim::Int=defaultobs)
+ return kernelmatrix(κ.kernel, X; obsdim=obsdim)
 end
 
 function kernelmatrix(
- κ::GaborKernel,
- X::AbstractMatrix,
- Y::AbstractMatrix;
- obsdim::Int=defaultobs)
- kernelmatrix(κ.kernel, X, Y; obsdim=obsdim)
+ κ::GaborKernel, X::AbstractMatrix, Y::AbstractMatrix;
+ obsdim::Int=defaultobs,
+)
+ return kernelmatrix(κ.kernel, X, Y; obsdim=obsdim)
 end
 
-function kerneldiagmatrix(
- κ::GaborKernel,
- X::AbstractMatrix;
- obsdim::Int=defaultobs) #TODO Add test
- kerneldiagmatrix(κ.kernel, X; obsdim=obsdim)
+function kerneldiagmatrix(κ::GaborKernel, X::AbstractMatrix; obsdim::Int=defaultobs) #TODO Add test
+ return kerneldiagmatrix(κ.kernel, X; obsdim=obsdim)
 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/matrix/kernelkroneckermat.jl

@@ -2,22 +2,16 @@ using .Kronecker
 
 export kernelkronmat
 
-function kernelkronmat(
- κ::Kernel,
- X::AbstractVector,
- dims::Int
- )
+function kernelkronmat(κ::Kernel, X::AbstractVector, dims::Int)
  @assert iskroncompatible(κ) "The chosen kernel is not compatible for kroenecker matrices (see [`iskroncompatible`](@ref))"
  k = kernelmatrix(κ, X)
  kronecker(k, dims)
 end
 
 function kernelkronmat(
- κ::Kernel,
- X::AbstractVector{<:AbstractVector};
- obsdim::Int=defaultobs
- )
- @assert iskroncompatible(κ) "The chosen kernel is not compatible for kroenecker matrices"
+ κ::Kernel, X::AbstractVector{<:AbstractVector}; obsdim::Int=defaultobs,
+)
+ @assert iskroncompatible(κ) "The chosen kernel is not compatible for Kronecker matrices"
  Ks = kernelmatrix.(κ, X)
  K = reduce(⊗, Ks)
 end

src/matrix/kernelmatrix.jl

@@ -74,10 +74,7 @@ function kernelmatrix!(K::AbstractMatrix, κ::SimpleKernel, x::AbstractVector)
 end
 
 function kernelmatrix!(
- K::AbstractMatrix,
- κ::SimpleKernel,
- x::AbstractVector,
- y::AbstractVector,
+ K::AbstractMatrix, κ::SimpleKernel, x::AbstractVector, y::AbstractVector,
 )
  validate_inplace_dims(K, x, y)
  pairwise!(K, metric(κ), x, y)
@@ -102,19 +99,13 @@ end
 const defaultobs = 2
 
 function kernelmatrix!(
- K::AbstractMatrix,
- κ::Kernel,
- X::AbstractMatrix;
- obsdim::Int = defaultobs
+ K::AbstractMatrix, κ::Kernel, X::AbstractMatrix; obsdim::Int = defaultobs
 )
  return kernelmatrix!(K, κ, vec_of_vecs(X; obsdim=obsdim))
 end
 
 function kernelmatrix!(
- K::AbstractMatrix,
- κ::Kernel,
- X::AbstractMatrix,
- Y::AbstractMatrix;
+ K::AbstractMatrix, κ::Kernel, X::AbstractMatrix, Y::AbstractMatrix;
  obsdim::Int = defaultobs
 )
  x = vec_of_vecs(X; obsdim=obsdim)
@@ -133,10 +124,7 @@ function kernelmatrix(κ::Kernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim=d
 end
 
 function kerneldiagmatrix!(
- K::AbstractVector,
- κ::Kernel,
- X::AbstractMatrix;
- obsdim::Int = defaultobs
+ K::AbstractVector, κ::Kernel, X::AbstractMatrix; obsdim::Int = defaultobs
 )
  return kerneldiagmatrix!(K, κ, vec_of_vecs(X; obsdim=obsdim))
 end

src/matrix/kernelpdmat.jl

@@ -3,24 +3,28 @@ using .PDMats: PDMat
 export kernelpdmat
 
 """
- Compute a positive-definite matrix in the form of a `PDMat` matrix see [PDMats.jl]() with the cholesky decomposition precomputed
- The algorithm recursively tries to add recursively a diagonal nugget until positive definiteness is achieved or that the noise is too big
+ Compute a positive-definite matrix in the form of a `PDMat` matrix see [PDMats.jl](https://github.com/JuliaStats/PDMats.jl)
+ with the cholesky decomposition precomputed.
+ The algorithm recursively tries to add recursively a diagonal nugget until positive
+ definiteness is achieved or that the noise is too big.
 """
-function kernelpdmat(
- κ::Kernel,
- X::AbstractMatrix;
- obsdim::Int = defaultobs
- )
- K = kernelmatrix(κ,X,obsdim=obsdim)
- Kmax =maximum(K)
+function kernelpdmat(κ::Kernel, X::AbstractMatrix; obsdim::Int=defaultobs)
+ 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
- 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"))
+ if α >= 0.01 * Kmax
+ error(
+ "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)
+function kernelpdmat(κ::Kernel, X::AbstractVector{<:Real}; obsdim=defaultobs)
+ return kernelpdmat(κ, reshape(X, 1, :); obsdim=2)
+end