Merge pull request #35 from DanielVandH/instab

Fix one instability

Project.toml

@@ -1,21 +1,19 @@
 name = "NaturalNeighbours"
 uuid = "f16ad982-4edb-46b1-8125-78e5a8b5a9e6"
 authors = ["Daniel VandenHeuvel <[email protected]>"]
-version = "1.3.2"
+version = "1.3.3"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
 DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
 ChunkSplitters = "1.0, 2.0"
 DelaunayTriangulation = "1.0"
 ElasticArrays = "1.2"
-PrecompileTools = "1.2"
 julia = "1.9"
 
 [extras]

README.md

@@ -42,12 +42,12 @@ hiyoshi_vals = itp(_x, _y; method=Hiyoshi(2))
 ## Plot 
 function plot_2d(fig, i, j, title, vals, xg, yg, x, y, show_scatter=true)
     ax = Axis(fig[i, j], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
-    contourf!(ax, xg, yg, reshape(vals, (length(xg), length(yg))), color=vals, colormap=:viridis, levels=-1:0.05:0, extendlow=:auto, extendhigh=:auto)
+    contourf!(ax, xg, yg, reshape(vals, (length(xg), length(yg))), colormap=:viridis, levels=-1:0.05:0, extendlow=:auto, extendhigh=:auto)
     show_scatter && scatter!(ax, x, y, color=:red, markersize=14)
 end
 function plot_3d(fig, i, j, title, vals, xg, yg)
     ax = Axis3(fig[i, j], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
-    surface!(ax, xg, yg, reshape(vals, (length(xg), length(yg))), color=vals, colormap=:viridis, levels=-1:0.05:0, extendlow=:auto, extendhigh=:auto)
+    surface!(ax, xg, yg, reshape(vals, (length(xg), length(yg))), color=vals, colormap=:viridis)
 end
 
 all_vals = (sibson_vals, triangle_vals, laplace_vals, sibson_1_vals, nearest_vals, farin_vals, hiyoshi_vals, exact)

docs/src/differentiation.md

@@ -33,7 +33,7 @@ Here is the surface of $f$ along with its derivatives.
 using CairoMakie
 function plot_f(fig, x, y, vals, title, i, show_3d=true, zlabel="z")
     ax = Axis(fig[1, i], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
-    c = contourf!(ax, x, y, vals, color=vals, colormap=:viridis,  extendhigh=:auto)
+    c = contourf!(ax, x, y, vals, colormap=:viridis,  extendhigh=:auto)
     if show_3d
         ax = Axis3(fig[2, i], xlabel="x", ylabel="y", zlabel=zlabel, width=600, height=600, title=" ", titlealign=:left, azimuth=0.49)
         surface!(ax, x, y, vals, color=vals, colormap=:viridis)
@@ -85,7 +85,7 @@ vorn = voronoi(tri)
 fig = Figure(fontsize=36, size=(1800, 600))
 ax = Axis(fig[1, 1], xlabel="x", ylabel="y", width=600, height=600, title="(a): Data and triangulation", titlealign=:left)
 scatter!(ax, x, y, color=:black, markersize=9)
-triplot!(ax, tri, color=:black, linewidth=2, show_convex_hull=false)
+triplot!(ax, tri, strokecolor=:black, strokewidth=2, show_convex_hull=false)
 voronoiplot!(ax, vorn, strokecolor=:blue)
 xlims!(ax, 0, 1)
 ylims!(ax, 0, 1)
@@ -123,7 +123,7 @@ using LinearAlgebra
 function plot_f2(fig, x, y, vals, title, i, tri, levels, show_3d=true, zlabel="z")
     triangles = [T[j] for T in each_solid_triangle(tri), j in 1:3]
     ax = Axis(fig[1, i], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
-    c = tricontourf!(ax, x, y, vals, color=vals, triangulation=triangles', colormap=:viridis, extendhigh=:auto, levels=levels)
+    c = tricontourf!(ax, x, y, vals, triangulation=triangles', colormap=:viridis, extendhigh=:auto, levels=levels)
     if show_3d
         ax = Axis3(fig[2, i], xlabel="x", ylabel="y", zlabel=zlabel, width=600, height=600, title=" ", titlealign=:left, azimuth=0.49)
         mesh!(ax, hcat(x, y, vals), triangles, color=vals, colormap=:viridis, colorrange=extrema(levels))
@@ -306,7 +306,7 @@ function rrmserr(z, ẑ)
 end
 function plot_f2(fig, x, y, vals, title, i, levels, show_3d=true, zlabel="z")
     ax = Axis(fig[1, i], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
-    c = contourf!(ax, x, y, vals, color=vals, colormap=:viridis, extendhigh=:auto, levels=levels)
+    c = contourf!(ax, x, y, vals, colormap=:viridis, extendhigh=:auto, levels=levels)
     if show_3d
         ax = Axis3(fig[2, i], xlabel="x", ylabel="y", zlabel=zlabel, width=600, height=600, title=" ", titlealign=:left, azimuth=0.49)
         surface!(ax, x, y, vals, color=vals, colormap=:viridis, colorrange=extrema(levels))

docs/src/interpolation.md

@@ -151,7 +151,7 @@ Now we can plot.
 ```julia
 function plot_itp(fig, x, y, vals, title, i, show_data_sites, itp, xd=nothing, yd=nothing, show_3d=true, levels=-0.1:0.05:0.3)
     ax = Axis(fig[1, i], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
-    c = contourf!(ax, x, y, vals, color=vals, colormap=:viridis, levels=levels, extendhigh=:auto)
+    c = contourf!(ax, x, y, vals, colormap=:viridis, levels=levels, extendhigh=:auto)
     show_data_sites && scatter!(ax, xd, yd, color=:red, markersize=9)
     tri = itp.triangulation
     ch_idx = get_convex_hull_vertices(tri)

docs/src/swiss.md

@@ -108,8 +108,6 @@ function update_boundary(boundary_points, data_sites)
     return boundary_points, data_sites, boundary_nodes
 end
 boundary_points, data_sites, boundary_nodes = update_boundary(boundary_points, data_sites)
-tri = triangulate(data_sites, boundary_nodes=boundary_nodes)
-triangles = [T[j] for T in each_solid_triangle(tri), j in 1:3]
 ```
 
 Next, before, we plot, let us downsample the data. We do this since the data set is quite large, so when we interpolate it'll be useful to have fewer data points for the purpose of this example.
@@ -134,15 +132,10 @@ Now let's look at the data.
 colorrange = (0, 4)
 levels = LinRange(colorrange..., 40)
 fig = Figure(fontsize=24)
-ax1 = Axis3(fig[1, 1], xlabel="Longitude", ylabel="Latitude", zlabel="Elevation (km)", width=600, height=400, azimuth=0.9, title="(a): Original height data (n = $(length(elevation_data)))", titlealign=:left)
-mesh!(ax1, data, triangles, color=elevation_data, colorrange=colorrange)
-ax2 = Axis(fig[1, 2], xlabel="Longitude", ylabel="Latitude", width=600, height=400, title="(b): Original height data (n = $(length(elevation_data)))", titlealign=:left)
-tf = tricontourf!(ax2, tri, elevation_data, color=elevation_data, levels=levels)
-
-ax3 = Axis3(fig[2, 1], xlabel="Longitude", ylabel="Latitude", zlabel="Elevation (km)", width=600, height=400, azimuth=0.9, title="(c): Downsampled height data (n = $(length(ds_elevation_data)))", titlealign=:left)
-mesh!(ax3, ds_data, ds_triangles, color=ds_elevation_data, colorrange=colorrange)
-ax4 = Axis(fig[2, 2], xlabel="Longitude", ylabel="Latitude", width=600, height=400, title="(d): Downsampled height data (n = $(length(ds_elevation_data)))", titlealign=:left)
-tricontourf!(ax4, ds_tri, ds_elevation_data, color=ds_elevation_data, levels=levels)
+ax1 = Axis3(fig[1, 1], xlabel="Longitude", ylabel="Latitude", zlabel="Elevation (km)", width=600, height=400, azimuth=0.9, title="(c): Downsampled height data (n = $(length(ds_elevation_data)))", titlealign=:left)
+mesh!(ax1, ds_data, ds_triangles, color=ds_elevation_data, colorrange=colorrange)
+ax2 = Axis(fig[1, 2], xlabel="Longitude", ylabel="Latitude", width=600, height=400, title="(d): Downsampled height data (n = $(length(ds_elevation_data)))", titlealign=:left)
+tf = tricontourf!(ax2, ds_tri, ds_elevation_data, levels=levels)
 Colorbar(fig[1:2, 3], tf)
 resize_to_layout!(fig)
 ```
@@ -184,7 +177,7 @@ hiyoshi_vals = interpolant(x, y; method=Hiyoshi(2), parallel=true)
 Let's look at our results for each of these methods.
 
 ```julia
-query_tri = triangulate([x'; y'])
+query_tri = triangulate([x'; y']; randomise=false)
 query_triangles = [T[j] for T in each_solid_triangle(query_tri), j in 1:3]
 function plot_results!(fig, i1, j1, i2, j2, x, y, xg, yg, vals, title1, title2, query_triangles, query_tri, a, b, c, d, e, f, nx, ny)
     ax = Axis3(fig[i1, j1], xlabel="Longitude", ylabel="Latitude", zlabel="Elevation (km)", width=600, height=400, azimuth=0.9, title=title1, titlealign=:left)
@@ -193,7 +186,7 @@ function plot_results!(fig, i1, j1, i2, j2, x, y, xg, yg, vals, title1, title2,
     ylims!(ax, c, d)
     zlims!(ax, e, f)
     ax = Axis(fig[i2, j2], xlabel="Longitude", ylabel="Latitude", width=600, height=400, title=title2, titlealign=:left)
-    contourf!(ax, xg, yg, reshape(vals, (nx, ny)), color=vals, levels=levels)
+    contourf!(ax, xg, yg, reshape(vals, (nx, ny)), levels=levels)
     lines!(ax, [get_point(query_tri, i) for i in get_convex_hull_vertices(query_tri)], color=:red, linewidth=4, linestyle=:dash)
     lines!(ax, ds_boundary_points, color=:white, linewidth=4)
     xlims!(ax, a, b)
@@ -215,15 +208,15 @@ function plot_results(sibson_vals, sibson_1_vals, laplace_vals, triangle_vals, n
     ylims!(ax, c, d)
     zlims!(ax, e, f)
     ax = Axis(fig[4, 4], xlabel="Longitude", ylabel="Latitude", width=600, height=400, title="(o): Original height data", titlealign=:left)
-    tricontourf!(ax, tri, elevation_data, color=elevation_data, levels=levels)
+    tricontourf!(ax, tri, elevation_data, levels=levels)
     xlims!(ax, a, b)
     ylims!(ax, c, d)
     Colorbar(fig[1:4, 5], m1)
     resize_to_layout!(fig)
     return fig
 end
 e, f = 0.0, 4.5
-fig = plot_results(sibson_vals, sibson_1_vals, laplace_vals, triangle_vals, nearest_vals, farin_vals, hiyoshi_vals, query_triangles, interpolant, a, b, c, d, e, f, nx, ny, data, triangles, elevation_data, tri)
+fig = plot_results(sibson_vals, sibson_1_vals, laplace_vals, triangle_vals, nearest_vals, farin_vals, hiyoshi_vals, query_triangles, interpolant, a, b, c, d, e, f, nx, ny, ds_data, ds_triangles, ds_elevation_data, ds_tri)
 ```
 
 ```@raw html
@@ -246,7 +239,7 @@ triangle_vals_p = interpolant(x, y; method=Triangle(), parallel=true, project=fa
 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)
-fig = plot_results(sibson_vals_p, sibson_1_vals_p, laplace_vals_p, triangle_vals_p, nearest_vals_p, farin_vals_p, hiyoshi_vals_p, query_triangles, interpolant, a, b, c, d, e, f, nx, ny, data, triangles, elevation_data, tri)
+fig = plot_results(sibson_vals_p, sibson_1_vals_p, laplace_vals_p, triangle_vals_p, nearest_vals_p, farin_vals_p, hiyoshi_vals_p, query_triangles, interpolant, a, b, c, d, e, f, nx, ny, ds_data, ds_triangles, ds_elevation_data, ds_tri)
 ```
 
 ```@raw html
@@ -270,7 +263,7 @@ triangle_vals_p[exterior_idx] .= Inf
 nearest_vals_p[exterior_idx] .= Inf
 farin_vals_p[exterior_idx] .= Inf
 hiyoshi_vals_p[exterior_idx] .= Inf
-fig = plot_results(sibson_vals_p, sibson_1_vals_p, laplace_vals_p, triangle_vals_p, nearest_vals_p, farin_vals_p, hiyoshi_vals_p, query_triangles, interpolant, a, b, c, d, e, f, nx, ny, data, triangles, elevation_data, tri)
+fig = plot_results(sibson_vals_p, sibson_1_vals_p, laplace_vals_p, triangle_vals_p, nearest_vals_p, farin_vals_p, hiyoshi_vals_p, query_triangles, interpolant, a, b, c, d, e, f, nx, ny, ds_data, ds_triangles, ds_elevation_data, ds_tri)
 ```
 
 ```@raw html

src/NaturalNeighbours.jl

@@ -58,7 +58,7 @@ import DelaunayTriangulation: DelaunayTriangulation,
     get_convex_hull_vertices,
     get_boundary_nodes,
     triangle_vertices
-    
+
 import ChunkSplitters: chunks
 using ElasticArrays
 using LinearAlgebra
@@ -98,31 +98,4 @@ include("differentiation/utils.jl")
 
 include("utils.jl")
 
-#=
-using PrecompileTools
-@setup_workload begin
-    @compile_workload begin
-        points = [(rand(), rand()) for _ in 1:10]
-        tri = triangulate(points)
-        f = (x, y) -> sin(x * y) - cos(x - y) * exp(-(x - y)^2)
-        z = [f(x, y) for (x, y) in points]
-        itp = interpolate(tri, z; derivatives=true, parallel=false)
-        for method in (Sibson(1), Sibson(), Laplace(), Farin(1), Hiyoshi(2), Triangle(), Nearest())
-            itp(rand(), rand(); method=method)
-            itp(rand(), rand(); method=method, project=false)
-            itp([rand()], [rand()]; method=method)
-            itp([rand()], [rand()]; method=method, project=false)
-        end
-        ∂1 = differentiate(itp, 1)
-        ∂2 = differentiate(itp, 2)
-        for method in (Iterative(), Direct())
-                for ∂ in (∂1, ∂2)
-                    ∂(rand(), rand(), method=method)
-                    ∂([rand()], [rand()], method=method)
-                end
-        end
-    end
-end
-=#
-
 end # module NaturalNeighbours

src/interpolation/interpolate.jl

@@ -35,6 +35,11 @@ all the points. The `kwargs...` argument is passed into `add_point!` from Delaun
 the `project` argument determines whether extrapolation is performed by projecting any exterior points onto the boundary of the convex hull 
 of the data sites and performing two-point interpolation, or to simply replace any extrapolated values with `Inf`.
 
+!!! performance 
+
+    For the best performance when evaluating the interpolant at many points, either of the second or 
+    third methods are preferred over repeatedly calling the first.
+
 !!! warning
 
     Until we implement ghost point extrapolation, behaviour near the convex hull of your data sites may in some cases be undesirable,
@@ -65,7 +70,8 @@ abstract type AbstractInterpolator{D} end
 
 Interpolate using Sibson's coordinates with `C(d)` continuity at the data sites.
 """ struct Sibson{D} <: AbstractInterpolator{D}
-    Sibson(d=0) = d ∈ (0, 1) ? new{d}() : throw(ArgumentError("The Sibson interpolant is only defined for d ∈ (0, 1)."))
+    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 Nearest{D} <: AbstractInterpolator{D} end
@@ -76,33 +82,38 @@ struct Hiyoshi{D} <: AbstractInterpolator{D} end
     Triangle()
 
 Interpolate using a piecewise linear interpolant over the underlying triangulation.
-""" Triangle(d=0) = Triangle{0}()
+""" Triangle() = Triangle{0}()
+Triangle(d) = Triangle()
 @doc """
     Nearest()
 
 Interpolate by taking the function value at the nearest data site.
-""" Nearest(d=0) = Nearest{0}()
+""" Nearest() = Nearest{0}()
+Nearest(d) = Nearest{d}()
 @doc """
     Laplace()
 
 Interpolate using Laplace's coordinates.
-""" Laplace(d=0) = Laplace{0}()
+""" Laplace() = Laplace{0}()
+Laplace(d) = Laplace()
 @doc """
     Farin()
 
 Interpolate using Farin's C(1) interpolant.
-""" Farin(d=0) = Farin{1}()
+""" Farin() = Farin{1}()
+Farin(d) = Farin()
 @doc """
-    Hiyoshi(d=2)
+    Hiyoshi(d)
 
-Interpolate using Hiyoshi's C(2) interpolant. Hiyoshi's interpolant C(0) is not yet implemented, 
+Interpolate using Hiyoshi's C(d) interpolant. Hiyoshi's interpolant C(0) is not yet implemented, 
 but we do not make any conversions to C(2) like in e.g. `Farin()`, e.g. `Farin()` gets 
 converted to `Farin(1)` but, to support possible later versions, `Hiyoshi()` does not get 
 converted to `Hiyoshi(2)`.
-""" Hiyoshi(d=0) = Hiyoshi{d}()
+""" Hiyoshi() = Hiyoshi{0}()
+Hiyoshi(d) = Hiyoshi{d}()
 
 @inline iwrap(s::AbstractInterpolator) = s
-@inline function iwrap(s::Symbol)
+@inline function iwrap(s::Symbol) # this is bad design, should just prohibit symbols
     if s == :sibson
         return Sibson()
     elseif s == :sibson_1

test/doc_examples/differentiation.jl

@@ -13,7 +13,7 @@ f′′ = (x, y) -> [(162*exp(-(9 * x - 4)^2 - (9 * y - 7)^2))/5-(243*exp(-(9 *
 
 function plot_f(fig, x, y, vals, title, i, show_3d=true, zlabel="z")
     ax = Axis(fig[1, i], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
-    c = contourf!(ax, x, y, vals, color=vals, colormap=:viridis, extendhigh=:auto)
+    c = contourf!(ax, x, y, vals, colormap=:viridis, extendhigh=:auto)
     if show_3d
         ax = Axis3(fig[2, i], xlabel="x", ylabel="y", zlabel=zlabel, width=600, height=600, title=" ", titlealign=:left, azimuth=0.49)
         surface!(ax, x, y, vals, color=vals, colormap=:viridis)
@@ -54,7 +54,7 @@ vorn = voronoi(tri)
 fig = Figure(fontsize=50, size=(1800, 600))
 ax = Axis(fig[1, 1], xlabel="x", ylabel="y", width=600, height=600, title="(a): Data and triangulation", titlealign=:left)
 scatter!(ax, x, y, color=:black, markersize=9)
-triplot!(ax, tri, color=:black, linewidth=2)
+triplot!(ax, tri, strokecolor=:black, strokewidth=2)
 voronoiplot!(ax, vorn, strokecolor=:blue, color=(:white, 0.0))
 xlims!(ax, 0, 1)
 ylims!(ax, 0, 1)
@@ -74,7 +74,7 @@ fig
 function plot_f2(fig, x, y, vals, title, i, tri, levels, show_3d=true, zlabel="z")
     triangles = [T[j] for T in each_solid_triangle(tri), j in 1:3]
     ax = Axis(fig[1, i], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
-    c = tricontourf!(ax, x, y, vals, color=vals, triangulation=triangles', colormap=:viridis, extendhigh=:auto, levels=levels)
+    c = tricontourf!(ax, x, y, vals, triangulation=triangles', colormap=:viridis, extendhigh=:auto, levels=levels)
     if show_3d
         ax = Axis3(fig[2, i], xlabel="x", ylabel="y", zlabel=zlabel, width=600, height=600, title=" ", titlealign=:left, azimuth=0.49)
         mesh!(ax, hcat(x, y, vals), triangles, color=vals, colormap=:viridis, colorrange=extrema(levels))
@@ -166,7 +166,7 @@ function rrmserr(z, ẑ)
 end
 function plot_f2(fig, x, y, vals, title, i, levels, show_3d=true, zlabel="z")
     ax = Axis(fig[1, i], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
-    c = contourf!(ax, x, y, vals, color=vals, colormap=:viridis, extendhigh=:auto, levels=levels)
+    c = contourf!(ax, x, y, vals, colormap=:viridis, extendhigh=:auto, levels=levels)
     if show_3d
         ax = Axis3(fig[2, i], xlabel="x", ylabel="y", zlabel=zlabel, width=600, height=600, title=" ", titlealign=:left, azimuth=0.49)
         surface!(ax, x, y, vals, color=vals, colormap=:viridis, colorrange=extrema(levels))

test/doc_examples/interpolation.jl

@@ -23,7 +23,7 @@ sib1_vals = itp(_x, _y; method=Sibson(1))
 
 function plot_itp(fig, x, y, vals, title, i, show_data_sites, itp, xd=nothing, yd=nothing, show_3d=true, levels=-0.1:0.05:0.3)
     ax = Axis(fig[1, i], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
-    c = contourf!(ax, x, y, vals, color=vals, colormap=:viridis, levels=levels, extendhigh=:auto)
+    c = contourf!(ax, x, y, vals, colormap=:viridis, levels=levels, extendhigh=:auto)
     show_data_sites && scatter!(ax, xd, yd, color=:red, markersize=9)
     tri = itp.triangulation
     ch_idx = get_convex_hull_vertices(tri)

test/doc_examples/interpolation_math.jl

@@ -17,7 +17,7 @@ points = [
 ]
 tri = triangulate(points)
 vorn = voronoi(tri)
-fig, ax, sc = voronoiplot(vorn, axis=(width=400, height=400), generator_color=:red)
+fig, ax, sc = voronoiplot(vorn, axis=(width=400, height=400), markercolor=:red)
 resize_to_layout!(fig)
 fig
 
@@ -62,7 +62,7 @@ vorn2 = voronoi(tri2)
 V = get_polygon(vorn2, DelaunayTriangulation.num_points(tri2))
 AX2 = get_area(vorn2, DelaunayTriangulation.num_points(tri2))
 
-fig, ax, sc = voronoiplot(vorn, axis=(width=400, height=400), generator_color=:red, markersize=7, color=:white)
+fig, ax, sc = voronoiplot(vorn, axis=(width=400, height=400), markercolor=:red, markersize=7, color=:white)
 xlims!(ax, 3, 9)
 ylims!(ax, 1.5, 7)
 Vcoords = [get_polygon_point(vorn2, i) for i in V]