Is the mutating code the problem here? How could I debug this?

[FelixBenning]

using KernelFunctions
using ForwardDiff: derivative

struct Pt{Dim}
	pos::AbstractArray
	partial
end

Pt(x;partial=()) = Pt{length(x)}(x, partial)

struct TaylorKernel <: KernelFunctions.Kernel
	k::KernelFunctions.Kernel
end

function (tk::TaylorKernel)(x::Pt{Dim}, y::Pt{Dim}) where Dim
	k = tk.k
	for ii in x.partial
		k = (x₁,x₂) -> derivative(0) do Δx
			x₁[ii] += Δx
			return k(x₁,x₂)
		end
	end
	for jj in y.partial
		k = (x₁,x₂) -> derivative(0) do Δx
			x₂[jj] += Δx
			return k(x₁,x₂)
		end
	end
	k(x.pos, y.pos)
end

k = TaylorKernel(MaternKernel())

k(Pt([1]), Pt([2])) # k(x,y)  with x=1, y=2
k(Pt([1], partial=(1,)), Pt([2])) # ∂ₓk(x,y) with x=1, y=2

The first call works, the second call results in

ERROR: InexactError: Int(Int64, Dual{ForwardDiff.Tag{var"#15#19"{Vector{Int64}, Vector{Int64}, Int64}, Int64}}(1,1))
Stacktrace:
 [1] Int64
   @ ~/.julia/packages/ForwardDiff/vXysl/src/dual.jl:364 [inlined]
 [2] convert
   @ ./number.jl:7 [inlined]
 [3] setindex!
   @ ./array.jl:966 [inlined]
 [4] (::var"#15#19"{Vector{Int64}, Vector{Int64}, Int64})(Δx::ForwardDiff.Dual{ForwardDiff.Tag{var"#15#19"{Vector{Int64}, Vector{Int64}, Int64}, Int64}, Int64, 1})
   @ Main ~/code/incrementalGP/scripts/experiments.jl:19
 [5] derivative
   @ ~/.julia/packages/ForwardDiff/vXysl/src/derivative.jl:14 [inlined]
 [6] (::var"#14#18"{Int64})(x₁::Vector{Int64}, x₂::Vector{Int64})
   @ Main ~/code/incrementalGP/scripts/experiments.jl:18
 [7] (::TaylorKernel)(x::Pt{1}, y::Pt{1})
   @ Main ~/code/incrementalGP/scripts/experiments.jl:29
 [8] top-level scope
   @ ~/code/incrementalGP/scripts/experiments.jl:35

[FelixBenning]

The following code (with less mutation) does not yield any output and (after spamming ^C) finally prints segfaults:

using KernelFunctions
using OneHotArrays: OneHotVector
using ForwardDiff: derivative

struct Pt{Dim}
	pos::AbstractArray
	partial
end

Pt(x;partial=()) = Pt{length(x)}(x, partial)

struct TaylorKernel <: KernelFunctions.Kernel
	k::KernelFunctions.Kernel
end

function (tk::TaylorKernel)(x::Pt{Dim}, y::Pt{Dim}) where Dim
	k = tk.k
	for ii in x.partial
		k = (x₁,x₂) -> derivative(0) do Δx
			return k(x₁ + Δx * OneHotVector(ii, Dim), x₂)
		end
	end
	for jj in y.partial
		k = (x₁,x₂) -> derivative(0) do Δx
			return k(x₁,x₂ + Δx * OneHotVector(jj, Dim))
		end
	end
	k(x.pos, y.pos)
end

k = TaylorKernel(MaternKernel())

k(Pt([1]), Pt([2])) # k(x,y)  with x=1, y=2
k(Pt([1], partial=(1,)), Pt([2])) # ∂ₓk(x,y)

[FelixBenning]

The second code snippet does not work, because unlike for variables you can not use the variable in the reassignment like this. This ends up in an infinite recursion which explains the segmentation fault I think.

The following seems to work

function (tk::TaylorKernel)(x::Pt{Dim}, y::Pt{Dim}) where Dim
	if !isnothing(local next = iterate(x.partial))
		ii, state = next # take partial derivative in direction ii
		return FD.derivative(0) do dx
			tk( # recursion
				Pt(
					x.pos + dx * OneHotVector(ii, Dim), # directional variation
					partial=Base.rest(x.partial, state) # remaining partial derivatives
				),
				y
			)
		end
	end
	if !isnothing(local next = iterate(y.partial))
		jj, state = next # take partial derivative in direction jj
		return FD.derivative(0) do dy
			tk( # recursion
				x,
				Pt(
					y.pos + dy * OneHotVector(jj, Dim), # directional variation
					partial=Base.rest(x.partial, state) # remaining partial derivatives
				)
			)
		end
	end
	tk.k(x.pos, y.pos)
end

this leads me to believe, that the mutation in the first snippet was actually the culprit. I'll close the issue