Interpolation on non-rectangular domain

Dear all,

I would like to ask for advice, how to interpolate function on a non-rectangular domain. I am working on a problem, where the domain of my function is 2 or 3 dimensional, where all variables are restricted to be in [0,1] and their sum has to be in [0,1], hence this problem produces a non-rectangular domain. Is there some convenient package for this type of problem? As far as I know (maybe I am wrong), Interpolations.jl works with rectangular domains.

Best,
Honza

Sounds like interpolation on triangles and tetrahedral domains. Do you want piecewise linear or high order polynomial interpolation?

You can perform such interpolations with GeoStats.jl. If you are looking for a lightweight dependency, other members can suggest alternatives.

2 Likes

Higher-order would be nice, but linear is also perfectly ok.

Is there some tutorial about interpolations in your package? What I need is to interpolate values of my function on a non-rectangular grid, and then use this interpolation as a representation of this function in further computations (I am using it in an iterative algorithm that solves for the unknown function). Compile time is ok for me, I am more concerned about runtime.

Say you have arbitrary locations with measurements of a variable z:

julia> using GeoStats
        
julia> data = georef((z=rand(100),), rand(2,100))
100 PointSet{Float64,2}
  variables
    └─z (Float64)

And another set of arbitrary locations where you want to perform the interpolation:

julia> locs = PointSet(rand(2,100))
100 PointSet{Float64,2}
 0.8877317040268236  0.5228607963917429  …  0.7112810163212397
 0.5372754529043251  0.7050593626184458     0.6308336193731663

You can define the “interpolation” problem:

julia> problem = EstimationProblem(data, locs, :z)
2D EstimationProblem
  data:      100 SpatialData{Float64,2}
  domain:    100 PointSet{Float64,2}
  variables: z (Float64)

and pick one of the solvers like for example inverse distance weighting (IDW):

julia> solution = solve(problem, IDW())
2D EstimationSolution
  domain: 100 PointSet{Float64,2}
  variables: z

These implementations should have a decent performance as they are using KDTrees from NearestNeighbors.jl and other tricks. But I never used them in the context of a tight loop like it is usually done in function approximation inside finite element codes for example.

You may also find ApproxFun.jl useful depending on what you want to do.

3 Likes

@juliohm Thank you very much! Would this be efficient for scalar operations? Systems that I am solving are typically implicit systems (macroeconomic discrete time functional equations), and hence I need to use at each grid point nonlinear solver…

If I would wrap EstimationProblem and solve as a scalar function, would it be efficient, or would the overhead be too high?

Best,
Honza

You are welcome @Honza9723. I’d benchmark the code with some typical setup for your use case. You can easily quantify possible bottlenecks if they exist. Take a look at BenchmarkTools.jl. The problem setup is cost-less (it is just a wrapper pointing to the other variables). The solvers have a compilation time, but runtime should be fine. The major downside I see is that GeoStats.jl is a big project to have as a dependency, and we are constantly evolving the API.

That sounds nice. Yeah, the compilation time of GeoStats.jl is quite large, but I can live with it. Thanks again!

1 Like

Sorry for the later reply, was actually working on the package I’m using for this answer. I think @juliohm’s approach is probably better suited (esp if you can’t choose the points at which you interpolate), but I figure I’d include this in case it’s helpful to anyone.

I have a package NodesAndModes.jl for high order polynomials and interpolation points on several shapes. A warning: I found a bug in the main branch, so I’m using a branch where I fixed this by making it more “Julian”.

You can get high order interpolation nodes for a Tri() or Tet() domain via

pkg> add NodesAndModes#dev_elem_types
julia> using NodesAndModes
julia> using Plots
julia> x,y = nodes(Tri(),25) # degree 25 interpolation nodes

This gives some “decent” interpolation nodes.

I do the actual interpolation using linear algebra (implicitly defined Lagrange basis functions). If I want to interpolate, say, a version of Runge’s function in 2D


with a degree N=20 interpolant, I can do this using

using NodesAndModes
N = 20 # polynomial degree
x,y = nodes(Tri(),N)

xp,yp = equi_nodes(Tri(),100) # plotting nodes
f(x,y) = 1/(1 + 5*((x+1/3)^2+(y+1/3)^2)) # Runge-ish function

function interp(elem,N,x_eval,y_eval,x,y,f_vals)
    coeffs = vandermonde(elem,N,x,y)\f_vals
    return vandermonde(elem,N,x_eval,y_eval)*coeffs
end

f_interp = interp(Tri(),N,xp,yp,x,y,f.(x,y))

err = maximum(abs.(f.(xp,yp)-f_interp))
@show err

which gives an error of 0.015041654247017144.

Runge’s function is sort of the worst-case scenario. If you interpolate a smoother function like f(x,y) = sin(pi*x)*sin(pi*y), then the error is 1.167277052793736e-9.

A few notes

  • if the points at which you want to interpolate are fixed, you can just replace the interpolation x,y in the above with your node locations. However, due to the nature of high order interpolation, this can lead to pretty poor results.
  • degree N > 25 starts to get iffy since some of the polynomials get really large and roundoff effects come into play. For many “nice” functions, N between 15-20 is more than enough though.
  • if your evaluation and interpolation points don’t change between iterations, you can reduce the interpolation procedure to just a matrix-vector product, which should be pretty fast.

Finally, for completeness, the plots were generated using
scatter(xp,yp,f_interp,zcolor=f_interp,msw=0,leg=false,title="Interpolant") and scatter(xp,yp,f.(xp,yp),zcolor=f.(xp,yp),msw=0,leg=false,title="Exact")

4 Likes

@jlchan Thank you! I will definitely try both approaches!

FastTransforms.jl has a fast transform for triangle OPs:

Note it’s only technically interpolation if your function is a polynomial (the number of sample points is double the degrees of freedom), but for function approximation it works great even with millions of unknowns

5 Likes

Note that docs are generated for examples via Literate.jl (for prettier output) https://juliaapproximation.github.io/FastTransforms.jl/dev/generated/triangle/

3 Likes