Fix 3d plots in examples
[linpy.git] / linpy / polyhedra.py
1 # Copyright 2014 MINES ParisTech
2 #
3 # This file is part of LinPy.
4 #
5 # LinPy is free software: you can redistribute it and/or modify
6 # it under the terms of the GNU General Public License as published by
7 # the Free Software Foundation, either version 3 of the License, or
8 # (at your option) any later version.
9 #
10 # LinPy is distributed in the hope that it will be useful,
11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 # GNU General Public License for more details.
14 #
15 # You should have received a copy of the GNU General Public License
16 # along with LinPy. If not, see <http://www.gnu.org/licenses/>.
17
18 import functools
19 import numbers
20
21 from . import islhelper
22
23 from .domains import Domain
24 from .geometry import GeometricObject, Point
25 from .islhelper import libisl, mainctx
26 from .linexprs import LinExpr, Rational
27
28
29 __all__ = [
30 'Empty',
31 'Eq',
32 'Ge',
33 'Gt',
34 'Le',
35 'Lt',
36 'Ne',
37 'Polyhedron',
38 'Universe',
39 ]
40
41
42 class Polyhedron(Domain):
43 """
44 A convex polyhedron (or simply "polyhedron") is the space defined by a
45 system of linear equalities and inequalities. This space can be unbounded.
46 A Z-polyhedron (simply called "polyhedron" in LinPy) is the set of integer
47 points in a convex polyhedron.
48 """
49
50 __slots__ = (
51 '_equalities',
52 '_inequalities',
53 '_symbols',
54 '_dimension',
55 )
56
57 def __new__(cls, equalities=None, inequalities=None):
58 """
59 Return a polyhedron from two sequences of linear expressions:
60 equalities is a list of expressions equal to 0, and inequalities is a
61 list of expressions greater or equal to 0. For example, the polyhedron
62 0 <= x <= 2, 0 <= y <= 2 can be constructed with:
63
64 >>> x, y = symbols('x y')
65 >>> square1 = Polyhedron([], [x, 2 - x, y, 2 - y])
66 >>> square1
67 And(0 <= x, x <= 2, 0 <= y, y <= 2)
68
69 It may be easier to use comparison operators LinExpr.__lt__(),
70 LinExpr.__le__(), LinExpr.__ge__(), LinExpr.__gt__(), or
71 functions Lt(), Le(), Eq(), Ge() and Gt(), using one of the following
72 instructions:
73
74 >>> x, y = symbols('x y')
75 >>> square1 = (0 <= x) & (x <= 2) & (0 <= y) & (y <= 2)
76 >>> square1 = Le(0, x, 2) & Le(0, y, 2)
77
78 It is also possible to build a polyhedron from a string.
79
80 >>> square1 = Polyhedron('0 <= x <= 2, 0 <= y <= 2')
81
82 Finally, a polyhedron can be constructed from a GeometricObject
83 instance, calling the GeometricObject.aspolyedron() method. This way,
84 it is possible to compute the polyhedral hull of a Domain instance,
85 i.e., the convex hull of two polyhedra:
86
87 >>> square1 = Polyhedron('0 <= x <= 2, 0 <= y <= 2')
88 >>> square2 = Polyhedron('1 <= x <= 3, 1 <= y <= 3')
89 >>> Polyhedron(square1 | square2)
90 And(0 <= x, 0 <= y, x <= y + 2, y <= x + 2, x <= 3, y <= 3)
91 """
92 if isinstance(equalities, str):
93 if inequalities is not None:
94 raise TypeError('too many arguments')
95 return cls.fromstring(equalities)
96 elif isinstance(equalities, GeometricObject):
97 if inequalities is not None:
98 raise TypeError('too many arguments')
99 return equalities.aspolyhedron()
100 sc_equalities = []
101 if equalities is not None:
102 for equality in equalities:
103 if isinstance(equality, LinExpr):
104 sc_equalities.append(equality.scaleint())
105 elif isinstance(equality, numbers.Rational):
106 sc_equalities.append(Rational(equality).scaleint())
107 else:
108 raise TypeError('equalities must be linear expressions '
109 'or rational numbers')
110 sc_inequalities = []
111 if inequalities is not None:
112 for inequality in inequalities:
113 if isinstance(inequality, LinExpr):
114 sc_inequalities.append(inequality.scaleint())
115 elif isinstance(inequality, numbers.Rational):
116 sc_inequalities.append(Rational(inequality).scaleint())
117 else:
118 raise TypeError('inequalities must be linear expressions '
119 'or rational numbers')
120 symbols = cls._xsymbols(sc_equalities + sc_inequalities)
121 islbset = cls._toislbasicset(sc_equalities, sc_inequalities, symbols)
122 return cls._fromislbasicset(islbset, symbols)
123
124 @property
125 def equalities(self):
126 """
127 The tuple of equalities. This is a list of LinExpr instances that are
128 equal to 0 in the polyhedron.
129 """
130 return self._equalities
131
132 @property
133 def inequalities(self):
134 """
135 The tuple of inequalities. This is a list of LinExpr instances that are
136 greater or equal to 0 in the polyhedron.
137 """
138 return self._inequalities
139
140 @property
141 def constraints(self):
142 """
143 The tuple of constraints, i.e., equalities and inequalities. This is
144 semantically equivalent to: equalities + inequalities.
145 """
146 return self._equalities + self._inequalities
147
148 @property
149 def polyhedra(self):
150 return self,
151
152 def make_disjoint(self):
153 return self
154
155 def isuniverse(self):
156 islbset = self._toislbasicset(self.equalities, self.inequalities,
157 self.symbols)
158 universe = bool(libisl.isl_basic_set_is_universe(islbset))
159 libisl.isl_basic_set_free(islbset)
160 return universe
161
162 def aspolyhedron(self):
163 return self
164
165 def convex_union(self, *others):
166 """
167 Return the convex union of two or more polyhedra.
168 """
169 for other in others:
170 if not isinstance(other, Polyhedron):
171 raise TypeError('arguments must be Polyhedron instances')
172 return Polyhedron(self.union(*others))
173
174 def __contains__(self, point):
175 if not isinstance(point, Point):
176 raise TypeError('point must be a Point instance')
177 if self.symbols != point.symbols:
178 raise ValueError('arguments must belong to the same space')
179 for equality in self.equalities:
180 if equality.subs(point.coordinates()) != 0:
181 return False
182 for inequality in self.inequalities:
183 if inequality.subs(point.coordinates()) < 0:
184 return False
185 return True
186
187 def subs(self, symbol, expression=None):
188 equalities = [equality.subs(symbol, expression)
189 for equality in self.equalities]
190 inequalities = [inequality.subs(symbol, expression)
191 for inequality in self.inequalities]
192 return Polyhedron(equalities, inequalities)
193
194 def asinequalities(self):
195 """
196 Express the polyhedron using inequalities, given as a list of
197 expressions greater or equal to 0.
198 """
199 inequalities = list(self.equalities)
200 inequalities.extend([-expression for expression in self.equalities])
201 inequalities.extend(self.inequalities)
202 return inequalities
203
204 def widen(self, other):
205 """
206 Compute the standard widening of two polyhedra, à la Halbwachs.
207
208 In its current implementation, this method is slow and should not be
209 used on large polyhedra.
210 """
211 if not isinstance(other, Polyhedron):
212 raise TypeError('argument must be a Polyhedron instance')
213 inequalities1 = self.asinequalities()
214 inequalities2 = other.asinequalities()
215 inequalities = []
216 for inequality1 in inequalities1:
217 if other <= Polyhedron(inequalities=[inequality1]):
218 inequalities.append(inequality1)
219 for inequality2 in inequalities2:
220 for i in range(len(inequalities1)):
221 inequalities3 = inequalities1[:i] + inequalities[i + 1:]
222 inequalities3.append(inequality2)
223 polyhedron3 = Polyhedron(inequalities=inequalities3)
224 if self == polyhedron3:
225 inequalities.append(inequality2)
226 break
227 return Polyhedron(inequalities=inequalities)
228
229 @classmethod
230 def _fromislbasicset(cls, islbset, symbols):
231 if bool(libisl.isl_basic_set_is_empty(islbset)):
232 return Empty
233 if bool(libisl.isl_basic_set_is_universe(islbset)):
234 return Universe
235 islconstraints = islhelper.isl_basic_set_constraints(islbset)
236 equalities = []
237 inequalities = []
238 for islconstraint in islconstraints:
239 constant = libisl.isl_constraint_get_constant_val(islconstraint)
240 constant = islhelper.isl_val_to_int(constant)
241 coefficients = {}
242 for index, symbol in enumerate(symbols):
243 coefficient = libisl.isl_constraint_get_coefficient_val(
244 islconstraint, libisl.isl_dim_set, index)
245 coefficient = islhelper.isl_val_to_int(coefficient)
246 if coefficient != 0:
247 coefficients[symbol] = coefficient
248 expression = LinExpr(coefficients, constant)
249 if libisl.isl_constraint_is_equality(islconstraint):
250 equalities.append(expression)
251 else:
252 inequalities.append(expression)
253 libisl.isl_basic_set_free(islbset)
254 self = object().__new__(Polyhedron)
255 self._equalities = tuple(equalities)
256 self._inequalities = tuple(inequalities)
257 self._symbols = cls._xsymbols(self.constraints)
258 self._dimension = len(self._symbols)
259 return self
260
261 @classmethod
262 def _toislbasicset(cls, equalities, inequalities, symbols):
263 dimension = len(symbols)
264 indices = {symbol: index for index, symbol in enumerate(symbols)}
265 islsp = libisl.isl_space_set_alloc(mainctx, 0, dimension)
266 islbset = libisl.isl_basic_set_universe(libisl.isl_space_copy(islsp))
267 islls = libisl.isl_local_space_from_space(islsp)
268 for equality in equalities:
269 isleq = libisl.isl_equality_alloc(
270 libisl.isl_local_space_copy(islls))
271 for symbol, coefficient in equality.coefficients():
272 islval = str(coefficient).encode()
273 islval = libisl.isl_val_read_from_str(mainctx, islval)
274 index = indices[symbol]
275 isleq = libisl.isl_constraint_set_coefficient_val(
276 isleq, libisl.isl_dim_set, index, islval)
277 if equality.constant != 0:
278 islval = str(equality.constant).encode()
279 islval = libisl.isl_val_read_from_str(mainctx, islval)
280 isleq = libisl.isl_constraint_set_constant_val(isleq, islval)
281 islbset = libisl.isl_basic_set_add_constraint(islbset, isleq)
282 for inequality in inequalities:
283 islin = libisl.isl_inequality_alloc(
284 libisl.isl_local_space_copy(islls))
285 for symbol, coefficient in inequality.coefficients():
286 islval = str(coefficient).encode()
287 islval = libisl.isl_val_read_from_str(mainctx, islval)
288 index = indices[symbol]
289 islin = libisl.isl_constraint_set_coefficient_val(
290 islin, libisl.isl_dim_set, index, islval)
291 if inequality.constant != 0:
292 islval = str(inequality.constant).encode()
293 islval = libisl.isl_val_read_from_str(mainctx, islval)
294 islin = libisl.isl_constraint_set_constant_val(islin, islval)
295 islbset = libisl.isl_basic_set_add_constraint(islbset, islin)
296 return islbset
297
298 @classmethod
299 def fromstring(cls, string):
300 domain = Domain.fromstring(string)
301 if not isinstance(domain, Polyhedron):
302 raise ValueError('non-polyhedral expression: {!r}'.format(string))
303 return domain
304
305 def __repr__(self):
306 strings = []
307 for equality in self.equalities:
308 left, right, swap = 0, 0, False
309 for i, (symbol, coefficient) in enumerate(equality.coefficients()):
310 if coefficient > 0:
311 left += coefficient * symbol
312 else:
313 right -= coefficient * symbol
314 if i == 0:
315 swap = True
316 if equality.constant > 0:
317 left += equality.constant
318 else:
319 right -= equality.constant
320 if swap:
321 left, right = right, left
322 strings.append('{} == {}'.format(left, right))
323 for inequality in self.inequalities:
324 left, right = 0, 0
325 for symbol, coefficient in inequality.coefficients():
326 if coefficient < 0:
327 left -= coefficient * symbol
328 else:
329 right += coefficient * symbol
330 if inequality.constant < 0:
331 left -= inequality.constant
332 else:
333 right += inequality.constant
334 strings.append('{} <= {}'.format(left, right))
335 if len(strings) == 1:
336 return strings[0]
337 else:
338 return 'And({})'.format(', '.join(strings))
339
340 @classmethod
341 def fromsympy(cls, expression):
342 domain = Domain.fromsympy(expression)
343 if not isinstance(domain, Polyhedron):
344 raise ValueError('non-polyhedral expression: {!r}'.format(
345 expression))
346 return domain
347
348 def tosympy(self):
349 import sympy
350 constraints = []
351 for equality in self.equalities:
352 constraints.append(sympy.Eq(equality.tosympy(), 0))
353 for inequality in self.inequalities:
354 constraints.append(sympy.Ge(inequality.tosympy(), 0))
355 return sympy.And(*constraints)
356
357
358 class EmptyType(Polyhedron):
359 """
360 The empty polyhedron, whose set of constraints is not satisfiable.
361 """
362
363 def __new__(cls):
364 self = object().__new__(cls)
365 self._equalities = (Rational(1),)
366 self._inequalities = ()
367 self._symbols = ()
368 self._dimension = 0
369 return self
370
371 def widen(self, other):
372 if not isinstance(other, Polyhedron):
373 raise ValueError('argument must be a Polyhedron instance')
374 return other
375
376 def __repr__(self):
377 return 'Empty'
378
379 Empty = EmptyType()
380
381
382 class UniverseType(Polyhedron):
383 """
384 The universe polyhedron, whose set of constraints is always satisfiable,
385 i.e. is empty.
386 """
387
388 def __new__(cls):
389 self = object().__new__(cls)
390 self._equalities = ()
391 self._inequalities = ()
392 self._symbols = ()
393 self._dimension = ()
394 return self
395
396 def __repr__(self):
397 return 'Universe'
398
399 Universe = UniverseType()
400
401
402 def _pseudoconstructor(func):
403 @functools.wraps(func)
404 def wrapper(expression1, expression2, *expressions):
405 expressions = (expression1, expression2) + expressions
406 for expression in expressions:
407 if not isinstance(expression, LinExpr):
408 if isinstance(expression, numbers.Rational):
409 expression = Rational(expression)
410 else:
411 raise TypeError('arguments must be rational numbers '
412 'or linear expressions')
413 return func(*expressions)
414 return wrapper
415
416
417 @_pseudoconstructor
418 def Lt(*expressions):
419 """
420 Create the polyhedron with constraints expr1 < expr2 < expr3 ...
421 """
422 inequalities = []
423 for left, right in zip(expressions, expressions[1:]):
424 inequalities.append(right - left - 1)
425 return Polyhedron([], inequalities)
426
427
428 @_pseudoconstructor
429 def Le(*expressions):
430 """
431 Create the polyhedron with constraints expr1 <= expr2 <= expr3 ...
432 """
433 inequalities = []
434 for left, right in zip(expressions, expressions[1:]):
435 inequalities.append(right - left)
436 return Polyhedron([], inequalities)
437
438
439 @_pseudoconstructor
440 def Eq(*expressions):
441 """
442 Create the polyhedron with constraints expr1 == expr2 == expr3 ...
443 """
444 equalities = []
445 for left, right in zip(expressions, expressions[1:]):
446 equalities.append(left - right)
447 return Polyhedron(equalities, [])
448
449
450 @_pseudoconstructor
451 def Ne(*expressions):
452 """
453 Create the domain such that expr1 != expr2 != expr3 ... The result is a
454 Domain object, not a Polyhedron.
455 """
456 domain = Universe
457 for left, right in zip(expressions, expressions[1:]):
458 domain &= ~Eq(left, right)
459 return domain
460
461
462 @_pseudoconstructor
463 def Ge(*expressions):
464 """
465 Create the polyhedron with constraints expr1 >= expr2 >= expr3 ...
466 """
467 inequalities = []
468 for left, right in zip(expressions, expressions[1:]):
469 inequalities.append(left - right)
470 return Polyhedron([], inequalities)
471
472
473 @_pseudoconstructor
474 def Gt(*expressions):
475 """
476 Create the polyhedron with constraints expr1 > expr2 > expr3 ...
477 """
478 inequalities = []
479 for left, right in zip(expressions, expressions[1:]):
480 inequalities.append(left - right - 1)
481 return Polyhedron([], inequalities)