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