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