[WIP] Refactoring of KernelSum and KernelProduct

src/kernels/kernelproduct.jl

@@ -1,61 +1,54 @@
 """
-    KernelProduct(kernels::Array{Kernel})
+    KernelProduct(k1::Kernel, k2::Kernel)
 
 Create a product of kernels.
 One can also use the operator `*` :
 ```
     k1 = SqExponentialKernel()
     k2 = LinearKernel()
-    k = KernelProduct([k1, k2]) == k1 * k2
+    k = KernelProduct(k1, k2) == k1 * k2
     kernelmatrix(k, X) == kernelmatrix(k1, X) .* kernelmatrix(k2, X)
     kernelmatrix(k, X) == kernelmatrix(k1 * k2, X)
 ```
 """
-struct KernelProduct <: Kernel
-    kernels::Vector{Kernel}
+struct KernelProduct{K₁<:Kernel, K₂<:Kernel} <: Kernel
+    κ₁::K₁
+    κ₂::K₂
 end
 
-Base.:*(k1::Kernel,k2::Kernel) = KernelProduct([k1,k2])
-Base.:*(k1::KernelProduct,k2::KernelProduct) = KernelProduct(vcat(k1.kernels,k2.kernels)) #TODO Add test
-Base.:*(k::Kernel,kp::KernelProduct) = KernelProduct(vcat(k,kp.kernels))
-Base.:*(kp::KernelProduct,k::Kernel) = KernelProduct(vcat(kp.kernels,k))
+Base.:*(k1::Kernel, k2::Kernel) = KernelProduct(k1, k2)
 
-Base.length(k::KernelProduct) = length(k.kernels)
-
-kappa(κ::KernelProduct, x ,y) = prod(kappa(k, x, y) for k in κ.kernels)
+kappa(κ::KernelProduct, x ,y) = kappa(κ.κ₁, x, y) * kappa(κ.κ₁, x, y)
 
 hadamard(x,y) = x.*y
 
-function kernelmatrix(
-    κ::KernelProduct,
-    X::AbstractMatrix;
-    obsdim::Int=defaultobs)
-    reduce(hadamard,kernelmatrix(κ.kernels[i],X,obsdim=obsdim) for i in 1:length(κ))
+function kernelmatrix(κ::KernelProduct, X::AbstractMatrix; obsdim::Int = defaultobs)
+    kernelmatrix(κ.κ₁, X, obsdim = obsdim) .* kernelmatrix(κ.κ₂, X, obsdim = obsdim)
 end
 
 function kernelmatrix(
     κ::KernelProduct,
     X::AbstractMatrix,
     Y::AbstractMatrix;
     obsdim::Int=defaultobs)
-    reduce(hadamard,_kernelmatrix(κ.kernels[i],X,Y,obsdim) for i in 1:length(κ))
+    kernelmatrix(κ.κ₁, X, Y, obsdim = obsdim) .* kernelmatrix(κ.κ₂, X, Y, obsdim = obsdim)
 end
 
 function kerneldiagmatrix(
     κ::KernelProduct,
     X::AbstractMatrix;
     obsdim::Int=defaultobs) #TODO Add test
-    reduce(hadamard,kerneldiagmatrix(κ.kernels[i],X,obsdim=obsdim) for i in 1:length(κ))
+    kerneldiagmatrix(κ.κ₁, X, obsdim = obsdim) .* kerneldiagmatrix(κ.κ₂, X, obsdim = obsdim)
 end
 
 function Base.show(io::IO, κ::KernelProduct)
     printshifted(io, κ, 0)
 end
 
 function printshifted(io::IO, κ::KernelProduct, shift::Int)
-    print(io, "Product of $(length(κ)) kernels:")
-    for i in 1:length(κ)
-        print(io, "\n" * ("\t" ^ (shift + 1))* "- ")
-        printshifted(io, κ.kernels[i], shift + 2)
-    end
+    print(io, "Kernel Product :")
+    print(io, "\n" * ("\t" ^ (shift + 1)))
+    printshifted(io, κ.κ₁, shift + 2)
+    print(io, "\n" * ("\t" ^ (shift + 1)))
+    printshifted(io, κ.κ₂, shift + 2)
 end

