Merge remote-tracking branch 'origin/master'

Author: theogf

examples/deepkernellearning.jl

@@ -0,0 +1,35 @@
+using KernelFunctions
+using MLDataUtils
+using Zygote
+using Flux
+using Distributions, LinearAlgebra
+using Plots
+
+Flux.@functor SqExponentialKernel
+Flux.@functor KernelSum
+Flux.@functor Matern32Kernel
+Flux.@functor FunctionTransform
+
+neuralnet = Chain(Dense(1,3),Dense(3,2))
+k = SqExponentialKernel(FunctionTransform(neuralnet))
+xmin = -3; xmax = 3
+x = range(xmin,xmax,length=100)
+x_test = rand(Uniform(xmin,xmax),200)
+x,y = noisy_function(sinc,x;noise=0.1)
+X = reshape(x,:,1)
+λ = [0.1]
+f(x,k,λ) = kernelmatrix(k,X,x,obsdim=1)*inv(kernelmatrix(k,X,obsdim=1)+exp(λ[1])*I)*y
+f(X,k,1.0)
+loss(k,λ) = f(X,k,λ) |> ŷ ->sum(y-ŷ)/length(y)+exp(λ[1])*norm(ŷ)
+loss(k,λ)
+ps = Flux.params(k)
+# push!(ps,λ)
+opt = Flux.Momentum(1.0)
+##
+for i in 1:10
+    grads = Zygote.gradient(()->loss(k,λ),ps)
+    Flux.Optimise.update!(opt,ps,grads)
+    p = Plots.scatter(x,y,lab="data",title="Loss = $(loss(k,λ))")
+    Plots.plot!(x,f(X,k,λ),lab="Prediction",lw=3.0)
+    display(p)
+end

examples/kernelridgeregression.jl

@@ -0,0 +1,35 @@
+using KernelFunctions
+using MLDataUtils
+using Zygote
+using Flux
+using Distributions, LinearAlgebra
+using Plots
+
+Flux.@functor SqExponentialKernel
+Flux.@functor ScaleTransform
+Flux.@functor KernelSum
+Flux.@functor Matern32Kernel
+
+xmin = -3; xmax = 3
+x = range(xmin,xmax,length=100)
+x_test = range(xmin,xmax,length=300)
+x,y = noisy_function(sinc,x;noise=0.1)
+X = reshape(x,:,1)
+X_test = reshape(x_test,:,1)
+k = SqExponentialKernel(1.0)#+Matern32Kernel(2.0)
+λ = [-1.0]
+f(x,k,λ) = kernelmatrix(k,x,X,obsdim=1)*inv(kernelmatrix(k,X,obsdim=1)+exp(λ[1])*I)*y
+f(X,k,1.0)
+loss(k,λ) = f(X,k,λ) |> ŷ ->sum(y-ŷ)/length(y)+exp(λ[1])*norm(ŷ)
+loss(k,λ)
+ps = Flux.params(k)
+push!(ps,λ)
+opt = Flux.Momentum(0.1)
+##
+for i in 1:10
+    grads = Zygote.gradient(()->loss(k,λ),ps)
+    Flux.Optimise.update!(opt,ps,grads)
+    p = Plots.scatter(x,y,lab="data",title="Loss = $(loss(k,λ))")
+    Plots.plot!(x_test,f(X_test,k,λ),lab="Prediction",lw=3.0)
+    display(p)
+end

examples/svm.jl

@@ -0,0 +1,23 @@
+using KernelFunctions
+using MLDataUtils
+using Zygote
+using Flux
+using Distributions, LinearAlgebra
+using Plots
+
+N = 100 #Number of samples
+μ = randn(2,2) # Random Centers
+xgrid = range(-3,3,length=100) # Create a grid
+Xgrid = hcat(collect.(Iterators.product(xgrid,xgrid))...)' #Combine into a 2D grid
+y = rand([-1,1],N) # Select randomly between the two classes
+X = zeros(N,2)
+X[y.==1,:] = rand(MvNormal(μ[:,1],I),count(y.==1))' #Attribute samples from class 1
+X[y.==-1,:] = rand(MvNormal(μ[:,2],I),count(y.==-1))' # Attribute samples from class 2
+
+
+k = SqExponentialKernel(2.0) # Create kernel function
+f(x,k,λ) = kernelmatrix(k,x,X,obsdim=1)*inv(kernelmatrix(k,X,obsdim=1)+exp(λ[1])*I)*y # Optimal prediction f
+svmloss(y,ŷ)= f(X,k,λ) |> ŷ -> sum(maximum.(0.0,1-y*ŷ)) - λ*norm(ŷ) # Total svm loss with regularisation
+pred = f(Xgrid,k,λ) #Compute prediction on a grid
+contourf(xgrid,xgrid,pred)
+scatter!(eachcol(X)...,color=y,lab="data")

