# Faces & Vertices of N Linear Constraints
**URL:** https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347
**Category:** Numerics
**Created:** [October 31, 2020, 1:08am UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347 "2020-10-31T01:08:43Z")
**Posts on this page:** 16
**Page:** 1
### Author:  [@anon92994695](https://discourse.julialang.org/u/anon92994695)
#### Post date: [October 31, 2020, 1:08am UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/1 "2020-10-31T01:08:43Z")
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.
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [October 31, 2020, 2:26am UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/2 "2020-10-31T02:26:44Z")
To get the vertices and rays of the polytope from the inequality representation, you can use [Polyhedra.jl](https://github.com/JuliaPolyhedra/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](https://github.com/JuliaGeometry/Meshes.jl). I am actively working there.
---
### Author:  [@anon92994695](https://discourse.julialang.org/u/anon92994695)
#### Post date: [October 31, 2020, 2:34am UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/3 "2020-10-31T02:34:28Z")
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.
I just found this paper, didn’t get a chance to read it yet: [http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.579.9132&rep=rep1&type=pdf](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.579.9132&rep=rep1&type=pdf) “Analysis of backtrack algorithms for listing all vertices and all faces of a convex polyhedron”. [Computational Geometry](https://www.sciencedirect.com/science/journal/09257721) [Volume 8, Issue 1](https://www.sciencedirect.com/science/journal/09257721/8/1), June 1997, Pages 1-12.
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.
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [October 31, 2020, 2:42am UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/4 "2020-10-31T02:42:07Z")
Thanks for sharing the papers. My research is requiring a lot of computational geometry lately, so I have no excuse to find the time. 🙂
---
### Author:  [@anon92994695](https://discourse.julialang.org/u/anon92994695)
#### Post date: [October 31, 2020, 11:26am UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/5 "2020-10-31T11:26:40Z")
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:
```nohighlight
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))
```
---
### Author:  [@anon92994695](https://discourse.julialang.org/u/anon92994695)
#### Post date: [October 31, 2020, 12:54pm UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/6 "2020-10-31T12:54:02Z")
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.
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [October 31, 2020, 1:41pm UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/7 "2020-10-31T13:41:27Z")
I am not using Polyhedra.jl, I think you will have a better chance opening an issue there.
---
### Author:  [@blegat](https://discourse.julialang.org/u/blegat)
#### Post date: [October 31, 2020, 5:20pm UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/8 "2020-10-31T17:20:37Z")
To get each facet, you can do
```julia
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](https://github.com/JuliaPolyhedra/LRSLib.jl)
---
### Author:  [@anon92994695](https://discourse.julialang.org/u/anon92994695)
#### Post date: [October 31, 2020, 5:36pm UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/9 "2020-10-31T17:36:54Z")
I just tried the chebryshevcenter and it gave me an error but - that’s okay. I just whipped up a little script:
```nohighlight
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.
---
### Author:  [@anon92994695](https://discourse.julialang.org/u/anon92994695)
#### Post date: [October 31, 2020, 9:58pm UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/10 "2020-10-31T21:58:08Z")
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.
```nohighlight
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)
```
---
### Author:  [@blegat](https://discourse.julialang.org/u/blegat)
#### Post date: [November 1, 2020, 9:53am UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/11 "2020-11-01T09:53:01Z")
You can use [`incidentpoints`](https://juliapolyhedra.github.io/Polyhedra.jl/stable/polyhedron/#Polyhedra.incidentpoints) to simplify the code.
CDD computes the incidence information and `incidentpoints` just retrieve this information. In your code sample, you recompute it.
```julia
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
```julia
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 !
---
### Author:  [@anon92994695](https://discourse.julialang.org/u/anon92994695)
#### Post date: [November 1, 2020, 11:15am UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/12 "2020-11-01T11:15:04Z")
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.
---
### Author:  [@anon92994695](https://discourse.julialang.org/u/anon92994695)
#### Post date: [November 1, 2020, 7:34pm UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/13 "2020-11-01T19:34:41Z")
Here’s what this has helped me do. [https://github.com/caseykneale/ExtremeVertexDesigns.jl](https://github.com/caseykneale/ExtremeVertexDesigns.jl)
Not a released package - yet, but, it’ll be a workspace for some design of experiments things to accompany: [https://github.com/phrb/ExperimentalDesign.jl](https://github.com/phrb/ExperimentalDesign.jl)
---
### Author:  [@blegat](https://discourse.julialang.org/u/blegat)
#### Post date: [November 2, 2020, 1:31pm UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/14 "2020-11-02T13:31:09Z")
Nice package, good to see the `incidentpoints` features in action! Thanks for the credits in the README 😉
---
### Author:  [@healyp](https://discourse.julialang.org/u/healyp)
#### Post date: [November 3, 2020, 8:25pm UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/15 "2020-11-03T20:25:21Z")
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)?
Míle buíochas.
---
### Author:  [@blegat](https://discourse.julialang.org/u/blegat)
#### Post date: [November 3, 2020, 11:07pm UTC](https://discourse.julialang.org/t/faces-vertices-of-n-linear-constraints/49347/16 "2020-11-03T23:07:29Z")
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](http://cgm.cs.mcgill.ca/~avis/doc/avis/ABS96a.ps) for a discussion on this