Interpolating vector quantities

[DanielVandH]

I wonder if it would be nice to have a way to interpolate data of the form $\mathbf z = f(x, y)$. One application would be solutions from a PDE, in which case we might want to interpolate at the same point but at many times. Being able to do this just a single time would save a lot of time.

[asinghvi17]

I guess this is kind of like barycentric interpolation? The way I did it in GeometryOps was to have two APIs, one that returns weights on each coordinate and one that returns a value. It causes a bit of code duplication, though.

You could potentially also return a namedtuple of (; coord1 = weight1, ...), or a vector of pairs.

[DanielVandH]

All the interpolants take the form (edit 13/05/2024: not true, the smooth interpolants are $z$-dependent obviously... woops)

$$f(\boldsymbol x_0) = \sum_{i \in N_0} \lambda_iz_i,$$

where $N_0$ are the natural neighbours of $\boldsymbol x_0$ and the natural coordinates $\boldsymbol\lambda$ depend only on $\boldsymbol x_0$ and its natural neighbours. In particular, the weights are independent of the data. So for this package, implementing this is really just a matter of replacing $z_i$ with the vector $\boldsymbol z_i$ of data points at $\boldsymbol x_i$. I think it'll literally be the same code otherwise, but some code duplication would be needed to separate the scalar and vector cases.

[asinghvi17]

