Merge pull request #36 from DanielVandH/triangle

Triangle caching

Project.toml

@@ -1,7 +1,7 @@
 name = "NaturalNeighbours"
 uuid = "f16ad982-4edb-46b1-8125-78e5a8b5a9e6"
 authors = ["Daniel VandenHeuvel <danj.vandenheuvel@gmail.com>"]
-version = "1.3.3"
+version = "1.3.4"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"

docs/src/compare.md

@@ -29,7 +29,7 @@ We also define the following constants and other useful variables:
 ```julia
 const itp_methods = (
     Sibson(0),
-    Triangle(),
+    Triangle(), 
     Nearest(),
     Laplace(),
     Sibson(1),

src/differentiation/differentiate.jl

@@ -82,6 +82,22 @@ end
 
 @inline default_diff_method(∂) = Direct() # isnothing(get_gradient(get_interpolant(∂))) ? Direct() : Iterative()
 
+function _get_nc_and_z(method::AbstractInterpolator{D}, p, z, gradients, hessians, tri, cache=NaturalNeighboursCache(tri); rng=Random.default_rng(), project=true) where {D}
+    if (method isa Triangle) || method == Nearest() # coordinates need to be the natural neighbours
+        nc = compute_natural_coordinates(Sibson(), tri, p, cache; rng, project)
+    else
+        nc = compute_natural_coordinates(method, tri, p, cache; rng, project)
+    end
+    if D == 0
+        zᵢ = _eval_natural_coordinates(nc, z)
+    elseif D == 1
+        zᵢ = _eval_natural_coordinates(method, nc, z, gradients, tri)
+    else # D == 2 
+        zᵢ = _eval_natural_coordinates(method, nc, z, gradients, hessians, tri)
+    end
+    return nc, zᵢ
+end
+
 @inline function (∂::NaturalNeighboursDifferentiator)(x, y, zᵢ, nc, id::Integer=1; parallel=false, method=default_diff_method(∂), kwargs...)
     method = dwrap(method)
     F = (number_type ∘ get_triangulation ∘ get_interpolant)(∂)
@@ -106,6 +122,7 @@ function (∂::NaturalNeighboursDifferentiator)(vals::AbstractVector, x::Abstrac
     @assert length(x) == length(y) == length(vals) "x, y, and vals must have the same length."
     method = dwrap(method)
     interpolant_method = iwrap(interpolant_method)
+    interpolant_method isa Triangle && populate_cache!(interpolant_method, get_triangulation(get_interpolant(∂)))
     if !parallel
         for i in eachindex(x, y)
             vals[i] = ∂(x[i], y[i], 1; method, interpolant_method, kwargs...)
@@ -134,6 +151,7 @@ function (∂::NaturalNeighboursDifferentiator{I,O})(x::AbstractVector, y::Abstr
     end
     method = dwrap(method)
     interpolant_method = iwrap(interpolant_method)
+    interpolant_method isa Triangle && populate_cache!(interpolant_method, tri)
     ∂(vals, x, y; method, interpolant_method, parallel, kwargs...)
     return vals
 end

src/interpolation/coordinates/triangle.jl

@@ -1,42 +1,18 @@
-function _compute_triangle_coordinates(
-    tri::Triangulation{P,T,BN,W,I,E,Es,BC,BEM,GVM,GVR,BPL,C,BE},
-    interpolation_point,
-    cache::NaturalNeighboursCache{F}=NaturalNeighboursCache(tri);
-    project=true,
-    kwargs...
-) where {P,T,BN,W,I,E,Es,BC,BEM,GVM,GVR,BPL,C,BE,F}
-    coordinates = get_coordinates(cache)
-    envelope = get_envelope(cache)
-    last_triangle = get_last_triangle(cache)
-    V = jump_and_march(tri, interpolation_point; try_points=last_triangle[], kwargs...)
-    i, j, return_flag = check_for_extrapolation(tri, V, interpolation_point, last_triangle)
-    return_flag && return two_point_interpolate!(coordinates, envelope, tri, i, j, interpolation_point, project)
-    i, j, k = triangle_vertices(V)
-    resize!(coordinates, 3)
-    resize!(envelope, 3)
-    λ₁, λ₂, λ₃ = _compute_triangle_barycentric_coordinates(tri, interpolation_point, i, j, k)
-    coordinates[1] = λ₁
-    coordinates[2] = λ₂
-    coordinates[3] = λ₃
-    envelope[1] = i
-    envelope[2] = j
-    envelope[3] = k
-    return NaturalCoordinates(coordinates, envelope, interpolation_point, tri)
-end
-
-function _compute_triangle_barycentric_coordinates(tri, interpolation_point, i, j, k)
+function _compute_triangle_shape_coefficients(tri, i, j, k)
     p, q, r = get_point(tri, i, j, k)
-    x₁, y₁ = getxy(p)
-    x₂, y₂ = getxy(q)
-    x₃, y₃ = getxy(r)
-    x, y = getxy(interpolation_point)
-    Δ = (y₂ - y₃) * (x₁ - x₃) + (x₃ - x₂) * (y₁ - y₃)
-    λ₁ = ((y₂ - y₃) * (x - x₃) + (x₃ - x₂) * (y - y₃)) / Δ
-    λ₂ = ((y₃ - y₁) * (x - x₃) + (x₁ - x₃) * (y - y₃)) / Δ
-    λ₃ = one(λ₁) - λ₁ - λ₂
-    return λ₁, λ₂, λ₃
-end
-
-function compute_natural_coordinates(::Triangle, tri, interpolation_point, cache=NaturalNeighboursCache(tri); kwargs...)
-    return _compute_triangle_coordinates(tri, interpolation_point, cache; kwargs...)
+    px, py = getxy(p)
+    qx, qy = getxy(q)
+    rx, ry = getxy(r)
+    Δ = qx * ry - qy * rx - px * ry + rx * py + px * qy - qx * py
+    s₁ = (qy - ry) / Δ
+    s₂ = (ry - py) / Δ
+    s₃ = (py - qy) / Δ
+    s₄ = (rx - qx) / Δ
+    s₅ = (px - rx) / Δ
+    s₆ = (qx - px) / Δ
+    s₇ = (qx * ry - rx * qy) / Δ
+    s₈ = (rx * py - px * ry) / Δ
+    s₉ = (px * qy - qx * py) / Δ
+    shape_function_coefficients = (s₁, s₂, s₃, s₄, s₅, s₆, s₇, s₈, s₉)
+    return shape_function_coefficients
 end

src/interpolation/eval.jl

@@ -44,6 +44,36 @@ end
     return _eval_natural_coordinates(nc, z)
 end
 
+@inline function _eval_interp(method::Triangle, itp::NaturalNeighboursInterpolant, p, cache; project=true, kwargs...)
+    tri = get_triangulation(itp)
+    z = get_z(itp)
+    last_triangle = get_last_triangle(cache)
+    V = jump_and_march(tri, p; try_points=last_triangle[], kwargs...)
+    i, j, return_flag = check_for_extrapolation(tri, V, p, last_triangle)
+    if return_flag
+        F = number_type(tri)
+        if project
+            t = two_point_interpolate(tri, i, j, p)
+            return z[i] * (1 - t) + z[j] * t
+        else
+            return typemax(F)
+        end
+    else        
+        F = number_type(tri)
+        i, j, k = triangle_vertices(sort_triangle(V))
+        if method.allow_cache && !isempty(method.s)
+            s₁, s₂, s₃, s₄, s₅, s₆, s₇, s₈, s₉ = method.s[(i, j, k)]
+        else
+            s₁, s₂, s₃, s₄, s₅, s₆, s₇, s₈, s₉ = _compute_triangle_shape_coefficients(tri, i, j, k)
+        end
+        α = s₁ * z[i] + s₂ * z[j] + s₃ * z[k]
+        β = s₄ * z[i] + s₅ * z[j] + s₆ * z[k]
+        γ = s₇ * z[i] + s₈ * z[j] + s₉ * z[k]
+        x, y = getxy(p)
+        return F(α * x + β * y + γ)
+    end
+end
+
 @inline function _eval_interp(method::Union{<:Farin,Sibson{1}}, itp::NaturalNeighboursInterpolant, p, cache; kwargs...)
     gradients = get_gradient(itp)
     if isnothing(gradients)
@@ -73,20 +103,4 @@ end
         return _eval_natural_coordinates(nc, z)
     end
     return _eval_natural_coordinates(method, nc, z, gradients, hessians, tri)
-end
-
-function _get_nc_and_z(method::AbstractInterpolator{D}, p, z, gradients, hessians, tri, cache=NaturalNeighboursCache(tri); rng=Random.default_rng(), project=true) where {D}
-    if method == Triangle() || method == Nearest() # coordinates need to be the natural neighbours
-        nc = compute_natural_coordinates(Sibson(), tri, p, cache; rng, project)
-    else
-        nc = compute_natural_coordinates(method, tri, p, cache; rng, project)
-    end
-    if D == 0
-        zᵢ = _eval_natural_coordinates(nc, z)
-    elseif D == 1
-        zᵢ = _eval_natural_coordinates(method, nc, z, gradients, tri)
-    else # D == 2 
-        zᵢ = _eval_natural_coordinates(method, nc, z, gradients, hessians, tri)
-    end
-    return nc, zᵢ
 end

src/interpolation/extrapolation.jl

@@ -1,4 +1,4 @@
-function two_point_interpolate!(tri, i, j, c)
+function two_point_interpolate(tri, i, j, c)
     # Project c onto the line through (a, b): https://stackoverflow.com/a/15187473
     # The orthogonal projection is not guaranteed to be on the line segment (corresponding 
     # to t < 0 or t > 1), in which case the weights are no longer a convex combination.
@@ -14,7 +14,7 @@ end
 
 function two_point_interpolate!(coordinates::AbstractVector{F}, envelope, tri, i, j, r, project=true) where {F} #interpolate r using two points i, j
     if project
-        t = two_point_interpolate!(tri, i, j, r)
+        t = two_point_interpolate(tri, i, j, r)
         resize!(coordinates, 2)
         resize!(envelope, 2)
         coordinates[1] = one(t) - t

src/interpolation/interpolate.jl

@@ -73,17 +73,51 @@ Interpolate using Sibson's coordinates with `C(d)` continuity at the data sites.
     Sibson(d) = d ∈ (0, 1) ? new{d}() : throw(ArgumentError("The Sibson interpolant is only defined for d ∈ (0, 1)."))
     Sibson() = new{0}()
 end
-struct Triangle{D} <: AbstractInterpolator{D} end
+struct Triangle{D} <: AbstractInterpolator{D}
+    allow_cache::Bool
+    s::Dict{NTuple{3,Int},NTuple{9,Float64}}
+end
+Triangle{D}(; allow_cache=false) where {D} = Triangle{D}(allow_cache, Dict{NTuple{3,Int},Float64}())
+Base.empty!(method::Triangle) = empty!(method.s)
+
+function populate_cache!(method::Triangle, tri::Triangulation)
+    method.allow_cache || return method
+    if length(method.s) == DelaunayTriangulation.num_solid_triangles(tri)
+        return method
+    elseif !isempty(method.s) # user is using a new triangulation
+        empty!(method)
+    end
+    for T in each_solid_triangle(tri)
+        V = sort_triangle(T)
+        i, j, k = triangle_vertices(V)
+        method.s[(i, j, k)] = _compute_triangle_shape_coefficients(tri, i, j, k)
+    end
+    return method
+end
+
 struct Nearest{D} <: AbstractInterpolator{D} end
 struct Laplace{D} <: AbstractInterpolator{D} end
 struct Farin{D} <: AbstractInterpolator{D} end
 struct Hiyoshi{D} <: AbstractInterpolator{D} end
 @doc """
-    Triangle()
+    Triangle(; allow_cache = true)
 
 Interpolate using a piecewise linear interpolant over the underlying triangulation.
-""" Triangle() = Triangle{0}()
-Triangle(d) = Triangle()
+
+!!! note "Cached coordinates with `allow_cache=true`"
+
+    The `Triangle()` interpolator is special as it will cache the coordinates used 
+    for each triangle. In particular, when an interpolator is evaluated with the 
+    `Triangle()` method, the object returned from `Triangle()` will store all 
+    the coordinates. For this reason, if you want to reuse `Triangle()` for different 
+    evaluations of the interpolant, you should be sure to reuse the same instance rather 
+    than reinstantiating it every single time. If you do not want this behaviour, set 
+    `allow_cache = false`.
+
+    If you only ever call the scalar-argument version of the interpolant, no caching will 
+    be done even with `allow_cache = true`.
+""" Triangle(; allow_cache=false) = Triangle{0}(; allow_cache)
+Triangle(d; allow_cache=false) = Triangle(; allow_cache)
 @doc """
     Nearest()
 
@@ -149,12 +183,14 @@ function (itp::NaturalNeighboursInterpolant)(x, y, id::Integer=1; parallel=false
     p = (F(x), F(y))
     cache = get_neighbour_cache(itp, id)
     method = iwrap(method)
+    # method isa Triangle && populate_cache!(method, tri)
     return _eval_interp(method, itp, p, cache; kwargs...)
 end
 
 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)
+    method isa Triangle && populate_cache!(method, get_triangulation(itp))
     if !parallel
         for i in eachindex(x, y)
             vals[i] = itp(x[i], y[i], 1; method, kwargs...)
@@ -178,6 +214,7 @@ end
     F = number_type(tri)
     vals = zeros(F, n)
     method = iwrap(method)
+    method isa Triangle && populate_cache!(method, tri)
     itp(vals, x, y; method, parallel, kwargs...)
     return vals
 end

test/doc_examples/interpolation_math.jl

@@ -141,7 +141,7 @@ x = vec([x for x in xg, _ in yg])
 y = vec([y for _ in xg, y in yg])
 ze = f.(x, y) |> x -> reshape(x, length(xg), length(yg))
 itp = interpolate(x1, y1, z1, derivatives=true)
-vals = itp(x, y, method=Triangle()) |> x -> reshape(x, length(xg), length(yg))
+vals = itp(x, y, method=Triangle(; allow_cache = false)) |> x -> reshape(x, length(xg), length(yg))
 fig = Figure(fontsize=36, size=(1700, 600))
 ax = Axis3(fig[1, 1], xlabel=L"x", ylabel=L"y", zlabel=L"z", title=L"(a): $f^{\text{TRI}}$", titlealign=:left)
 surface!(ax, xg, yg, vals)

test/doc_examples/swiss.jl

@@ -139,7 +139,7 @@ fig = plot_results(sibson_vals, sibson_1_vals, laplace_vals, triangle_vals, near
 sibson_vals_p = interpolant(x, y; method=Sibson(), parallel=true, project=false)
 sibson_1_vals_p = interpolant(x, y; method=Sibson(1), parallel=true, project=false)
 laplace_vals_p = interpolant(x, y; method=Laplace(), parallel=true, project=false)
-triangle_vals_p = interpolant(x, y; method=Triangle(), parallel=true, project=false)
+triangle_vals_p = interpolant(x, y; method=Triangle(; allow_cache = false), parallel=true, project=false)
 nearest_vals_p = interpolant(x, y; method=Nearest(), parallel=true, project=false)
 farin_vals_p = interpolant(x, y; method=Farin(), parallel=true, project=false)
 hiyoshi_vals_p = interpolant(x, y; method=Hiyoshi(2), parallel=true, project=false)

test/interpolation/basic_tests.jl

@@ -87,7 +87,6 @@ end
     @test_throws ArgumentError("Gradients and Hessians must be provided for Hiyoshi-2 interpolation. Consider using e.g. interpolate(tri, z; derivatives = true).") itp(0.5, 0.5; method=Hiyoshi(2))
 end
 
-
 @testset "Test Float32" begin
     rng = StableRNG(123)
     xs = randn(rng, 100)
@@ -98,7 +97,7 @@ end
     itp1 = interpolate(tri1, Float32.(zs); derivatives=true)
     itp2 = interpolate(tri2, zs; derivatives=true)
     for itp in (itp1, itp2)
-        for method in (Sibson(1), Sibson(), Laplace(), Farin(1), Hiyoshi(2), Triangle(), Nearest())
+        for method in (Sibson(1), Sibson(), Laplace(), Farin(1), Hiyoshi(2), Triangle(), Triangle(; allow_cache = false), Nearest())
             @inferred itp(rand(), rand(); method=method)
             @inferred itp(rand(), rand(); method=method, project=false)
             @inferred itp(rand(Float32), rand(Float32); method=method)
@@ -167,4 +166,35 @@ end
         H = NaturalNeighbours.get_hessian(itp, i)
         @test all(iszero, H)
     end
+end
+
+@testset "Testing Triangle()'s cache" begin
+    tri = triangulate(rand(2, 50))
+    method = Triangle()
+    method2 = Triangle(; allow_cache=true)
+    s = Dict{NTuple{3,Int},NTuple{9,Float64}}()
+    for T in each_solid_triangle(tri)
+        V = DelaunayTriangulation.sort_triangle(T)
+        s[V] = NaturalNeighbours._compute_triangle_shape_coefficients(tri, V...)
+        @test sum(s[V]) ≈ 1.0
+    end
+    itp = interpolate(tri, rand(50))
+    itp(rand(50), rand(50); method)
+    @test isempty(method.s)
+    itp(1 / 2, 1 / 2; method=method2)
+    @test isempty(method2.s)
+    itp(rand(50), rand(50); method=method2)
+    @test !isempty(method2.s)
+    @test method2.s == s
+    tri2 = triangulate(rand(2, 68))
+    itp = interpolate(tri2, rand(68))
+    itp(rand(50), rand(50); method=method2)
+    s = Dict{NTuple{3,Int},NTuple{9,Float64}}()
+    for T in each_solid_triangle(tri2)
+        V = DelaunayTriangulation.sort_triangle(T)
+        s[V] = NaturalNeighbours._compute_triangle_shape_coefficients(tri2, V...)
+        @inferred NaturalNeighbours._compute_triangle_shape_coefficients(tri2, V...)
+        @test sum(s[V]) ≈ 1.0
+    end
+    @test method2.s == s # make sure that method2 knows to regenerate the cache
 end

test/interpolation/constrained.jl

@@ -28,7 +28,7 @@ end
     _y = vec([y for _ in xg, y in yg])
     exact = f.(_x, _y)
     sibson_vals = itp(_x, _y; method=Sibson())
-    triangle_vals = itp(_x, _y; method=Triangle())
+    triangle_vals = itp(_x, _y; method=Triangle(; allow_cache = false))
     laplace_vals = itp(_x, _y; method=Laplace())
     sibson_1_vals = itp(_x, _y; method=Sibson(1))
     nearest_vals = itp(_x, _y; method=Nearest())

test/interpolation/coordinates/natural_coordinates.jl

@@ -10,7 +10,7 @@ include(normpath(@__DIR__, "../..", "helper_functions", "point_generator.jl"))
 
 to_mat(H) = [H[1] H[3]; H[3] H[2]]
 @testset "Natural coordinates" begin
-    for method in (Sibson(), Triangle(), Nearest(), Laplace())
+    for method in (Sibson(), Nearest(), Laplace())
         method = NNI.iwrap(method)
         pts = [(0.0, 8.0), (0.0, 0.0), (14.0, 0.0), (14.0, 8.0), (4.0, 4.0), (10.0, 6.0), (6.0, 2.0), (12.0, 4.0), (0.0, 4.0)]
         tri = triangulate(pts, randomise=false, delete_ghosts=false)
@@ -96,12 +96,6 @@ end
         @test _circular_equality(nc.indices, [23, 12, 5, 24, 7, 11])
         @test _circular_equality(nc.coordinates, [K1I1J1, F1G1J1K1, AF1G1D1B1A1, A1B1C1, B1D1E1C1, H1I1E1D1G1] ./ AX, ≈, rtol = 1e-2)
 
-        # Triangle 
-        nc = NNI.compute_natural_coordinates(NNI.Triangle(), tri, q; rng=rng)
-        V = jump_and_march(tri, q; rng)
-        @test _circular_equality(nc.indices, [5, 11, 12])
-        @test _circular_equality(nc.coordinates, [0.52, 0.3, 0.18], ≈, rtol = 1e-2)
-
         # Laplace 
         nc = NNI.compute_natural_coordinates(NNI.Laplace(), tri, q; rng=rng)
         dqw = 4.0311288741493

test/interpolation/structs.jl

@@ -4,17 +4,17 @@ const NNI = NaturalNeighbours
 
 @testset "iwrap" begin
     @test NNI.iwrap(NNI.Sibson()) == NNI.Sibson()
-    @test NNI.iwrap(NNI.Triangle()) == NNI.Triangle()
+    @test NNI.iwrap(NNI.Triangle()) isa NNI.Triangle
     @test NNI.iwrap(NNI.Nearest()) == NNI.Nearest()
     @test NNI.iwrap(NNI.Laplace()) == NNI.Laplace()
     @test NNI.iwrap(:sibson) == NNI.Sibson()
-    @test NNI.iwrap(:triangle) == NNI.Triangle()
+    @test NNI.iwrap(:triangle) isa NNI.Triangle
     @test NNI.iwrap(:nearest) == NNI.Nearest()
     @test NNI.iwrap(:laplace) == NNI.Laplace()
     @test_throws ArgumentError NNI.iwrap(:lap)
 
     @test NNI.iwrap(NNI.Sibson(1)) == NNI.Sibson(1)
-    @test NNI.iwrap(NNI.Triangle(1)) == NNI.Triangle(0)
+    @test NNI.iwrap(NNI.Triangle(1)) isa NNI.Triangle{0}
     @test_throws ArgumentError NNI.Sibson(5)
     @test NNI.iwrap(NNI.Laplace(1)) == NNI.Laplace(0)
     @test NNI.iwrap(NNI.Hiyoshi(2)) == NNI.Hiyoshi(2)