Is the mutating code the problem here? How could I debug this? about forwarddiff.jl HOT 2 CLOSED

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)

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