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