Merge branch 'master' into compathelper/new_version/2023-09-16-00-52-…

…07-949-02230791788

.github/dependabot.yml

@@ -0,0 +1,7 @@
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"

.github/workflows/CompatHelper.yml

@@ -5,11 +5,30 @@ on:
     - cron: 0 0 * * *
   workflow_dispatch:
 
+permissions:
+  contents: write
+  pull-requests: write
+
 jobs:
   CompatHelper:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
+      - name: Check if Julia is already available in the PATH
+        id: julia_in_path
+        run: which julia
+        continue-on-error: true
+      - name: Install Julia, but only if it is not already available in the PATH
+        uses: julia-actions/setup-julia@v1
+        with:
+          version: '1'
+          arch: ${{ runner.arch }}
+        if: steps.julia_in_path.outcome != 'success'
+      - name: "Add the General registry via Git"
+        run: |
+          import Pkg
+          ENV["JULIA_PKG_SERVER"] = ""
+          Pkg.Registry.add("General")
+        shell: julia --color=yes {0}
       - name: "Install CompatHelper"
         run: |
           import Pkg

.github/workflows/benchmark-comment.yml

@@ -24,7 +24,7 @@ jobs:
       ${{ github.event.workflow_run.event == 'pull_request' &&
           github.event.workflow_run.conclusion == 'success' }}
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
       # restore records from the artifacts
       - uses: dawidd6/action-download-artifact@v2
         with:
@@ -41,7 +41,7 @@ jobs:
           echo ::set-output name=body::$(cat ./pull-request-number.artifact)
       # check if the previous comment exists
       - name: find comment
-        uses: peter-evans/find-comment@v1
+        uses: peter-evans/find-comment@v3
         id: fc
         with:
           issue-number: ${{ steps.output-pull-request-number.outputs.body }}

.github/workflows/benchmark.yml

@@ -16,11 +16,11 @@ jobs:
     if: contains(github.event.pull_request.labels.*.name, 'performance critical')
     steps:
       # setup
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@latest
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v1
         with:
           version: '1.7'
-      - uses: julia-actions/julia-buildpkg@latest
+      - uses: julia-actions/julia-buildpkg@v1
       - name: install dependencies
         run: julia -e 'using Pkg; pkg"add PkgBenchmark [email protected]"'
       # run the benchmark suite
@@ -47,7 +47,7 @@ jobs:
       - name: record pull request number
         run: echo ${{ github.event.pull_request.number }} > ./pull-request-number.artifact
       # save as artifacts (performance tracking (comment) workflow will use it)
