Here’s a fun one, so instead of doing what I normally do(hacking away at code), I figured I’d ask if people had tools to recommend for solving a problem. I know you people working in optimization and geometry live in this world, so my guess is there’s a bunch of tools to make life easier.
Problem statement. I have M linear inequalities pertaining to a finite volume of an N dimensional space. How do I find the vertices of the shape, and the barycenters of the faces?
For example:
0 < X < 1.
0 < Y < 1.
Makes a square. (4 vertices, 4 faces)
But we could
add a constraint…
0 < X < Y
Now we have a traingle (3 vertices, 3 faces)
But now what if I have 20 dimensions and 30 constraints, where are my vertices?
I can imagine a naive approach to doing this… But, it would likely be pretty slow, and I’d like to hear how others could approach the problem (as easily as possible) using the julia ecosystem. Turns out this is a useful thing for a missing component to our ecosystem, but also - a fun algorithm problem.
To get the vertices and rays of the polytope from the inequality representation, you can use Polyhedra.jl. For the barycenters, please let me know if you figure something out. I know that the Chebyshev center of the polytope can be obtained via convex optimization for example. I would first retrieve the facets of the polytope, and then try to formulate the barycenter as an optimization.
If you are interested in computational geometry, take a look at Meshes.jl. I am actively working there.
Cool! Never seen polyhedra.jl. I’ll take a crack at it using it - might ask some questions along the way. Also will likely lead to a PR to another package.
Yes I think that the naive thing is best for the barycenter calculations. Quickly determining the facets is a neat problem too.
I am very interested in computational geometry, my biggest issue is - finding time to formally learn a lot of it. I’ll star Meshes.jl though and keep track of it.
Hey the good thing about comp. geom. is it get’s used everywhere. Graphics, physics, data stuff, stat’s, pure math/opt, etc. So having a backbone in it will only benefit you.
So I got a MWE for how to get the vertices, really happy with how this comes together:
using Polyhedra, CDDLib, JuMP
model = Model()
@variables(model, begin
x
y
end)
con = @constraints(model, begin
0. <= x <= 1.0
0. <= y <= 1.0
x <= y
end)
poly = polyhedron(model, CDDLib.Library(:exact))
points(vrep(poly))
any idea why I can’t get rays/lines for the example above? I believe it’s due to the use of <= and >= but if you swap them for just the </> oprators everything crashes.
for hs in halfspaces(poly)
hp = Polyhedra.hyperplane(hs)
facet = intersect(poly, hp)
end
You can get the Chebyshev center with chebyshevcenter(facet) but the barycenter is not implemented.
The backtracking algorithm you are mentioning is implemented in https://github.com/JuliaPolyhedra/LRSLib.jl
I just tried the chebryshevcenter and it gave me an error but - that’s okay. I just whipped up a little script:
plane_point_map = Dict()
for (h,halfspace) in enumerate(halfspaces(hrep(poly))), (p,point) in enumerate(points(vrep(poly) ) )
if ( in(point,Polyhedra.hyperplane(halfspace)) )
if haskey(plane_point_map, h)
push!(plane_point_map[h], p)
else
plane_point_map[h] = [p]
end
end
end
gives the points on each plane/facet. So basically you could take the column means of the coordinates(with > 1 point) and get the barycenter of each face? Not a “fast” algorithm but - speed isn’t a huge concern yet and it could be put into embarassingly parallel no problem.
I’ll try to switch to the LRSLib backend and see if I notice anything different. The CDDLib backend is working great so far, but is a little slow to load on the first go.
Looks like I was able to get what I wanted. Code quality isn’t top notch but I’m just playing with an idea for now.
using Polyhedra, CDDLib, JuMP, Statistics
model = Model()
@variables(model, begin
x
y
end)
@constraints(model, begin
0 <= x <= 1.0
0 <= y <= 1.0
x <= y
end)
print(model)
poly = polyhedron(model, CDDLib.Library())
plane_point_map = Dict()
for (h,halfspace) in enumerate(halfspaces(hrep(poly))), (p,point) in enumerate(points(vrep(poly) ) )
if ( in(point,Polyhedra.hyperplane(halfspace)) )
if haskey(plane_point_map, h)
push!(plane_point_map[h], p)
else
plane_point_map[h] = [p]
end
end
end
plane_point_map
pnts = hcat(points( vrep( poly ) )...)'
vars = last(size(pnts))
barycenters = []
for (_,v) in plane_point_map
if length( v ) > 1
push!( barycenters, mean(pnts[v,:], dims = 1) )
end
end
barycenters = vcat( barycenters... )
pnts
using Plots
scatter(pnts[:,1], pnts[:,2], legend = false)
scatter!(barycenters[:,1], barycenters[:,2], legend = false)
You can use incidentpoints to simplify the code.
CDD computes the incidence information and incidentpoints just retrieve this information. In your code sample, you recompute it.
using Polyhedra, CDDLib, JuMP, Statistics
model = Model()
@variables(model, begin
x
y
end)
@constraints(model, begin
0 <= x <= 1.0
0 <= y <= 1.0
x <= y
end)
print(model)
poly = polyhedron(model, CDDLib.Library())
removehredundancy!(poly)
barycenters = [mean(incidentpoints(poly, hidx)) for hidx in eachindex(halfspaces(poly))]
I get the output
Feasibility
Subject to
x - y ≤ 0.0
x ∈ [0.0, 1.0]
y ∈ [0.0, 1.0]
3-element Array{Array{Float64,1},1}:
[0.5, 0.5]
[0.0, 0.5]
[0.5, 1.0]
Also: try to avoid things like hcat(points( vrep( poly ) )...) or vcat( barycenters... ). If you have millions of points, it will compile a hcat or vcat with millions of arguments which will be quite slow to compile !
I did try to dereference/unsafe_load the pointers from the representation to get a 2 array, but there’s some issues with that. I think I can use them just fine as a list of 1-arrays, I do like what you did for the barycenters.
Thank you for sharing this snippit this is a really nice example of how to use the outputs from the types in the package. I think I’m going to use this code to build the skeleton of a package today. If you replied in this thread, please message me if you want a special form of credit otherwise I will give a shout out in the readme. I’ll let you all know when the git is made and all of that.
A curiosity: I have never used CDDLib before and I’m wondering how far I could push it as a teaching aid. For example could it handle 50 or 100 or a couple of hundred constraints and demonstrate to students the task that LP solvers must confront (in terms of number of vertices of the polyhedron)?
It really depends on the dimension, you can handle billions of constraints in 2D as there is an n log n algorithm in that case. For higher dimension, it depends also on the number of vertices you expect to have. See http://cgm.cs.mcgill.ca/~avis/doc/avis/ABS96a.ps for a discussion on this
