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.
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.
@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!
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
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
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.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")
FastTransforms.jl has a fast transform for triangle OPs:
https://github.com/JuliaApproximation/FastTransforms.jl/blob/master/examples/triangle.jl
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
Note that docs are generated for examples via Literate.jl (for prettier output) Calculus on the reference triangle · FastTransforms.jl