Gotcha. I guess this could be built off of compute_natural_coordinates? (which should probably be documented as well, unless it's meant to be non-public API)

[DanielVandH]

I haven't really intended to make much public API outside of the interpolants themselves. The natural coordinates stuff isn't really public, but it could be. Everything that is public is on the first page of the docs, but if there are functions in the package you find useful feel free to make a PR putting them in there so that they are a part of the API.

There'd be a couple steps to implement this issue:

1. Change z::Vector{F} here

   NaturalNeighbours.jl/src/data_structures/interpolant.jl

   Lines 1 to 8 in 95a1f50

   	struct NaturalNeighboursInterpolant{T<:Triangulation,F,G,H,N,D}
   	triangulation::T
   	z::Vector{F}
   	gradient::G # (∂ˣf, ∂ʸf)
   	hessian::H # (∂ˣˣf, ∂ʸʸf, ∂ˣʸf)
   	neighbour_cache::N
   	derivative_cache::D
   	function NaturalNeighboursInterpolant(

   to z::Z, specifying Z <: Union{Vector{F}, Matrix{F}}.
2. Define is_scalar(itp::NaturalNeighborusInterpolant) = get_z(itp) isa Vector
3. I think then we can change

   NaturalNeighbours.jl/src/interpolation/eval.jl

   Lines 1 to 8 in 95a1f50

   	@inline function _eval_natural_coordinates(coordinates, indices, z)
   	val = zero(eltype(z))
   	for (λ, k) in zip(coordinates, indices)
   	zₖ = z[k]
   	val += λ * zₖ
   	end
   	return val
   	end

   into two functions:

@inline function _eval_natural_coordinates(coordinates, indices, z::Vector)
    val = zero(eltype(z))
    for (λ, k) in zip(coordinates, indices)
        zₖ = z[k]
        val += λ * zₖ
    end
    return val
end
@inline function _eval_natural_coordinates(coordinates, indices, z::Matrix)
    val = zeros(eltype(z), size(z, 1))
    for (λ, k) in zip(coordinates, indices)
        zₖ = @views z[:, k]
        @. val += λ * zₖ
    end
    return val
end

4. Inside https://github.com/DanielVandH/NaturalNeighbours.jl/blob/95a1f50d01be75026375304f45cd1b480c4c5571/src/interpolation/interpolate.jl, I think I'd want to define

function _val_cache(itp::NaturalNeighboursInterpolant, n)
	tri = get_triangulation(itp)
	F = number_type(tri)
	if is_scalar(itp)
		return zeros(F, n)
	else
		z = get_z(itp)
	    m = size(z, 1)
		return zeros(F, m, n)
   end 
end

so that vals = _val_cache(itp, n) in

NaturalNeighbours.jl/src/interpolation/interpolate.jl

Line 162 in 95a1f50

vals = zeros(F, n)

. Then inside

NaturalNeighbours.jl/src/interpolation/interpolate.jl

Lines 138 to 156 in 95a1f50

	function (itp::NaturalNeighboursInterpolant)(vals::AbstractVector, x::AbstractVector, y::AbstractVector; parallel=true, method=Sibson(), kwargs...)
	@assert length(x) == length(y) == length(vals) "x, y, and vals must have the same length."
	method = iwrap(method)
	if !parallel
	for i in eachindex(x, y)
	vals[i] = itp(x[i], y[i], 1; method, kwargs...)
	end
	else
	caches = get_neighbour_cache(itp)
	nt = length(caches)
	chunked_iterator = chunks(vals, nt)
	Threads.@threads for (xrange, chunk_id) in chunked_iterator
	for i in xrange
	vals[i] = itp(x[i], y[i], chunk_id; method, kwargs...)
	end
	end
	end
	return nothing
	end

we'd probably have to replace vals[i] = itp(x[i], y[i], chunk_id; method, kwargs...) (and similarly for Line 143) with something like

_set_vals!(vals, itp::NaturalNeighboursInterpolant, idx, val) = is_scalar(itp) ? (vals[idx] = val) : (vals[:, idx] .= val)
val = itp(x[i], y[i], chunk_id; method, kwargs...)
_set_vals!(vals, itp, i, val)

or

if is_scalar(itp)
	vals[i] = itp(x[i], y[i], chunk_id; method, kwargs...)
else
	vals[:, i] .= itp(x[i], y[i], chunk_id; method, kwargs...)
end

The vector case still wouldn't be so optimised at this point since itp(x[i], y[i], chunk_id; method, kwargs...) would be allocating an entire vector of length size(z, 1) each time. You'd probably need to add an extra field to NaturalNeighboursCache here

NaturalNeighbours.jl/src/data_structures/neighbour_cache.jl

Lines 1 to 9 in 95a1f50

	struct NaturalNeighboursCache{F,I,H,E,R}
	coordinates::Vector{F}
	envelope::Vector{I}
	insertion_event_history::H
	poly_points::Vector{NTuple{2,F}}
	temp_adjacent::Adjacent{I,E}
	last_triangle::R
	end
	get_coordinates(cache::NaturalNeighboursCache) = cache.coordinates

e.g. data_cache::C which is nothing if is_scalar(itp) and a vector of length size(z, 1) otherwise. Then, an in-place version of _eval_natural_coordinates(coordinates, indices, z::Vector) (above) would need to be defined to make use of this cache appropriately to avoid these extra allocations.

I think that would be enough. I haven't run any of the above code so I've not seen if I've overlooked anything, but that'd be a good start. If we wanted to go a step further and allow for differentiating such data, I think a similar amount of work might have to be done - I think the linear system is just of the form $\boldsymbol A\boldsymbol X = \boldsymbol B$ (matrix RHS) instead of a vector RHS.

[DanielVandH]

I haven't run any of the above code so I've not seen if I've overlooked anything, but that'd be a good start. If we wanted to go a step further and allow for differentiating such data, I think a similar amount of work might have to be done - I think the linear system is just of the form (matrix RHS) instead of a vector RHS.

Differentiating might be one of the first things needed actually, at least the generation of the derivatives, so that methods like Sibson(1) could actually be used also. In that case lines like

NaturalNeighbours.jl/src/differentiation/generate.jl

Lines 76 to 77 in 0d79b4d

	∇ = Vector{NTuple{2,F}}(undef, n)
	ℋ = Vector{NTuple{3,F}}(undef, n)

would need to be changed, in the matrix case, to

∇ = Matrix{NTuple{2,F}}(undef, size(z)...)
ℋ = Matrix{NTuple{3,F}}(undef, size(z)...)

[DanielVandH]

Started an implementation in #32, but just need to get around to finishing it off testing and documenting it. Not sure when I'll do that.