9bfc64b29c1f696f2a10acfba5cf8baa5a5ebd67
[linpy.git] / pypol / 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 math
20 import numbers
21
22 from . import islhelper
23
24 from .islhelper import mainctx, libisl
25 from .geometry import GeometricObject, Point
26 from .linexprs import Expression, Rational
27 from .domains import Domain
28
29
30 __all__ = [
31 'Polyhedron',
32 'Lt', 'Le', 'Eq', 'Ne', 'Ge', 'Gt',
33 'Empty', 'Universe',
34 ]
35
36
37 class Polyhedron(Domain):
38
39 __slots__ = (
40 '_equalities',
41 '_inequalities',
42 '_constraints',
43 '_symbols',
44 '_dimension',
45 )
46
47 def __new__(cls, equalities=None, inequalities=None):
48 if isinstance(equalities, str):
49 if inequalities is not None:
50 raise TypeError('too many arguments')
51 return cls.fromstring(equalities)
52 elif isinstance(equalities, GeometricObject):
53 if inequalities is not None:
54 raise TypeError('too many arguments')
55 return equalities.aspolyhedron()
56 if equalities is None:
57 equalities = []
58 else:
59 for i, equality in enumerate(equalities):
60 if not isinstance(equality, Expression):
61 raise TypeError('equalities must be linear expressions')
62 equalities[i] = equality.scaleint()
63 if inequalities is None:
64 inequalities = []
65 else:
66 for i, inequality in enumerate(inequalities):
67 if not isinstance(inequality, Expression):
68 raise TypeError('inequalities must be linear expressions')
69 inequalities[i] = inequality.scaleint()
70 symbols = cls._xsymbols(equalities + inequalities)
71 islbset = cls._toislbasicset(equalities, inequalities, symbols)
72 return cls._fromislbasicset(islbset, symbols)
73
74 @property
75 def equalities(self):
76 """
77 Return a list of the equalities in a set.
78 """
79 return self._equalities
80
81 @property
82 def inequalities(self):
83 """
84 Return a list of the inequalities in a set.
85 """
86 return self._inequalities
87
88 @property
89 def constraints(self):
90 """
91 Return ta list of the constraints of a set.
92 """
93 return self._constraints
94
95 @property
96 def polyhedra(self):
97 return self,
98
99 def disjoint(self):
100 """
101 Return a set as disjoint.
102 """
103 return self
104
105 def isuniverse(self):
106 """
107 Return true if a set is the Universe set.
108 """
109 islbset = self._toislbasicset(self.equalities, self.inequalities,
110 self.symbols)
111 universe = bool(libisl.isl_basic_set_is_universe(islbset))
112 libisl.isl_basic_set_free(islbset)
113 return universe
114
115 def aspolyhedron(self):
116 """
117 Return polyhedral hull of a set.
118 """
119 return self
120
121 def __contains__(self, point):
122 if not isinstance(point, Point):
123 raise TypeError('point must be a Point instance')
124 if self.symbols != point.symbols:
125 raise ValueError('arguments must belong to the same space')
126 for equality in self.equalities:
127 if equality.subs(point.coordinates()) != 0:
128 return False
129 for inequality in self.inequalities:
130 if inequality.subs(point.coordinates()) < 0:
131 return False
132 return True
133
134 def subs(self, symbol, expression=None):
135 """
136 Subsitute the given value into an expression and return the resulting
137 expression.
138 """
139 equalities = [equality.subs(symbol, expression)
140 for equality in self.equalities]
141 inequalities = [inequality.subs(symbol, expression)
142 for inequality in self.inequalities]
143 return Polyhedron(equalities, inequalities)
144
145 def _asinequalities(self):
146 inequalities = list(self.equalities)
147 inequalities.extend([-expression for expression in self.equalities])
148 inequalities.extend(self.inequalities)
149 return inequalities
150
151 def widen(self, other):
152 if not isinstance(other, Polyhedron):
153 raise ValueError('argument must be a Polyhedron instance')
154 inequalities1 = self._asinequalities()
155 inequalities2 = other._asinequalities()
156 inequalities = []
157 for inequality1 in inequalities1:
158 if other <= Polyhedron(inequalities=[inequality1]):
159 inequalities.append(inequality1)
160 for inequality2 in inequalities2:
161 for i in range(len(inequalities1)):
162 inequalities3 = inequalities1[:i] + inequalities[i + 1:]
163 inequalities3.append(inequality2)
164 polyhedron3 = Polyhedron(inequalities=inequalities3)
165 if self == polyhedron3:
166 inequalities.append(inequality2)
167 break
168 return Polyhedron(inequalities=inequalities)
169
170 @classmethod
171 def _fromislbasicset(cls, islbset, symbols):
172 islconstraints = islhelper.isl_basic_set_constraints(islbset)
173 equalities = []
174 inequalities = []
175 for islconstraint in islconstraints:
176 constant = libisl.isl_constraint_get_constant_val(islconstraint)
177 constant = islhelper.isl_val_to_int(constant)
178 coefficients = {}
179 for index, symbol in enumerate(symbols):
180 coefficient = libisl.isl_constraint_get_coefficient_val(islconstraint,
181 libisl.isl_dim_set, index)
182 coefficient = islhelper.isl_val_to_int(coefficient)
183 if coefficient != 0:
184 coefficients[symbol] = coefficient
185 expression = Expression(coefficients, constant)
186 if libisl.isl_constraint_is_equality(islconstraint):
187 equalities.append(expression)
188 else:
189 inequalities.append(expression)
190 libisl.isl_basic_set_free(islbset)
191 self = object().__new__(Polyhedron)
192 self._equalities = tuple(equalities)
193 self._inequalities = tuple(inequalities)
194 self._constraints = tuple(equalities + inequalities)
195 self._symbols = cls._xsymbols(self._constraints)
196 self._dimension = len(self._symbols)
197 return self
198
199 @classmethod
200 def _toislbasicset(cls, equalities, inequalities, symbols):
201 dimension = len(symbols)
202 indices = {symbol: index for index, symbol in enumerate(symbols)}
203 islsp = libisl.isl_space_set_alloc(mainctx, 0, dimension)
204 islbset = libisl.isl_basic_set_universe(libisl.isl_space_copy(islsp))
205 islls = libisl.isl_local_space_from_space(islsp)
206 for equality in equalities:
207 isleq = libisl.isl_equality_alloc(libisl.isl_local_space_copy(islls))
208 for symbol, coefficient in equality.coefficients():
209 islval = str(coefficient).encode()
210 islval = libisl.isl_val_read_from_str(mainctx, islval)
211 index = indices[symbol]
212 isleq = libisl.isl_constraint_set_coefficient_val(isleq,
213 libisl.isl_dim_set, index, islval)
214 if equality.constant != 0:
215 islval = str(equality.constant).encode()
216 islval = libisl.isl_val_read_from_str(mainctx, islval)
217 isleq = libisl.isl_constraint_set_constant_val(isleq, islval)
218 islbset = libisl.isl_basic_set_add_constraint(islbset, isleq)
219 for inequality in inequalities:
220 islin = libisl.isl_inequality_alloc(libisl.isl_local_space_copy(islls))
221 for symbol, coefficient in inequality.coefficients():
222 islval = str(coefficient).encode()
223 islval = libisl.isl_val_read_from_str(mainctx, islval)
224 index = indices[symbol]
225 islin = libisl.isl_constraint_set_coefficient_val(islin,
226 libisl.isl_dim_set, index, islval)
227 if inequality.constant != 0:
228 islval = str(inequality.constant).encode()
229 islval = libisl.isl_val_read_from_str(mainctx, islval)
230 islin = libisl.isl_constraint_set_constant_val(islin, islval)
231 islbset = libisl.isl_basic_set_add_constraint(islbset, islin)
232 return islbset
233
234 @classmethod
235 def fromstring(cls, string):
236 domain = Domain.fromstring(string)
237 if not isinstance(domain, Polyhedron):
238 raise ValueError('non-polyhedral expression: {!r}'.format(string))
239 return domain
240
241 def __repr__(self):
242 strings = []
243 for equality in self.equalities:
244 strings.append('Eq({}, 0)'.format(equality))
245 for inequality in self.inequalities:
246 strings.append('Ge({}, 0)'.format(inequality))
247 if len(strings) == 1:
248 return strings[0]
249 else:
250 return 'And({})'.format(', '.join(strings))
251
252
253 def _repr_latex_(self):
254 strings = []
255 for equality in self.equalities:
256 strings.append('{} = 0'.format(equality._repr_latex_().strip('$')))
257 for inequality in self.inequalities:
258 strings.append('{} \\ge 0'.format(inequality._repr_latex_().strip('$')))
259 return '$${}$$'.format(' \\wedge '.join(strings))
260
261 @classmethod
262 def fromsympy(cls, expr):
263 """
264 Convert a sympy object to an expression.
265 """
266 domain = Domain.fromsympy(expr)
267 if not isinstance(domain, Polyhedron):
268 raise ValueError('non-polyhedral expression: {!r}'.format(expr))
269 return domain
270
271 def tosympy(self):
272 """
273 Return an expression as a sympy object.
274 """
275 import sympy
276 constraints = []
277 for equality in self.equalities:
278 constraints.append(sympy.Eq(equality.tosympy(), 0))
279 for inequality in self.inequalities:
280 constraints.append(sympy.Ge(inequality.tosympy(), 0))
281 return sympy.And(*constraints)
282
283
284 class EmptyType(Polyhedron):
285
286 __slots__ = Polyhedron.__slots__
287
288 def __new__(cls):
289 self = object().__new__(cls)
290 self._equalities = (Rational(1),)
291 self._inequalities = ()
292 self._constraints = self._equalities
293 self._symbols = ()
294 self._dimension = 0
295 return self
296
297 def widen(self, other):
298 if not isinstance(other, Polyhedron):
299 raise ValueError('argument must be a Polyhedron instance')
300 return other
301
302 def __repr__(self):
303 return 'Empty'
304
305 def _repr_latex_(self):
306 return '$$\\emptyset$$'
307
308 Empty = EmptyType()
309
310
311 class UniverseType(Polyhedron):
312
313 __slots__ = Polyhedron.__slots__
314
315 def __new__(cls):
316 self = object().__new__(cls)
317 self._equalities = ()
318 self._inequalities = ()
319 self._constraints = ()
320 self._symbols = ()
321 self._dimension = ()
322 return self
323
324 def __repr__(self):
325 return 'Universe'
326
327 def _repr_latex_(self):
328 return '$$\\Omega$$'
329
330 Universe = UniverseType()
331
332
333 def _polymorphic(func):
334 @functools.wraps(func)
335 def wrapper(left, right):
336 if not isinstance(left, Expression):
337 if isinstance(left, numbers.Rational):
338 left = Rational(left)
339 else:
340 raise TypeError('left must be a a rational number '
341 'or a linear expression')
342 if not isinstance(right, Expression):
343 if isinstance(right, numbers.Rational):
344 right = Rational(right)
345 else:
346 raise TypeError('right must be a a rational number '
347 'or a linear expression')
348 return func(left, right)
349 return wrapper
350
351 @_polymorphic
352 def Lt(left, right):
353 """
354 Assert first set is less than the second set.
355 """
356 return Polyhedron([], [right - left - 1])
357
358 @_polymorphic
359 def Le(left, right):
360 """
361 Assert first set is less than or equal to the second set.
362 """
363 return Polyhedron([], [right - left])
364
365 @_polymorphic
366 def Eq(left, right):
367 """
368 Assert first set is equal to the second set.
369 """
370 return Polyhedron([left - right], [])
371
372 @_polymorphic
373 def Ne(left, right):
374 """
375 Assert first set is not equal to the second set.
376 """
377 return ~Eq(left, right)
378
379 @_polymorphic
380 def Gt(left, right):
381 """
382 Assert first set is greater than the second set.
383 """
384 return Polyhedron([], [left - right - 1])
385
386 @_polymorphic
387 def Ge(left, right):
388 """
389 Assert first set is greater than or equal to the second set.
390 """
391 return Polyhedron([], [left - right])