Merge branch 'master-dev'

Author: theogf

Project.toml

@@ -4,19 +4,23 @@ version = "0.2.0"
 
 [deps]
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+IRTools = "7869d1d1-7146-5819-86e3-90919afe41df"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
+ZygoteRules = "700de1a5-db45-46bc-99cf-38207098b444"
 
 [compat]
 FiniteDifferences = ">= 0.7.2"
 
 [extras]
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["FiniteDifferences", "Random", "Test"]
+test = ["FiniteDifferences", "Tracker", "ForwardDiff", "Random", "Test"]

README.md

@@ -1,7 +1,12 @@
 [![Build Status](https://travis-ci.org/theogf/KernelFunctions.jl.svg?branch=master)](https://travis-ci.org/theogf/AugmentedGaussianProcesses.jl)
+[![Coverage Status](https://coveralls.io/repos/github/theogf/KernelFunctions.jl/badge.svg?branch=master)](https://coveralls.io/github/theogf/KernelFunctions.jl?branch=master)
 [![Documentation](https://img.shields.io/badge/docs-dev-blue.svg)](https://theogf.github.io/KernelFunctions.jl/dev/)
-# KernelFunctions.jl (WIP)
-Julia Package for kernel functions for machine learning
+# KernelFunctions.jl
+## Kernel functions for machine learning
+
+KernelFunctions.jl provide a flexible and complete framework for kernel functions, pretransforming the input data.
+
+The aim is to make the API as model-agnostic as possible while still being user-friendly.
 
 ## Objectives (by priority)
 - ARD Kernels

dev/debugAD.jl

@@ -9,12 +9,12 @@ K = [zeros(dims[1],dims[1]),zeros(dims[2],dims[2])]
 l = 0.1
 vl = l*ones(dims[1])
 testfunction(k,A,B) = det(kernelmatrix(k,A,B))
-testfunction(k,A) = sum(kernelmatrix(k,A))
+testfunction(k,A) = sum(kernelmatrix(k,A,obsdim=2))
 k = MaternKernel(vl)
 KernelFunctions.kappa(k,3)
 testfunction(SqExponentialKernel(vl),A)
 testfunction(MaternKernel(vl),A)
-@which kernelmatrix(MaternKernel(vl),A,B)
+kernelmatrix(MaternKernel(vl),A)
 #For debugging
 @info "Running Zygote gradients"
 Zygote.refresh()
@@ -40,10 +40,10 @@ Zygote.gradient(x->kernelmatrix(MaternKernel(x,1.0),A)[1],l)
 
 @info "Running ForwardDiff gradients"
 ## ForwardDiff
-ForwardDiff.gradient(x->testfunction(SqExponentialKernel(x),A,B),vl) #✓
-ForwardDiff.gradient(x->testfunction(MaternKernel(x),A,B),vl) #✓
 ForwardDiff.gradient(x->testfunction(SqExponentialKernel(x),A),vl) #✓
 ForwardDiff.gradient(x->testfunction(MaternKernel(x),A),vl) #✓
+ForwardDiff.gradient(x->testfunction(SqExponentialKernel(x),A,B),vl) #✓
+ForwardDiff.gradient(x->testfunction(MaternKernel(x),A,B),vl) #✓
 ForwardDiff.gradient(x->testfunction(SqExponentialKernel(x[1]),A,B),[l])
 ForwardDiff.gradient(x->testfunction(MaternKernel(x[1]),A,B),[l])
 ForwardDiff.gradient(x->testfunction(SqExponentialKernel(x[1]),A),[l])

docs/create_kernel_plots.jl

@@ -0,0 +1,49 @@
+using Plots; pyplot();
+using Distributions
+using LinearAlgebra
+using KernelFunctions
+# Translational invariants kernels
+
+default(lw=3.0,titlefontsize=28,tickfontsize=18)
+
+x₀ = 0.0; l=0.1
+n_grid = 101
+fill(x₀,n_grid,1)
+xrange = reshape(collect(range(-3,3,length=n_grid)),:,1)
+for k in [SqExponentialKernel,ExponentialKernel]
+    K = kernelmatrix(k(),xrange,obsdim=1)
+    v = rand(MvNormal(K+1e-7I))
+    plot(xrange,v,lab="",title="f(x)",framestyle=:none) |> display
+    savefig(joinpath(@__DIR__,"src","assets","GP_sample_$(k).png"))
+    plot(xrange,kernel.(k(),x₀,xrange),lab="",ylims=(0,1.1),title="k(0,x)") |> display
+    savefig(joinpath(@__DIR__,"src","assets","kappa_function_$(k).png"))
+end
+
+for k in [GammaExponentialKernel(1.0,1.5)]
+    sparse =1
+    while !isposdef(kernelmatrix(k,xrange*sparse,obsdim=1) + 1e-5I); sparse += 1; end
+    v = rand(MvNormal(kernelmatrix(k,xrange*sparse,obsdim=1)+1e-7I))
+    plot(xrange,v,lab="",title="f(x)",framestyle=:none) |> display
+    savefig(joinpath(@__DIR__,"src","assets","GP_sample_GammaExponentialKernel.png"))
+    plot(xrange,kernel.(k,x₀,xrange),lab="",ylims=(0,1.1),title="k(0,x)") |> display
+    savefig(joinpath(@__DIR__,"src","assets","kappa_function_GammaExponentialKernel.png"))
+end
+
+for k in [MaternKernel,Matern32Kernel,Matern52Kernel]
+    K = kernelmatrix(k(),xrange,obsdim=1)
+    v = rand(MvNormal(K+1e-7I))
+    plot(xrange,v,lab="",title="f(x)",framestyle=:none) |> display
+    savefig(joinpath(@__DIR__,"src","assets","GP_sample_$(k).png"))
+    plot(xrange,kernel.(k(),x₀,xrange),lab="",ylims=(0,1.1),title="k(0,x)") |> display
+    savefig(joinpath(@__DIR__,"src","assets","kappa_function_$(k).png"))
+end
+
+
+for k in [RationalQuadraticKernel]
+    K = kernelmatrix(k(),xrange,obsdim=1)
+    v = rand(MvNormal(K+1e-7I))
+    plot(xrange,v,lab="",title="f(x)",framestyle=:none) |> display
+    savefig(joinpath(@__DIR__,"src","assets","GP_sample_$(k).png"))
+    plot(xrange,kernel.(k(),x₀,xrange),lab="",ylims=(0,1.1),title="k(0,x)") |> display
+    savefig(joinpath(@__DIR__,"src","assets","kappa_function_$(k).png"))
+end

docs/make.jl

@@ -4,7 +4,14 @@ using KernelFunctions
 makedocs(
     sitename = "KernelFunctions",
     format = Documenter.HTML(),
-    modules = [KernelFunctions]
+    modules = [KernelFunctions],
+    pages = ["Home"=>"index.md",
+            "User Guide" => "userguide.md",
+             "Kernel Functions"=>"kernels.md",
+             "Transform"=>"transform.md",
+             "Metrics"=>"metrics.md",
+             "Theory"=>"theory.md",
+             "API"=>"api.md"]
 )
 
 # Documenter can also automatically deploy documentation to gh-pages.

docs/src/api.md

@@ -0,0 +1,67 @@
+# API Library
+
+---
+```@contents
+Pages = ["api.md"]
+```
+
+```@meta
+CurrentModule = KernelFunctions
+```
+
+## Module
+```@docs
+KernelFunctions
+```
+
+## Kernel Functions
+
+```@docs
+SqExponentialKernel
+Exponential
+GammaExponentialKernel
+ExponentiatedKernel
+MaternKernel
+Matern32Kernel
+Matern52Kernel
+LinearKernel
+PolynomialKernel
+RationalQuadraticKernel
+GammaRationalQuadraticKernel
+```
+
+## Kernel Combinations
+
+```@docs
+KernelSum
+KernelProduct
+```
+
+## Transforms
+
+```@docs
+Transform
+IdentityTransform
+ScaleTransform
+LowRankTransform
+FunctionTransform
+```
+
+## Functions
+
+```@docs
+kernelmatrix
+kernelmatrix!
+kerneldiagmatrix
+kerneldiagmatrix!
+transform
+```
+
+
+## Index
+
+```@index
+Pages = ["api.md"]
+Module = ["KernelFunctions"]
+Order = [:type, :function]
+```

docs/src/index.md

@@ -1,5 +1,5 @@
 # KernelFunctions.jl
 
-Documentation for KernelFunctions.jl
+Model agnostic kernel functions compatible with automatic differentiation
 
 *** In Construction ***

docs/src/kernels.md

@@ -0,0 +1,34 @@
+```@meta
+CurrentModule = KernelFunctions
+```
+
+## Exponential Kernels
+
+```@docs
+ExponentialKernel
+SqExponentialKernel
+GammaExponentialKernel
+```
+
+## Matern Kernels
+
+```@docs
+MaternKernel
+Matern32Kernel
+Matern52Kernel
+```
+
+## Polynomial Kernels
+
+```@docs
+LinearKernel
+PolynomialKernel
+```
+
+## Constant Kernels
+
+```@docs
+ConstantKernel
+WhiteKernel
+ZeroKernel
+```

docs/src/metrics.md

@@ -0,0 +1,20 @@
+# Metrics
+
+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 :
+
+```julia
+struct Delta <: Distances.PreMetric
+end
+
+@inline function Distances._evaluate(::Delta,a::AbstractVector{T},b::AbstractVector{T}) where {T}
+    @boundscheck if length(a) != length(b)
+        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
+    end
+    return a==b
+end
+
+@inline (dist::Delta)(a::AbstractArray,b::AbstractArray) = Distances._evaluate(dist,a,b)
+@inline (dist::Delta)(a::Number,b::Number) = a==b
+```

docs/src/theory.md

@@ -0,0 +1,17 @@
+# 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.
+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)`
+
+## Transforms :
+
+```@docs
+  IdentityTransform
+  ScaleTransform
+  LowRankTransform
+  FunctionTransform
+  ChainTransform
+```

docs/src/transform.md

@@ -0,0 +1 @@
+# Building kernel and matrices easily!

docs/src/userguide.md

@@ -1,27 +1,40 @@
 module KernelFunctions
 
-export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa
-export Kernel, SqExponentialKernel, MaternKernel, Matern32Kernel, Matern52Kernel
+export kernel, kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa
+export Kernel
+export ConstantKernel, WhiteKernel, ZeroKernel
+export SqExponentialKernel, ExponentialKernel, GammaExponentialKernel
+export ExponentiatedKernel
+export MaternKernel, Matern32Kernel, Matern52Kernel
+export LinearKernel, PolynomialKernel
+export RationalQuadraticKernel, GammaRationalQuadraticKernel
+export KernelSum, KernelProduct
+
 
-export Transform, ScaleTransform
 
 using Distances, LinearAlgebra
 using Zygote: @adjoint
 using SpecialFunctions: lgamma, besselk
 using StatsFuns: logtwo
 
 const defaultobs = 2
-abstract type Kernel{T,Tr} end
 
-include("zygote_rules.jl")
+# include("zygote_rules.jl")
 include("utils.jl")
+include("distances/dotproduct.jl")
+include("distances/delta.jl")
 include("transform/transform.jl")
-include("kernelmatrix.jl")
 
-kernels = ["sqexponential","matern"]
+
+abstract type Kernel{T,Tr<:Transform} end
+
+kernels = ["exponential","matern","polynomial","constant","rationalquad","exponentiated"]
 for k in kernels
     include(joinpath("kernels",k*".jl"))
 end
+include("kernelmatrix.jl")
+include("kernels/kernelsum.jl")
+include("kernels/kernelproduct.jl")
 
 include("generic.jl")
 

src/KernelFunctions.jl

@@ -0,0 +1,12 @@
+struct Delta <: Distances.PreMetric
+end
+
+@inline function Distances._evaluate(::Delta,a::AbstractVector{T},b::AbstractVector{T}) where {T}
+    @boundscheck if length(a) != length(b)
+        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
+    end
+    return a==b
+end
+
+@inline (dist::Delta)(a::AbstractArray,b::AbstractArray) = Distances._evaluate(dist,a,b)
+@inline (dist::Delta)(a::Number,b::Number) = a==b

src/distances/delta.jl

@@ -0,0 +1,12 @@
+struct DotProduct <: Distances.PreMetric
+end
+
+@inline function Distances._evaluate(::DotProduct,a::AbstractVector{T},b::AbstractVector{T}) where {T}
+    @boundscheck if length(a) != length(b)
+        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
+    end
+    return dot(a,b)
+end
+
+@inline (dist::DotProduct)(a::AbstractArray,b::AbstractArray) = Distances._evaluate(dist,a,b)
+@inline (dist::DotProduct)(a::Number,b::Number) = a*b

src/distances/dotproduct.jl

@@ -1,16 +1,31 @@
-
 @inline metric(κ::Kernel) = κ.metric
-kernels =
-for k in [:SqExponentialKernel,:MaternKernel,:Matern32Kernel,:Matern52Kernel]
-    eval(quote
+
+## Allows to iterate over kernels
+Base.length(::Kernel) = 1
+
+Base.iterate(k::Kernel) = (k,nothing)
+Base.iterate(k::Kernel, ::Any) = nothing
+
+### 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]
+    @eval begin
         @inline (κ::$k)(d::Real) = kappa(κ,d)
-        @inline (κ::$k)(x::AbstractVector{T},y::AbstractVector{T}) where {T} = kernel(κ,evaluate(κ.(metric),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)
-    end)
+        @inline (κ::$k)(x::AbstractVector{<:Real},y::AbstractVector{<:Real}) = kappa(κ,evaluate(κ.metric,transform(κ,x),transform(κ,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)
+    end
 end
-### Transform generics
 
+### Transform generics
 @inline transform(κ::Kernel) = κ.transform
 @inline transform(κ::Kernel,x::AbstractVecOrMat) = transform(κ.transform,x)
 @inline transform(κ::Kernel,x::AbstractVecOrMat,obsdim::Int) = transform(κ.transform,x,obsdim)
+
+## Constructors for kernels without parameters
+for kernel in [:ExponentialKernel,:SqExponentialKernel,:Matern32Kernel,:Matern52Kernel,:ExponentiatedKernel]
+    @eval begin
+        $kernel(ρ::T=1.0) where {T<:Real} =   $kernel{T,ScaleTransform{T}}(ScaleTransform(ρ))
+        $kernel(ρ::A) where {A<:AbstractVector{<:Real}} = $kernel{eltype(A),ScaleTransform{A}}(ScaleTransform(ρ))
+        $kernel(t::Tr) where {Tr<:Transform} = $kernel{eltype(t),Tr}(t)
+    end
+end