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,
###
### 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