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