# How to transform a derived set of inequation in the good polyhedron format

This is a follow-up to:

Applying this code

def one_dimension_less(representationH):
x = list(var('x_%i' % i) for i in (0..2))
sol = [solve_ineq(ieq, [x[0]]) for ieq in Ineq]
sol_r = flatten([x[1] for x in sol])
ineq_ge = [z for z in sol_r if z.rhs() == x[0]]
ineq_le = [z for z in sol_r if z.lhs() == x[0]]
ineq_a = [z for z in sol_r if z.lhs() != x[0] and z.rhs() != x[0]]
ineq_aa = [(z.lhs()).full_simplify() >= 0 for z in ineq_a]
result = flatten(ineq_aa + [
[(ineq_le[j].rhs()-ineq_ge[i].lhs()).full_simplify() >= 0
for i in range(len(ineq_ge))]
for  j in range(len(ineq_le))])
# result1 = [result[i].factor() for i in range(len(result))]
return result


to

D = polytopes.dodecahedron()
DH = D.Hrepresentation()


with

PD1 = one_dimension_less(DH)
show(PD1)


I can eliminate one variable and obtain a set of inequalities defining a polyhedron in a lower dimension (here I go from 3 to 2 but it can be higher.

For each x in PD1, I can isolate its right hand side but I don't know how to isolate the parameters in such a way to write it as

sage: Polyhedron(ieqs=[(0, 1, 0), (0, 0, 1), (1, -1, -1)]).Hrepresentation()
(An inequality (-1, -1) x + 1 >= 0,
An inequality (1, 0) x + 0 >= 0,
An inequality (0, 1) x + 0 >= 0)


the lists needed to feed ieqs in Polyhedron. How to do that?

edit retag close merge delete

1

Please, when asking a question, check that copying and pasting the provided code in a fresh Sage session allows to reproduce the behaviour you are asking about.

When running the code provided, I get

sage: PD1 = one_dimension_less(DH)
Traceback (most recent call last)
...
NameError: name 'Ineq' is not defined

( 2021-05-31 13:56:56 +0200 )edit

For convenience, I added links to the stream of questions this is part of.

( 2021-05-31 14:08:35 +0200 )edit

Maybe you meant this:

def one_dimension_less(representationH):
x = vector(SR, SR.var('x_', 3))
ieqs = [[el.A()*x + el.b() >= 0] for el in DH]
sol = [solve_ineq(ieq, [x[0]]) for ieq in ieqs]
sol_r = flatten([x[1] for x in sol])
ineq_ge = [z for z in sol_r if z.rhs() == x[0]]
ineq_le = [z for z in sol_r if z.lhs() == x[0]]
ineq_a = [z for z in sol_r if z.lhs() != x[0] and z.rhs() != x[0]]
ineq_aa = [(z.lhs()).full_simplify() >= 0 for z in ineq_a]
result = flatten(ineq_aa + [
[(ineq_le[j].rhs() - ineq_ge[i].lhs()).full_simplify() >= 0
for i in range(len(ineq_ge))]
for  j in range(len(ineq_le))])
return result

( 2021-05-31 14:24:52 +0200 )edit

Slelievre I am sincerely sorry for Ineq. Jupyter keep in memory the first tries and as the command was working I have not been enough vigilent. From now on, I will try to restard the kernel before sending a question.

( 2021-05-31 17:03:52 +0200 )edit

Sort by » oldest newest most voted

Not sure I completely understand the question.

The following should at least be a step in the right direction.

Define one_dimension_less as intended in the question.

def one_dimension_less(representationH):
x = vector(SR, SR.var('x_', 3))
ieqs = [[el.A()*x + el.b() >= 0] for el in DH]
sol = [solve_ineq(ieq, [x[0]]) for ieq in ieqs]
sol_r = flatten([x[1] for x in sol])
ineq_ge = [z for z in sol_r if z.rhs() == x[0]]
ineq_le = [z for z in sol_r if z.lhs() == x[0]]
ineq_a = [z for z in sol_r if z.lhs() != x[0] and z.rhs() != x[0]]
ineq_aa = [(z.lhs()).full_simplify() >= 0 for z in ineq_a]
result = flatten(ineq_aa + [
[(ineq_le[j].rhs() - ineq_ge[i].lhs()).full_simplify() >= 0
for i in range(len(ineq_ge))]
for  j in range(len(ineq_le))])
return result


Use it:

sage: D = polytopes.dodecahedron()
sage: DH = D.Hrepresentation()
sage: PD1 = one_dimension_less(DH)
sage: PD1
[-x_2*(sqrt(5) + 1) - 2*x_1 + 4 >= 0,
-x_2*(sqrt(5) + 1) + 2*x_1 + 4 >= 0,
x_2*(sqrt(5) + 1) - 2*x_1 + 4 >= 0,
x_2*(sqrt(5) + 1) + 2*x_1 + 4 >= 0,
-x_1*(sqrt(5) + 1) + 4 >= 0,
4 >= 0,
-(2*x_1*(sqrt(5) + 2) + x_2*(sqrt(5) + 1) - 4*sqrt(5) - 8)/(sqrt(5) + 3) >= 0,
-(2*x_1*(sqrt(5) + 2) - x_2*(sqrt(5) + 1) - 4*sqrt(5) - 8)/(sqrt(5) + 3) >= 0,
4 >= 0,
x_1*(sqrt(5) + 1) + 4 >= 0,
(2*x_1*(sqrt(5) + 2) - x_2*(sqrt(5) + 1) + 4*sqrt(5) + 8)/(sqrt(5) + 3) >= 0,
(2*x_1*(sqrt(5) + 2) + x_2*(sqrt(5) + 1) + 4*sqrt(5) + 8)/(sqrt(5) + 3) >= 0,
-(2*x_1*(sqrt(5) + 2) + x_2*(sqrt(5) + 1) - 4*sqrt(5) - 8)/(sqrt(5) + 3) >= 0,
(2*x_1*(sqrt(5) + 2) - x_2*(sqrt(5) + 1) + 4*sqrt(5) + 8)/(sqrt(5) + 3) >= 0,
-2*(x_2*(sqrt(5) + 1) - 2*sqrt(5) - 2)/(sqrt(5) + 3) >= 0,
8/(sqrt(5) + 1) >= 0,
-(2*x_1*(sqrt(5) + 2) - x_2*(sqrt(5) + 1) - 4*sqrt(5) - 8)/(sqrt(5) + 3) >= 0,
(2*x_1*(sqrt(5) + 2) + x_2*(sqrt(5) + 1) + 4*sqrt(5) + 8)/(sqrt(5) + 3) >= 0,
8/(sqrt(5) + 1) >= 0,
2*(x_2*(sqrt(5) + 1) + 2*sqrt(5) + 2)/(sqrt(5) + 3) >= 0]


Examine one of the inequalities:

sage: ieq = PD1[0]
sage: ieq
-x_2*(sqrt(5) + 1) - 2*x_1 + 4 >= 0


Grab its left hand side:

sage: lhs = ieq.lhs()
sage: lhs
-x_2*(sqrt(5) + 1) - 2*x_1 + 4


Get the constant coefficient:

sage: lhs(x_0=0, x_1=0, x_2=0)
4


Same with a substitution dictionary:

sage: lhs({x_0: 0, x_1: 0, x_2: 0})
4


Get the coefficients of the variables:

sage: lhs.coefficient(x_0)
0
sage: lhs.coefficient(x_1)
-2
sage: lhs.coefficient(x_2)
-sqrt(5) - 1


All in one go:

sage: x = vector(SR.var('x_', 3))
sage: A = vector(lhs.coefficient(xi) for xi in x)
sage: A
(0, -2, -sqrt(5) - 1)

sage: b = lhs({xi: 0 for xi in x})
sage: b
4


Reform the inequation:

sage: ieq_tuple = (b,) + tuple(A)
sage: ieq_tuple
(4, 0, -2, -sqrt(5) - 1)


Make that into a function:

def tuple_from_ieq(ieq, x):
r"""
Return the tuple for an inequation.

INPUT:

- ieq -- an inequation, assumed to be of the form
a_0 * x_0 + ... + a_n * x_n + b >= 0
- x -- a multidimensional variable
as a list, tuple or vector (x_0, ..., x_n)

OUTPUT: the tuple (b, a_0, ..., a_n).
"""
b = lhs({xi: 0 for xi in x})
A = vector(lhs.coefficient(xi) for xi in x)
return (b,) + tuple(A)


Use it:

sage: x = vector(SR.var('x_', 3))
sage: t = tuple_from_ieq(PD1[0], x)
sage: t
(4, 0, -2, -sqrt(5) - 1)


Convert to algebraic numbers:

sage: tt = vector(AA, t)
sage: tt
(4, 0, -2, -3.236067977499790?)

sage: Polyhedron(ieqs=[tt])
A 3-dimensional polyhedron in AA^3 defined as
the convex hull of 1 vertex, 1 ray, 2 lines

more

Let me change the perspective: to project a polytope (or apply any linear map to it), you can directly apply a matrix to your polytope, for example:

sage: D = polytopes.dodecahedron()
sage: M = matrix([[0,1,0],[0,0,1]])
sage: M * D
A 2-dimensional polyhedron
in (Number Field in sqrt5 with defining polynomial x^2 - 5 with sqrt5 = 2.236...?)^2
defined as the convex hull of 6 vertices
sage: (M * D).plot()
`
more