# Package/algorithm for barycentric grid interpolation
**URL:** https://discourse.julialang.org/t/package-algorithm-for-barycentric-grid-interpolation/50803
**Category:** Numerics
**Tags:** question
**Created:** [November 26, 2020, 8:58am UTC](https://discourse.julialang.org/t/package-algorithm-for-barycentric-grid-interpolation/50803 "2020-11-26T08:58:29Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 26, 2020, 8:58am UTC](https://discourse.julialang.org/t/package-algorithm-for-barycentric-grid-interpolation/50803/1 "2020-11-26T08:58:29Z")
I am looking for a package or an algorithm (which I am happy to code into a package and make it available) that expresses a coordinate inside a simplex as a convex combination “near” of points with integer coordinates.
Formally, for a given integer K \in \mathbb{N}, consider the simplex
S = \{ x \in \mathbb{R}^M \mid x\_i \ge 0\ \forall i, \sum\_{i=0}^M x\_i \le K \}
defined by the sum of coordinates being \le K, and the points with integer coordinates
I = S \cap \mathbb{N}^M
For an z \in S, I am looking for vertices v\_i \in I, i = 0, \dots, M and weights \sum\_{i = 0}^M \alpha\_i = 1 such that
z = \sum\_{i = 0}^M \alpha\_i v\_i
Ideally the v\_i should be “near” z, but not necessarily the nearest if that’s difficult.
The problem is very easy for M = 2, but the naive algorithms I came up with fail for M = 3 — 2D is very easy to partition with self-similar simplexes like this:
But, importantly, for this algorithm the simplexes formed by the v\_i for each z do not have to form a partition of S.
Here is some test code I would like to pass, if anyone wants to provide concrete code:
```julia
using LinearAlgebra, Test
f(z, K) = IMPLEMENTATION # returns a vector of vs and αs
ϵ = 1e-8 # just picked some arbitrary tolerance
for _ in 1:100
M = rand(3:5)
K = rand(2:10)
z = normalize(rand(M), 1) .* K
vs, αs = f(z, K)
ẑ = mapreduce((v, α) -> v .* α, (v1, v2) -> v1 .+ v2, vs, αs)
@test ẑ ≈ z atol = ϵ
@test sum(αs) ≈ 1 atol = ϵ
@test all(αs .≥ 0)
for v in vs
@test eltype(v) == Int # or all v are integers, can be stored as floats
@test all(v .≥ 0)
@test sum(v) ≤ K
end
end
```
But, again, I appreciate any kind of suggestions.
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [November 26, 2020, 9:32am UTC](https://discourse.julialang.org/t/package-algorithm-for-barycentric-grid-interpolation/50803/2 "2020-11-26T09:32:16Z")
Are you looking for something like
“Wasserstein Barycentric Coordinates:
Histogram Regression Using Optimal Transport”
Or Are you looking for the barycenter itself?
Barycenters under and optimal transport cost are relatively straightforward to calculate, especially using the sinkhorn algorithm.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 26, 2020, 11:45am UTC](https://discourse.julialang.org/t/package-algorithm-for-barycentric-grid-interpolation/50803/3 "2020-11-26T11:45:10Z")
I am looking for the solution of the problem above. I read the article, but I am not sure how it is related (sorry for being dense, this is not my area of expertise).
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 26, 2020, 4:28pm UTC](https://discourse.julialang.org/t/package-algorithm-for-barycentric-grid-interpolation/50803/5 "2020-11-26T16:28:08Z")
I found a nice thesis about these topics, [Moore (1992)](https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.73.7084&rep=rep1&type=pdf). I may be able to employ the symmetric subdivision algorithm (Section 2.2.3) using the Kuhn triangulation, I just need to understand it first.
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [November 26, 2020, 4:37pm UTC](https://discourse.julialang.org/t/package-algorithm-for-barycentric-grid-interpolation/50803/6 "2020-11-26T16:37:47Z")
If I understand your problem, the solution may be a simplicial triangulation of “points” with integer coordinates. The problem is tricky with the Delaunay triangulation, since the triangulation is not unique. However, a regular mesh like this may be readily generated in 2 and 3 dimensions using a cell template. You mentioned Kuhn triangulation: in 3d that is one possibility. [GitHub - PetrKryslUCSD/FinEtools.jl: Finite Element tools in Julia](https://github.com/PetrKryslUCSD/FinEtools.jl) can generate these kinds of meshes, with different templates for the cells. In particular, refer to [the lines below](https://github.com/PetrKryslUCSD/FinEtools.jl/blob/ce81150e3176ec3d70ae25586813111380218bd3/src/MeshTetrahedronModule.jl#L76).
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 27, 2020, 8:15am UTC](https://discourse.julialang.org/t/package-algorithm-for-barycentric-grid-interpolation/50803/7 "2020-11-27T08:15:42Z")
> [@PetrKryslUCSD](#):
>
> In particular, refer to [the lines below](https://github.com/PetrKryslUCSD/FinEtools.jl/blob/ce81150e3176ec3d70ae25586813111380218bd3/src/MeshTetrahedronModule.jl#L76).
Thanks — I am may be able to generalize this to \mathbb{R}^N.
Is there a book/introductory article you would recommend about this topic? (symmetric subdivision of simplexes).
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 27, 2020, 3:16pm UTC](https://discourse.julialang.org/t/package-algorithm-for-barycentric-grid-interpolation/50803/8 "2020-11-27T15:16:35Z")
I think I figured it out. The key is to transform to what Coxeter calls R-simplexes, where the Kuhn triangulation will always be in the large simplex by construction.
I include a _pedagogical_ implementation below if anyone is interested — a much nicer, optimized, non-allocating, and efficient version will be made available in a package soon (which accompanies a paper that motivated this whole thing) — the `vs` of course do not need to be constructed, the permutation is sufficient. I will of course link in both here. In the meantime,
```julia
###
### horribly inefficient pedagogical implementation
###
StoR(z) = reverse(cumsum(reverse(z))) # did I say horribily inefficient?
RtoS(x) = vcat(.- diff(x), x[end:end])
function f(z, K)
# go to canonical R-simplex
c = StoR(z)
# fractional and integer parts
fv = map(x -> ((f, i) = modf(x); (f, Int(i))), c)
# “lower” corner
v0 = last.(fv)
# fractional part as a convex combination, Kuhn triangulation
fi = sort!(collect(enumerate(first.(fv))); by = last, rev = true)
# gridpoints mutated buffer for walking the vertices
vv = copy(v0)
vs = [(vv[fi[1]] += 1; copy(vv)) for fi in fi]
push!(vs, v0)
# weights
αs = RtoS((last.(fi)))
push!(αs, 1 - sum(αs))
# back to inflated standard simplex
map(RtoS, vs), αs
end
ϵ = 1e-8 # just picked some arbitrary tolerance
for _ in 1:1000
M = rand(3:5)
K = rand(2:10)
z = normalize(rand(M), 1) .* (K * rand())
vs, αs = f(z, K)
ẑ = mapreduce((v, α) -> v .* α, (v1, v2) -> v1 .+ v2, vs, αs)
@test ẑ ≈ z atol = ϵ
@test sum(αs) ≈ 1 atol = ϵ
@test all(αs .≥ 0)
@test length(αs) == length(vs) == M + 1
for v in vs
@test eltype(v) == Int
@test length(v) == M
@test all(v .≥ 0)
@test sum(v) ≤ K
end
end
```
BTW, the book I have been looking for is
> **[Regular Polytopes (book)](https://en.wikipedia.org/wiki/Regular_Polytopes_(book))**
>
> Regular Polytopes is a geometry book on regular polytopes written by Harold Scott MacDonald Coxeter. It was originally published by Methuen in 1947 and by Pitman Publishing in 1948, with a second edition published by Macmillan in 1963 and a third edition by Dover Publications in 1973.
> The Basic Library List Committee of the Mathematical Association of America has recommended that it be included in undergraduate mathematics libraries.
> The main topics of the book are the Platonic solids (regular c...