Merge pull request #32 from theogf/remove-transform

Removing transform field and creating TransformedKernel (and ScaledKernel)

Author: theogf

.gitignore

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

Project.toml

@@ -5,6 +5,7 @@ version = "0.2.4"
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -13,7 +14,7 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 ZygoteRules = "700de1a5-db45-46bc-99cf-38207098b444"
 
 [compat]
-Compat = "2.2, 3.2"
+Compat = "2.2, 3"
 Distances = "0.8"
 Requires = "1.0.1"
 SpecialFunctions = "0.8, 0.9, 0.10"

README.md

@@ -13,19 +13,19 @@ The aim is to make the API as model-agnostic as possible while still being user-
 ```julia
   X = reshape(collect(range(-3.0,3.0,length=100)),:,1)
   # Set simple scaling of the data
-  k₁ = SqExponentialKernel(1.0)
+  k₁ = SqExponentialKernel()
   K₁ = kernelmatrix(k₁,X,obsdim=1)
 
   # Set a function transformation on the data
-  k₂ = MaternKernel(FunctionTransform(x->sin.(x)))
+  k₂ = TransformedKernel(Matern32Kernel(),FunctionTransform(x->sin.(x)))
   K₂ = kernelmatrix(k₂,X,obsdim=1)
 
   # Set a matrix premultiplication on the data
-  k₃ = PolynomialKernel(LowRankTransform(randn(4,1)),2.0,0.0)
+  k₃ = transform(PolynomialKernel(c=2.0,d=2.0),LowRankTransform(randn(4,1)))
   K₃ = kernelmatrix(k₃,X,obsdim=1)
 
   # Add and sum kernels
-  k₄ = 0.5*SqExponentialKernel()*LinearKernel(0.5) + 0.4*k₂
+  k₄ = 0.5*SqExponentialKernel()*LinearKernel(c=0.5) + 0.4*k₂
   K₄ = kernelmatrix(k₄,X,obsdim=1)
 
   plot(heatmap.([K₁,K₂,K₃,K₄],yflip=true,colorbar=false)...,layout=(2,2),title=["K₁" "K₂" "K₃" "K₄"])

docs/src/api.md

@@ -14,7 +14,7 @@ CurrentModule = KernelFunctions
 KernelFunctions
 ```
 
-## Kernel Functions
+## Base Kernels
 
 ```@docs
 SqExponentialKernel
@@ -33,9 +33,11 @@ ConstantKernel
 WhiteKernel
 ```
 
-## Kernel Combinations
+## Composite Kernels
 
 ```@docs
+TransformedKernel
+ScaledKernel
 KernelSum
 KernelProduct
 ```

docs/src/kernels.md