src/kernels/kernelsum.jl

@@ -1,55 +1,29 @@
 """
-    KernelSum(kernels::Array{Kernel}; weights::Array{Real}=ones(length(kernels)))
+    KernelSum(k1::Kernel, k2::Kernel)
 
 Create a positive weighted sum of kernels. All weights should be positive.
 One can also use the operator `+`
 ```
     k1 = SqExponentialKernel()
     k2 = LinearKernel()
-    k = KernelSum([k1, k2]) == k1 + k2
+    k = KernelSum(k1, k2) == k1 + k2
     kernelmatrix(k, X) == kernelmatrix(k1, X) .+ kernelmatrix(k2, X)
     kernelmatrix(k, X) == kernelmatrix(k1 + k2, X)
-    kweighted = 0.5* k1 + 2.0*k2 == KernelSum([k1, k2], weights = [0.5, 2.0])
 ```
 """
-struct KernelSum <: Kernel
-    kernels::Vector{Kernel}
-    weights::Vector{Real}
+struct KernelSum{K₁<:Kernel, K₂<:Kernel} <: Kernel
+    κ₁::K₁
+    κ₂::K₂
 end
 
-function KernelSum(
-    kernels::AbstractVector{<:Kernel};
-    weights::AbstractVector{<:Real} = ones(Float64, length(kernels)),
-)
-    @assert length(kernels) == length(weights) "Weights and kernel vector should be of the same length"
-    @assert all(weights .>= 0) "All weights should be positive"
-    return KernelSum(kernels, weights)
-end
-
-Base.:+(k1::Kernel, k2::Kernel) = KernelSum([k1, k2], weights = [1.0, 1.0])
-Base.:+(k1::ScaledKernel, k2::ScaledKernel) = KernelSum([kernel(k1), kernel(k2)], weights = [first(k1.σ²), first(k2.σ²)])
-Base.:+(k1::KernelSum, k2::KernelSum) =
-    KernelSum(vcat(k1.kernels, k2.kernels), weights = vcat(k1.weights, k2.weights))
-Base.:+(k::Kernel, ks::KernelSum) =
-    KernelSum(vcat(k, ks.kernels), weights = vcat(1.0, ks.weights))
-Base.:+(k::ScaledKernel, ks::KernelSum) =
-        KernelSum(vcat(kernel(k), ks.kernels), weights = vcat(first(k.σ²), ks.weights))
-Base.:+(k::ScaledKernel, ks::Kernel) =
-        KernelSum(vcat(kernel(k), ks), weights = vcat(first(k.σ²), 1.0))
-Base.:+(ks::KernelSum, k::Kernel) =
-    KernelSum(vcat(ks.kernels, k), weights = vcat(ks.weights, 1.0))
-Base.:+(ks::KernelSum, k::ScaledKernel) =
-        KernelSum(vcat(ks.kernels, kernel(k)), weights = vcat(ks.weights, first(k.σ²)))
-Base.:+(ks::Kernel, k::ScaledKernel) =
-        KernelSum(vcat(ks, kernel(k)), weights = vcat(1.0, first(k.σ²)))
-Base.:*(w::Real, k::KernelSum) = KernelSum(k.kernels, weights = w * k.weights) #TODO add tests
+Base.:+(k1::Kernel, k2::Kernel) = KernelSum(k1, k2)
 
-Base.length(k::KernelSum) = length(k.kernels)
+nmetrics(κ::KernelSum) = metric(κ.κ₁) == metric(κ.κ₂) ? 2 : 1
 
-kappa(κ::KernelSum, x, y) = sum(κ.weights[i] * kappa(κ.kernels[i], x, y) for i in 1:length(κ))
+kappa(κ::KernelSum, x, y) = kappa(κ.κ₁, x, y) + kappa(κ.κ₂, x, y)
 
 function kernelmatrix(κ::KernelSum, X::AbstractMatrix; obsdim::Int = defaultobs)