src/KernelFunctions.jl

@@ -1,7 +1,7 @@
 module KernelFunctions
 
-export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa, kernelpdmat
-export get_params, set_params!
+export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa, kernelpdmat # Main matrix functions
+export params, duplicate, set! # Helpers
 
 export Kernel
 export ConstantKernel, WhiteKernel, ZeroKernel
@@ -12,8 +12,7 @@ export LinearKernel, PolynomialKernel
 export RationalQuadraticKernel, GammaRationalQuadraticKernel
 export KernelSum, KernelProduct
 
-export SelectTransform, ChainTransform, ScaleTransform, LowRankTransform, IdentityTransform, FunctionTransform
-
+export Transform, SelectTransform, ChainTransform, ScaleTransform, LowRankTransform, IdentityTransform, FunctionTransform
 
 using Distances, LinearAlgebra
 using SpecialFunctions: lgamma, besselk
@@ -42,6 +41,5 @@ include("kernels/kernelsum.jl")
 include("kernels/kernelproduct.jl")
 
 include("generic.jl")
-include("squeeze.jl")
 
 end

src/generic.jl

@@ -28,11 +28,3 @@ for kernel in [:ExponentialKernel,:SqExponentialKernel,:Matern32Kernel,:Matern52
         $kernel(t::Tr) where {Tr<:Transform} = $kernel{eltype(t),Tr}(t)
     end
 end
-
-function set_params!(k::Kernel,x)
-    set!(k.transform,first(x))
-end
-
-
-params(k::Kernel) = (params(k.transform),)
-opt_params(k::Kernel) = (opt_params(k.transform),)

src/kernels/kernelproduct.jl

@@ -20,6 +20,7 @@ end
 
 params(k::KernelProduct) = params.(k.kernels)
 opt_params(k::KernelProduct) = opt_params.(k.kernels)
+duplicate(k::KernelProduct,θ) = KernelProduct(duplicate.(k.kernels,θ))
 
 Base.:*(k1::Kernel,k2::Kernel) = KernelProduct([k1,k2])
 Base.:*(k1::KernelProduct,k2::KernelProduct) = KernelProduct(vcat(k1.kernels,k2.kernels)) #TODO Add test

src/kernels/kernelsum.jl

@@ -27,6 +27,7 @@ end
 
 params(k::KernelSum) = (k.weights,params.(k.kernels))
 opt_params(k::KernelSum) = (k.weights,opt_params.(k.kernels))
+duplicate(k::KernelSum,θ) = KernelSum(duplicate.(k.kernels,θ[end]),weights=first(θ))
 
 Base.:+(k1::Kernel,k2::Kernel) = KernelSum([k1,k2],weights=[1.0,1.0])
 Base.:+(k1::KernelSum,k2::KernelSum) = KernelSum(vcat(k1.kernels,k2.kernels),weights=vcat(k1.weights,k2.weights))

src/squeeze.jl

@@ -26,7 +26,6 @@ function set!(t::ARDTransform{T},ρ::AbstractVector{T}) where {T<:Real}
 end
 
 params(t::ARDTransform) = t.v
-opt_params(t::ARDTransform) = t.v
 dim(t::ARDTransform) = length(t.v)
 
 function transform(t::ARDTransform,X::AbstractMatrix{<:Real},obsdim::Int)

src/transform/ardtransform.jl

@@ -32,7 +32,8 @@ end
 
 set!(t::ChainTransform,θ) = set!.(t.transforms,θ)
 params(t::ChainTransform) = (params.(t.transforms))
-opt_params(t::ChainTransform) = (opt_params.(t.transforms))
+duplicate(t::ChainTransform,θ) = ChainTransform(duplicate.(t.transforms,θ))
+
 
 Base.:∘(t₁::Transform,t₂::Transform) = ChainTransform([t₂,t₁])
 Base.:∘(t::Transform,tc::ChainTransform) = ChainTransform(vcat(tc.transforms,t)) #TODO add test