@@ -2,6 +2,10 @@
   CurrentModule = KernelFunctions
 ```
 
+# Base Kernels
+
+These are the basic kernels without any transformation of the data. They are the building blocks of KernelFunctions
+
 ## Exponential Kernels
 
 ### Exponential Kernel
@@ -13,7 +17,7 @@ The [Exponential Kernel](@ref ExponentialKernel) is defined as
 
 ### Square Exponential Kernel
 
-The [Square Exponential Kernel](@ref KernelFunctions.SqExponentialKernel) is defined as 
+The [Square Exponential Kernel](@ref KernelFunctions.SqExponentialKernel) is defined as
 ```math
   k(x,x') = \exp\left(-\|x-x'\|^2\right)
 ```
@@ -91,3 +95,13 @@ The [Square Exponential Kernel](@ref KernelFunctions.SqExponentialKernel) is def
 ```math
   k(x,x') = 0
 ```
+
+# Composite Kernels
+
+## TransformedKernel
+
+## ScaledKernel
+
+## KernelSum
+
+## KernelProduct

docs/src/metrics.md

@@ -2,7 +2,15 @@
 
 KernelFunctions.jl relies on [Distances.jl]() for computing the pairwise matrix.
 To do so a distance measure is needed for each kernel. Two very common ones can already be used : `SqEuclidean` and `Euclidean`.
-However all kernels do not rely on distances metrics respecting all the definitions. That's why two additional metrics come with the package : `DotProduct` (`<x,y>`) and `Delta` (`δ(x,y)`). If you want to create a new distance just implement the following :
+However all kernels do not rely on distances metrics respecting all the definitions. That's why two additional metrics come with the package : `DotProduct` (`<x,y>`) and `Delta` (`δ(x,y)`).
+Note that all base kernels must have a defined metric defined as :
+```julia
+    metric(::CustomKernel) = SqEuclidean()
+```
+
+## Adding a new metric
+
+If you want to create a new distance just implement the following :
 
 ```julia
 struct Delta <: Distances.PreMetric

docs/src/transform.md

@@ -1,10 +1,10 @@
 # Transform
 
-`Transform` is the object that takes care of transforming the input data before distances are being computed. It can be as standard as `IdentityTransform` returning the same input, can be a scalar with `ScaleTransform` multiplying the vectors by a scalar or a vector.
+`Transform` is the object that takes care of transforming the input data before distances are being computed. It can be as standard as `IdentityTransform` returning the same input, or multiplying the data by a scalar with `ScaleTransform` or by a vector with `ARDTransform`.
 There is a more general `Transform`: `FunctionTransform` that uses a function and apply it on each vector via `mapslices`.
 You can also create a pipeline of `Transform` via `TransformChain`. For example `LowRankTransform(rand(10,5))∘ScaleTransform(2.0)`.
 
-One apply a transformation on a matrix or a vector via `transform(t::Transform,v::AbstractVecOrMat)`
+One apply a transformation on a matrix or a vector via `KernelFunctions.apply(t::Transform,v::AbstractVecOrMat)`
 
 ## Transforms :
 ```@meta
@@ -14,6 +14,7 @@ CurrentModule = KernelFunctions
 ```@docs
   IdentityTransform
   ScaleTransform
+  ARDTransform
   LowRankTransform
   FunctionTransform
   ChainTransform

src/KernelFunctions.jl

@@ -1,6 +1,7 @@
 module KernelFunctions
 
-export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa # Main matrix functions
+export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa
+export transform
 export params, duplicate, set! # Helpers
 
 export Kernel
@@ -11,6 +12,7 @@ export MaternKernel, Matern32Kernel, Matern52Kernel
 export LinearKernel, PolynomialKernel
 export RationalQuadraticKernel, GammaRationalQuadraticKernel
 export KernelSum, KernelProduct
+export TransformedKernel, ScaledKernel
 
 export Transform, SelectTransform, ChainTransform, ScaleTransform, LowRankTransform, IdentityTransform, FunctionTransform
 
@@ -22,6 +24,7 @@ using Distances, LinearAlgebra
 using SpecialFunctions: logabsgamma, besselk
 using ZygoteRules: @adjoint
 using StatsFuns: logtwo
+using InteractiveUtils: subtypes
 using StatsBase
 
 const defaultobs = 2
@@ -30,7 +33,8 @@ const defaultobs = 2
 Abstract type defining a slice-wise transformation on an input matrix
 """
 abstract type Transform end
-abstract type Kernel{Tr<:Transform} end
+abstract type Kernel end
+abstract type BaseKernel <: Kernel end
 
 include("utils.jl")
 include("distances/dotproduct.jl")
@@ -40,6 +44,8 @@ include("transform/transform.jl")
 for k in ["exponential","matern","polynomial","constant","rationalquad","exponentiated"]
     include(joinpath("kernels",k*".jl"))
 end
+include("kernels/transformedkernel.jl")
+include("kernels/scaledkernel.jl")
 include("matrix/kernelmatrix.jl")
 include("kernels/kernelsum.jl")
 include("kernels/kernelproduct.jl")

src/generic.jl

@@ -1,33 +1,32 @@
-@inline metric(κ::Kernel) = κ.metric
-
 ## Allows to iterate over kernels
 Base.length(::Kernel) = 1
 Base.iterate(k::Kernel) = (k,nothing)
 Base.iterate(k::Kernel, ::Any) = nothing
 
 # default fallback for evaluating a kernel with two arguments (such as vectors etc)
-kappa(κ::Kernel, x, y) = kappa(κ, evaluate(metric(κ), transform(κ, x), transform(κ, y)))
+kappa(κ::Kernel, x, y) = kappa(κ, evaluate(metric(κ), x, y))
+kappa(κ::TransformedKernel, x, y) = kappa(kernel(κ), apply(κ.transform,x), apply(κ.transform,y))
+kappa(κ::TransformedKernel{<:BaseKernel,<:ScaleTransform}, x, y) = kappa(κ, _scale(κ.transform, metric(κ), x, y))
+_scale(t::ScaleTransform, metric::Euclidean, x, y) =  first(t.s) * evaluate(metric, x, y)
+_scale(t::ScaleTransform, metric::Union{SqEuclidean,DotProduct}, x, y) =  first(t.s)^2 * evaluate(metric, x, y)
+_scale(t::ScaleTransform, metric, x, y) = evaluate(metric, apply(t, x), apply(t, y))
+
+printshifted(io::IO,κ::Kernel,shift::Int) = print(io,"$κ")
+Base.show(io::IO,κ::Kernel) = print(io,nameof(typeof(κ)))
 
 ### Syntactic sugar for creating matrices and using kernel functions
-for k in [:ExponentialKernel,:SqExponentialKernel,:GammaExponentialKernel,:MaternKernel,:Matern32Kernel,:Matern52Kernel,:LinearKernel,:PolynomialKernel,:ExponentiatedKernel,:ZeroKernel,:WhiteKernel,:ConstantKernel,:RationalQuadraticKernel,:GammaRationalQuadraticKernel]
+for k in subtypes(BaseKernel)
     @eval begin
         @inline (κ::$k)(d::Real) = kappa(κ,d) #TODO Add test
         @inline (κ::$k)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = kappa(κ, x, y)
-        @inline (κ::$k)(X::AbstractMatrix{T},Y::AbstractMatrix{T};obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ,X,Y,obsdim=obsdim)
-        @inline (κ::$k)(X::AbstractMatrix{T};obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ,X,obsdim=obsdim)
+        @inline (κ::$k)(X::AbstractMatrix{T}, Y::AbstractMatrix{T}; obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ, X, Y, obsdim=obsdim)
+        @inline (κ::$k)(X::AbstractMatrix{T}; obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ, X, obsdim=obsdim)
     end
 end
 
-### Transform generics
-@inline transform(κ::Kernel) = κ.transform
-@inline transform(κ::Kernel, x) = transform(transform(κ), x)
-@inline transform(κ::Kernel, x, obsdim::Int) = transform(transform(κ), x, obsdim)
-
-## Constructors for kernels without parameters
-for kernel in [:ExponentialKernel,:SqExponentialKernel,:Matern32Kernel,:Matern52Kernel,:ExponentiatedKernel]
+for k in nameof.(subtypes(BaseKernel))
     @eval begin
-        $kernel() = $kernel(IdentityTransform())
-        $kernel(ρ::Real) = $kernel(ScaleTransform(ρ))
-        $kernel(ρ::AbstractVector{<:Real}) = $kernel(ARDTransform(ρ))
+        @deprecate($k(ρ::Real;args...),transform($k(args...),ρ))
+        @deprecate($k(ρ::AbstractVector{<:Real};args...),transform($k(args...),ρ))
     end
 end

src/kernels/constant.jl

@@ -1,55 +1,49 @@
 """
-ZeroKernel([tr=IdentityTransform()])
+ZeroKernel()
 
-Create a kernel always returning zero
+Create a kernel that always returning zero
+```
+    κ(x,y) = 0.0
+```
+The output type depends of `x` and `y`
 """
-struct ZeroKernel{Tr} <: Kernel{Tr}
-    transform::Tr
-end
-
-ZeroKernel() = ZeroKernel(IdentityTransform())
+struct ZeroKernel <: BaseKernel end
 
-@inline kappa(κ::ZeroKernel, d::T) where {T<:Real} = zero(T)
+kappa(κ::ZeroKernel, d::T) where {T<:Real} = zero(T)
 
 metric(::ZeroKernel) = Delta()
 
 """
-`WhiteKernel([tr=IdentityTransform()])`
+`WhiteKernel()`
 
 ```
     κ(x,y) = δ(x,y)
 ```
 Kernel function working as an equivalent to add white noise.
 """
-struct WhiteKernel{Tr} <: Kernel{Tr}
-    transform::Tr
-end
-
-WhiteKernel() = WhiteKernel(IdentityTransform())
+struct WhiteKernel <: BaseKernel end
 
-@inline kappa(κ::WhiteKernel,δₓₓ::Real) = δₓₓ
+kappa(κ::WhiteKernel,δₓₓ::Real) = δₓₓ
 
 metric(::WhiteKernel) = Delta()
 
 """
-`ConstantKernel([tr=IdentityTransform(),[c=1.0]])`
+`ConstantKernel(c=1.0)`
 ```
     κ(x,y) = c
 ```
 Kernel function always returning a constant value `c`
 """
-struct ConstantKernel{Tr, Tc<:Real} <: Kernel{Tr}
-    transform::Tr
+struct ConstantKernel{Tc<:Real} <: BaseKernel
     c::Tc
+    function ConstantKernel(;c::T=1.0) where {T<:Real}
+        new{T}(c)
+    end
 end
 
-params(k::ConstantKernel) = (params(k.transform),k.c)
-opt_params(k::ConstantKernel) = (opt_params(k.transform),k.c)
-
-ConstantKernel(c::Real=1.0) = ConstantKernel(IdentityTransform(),c)
-
-ConstantKernel(t::Tr,c::Tc=1.0) where {Tr<:Transform,Tc<:Real} = ConstantKernel{Tr,Tc}(t,c)
+params(k::ConstantKernel) = (k.c,)
+opt_params(k::ConstantKernel) = (k.c,)
 
-@inline kappa(κ::ConstantKernel,x::Real) = κ.c
+kappa(κ::ConstantKernel,x::Real) = κ.c*one(x)
 
 metric(::ConstantKernel) = Delta()

src/kernels/exponential.jl

@@ -1,22 +1,22 @@
 """
-`SqExponentialKernel([ρ=1.0])`
+`SqExponentialKernel()`
 
 The squared exponential kernel is an isotropic Mercer kernel given by the formula:
 ```
-    κ(x,y) = exp(-ρ²‖x-y‖²)
+    κ(x,y) = exp(-‖x-y‖²)
 ```
 See also [`ExponentialKernel`](@ref) for a
 related form of the kernel or [`GammaExponentialKernel`](@ref) for a generalization.
 """
-struct SqExponentialKernel{Tr} <: Kernel{Tr}
-    transform::Tr
-end
+struct SqExponentialKernel <: BaseKernel end
 
-@inline kappa(κ::SqExponentialKernel, d²::Real) = exp(-d²)
-@inline iskroncompatible(::SqExponentialKernel) = true
+kappa(κ::SqExponentialKernel, d²::Real) = exp(-d²)
+iskroncompatible(::SqExponentialKernel) = true
 
 metric(::SqExponentialKernel) = SqEuclidean()
 
+Base.show(io::IO,::SqExponentialKernel) = print(io,"Squared Exponential Kernel")
+
 ## Aliases ##
 const RBFKernel = SqExponentialKernel
 const GaussianKernel = SqExponentialKernel
@@ -28,14 +28,14 @@ The exponential kernel is an isotropic Mercer kernel given by the formula:
     κ(x,y) = exp(-ρ‖x-y‖)
 ```
 """
-struct ExponentialKernel{Tr} <: Kernel{Tr}
-    transform::Tr
-end
+struct ExponentialKernel <: BaseKernel end
 
-@inline kappa(κ::ExponentialKernel, d::Real) = exp(-d)
-@inline iskroncompatible(::ExponentialKernel) = true
+kappa(κ::ExponentialKernel, d::Real) = exp(-d)
+iskroncompatible(::ExponentialKernel) = true
 metric(::ExponentialKernel) = Euclidean()
 
+Base.show(io::IO,::ExponentialKernel) = print(io,"Exponential Kernel")
+
 ## Alias ##
 const LaplacianKernel = ExponentialKernel
 
@@ -46,30 +46,17 @@ The γ-exponential kernel is an isotropic Mercer kernel given by the formula:
     κ(x,y) = exp(-ρ^(2γ)‖x-y‖^(2γ))
 ```
 """
-struct GammaExponentialKernel{Tr, Tγ<:Real} <: Kernel{Tr}
-    transform::Tr
+struct GammaExponentialKernel{Tγ<:Real} <: BaseKernel
     γ::Tγ
-    function GammaExponentialKernel{Tr,Tγ}(t::Tr, γ::Tγ) where {Tr<:Transform,Tγ<:Real}
-        @check_args(GammaExponentialKernel, γ, γ >= zero(Tγ), "γ > 0")
-        return new{Tr, Tγ}(t, γ)
+    function GammaExponentialKernel(;γ::T=2.0) where {T<:Real}
+        @check_args(GammaExponentialKernel, γ, γ >= zero(T), "γ > 0")
+        return new{T}(γ)
     end
 end
 
-params(k::GammaExponentialKernel) = (params(transform),γ)
-opt_params(k::GammaExponentialKernel) = (opt_params(transform),γ)
-
-function GammaExponentialKernel(ρ::Real=1.0, γ::Real=2.0)
-    GammaExponentialKernel(ScaleTransform(ρ), γ)
-end
-
-function GammaExponentialKernel(ρ::AbstractVector{<:Real}, γ::Real=2.0)
-    GammaExponentialKernel(ARDTransform(ρ), γ)
-end
-
-function GammaExponentialKernel(t::Tr, γ::Tγ=2.0) where {Tr<:Transform, Tγ<:Real}
-    GammaExponentialKernel{Tr, Tγ}(t, γ)
-end
+params(k::GammaExponentialKernel) = (γ,)
+opt_params(k::GammaExponentialKernel) = (γ,)
 
-@inline kappa(κ::GammaExponentialKernel, d²::Real) = exp(-d²^κ.γ)
-@inline iskroncompatible(::GammaExponentialKernel) = true
+kappa(κ::GammaExponentialKernel, d²::Real) = exp(-d²^κ.γ)
+iskroncompatible(::GammaExponentialKernel) = true
 metric(::GammaExponentialKernel) = SqEuclidean()