-    sum(κ.weights[i] * kernelmatrix(κ.kernels[i], X, obsdim = obsdim) for i in 1:length(κ))
+    kernelmatrix(κ.κ₁, X, obsdim = obsdim) + kernelmatrix(κ.κ₂, X, obsdim = obsdim)
 end
 
 function kernelmatrix(
@@ -58,25 +32,25 @@ function kernelmatrix(
     Y::AbstractMatrix;
     obsdim::Int = defaultobs,
 )
-    sum(κ.weights[i] * _kernelmatrix(κ.kernels[i], X, Y, obsdim) for i in 1:length(κ))
+    kernelmatrix(κ.κ₁, X, Y, obsdim = obsdim) + kernelmatrix(κ.κ₂, X, Y, obsdim = obsdim)
 end
 
 function kerneldiagmatrix(
     κ::KernelSum,
     X::AbstractMatrix;
     obsdim::Int = defaultobs,
 )
-    sum(κ.weights[i] * kerneldiagmatrix(κ.kernels[i], X, obsdim = obsdim) for i in 1:length(κ))
+    kerneldiagmatrix(κ.κ₁, X, obsdim = obsdim) + kerneldiagmatrix(κ.κ₂, X, obsdim = obsdim)
 end
 
 function Base.show(io::IO, κ::KernelSum)
     printshifted(io, κ, 0)
 end
 
-function printshifted(io::IO,κ::KernelSum, shift::Int)
-    print(io,"Sum of $(length(κ)) kernels:")
-    for i in 1:length(κ)
-        print(io, "\n" * ("\t" ^ (shift + 1)) * "- (w = $(κ.weights[i])) ")
-        printshifted(io, κ.kernels[i], shift + 2)
-    end
+function printshifted(io::IO, κ::KernelSum, shift::Int)
+    print(io,"Kernel Sum :")
+    print(io, "\n" * ("\t" ^ (shift + 1)))
+    printshifted(io, κ.κ₁, shift + 2)
+    print(io, "\n" * ("\t" ^ (shift + 1)))
+    printshifted(io, κ.κ₂, shift + 2)
 end

test/kernels/kernelproduct.jl

@@ -7,8 +7,7 @@
     k2 = SqExponentialKernel()
     k3 = RationalQuadraticKernel()
 
-    k = KernelProduct([k1,k2])
-    @test length(k) == 2
+    k = KernelProduct(k1, k2)
     @test kappa(k,v1,v2) == kappa(k1*k2,v1,v2)
     @test kappa(k*k,v1,v2) ≈ kappa(k,v1,v2)^2
     @test kappa(k*k3,v1,v2) ≈ kappa(k3*k,v1,v2)

test/kernels/kernelsum.jl

@@ -7,16 +7,15 @@
     k1 = LinearKernel()
     k2 = SqExponentialKernel()
     k3 = RationalQuadraticKernel()
-    X = rand(rng, 2,2)
+    X = rand(rng, 2, 2)
 
     w = [2.0,0.5]
-    k = KernelSum([k1,k2],w)
-    ks1 = 2.0*k1
-    ks2 = 0.5*k2
-    @test length(k) == 2
-    @test kappa(k,v1,v2) == kappa(2.0*k1+0.5*k2,v1,v2)
-    @test kappa(k+k3,v1,v2) ≈ kappa(k3+k,v1,v2)
-    @test kappa(k1+k2,v1,v2) == kappa(KernelSum([k1,k2]),v1,v2)
-    @test kappa(k+ks1,v1,v2) ≈ kappa(ks1+k,v1,v2)
-    @test kappa(k+k,v1,v2) == kappa(KernelSum([k1,k2,k1,k2],vcat(w,w)),v1,v2)
+    k = KernelSum(k1, k2)
+    ks1 = 2.0 * k1
+    ks2 = 0.5 * k2
+    @test kappa(k, v1, v2) == kappa(2.0 * k1 + 0.5 * k2, v1, v2)
+    @test kappa(k + k3, v1, v2) ≈ kappa(k3 + k, v1, v2)
+    @test kappa(k1 + k2, v1, v2) == kappa(k, v1, v2)
+    @test kappa(k + ks1, v1, v2) ≈ kappa(ks1 + k, v1, v2)
+    # @test kappa(k+k,v1,v2) == kappa(KernelSum([k1,k2,k1,k2],vcat(w,w)),v1,v2)
 end