-      - uses: actions/upload-artifact@v2
+      - uses: actions/upload-artifact@v4
         with:
           name: Benchmarking
           path: ./*.artifact

.github/workflows/ci.yml

@@ -35,21 +35,12 @@ jobs:
           - 'MultiOutput'
           - 'Others'
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
+      - uses: julia-actions/cache@v1
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         with:

.github/workflows/doc_preview_cleanup.yml

@@ -9,7 +9,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout gh-pages branch
-        uses: actions/checkout@v2
+        uses: actions/checkout@v4
         with:
           ref: gh-pages
       - name: Delete preview and history + push changes

.github/workflows/docs.yml

@@ -17,7 +17,7 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1'

.github/workflows/format.yml

@@ -13,8 +13,8 @@ jobs:
   format:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@latest
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v1
         with:
           version: 1
       - run: |

Project.toml

@@ -1,6 +1,6 @@
 name = "KernelFunctions"
 uuid = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
-version = "0.10.56"
+version = "0.10.63"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -25,13 +25,14 @@ ZygoteRules = "700de1a5-db45-46bc-99cf-38207098b444"
 ChainRulesCore = "1"
 Compat = "3.7, 4"
 CompositionsBase = "0.1"
-Distances = "0.10"
+Distances = "0.10.9"
 FillArrays = "0.10, 0.11, 0.12, 0.13, 1"
 Functors = "0.1, 0.2, 0.3, 0.4"
 IrrationalConstants = "0.1, 0.2"
 LogExpFunctions = "0.2.1, 0.3"
 Requires = "1.0.1"
 SpecialFunctions = "0.8, 0.9, 0.10, 1, 2"
+Statistics = "1"
 StatsBase = "0.32, 0.33, 0.34"
 TensorCore = "0.1"
 ZygoteRules = "0.2"

README.md

@@ -1,6 +1,6 @@
 # KernelFunctions.jl
 
-![CI](https://github.com/JuliaGaussianProcesses/KernelFunctions.jl/workflows/CI/badge.svg?branch=master)
+[![CI](https://github.com/JuliaGaussianProcesses/KernelFunctions.jl/actions/workflows/ci.yml/badge.svg?branch=master)](https://github.com/JuliaGaussianProcesses/KernelFunctions.jl/actions/workflows/ci.yml?query=branch%3Amaster)
 [![codecov](https://codecov.io/gh/JuliaGaussianProcesses/KernelFunctions.jl/branch/master/graph/badge.svg?token=rmDh3gb7hN)](https://codecov.io/gh/JuliaGaussianProcesses/KernelFunctions.jl)
 [![Documentation (stable)](https://img.shields.io/badge/docs-stable-blue.svg)](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable)
 [![Documentation (latest)](https://img.shields.io/badge/docs-dev-blue.svg)](https://juliagaussianprocesses.github.io/KernelFunctions.jl/dev)

docs/Project.toml

@@ -12,4 +12,5 @@ Documenter = "0.27"
 KernelFunctions = "0.10"
 Kronecker = "0.4, 0.5"
 PDMats = "0.11"
+Statistics = "1"
 julia = "1.3"

src/KernelFunctions.jl

@@ -48,7 +48,7 @@ export tensor, ⊗, compose
 
 using Compat
 using ChainRulesCore: ChainRulesCore, Tangent, ZeroTangent, NoTangent
-using ChainRulesCore: @thunk, InplaceableThunk
+using ChainRulesCore: @thunk, InplaceableThunk, ProjectTo, unthunk
 using CompositionsBase
 using Distances
 using FillArrays

src/basekernels/fbm.jl

@@ -57,7 +57,7 @@ end
 
 function kernelmatrix!(K::AbstractMatrix, κ::FBMKernel, x::AbstractVector)
     modx = _mod(x)
-    pairwise!(K, SqEuclidean(), x)
+    pairwise!(SqEuclidean(), K, x)
     K .= _fbm.(modx, modx', K, κ.h)
     return K
 end
@@ -72,7 +72,7 @@ end
 function kernelmatrix!(
     K::AbstractMatrix, κ::FBMKernel, x::AbstractVector, y::AbstractVector
 )
-    pairwise!(K, SqEuclidean(), x, y)
+    pairwise!(SqEuclidean(), K, x, y)
     K .= _fbm.(_mod(x), _mod(y)', K, κ.h)
     return K
 end

src/basekernels/matern.jl

@@ -41,7 +41,8 @@ MaternKernel(; nu::Real=1.5, ν::Real=nu, metric=Euclidean()) = MaternKernel(ν,
 
 function _matern(ν::Real, d::Real)
     if iszero(d)
-        return one(d)
+        c = -ν / (ν - 1)
+        return one(d) + c * d^2 / 2
     else
         y = sqrt(2ν) * d
         b = log(besselk(ν, y))

src/chainrules.jl

@@ -111,36 +111,111 @@ end
 
 function ChainRulesCore.rrule(s::Sinus, x::AbstractVector, y::AbstractVector)
     d = x - y
-    sind = sinpi.(d)
-    abs2_sind_r = abs2.(sind) ./ s.r
+    abs2_sind_r = (sinpi.(d) ./ s.r) .^ 2
     val = sum(abs2_sind_r)
-    gradx = twoπ .* cospi.(d) .* sind ./ (s.r .^ 2)
+    gradx = π .* sinpi.(2 .* d) ./ s.r .^ 2
     function evaluate_pullback(Δ::Any)
-        return (r=-2Δ .* abs2_sind_r,), Δ * gradx, -Δ * gradx
+        r̄ = -2Δ .* abs2_sind_r ./ s.r
+        s̄ = ChainRulesCore.Tangent{typeof(s)}(; r=r̄)
+        return s̄, Δ * gradx, -Δ * gradx
     end
     return val, evaluate_pullback
 end
 
-## Reverse Rules SqMahalanobis
+function ChainRulesCore.rrule(
+    ::typeof(Distances.pairwise), d::Sinus, x::AbstractMatrix; dims=2
+)
+    project_x = ProjectTo(x)
+    function pairwise_pullback(z̄)
+        Δ = unthunk(z̄)
+        n = size(x, dims)
+        x̄ = collect(zero(x))
+        r̄ = zero(d.r)
+        if dims == 1
+            for j in 1:n, i in 1:n
+                xi = view(x, i, :)
+                xj = view(x, j, :)
+                ds = π .* Δ[i, j] .* sinpi.(2 .* (xi .- xj)) ./ d.r .^ 2
+                r̄ .-= 2 .* Δ[i, j] .* sinpi.(xi .- xj) .^ 2 ./ d.r .^ 3
+                x̄[i, :] += ds
+                x̄[j, :] -= ds
+            end
+        elseif dims == 2
+            for j in 1:n, i in 1:n
+                xi = view(x, :, i)
+                xj = view(x, :, j)
+                ds = twoπ .* Δ[i, j] .* sinpi.(xi .- xj) .* cospi.(xi .- xj) ./ d.r .^ 2
+                r̄ .-= 2 .* Δ[i, j] .* sinpi.(xi .- xj) .^ 2 ./ d.r .^ 3
+                x̄[:, i] .+= ds
+                x̄[:, j] .-= ds
+            end
+        end
+        d̄ = ChainRulesCore.Tangent{typeof(d)}(; r=r̄)
+        return NoTangent(), d̄, @thunk(project_x(x̄))
+    end
+    return Distances.pairwise(d, x; dims), pairwise_pullback
+end
 
 function ChainRulesCore.rrule(
-    dist::Distances.SqMahalanobis, a::AbstractVector, b::AbstractVector
+    ::typeof(Distances.pairwise), d::Sinus, x::AbstractMatrix, y::AbstractMatrix; dims=2
 )
-    d = dist(a, b)
-    function SqMahalanobis_pullback(Δ::Real)
-        a_b = a - b
-        ∂qmat = InplaceableThunk(
-            X̄ -> mul!(X̄, a_b, a_b', true, Δ), @thunk((a_b * a_b') * Δ)
-        )
-        ∂a = InplaceableThunk(
-            X̄ -> mul!(X̄, dist.qmat, a_b, true, 2 * Δ), @thunk((2 * Δ) * dist.qmat * a_b)
-        )
-        ∂b = InplaceableThunk(
-            X̄ -> mul!(X̄, dist.qmat, a_b, true, -2 * Δ), @thunk((-2 * Δ) * dist.qmat * a_b)
-        )
-        return Tangent{typeof(dist)}(; qmat=∂qmat), ∂a, ∂b
+    project_x = ProjectTo(x)
+    project_y = ProjectTo(y)
+    function pairwise_pullback(z̄)
+        Δ = unthunk(z̄)
+        n = size(x, dims)
+        m = size(y, dims)
+        x̄ = collect(zero(x))
+        ȳ = collect(zero(y))
+        r̄ = zero(d.r)
+        if dims == 1
+            for j in 1:m, i in 1:n
+                xi = view(x, i, :)
+                yj = view(y, j, :)
+                ds = π .* Δ[i, j] .* sinpi.(2 .* (xi .- yj)) ./ d.r .^ 2
+                r̄ .-= 2 .* Δ[i, j] .* sinpi.(xi .- yj) .^ 2 ./ d.r .^ 3
+                x̄[i, :] .+= ds
+                ȳ[j, :] .-= ds
+            end
+        elseif dims == 2
+            for j in 1:m, i in 1:n
+                xi = view(x, :, i)
+                yj = view(y, :, j)
+                ds = π .* Δ[i, j] .* sinpi.(2 .* (xi .- yj)) ./ d.r .^ 2
+                r̄ .-= 2 .* Δ[i, j] .* sinpi.(xi .- yj) .^ 2 ./ d.r .^ 3
+                x̄[:, i] .+= ds
+                ȳ[:, j] .-= ds
+            end
+        end
+        d̄ = ChainRulesCore.Tangent{typeof(d)}(; r=r̄)
+        return NoTangent(), d̄, @thunk(project_x(x̄)), @thunk(project_y(ȳ))
+    end
+    return Distances.pairwise(d, x, y; dims), pairwise_pullback
+end
+
+function ChainRulesCore.rrule(
+    ::typeof(Distances.colwise), d::Sinus, x::AbstractMatrix, y::AbstractMatrix
+)
+    project_x = ProjectTo(x)
+    project_y = ProjectTo(y)
+    function colwise_pullback(z̄)
+        Δ = unthunk(z̄)
+        n = size(x, 2)
+        x̄ = collect(zero(x))
+        ȳ = collect(zero(y))
+        r̄ = zero(d.r)
+        for i in 1:n
+            xi = view(x, :, i)
+            yi = view(y, :, i)
+            ds = π .* Δ[i] .* sinpi.(2 .* (xi .- yi)) ./ d.r .^ 2
+            r̄ .-= 2 .* Δ[i] .* sinpi.(xi .- yi) .^ 2 ./ d.r .^ 3
+            x̄[:, i] .+= ds
+            ȳ[:, i] .-= ds
+        end
+        d̄ = ChainRulesCore.Tangent{typeof(d)}(; r=r̄)
+        return NoTangent(), d̄, @thunk(project_x(x̄)), @thunk(project_y(ȳ))
     end
-    return d, SqMahalanobis_pullback
+    return Distances.colwise(d, x, y), colwise_pullback
 end
 
 ## Reverse Rules for matrix wrappers
@@ -150,8 +225,11 @@ function ChainRulesCore.rrule(::Type{<:ColVecs}, X::AbstractMatrix)
     function ColVecs_pullback(::AbstractVector{<:AbstractVector{<:Real}})
         return error(
             "Pullback on AbstractVector{<:AbstractVector}.\n" *
-            "This might happen if you try to use gradients on the generic `kernelmatrix` or `kernelmatrix_diag`.\n" *
-            "To solve this issue overload `kernelmatrix(_diag)` for your kernel for `ColVecs`",
+            "This might happen if you try to use gradients on the generic `kernelmatrix` or `kernelmatrix_diag`,\n" *
+            "or because some external computation has acted on `ColVecs` to produce a vector of vectors." *
+            "In the former case, to solve this issue overload `kernelmatrix(_diag)` for your kernel for `ColVecs`." *
+            "In the latter case, one needs to track down the `rrule` whose pullback returns a `Vector{Vector{T}}`," *
+            " rather than a `Tangent`, as the cotangent / gradient for `ColVecs` input, and circumvent it."
         )
     end
     return ColVecs(X), ColVecs_pullback
@@ -162,8 +240,9 @@ function ChainRulesCore.rrule(::Type{<:RowVecs}, X::AbstractMatrix)
     function RowVecs_pullback(::AbstractVector{<:AbstractVector{<:Real}})
         return error(
             "Pullback on AbstractVector{<:AbstractVector}.\n" *
-            "This might happen if you try to use gradients on the generic `kernelmatrix` or `kernelmatrix_diag`.\n" *
-            "To solve this issue overload `kernelmatrix(_diag)` for your kernel for `RowVecs`",
+            "This might happen if you try to use gradients on the generic `kernelmatrix` or `kernelmatrix_diag`,\n" *
+            "or because some external computation has acted on `RowVecs` to produce a vector of vectors." *
+            "If it is the former, to solve this issue overload `kernelmatrix(_diag)` for your kernel for `RowVecs`",
         )
     end
     return RowVecs(X), RowVecs_pullback