<jupyter_start><jupyter_text>**Objective**. This notebook contains illustrating examples for the utilities in the [polyhedron_tools](https://github.com/mforets/polyhedron_tools) module.<jupyter_code>%display typeset<jupyter_output><empty_output><jupyter_text>1. Modeling with Polyhedra: back and forth with half-space representation We present examples for creating Polyhedra from matrices and conversely to obtain matrices from Polyhedra.<jupyter_code>from polyhedron_tools.misc import polyhedron_from_Hrep, polyhedron_to_Hrep

A = matrix([[-1.0, 0.0,  0.0,  0.0,  0.0,  0.0],
[ 1.0,  0.0,  0.0,  0.0,  0.0,  0.0],
[ 0.0,  1.0,  0.0,  0.0,  0.0,  0.0],
[ 0.0, -1.0,  0.0,  0.0,  0.0,  0.0],
[ 0.0,  0.0, -1.0,  0.0,  0.0,  0.0],
[ 0.0,  0.0,  1.0,  0.0,  0.0,  0.0],
[ 0.0,  0.0,  0.0, -1.0,  0.0,  0.0],
[ 0.0,  0.0,  0.0,  1.0,  0.0,  0.0],
[ 0.0,  0.0,  0.0,  0.0,  1.0,  0.0],
[ 0.0,  0.0,  0.0,  0.0, -1.0,  0.0],
[ 0.0,  0.0,  0.0,  0.0,  0.0,  1.0],
[ 0.0,  0.0,  0.0,  0.0,  0.0, -1.0]])

b = vector([0.0, 10.0, 0.0, 0.0, 0.2, 0.2, 0.1, 0.1, 0.0, 0.0, 0.0, 0.0])

P = polyhedron_from_Hrep(A, b); 
P
P.inequalities()
P.equations()<jupyter_output><empty_output><jupyter_text>It is possible to obtain the matrices that represent the inequality and the equality constraints separately, using the keyword argument `separate_equality_constraints`. This type of information is somtimes useful for optimization solvers.<jupyter_code>[A, b, Aeq, beq] = polyhedron_to_Hrep(P, separate_equality_constraints = True)
A, b
Aeq, beq<jupyter_output><empty_output><jupyter_text>2. Generating polyhedra  Constructing hyperrectangles Let's construct a ball in the infinity norm, specifying the center and radius. We remark that the case of a hyperbox can be done in Sage's library `polytopes.hypercube(n)` where `n` is the dimension. However, as of Sage v7.6. there is no such hyperrectangle function (or n-orthotope, see the [wikipedia page](https://en.wikipedia.org/wiki/Hyperrectangle)), so we use `BoxInfty`.<jupyter_code>from polyhedron_tools.misc import BoxInfty
P = BoxInfty(center=[1,2,3], radius=0.1); P.plot(aspect_ratio=1)<jupyter_output><empty_output><jupyter_text>As a side note, the function also works when the arguments are not named, as in<jupyter_code>P = BoxInfty([1,2,3], 0.1); P<jupyter_output><empty_output><jupyter_text>Another use of `BoxInfty` is to specify the lengths of the sides. For example:<jupyter_code>P = BoxInfty([[0,1], [2,3]]); P<jupyter_output><empty_output><jupyter_text>Random polyhedra The `random_polygon_2d` receives the number of arguments as input and produces a polygon whose vertices are randomly sampled from the unit circle. See the docstring for more another options.<jupyter_code>from polyhedron_tools.polygons import random_polygon_2d

random_polygon_2d(5)<jupyter_output><empty_output><jupyter_text>Opposite polyhedron<jupyter_code>from polyhedron_tools.misc import BoxInfty, opposite_polyhedron

P = BoxInfty([1,1], 0.5);
mp = opposite_polyhedron(P);
P.plot(aspect_ratio=1) + mp.plot(color='red')<jupyter_output><empty_output><jupyter_text>3. Miscelaneous functions  Support function of a polytope<jupyter_code>from polyhedron_tools.misc import support_function

P = BoxInfty([1,2,3,4,5], 1); P
support_function(P, [1,-1,1,-1,1], verbose=1)<jupyter_output>**** Solve LP  ****
Maximization:
  x_0 - x_1 + x_2 - x_3 + x_4 

Constraints:
  x_0 <= 2.0
  x_1 <= 3.0
  x_2 <= 4.0
  x_3 <= 5.0
  x_4 <= 6.0
  - x_0 <= 0.0
  - x_4 <= -4.0
  - x_3 <= -3.0
  - x_2 <= -2.0
  - x_1 <= -1.0
Variables:
  x_0 is a continuous variable (min=-oo, max=+oo)
  x_1 is a continuous variable (min=-oo, max=+oo)
  x_2 is a continuous variable (min=-oo, max=+oo)
  x_3 is a continuous variable (min=-oo, max=+oo)
  x_4 is a continuous variable (min=-oo, max=+oo)
Objective Value: 8.0
x_0 = 2.000000
x_1 = 1.000000
x_2 = 4.000000
x_3 = 3.000000
x_4 = 6.000000<jupyter_text>It is also possible to input the polyhedron in matrix form, $[A, b]$. If this is possible, it is preferable, since it is often faster. Below is an example with $12$ variables. We get beteen 3x and 4x improvement in the second case.<jupyter_code>reset('P, A, b')
P = BoxInfty([1,2,3,4,5,6,7,8,9,10,11,12], 1); P
[A, b] = polyhedron_to_Hrep(P)
timeit('support_function(P, [1,-1,1,-1,1,-1,1,-1,1,-1,1,-1])')
support_function([A, b], [1,-1,1,-1,1,-1,1,-1,1,-1,1,-1])
timeit('support_function([A, b], [1,-1,1,-1,1,-1,1,-1,1,-1,1,-1])')<jupyter_output>625 loops, best of 3: 421 µs per loop<jupyter_text>Support function of an ellipsoid<jupyter_code>from polyhedron_tools.misc import support_function, support_function_ellipsoid 
import random

# Generate a random ellipsoid and check support function outer approximation.
# Define an ellipse as: x^T*Q*x <= 1
M = random_matrix(RR, 2, distribution="uniform")
Q = M.T*M
f = lambda x, y : Q[0,0]*x^2 + Q[1,1]*y^2 + (Q[0,1]+Q[1,0])*x*y-1
E = implicit_plot(f,(-5,5),(-3,3),fill=True,alpha=0.5,plot_points=600)

# generate at random k directions, and compute the overapproximation of E using support functions
# It works 'in average': we might get unbounded domains (random choice did not enclose the ellipsoid).
# It is recommended to use QQ as base_ring to avoid 'frozen set' issues.
k=15
A = matrix(RR,k,2); b = vector(RR,k)
for i in range(k):
    theta = random.uniform(0, 2*pi.n(digits=5))
    d = vector(RR,[cos(theta), sin(theta)])
    s_fun = support_function_ellipsoid(Q, d)
    A.set_row(i,d); b[i] = s_fun

OmegaApprox = polyhedron_from_Hrep(A, b, base_ring = QQ)
E + OmegaApprox.plot(fill=False, color='red')<jupyter_output><empty_output><jupyter_text>Supremum norm of a polyhedron<jupyter_code>from polyhedron_tools.misc import BoxInfty, radius
P = BoxInfty([-13,24,-51,18.54,309],27.04);
radius(P)
got_lengths, got_center_and_radius = False, False
got_lengths is not False
radius(polyhedron_to_Hrep(P))
8401/25.n()<jupyter_output><empty_output><jupyter_text>Consider a higher-dimensional system. We obtain almost a 200x improvement for a 15-dimensional set. This is because in the case we call ```poly_sup_norm``` with a polytope, the ```bounding_box()``` function consumes time.<jupyter_code>%%time 
P = BoxInfty([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1], 14.28);
radius(P)
%%time 
[A, b] = BoxInfty([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1], 14.28, base_ring = RDF, return_HSpaceRep = True)
radius([A, b])<jupyter_output>CPU times: user 20.4 ms, sys: 1.19 ms, total: 21.6 ms
Wall time: 20.9 ms<jupyter_text>The improvement in speed is quite interesting! Moreover, for a 20-dimensional set, the polyhedron construct does not finish.   Linear map of a polyhedron The operation `matrix x polyhedron` can be performed in Sage with the usual command, `*`.<jupyter_code>U = BoxInfty([[0,1],[0,1]])
B = matrix([[1,1],[1,-1]])
P = B * U
P.plot(color='red', alpha=0.5) + U.plot()<jupyter_output><empty_output><jupyter_text>Chebyshev center<jupyter_code>from polyhedron_tools.misc import chebyshev_center, BoxInfty
from polyhedron_tools.polygons import random_polygon_2d

P = random_polygon_2d(10, base_ring = QQ)
c = chebyshev_center(P); 

B = BoxInfty([[1,2],[0,1]])
b = chebyshev_center(B)

fig = point(c, color='blue') + P.plot(color='blue', alpha=0.2)
fig += point(b, color='red') + B.plot(color='red', alpha=0.2)

fig += point(P.center().n(), color='green',marker='x')
fig += point(B.center().n(), color='green',marker='x')

fig<jupyter_output><empty_output><jupyter_text>The method ```center()``` existent in the Polyhedra class, computes the average of the vertices. In contrast, the Chebyshev center is the center of the largest box enclosed by the polytope.<jupyter_code>B.bounding_box()
P.bounding_box()
e = [ (P.bounding_box()[0][0] + P.bounding_box()[1][0])/2, (P.bounding_box()[0][1] + P.bounding_box()[1][1])/2]
l = [[P.bounding_box()[0][0], P.bounding_box()[1][0]], [P.bounding_box()[0][1], P.bounding_box()[1][1]] ]

fig += point(e,color='black') + BoxInfty(lengths=l).plot(alpha=0.1,color='grey')

fig<jupyter_output><empty_output><jupyter_text>Here we have added in grey the bounding box that is obtained from the method ```bounding_box()```. To make the picture complete, we should also add the box of center Cheby center, and of maximal radius which is included in it.  4. Approximate projections<jupyter_code>from polyhedron_tools.projections import lotov_algo
from polyhedron_tools.misc import polyhedron_to_Hrep, polyhedron_from_Hrep 

A = matrix([[-1.0, 0.0,  0.0,  0.0,  0.0,  0.0],
[ 1.0,  0.0,  0.0,  0.0,  0.0,  0.0],
[ 0.0,  1.0,  0.0,  0.0,  0.0,  0.0],
[ 0.0, -1.0,  0.0,  0.0,  0.0,  0.0],
[ 0.0,  0.0, -1.0,  0.0,  0.0,  0.0],
[ 0.0,  0.0,  1.0,  0.0,  0.0,  0.0],
[ 0.0,  0.0,  0.0, -1.0,  0.0,  0.0],
[ 0.0,  0.0,  0.0,  1.0,  0.0,  0.0],
[ 0.0,  0.0,  0.0,  0.0,  1.0,  0.0],
[ 0.0,  0.0,  0.0,  0.0, -1.0,  0.0],
[ 0.0,  0.0,  0.0,  0.0,  0.0,  1.0],
[ 0.0,  0.0,  0.0,  0.0,  0.0, -1.0]])

b = vector([0.0, 10.0, 0.0, 0.0, 0.2, 0.2, 0.1, 0.1, 0.0, 0.0, 0.0, 0.0])

P = polyhedron_from_Hrep(A, b); P
lotov_algo(A, b, [1,0,0], [0,1,0], 0.5)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16954787234042554, 'fraction_numerical': 0.07568705673758866, 'mean_word_length': 4.644152595372107, 'pattern_counts': {'":': 0, '<': 77, '<?xml version=': 0, '>': 66, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '7525', 'license': 'MIT', 'n_tokens_mistral': 4064, 'n_tokens_neox': 3833, 'n_words': 945, 'path': 'examples/misc.ipynb', 'repo_name': 'mforets/polyhedron_tools', 'size': '9024'}