Added adjoint and kmatrix(X,Y)

Author: theogf

Project.toml

@@ -6,3 +6,4 @@ version = "0.1.0"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

src/KernelFunctions.jl

@@ -6,6 +6,7 @@ export Kernel, SquaredExponentialKernel
 export Transform, ScaleTransform
 
 using Distances, LinearAlgebra
+using Zygote: @adjoint
 
 const defaultobs = 2
 abstract type Kernel{T,Tr} end

src/kernelmatrix.jl

@@ -116,10 +116,18 @@ function kernelmatrix(
         Y::AbstractMatrix{T₂};
         obsdim=defaultobs
     ) where {T,T₁<:Real,T₂<:Real}
-    Tₖ = typeof(zero(eltype(X))*zero(eltype(Y))*zero(T))
-    m = size(X,obsdim)
-    n = size(Y,obsdim)
-    kernelmatrix!(Matrix{Tₖ}(undef,m,n),κ,X,Y,obsdim=obsdim)
+    # Tₖ = typeof(zero(eltype(X))*zero(T))
+    # m = size(X,obsdim)
+    K = map(x->kappa(κ,x),pairwise(metric(κ),transform(κ,X,obsdim),transform(κ,Y,obsdim),dims=obsdim))
+    # K = Matrix{Tₖ}(undef,m,m)
+    # for i in 1:m
+    #     tx = transform(κ,@view X[i,:])
+    #     for j in 1:i
+    #         K[i,j] = kappa(κ,kernel(κ,tx,transform(@view X[j,:])))
+    #     end
+    # end
+    return K
+    # return kernelmatrix!(Matrix{Tₖ}(undef,m,m),κ,X,obsdim=obsdim,symmetrize=symmetrize)
 end
 
 

src/transform/transform.jl

@@ -32,4 +32,12 @@ end
 transform(t::ScaleTransform{<:AbstractVector{<:Real}},x::AbstractVector{<:Real}) = t.s.*x
 transform(t::ScaleTransform{<:AbstractVector{<:Real}},X::AbstractMatrix{<:Real},obsdim::Int) = obsdim == 1 ? t.s'.*X : t.s.*X
 
-transform(t::ScaleTransform{<:Real},x::AbstractVecOrMat) = t.s*x
+transform(t::ScaleTransform{<:Real},x::AbstractVecOrMat,obsdim::Int) = t.s*x
+
+@adjoint transform(t::ScaleTransform{<:AbstractVector{<:Real}},x::AbstractVector{<:Real}) = transform(t,x),Δ->(Δ.*x,t.s.*Δ)
+@adjoint transform(t::ScaleTransform{<:AbstractVector{<:Real}},X::AbstractMatrix{<:Real},obsdim::Int) = transform(t,X,obsdim),Δ->begin
+@show Δ,size(Δ);
+return (obsdim == 1 ? Δ'.*X : Δ.*X,transform(t,Δ,obsdim),nothing)
+end
+
+@adjoint transform(t::ScaleTransform{<:Real},x::AbstractVecOrMat,obsdim::Int) = transform(t,x), Δ->(Δ.s.*x,t.s.*Δ)

test/testAD.jl

@@ -13,7 +13,7 @@ vl = l*ones(dims[1])
 testfunction(k,A,B) = sum(kernelmatrix(k,A,B))
 testfunction(k,A) = sum(kernelmatrix(k,A))
 
-
+testfunction(SquaredExponentialKernel(vl),A)
 #For debugging
 
 ## Zygote