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