Implement Polyhedron.fromstring
[linpy.git] / pypol / linear.py
1 import ast
2 import functools
3 import numbers
4 import re
5
6 from fractions import Fraction, gcd
7
8 from pypol import isl
9 from pypol.isl import libisl
10
11
12 __all__ = [
13 'Expression', 'Constant', 'Symbol', 'symbols',
14 'eq', 'le', 'lt', 'ge', 'gt',
15 'Polyhedron',
16 'Empty', 'Universe'
17 ]
18
19
20 def _polymorphic_method(func):
21 @functools.wraps(func)
22 def wrapper(a, b):
23 if isinstance(b, Expression):
24 return func(a, b)
25 if isinstance(b, numbers.Rational):
26 b = Constant(b)
27 return func(a, b)
28 return NotImplemented
29 return wrapper
30
31 def _polymorphic_operator(func):
32 # A polymorphic operator should call a polymorphic method, hence we just
33 # have to test the left operand.
34 @functools.wraps(func)
35 def wrapper(a, b):
36 if isinstance(a, numbers.Rational):
37 a = Constant(a)
38 return func(a, b)
39 elif isinstance(a, Expression):
40 return func(a, b)
41 raise TypeError('arguments must be linear expressions')
42 return wrapper
43
44
45 _main_ctx = isl.Context()
46
47
48 class Expression:
49 """
50 This class implements linear expressions.
51 """
52
53 def __new__(cls, coefficients=None, constant=0):
54 if isinstance(coefficients, str):
55 if constant:
56 raise TypeError('too many arguments')
57 return cls.fromstring(coefficients)
58 if isinstance(coefficients, dict):
59 coefficients = coefficients.items()
60 if coefficients is None:
61 return Constant(constant)
62 coefficients = [(symbol, coefficient)
63 for symbol, coefficient in coefficients if coefficient != 0]
64 if len(coefficients) == 0:
65 return Constant(constant)
66 elif len(coefficients) == 1 and constant == 0:
67 symbol, coefficient = coefficients[0]
68 if coefficient == 1:
69 return Symbol(symbol)
70 self = object().__new__(cls)
71 self._coefficients = {}
72 for symbol, coefficient in coefficients:
73 if isinstance(symbol, Symbol):
74 symbol = symbol.name
75 elif not isinstance(symbol, str):
76 raise TypeError('symbols must be strings or Symbol instances')
77 if isinstance(coefficient, Constant):
78 coefficient = coefficient.constant
79 if not isinstance(coefficient, numbers.Rational):
80 raise TypeError('coefficients must be rational numbers or Constant instances')
81 self._coefficients[symbol] = coefficient
82 if isinstance(constant, Constant):
83 constant = constant.constant
84 if not isinstance(constant, numbers.Rational):
85 raise TypeError('constant must be a rational number or a Constant instance')
86 self._constant = constant
87 self._symbols = tuple(sorted(self._coefficients))
88 self._dimension = len(self._symbols)
89 return self
90
91 @classmethod
92 def _fromast(cls, node):
93 if isinstance(node, ast.Module):
94 assert len(node.body) == 1
95 return cls._fromast(node.body[0])
96 elif isinstance(node, ast.Expr):
97 return cls._fromast(node.value)
98 elif isinstance(node, ast.Name):
99 return Symbol(node.id)
100 elif isinstance(node, ast.Num):
101 return Constant(node.n)
102 elif isinstance(node, ast.UnaryOp):
103 if isinstance(node.op, ast.USub):
104 return -cls._fromast(node.operand)
105 elif isinstance(node, ast.BinOp):
106 left = cls._fromast(node.left)
107 right = cls._fromast(node.right)
108 if isinstance(node.op, ast.Add):
109 return left + right
110 elif isinstance(node.op, ast.Sub):
111 return left - right
112 elif isinstance(node.op, ast.Mult):
113 return left * right
114 elif isinstance(node.op, ast.Div):
115 return left / right
116 raise SyntaxError('invalid syntax')
117
118 @classmethod
119 def fromstring(cls, string):
120 string = re.sub(r'(\d+|\))\s*([^\W\d_]\w*|\()', r'\1*\2', string)
121 tree = ast.parse(string, 'eval')
122 return cls._fromast(tree)
123
124 @property
125 def symbols(self):
126 return self._symbols
127
128 @property
129 def dimension(self):
130 return self._dimension
131
132 def coefficient(self, symbol):
133 if isinstance(symbol, Symbol):
134 symbol = str(symbol)
135 elif not isinstance(symbol, str):
136 raise TypeError('symbol must be a string or a Symbol instance')
137 try:
138 return self._coefficients[symbol]
139 except KeyError:
140 return 0
141
142 __getitem__ = coefficient
143
144 def coefficients(self):
145 for symbol in self.symbols:
146 yield symbol, self.coefficient(symbol)
147
148 @property
149 def constant(self):
150 return self._constant
151
152 def isconstant(self):
153 return False
154
155 def values(self):
156 for symbol in self.symbols:
157 yield self.coefficient(symbol)
158 yield self.constant
159
160 def issymbol(self):
161 return False
162
163 def __bool__(self):
164 return True
165
166 def __pos__(self):
167 return self
168
169 def __neg__(self):
170 return self * -1
171
172 @_polymorphic_method
173 def __add__(self, other):
174 coefficients = dict(self.coefficients())
175 for symbol, coefficient in other.coefficients():
176 if symbol in coefficients:
177 coefficients[symbol] += coefficient
178 else:
179 coefficients[symbol] = coefficient
180 constant = self.constant + other.constant
181 return Expression(coefficients, constant)
182
183 __radd__ = __add__
184
185 @_polymorphic_method
186 def __sub__(self, other):
187 coefficients = dict(self.coefficients())
188 for symbol, coefficient in other.coefficients():
189 if symbol in coefficients:
190 coefficients[symbol] -= coefficient
191 else:
192 coefficients[symbol] = -coefficient
193 constant = self.constant - other.constant
194 return Expression(coefficients, constant)
195
196 def __rsub__(self, other):
197 return -(self - other)
198
199 @_polymorphic_method
200 def __mul__(self, other):
201 if other.isconstant():
202 coefficients = dict(self.coefficients())
203 for symbol in coefficients:
204 coefficients[symbol] *= other.constant
205 constant = self.constant * other.constant
206 return Expression(coefficients, constant)
207 if isinstance(other, Expression) and not self.isconstant():
208 raise ValueError('non-linear expression: '
209 '{} * {}'.format(self._parenstr(), other._parenstr()))
210 return NotImplemented
211
212 __rmul__ = __mul__
213
214 @_polymorphic_method
215 def __truediv__(self, other):
216 if other.isconstant():
217 coefficients = dict(self.coefficients())
218 for symbol in coefficients:
219 coefficients[symbol] = \
220 Fraction(coefficients[symbol], other.constant)
221 constant = Fraction(self.constant, other.constant)
222 return Expression(coefficients, constant)
223 if isinstance(other, Expression):
224 raise ValueError('non-linear expression: '
225 '{} / {}'.format(self._parenstr(), other._parenstr()))
226 return NotImplemented
227
228 def __rtruediv__(self, other):
229 if isinstance(other, self):
230 if self.isconstant():
231 constant = Fraction(other, self.constant)
232 return Expression(constant=constant)
233 else:
234 raise ValueError('non-linear expression: '
235 '{} / {}'.format(other._parenstr(), self._parenstr()))
236 return NotImplemented
237
238 def __str__(self):
239 string = ''
240 i = 0
241 for symbol in self.symbols:
242 coefficient = self.coefficient(symbol)
243 if coefficient == 1:
244 if i == 0:
245 string += symbol
246 else:
247 string += ' + {}'.format(symbol)
248 elif coefficient == -1:
249 if i == 0:
250 string += '-{}'.format(symbol)
251 else:
252 string += ' - {}'.format(symbol)
253 else:
254 if i == 0:
255 string += '{}*{}'.format(coefficient, symbol)
256 elif coefficient > 0:
257 string += ' + {}*{}'.format(coefficient, symbol)
258 else:
259 assert coefficient < 0
260 coefficient *= -1
261 string += ' - {}*{}'.format(coefficient, symbol)
262 i += 1
263 constant = self.constant
264 if constant != 0 and i == 0:
265 string += '{}'.format(constant)
266 elif constant > 0:
267 string += ' + {}'.format(constant)
268 elif constant < 0:
269 constant *= -1
270 string += ' - {}'.format(constant)
271 if string == '':
272 string = '0'
273 return string
274
275 def _parenstr(self, always=False):
276 string = str(self)
277 if not always and (self.isconstant() or self.issymbol()):
278 return string
279 else:
280 return '({})'.format(string)
281
282 def __repr__(self):
283 string = '{}({{'.format(self.__class__.__name__)
284 for i, (symbol, coefficient) in enumerate(self.coefficients()):
285 if i != 0:
286 string += ', '
287 string += '{!r}: {!r}'.format(symbol, coefficient)
288 string += '}}, {!r})'.format(self.constant)
289 return string
290
291 @_polymorphic_method
292 def __eq__(self, other):
293 # "normal" equality
294 # see http://docs.sympy.org/dev/tutorial/gotchas.html#equals-signs
295 return isinstance(other, Expression) and \
296 self._coefficients == other._coefficients and \
297 self.constant == other.constant
298
299 def __hash__(self):
300 return hash((tuple(sorted(self._coefficients.items())), self._constant))
301
302 def _toint(self):
303 lcm = functools.reduce(lambda a, b: a*b // gcd(a, b),
304 [value.denominator for value in self.values()])
305 return self * lcm
306
307 @_polymorphic_method
308 def _eq(self, other):
309 return Polyhedron(equalities=[(self - other)._toint()])
310
311 @_polymorphic_method
312 def __le__(self, other):
313 return Polyhedron(inequalities=[(other - self)._toint()])
314
315 @_polymorphic_method
316 def __lt__(self, other):
317 return Polyhedron(inequalities=[(other - self)._toint() - 1])
318
319 @_polymorphic_method
320 def __ge__(self, other):
321 return Polyhedron(inequalities=[(self - other)._toint()])
322
323 @_polymorphic_method
324 def __gt__(self, other):
325 return Polyhedron(inequalities=[(self - other)._toint() - 1])
326
327
328 class Constant(Expression):
329
330 def __new__(cls, numerator=0, denominator=None):
331 self = object().__new__(cls)
332 if denominator is None:
333 if isinstance(numerator, numbers.Rational):
334 self._constant = numerator
335 elif isinstance(numerator, Constant):
336 self._constant = numerator.constant
337 else:
338 raise TypeError('constant must be a rational number or a Constant instance')
339 else:
340 self._constant = Fraction(numerator, denominator)
341 self._coefficients = {}
342 self._symbols = ()
343 self._dimension = 0
344 return self
345
346 def isconstant(self):
347 return True
348
349 def __bool__(self):
350 return bool(self.constant)
351
352 def __repr__(self):
353 return '{}({!r})'.format(self.__class__.__name__, self._constant)
354
355
356 class Symbol(Expression):
357
358 def __new__(cls, name):
359 if isinstance(name, Symbol):
360 name = name.name
361 elif not isinstance(name, str):
362 raise TypeError('name must be a string or a Symbol instance')
363 self = object().__new__(cls)
364 self._coefficients = {name: 1}
365 self._constant = 0
366 self._symbols = tuple(name)
367 self._name = name
368 self._dimension = 1
369 return self
370
371 @property
372 def name(self):
373 return self._name
374
375 def issymbol(self):
376 return True
377
378 def __repr__(self):
379 return '{}({!r})'.format(self.__class__.__name__, self._name)
380
381 def symbols(names):
382 if isinstance(names, str):
383 names = names.replace(',', ' ').split()
384 return (Symbol(name) for name in names)
385
386
387 @_polymorphic_operator
388 def eq(a, b):
389 return a.__eq__(b)
390
391 @_polymorphic_operator
392 def le(a, b):
393 return a.__le__(b)
394
395 @_polymorphic_operator
396 def lt(a, b):
397 return a.__lt__(b)
398
399 @_polymorphic_operator
400 def ge(a, b):
401 return a.__ge__(b)
402
403 @_polymorphic_operator
404 def gt(a, b):
405 return a.__gt__(b)
406
407
408 class Polyhedron:
409 """
410 This class implements polyhedrons.
411 """
412
413 def __new__(cls, equalities=None, inequalities=None):
414 if isinstance(equalities, str):
415 if inequalities is not None:
416 raise TypeError('too many arguments')
417 return cls.fromstring(equalities)
418 self = super().__new__(cls)
419 self._equalities = []
420 if equalities is not None:
421 for constraint in equalities:
422 for value in constraint.values():
423 if value.denominator != 1:
424 raise TypeError('non-integer constraint: '
425 '{} == 0'.format(constraint))
426 self._equalities.append(constraint)
427 self._equalities = tuple(self._equalities)
428 self._inequalities = []
429 if inequalities is not None:
430 for constraint in inequalities:
431 for value in constraint.values():
432 if value.denominator != 1:
433 raise TypeError('non-integer constraint: '
434 '{} <= 0'.format(constraint))
435 self._inequalities.append(constraint)
436 self._inequalities = tuple(self._inequalities)
437 self._constraints = self._equalities + self._inequalities
438 self._symbols = set()
439 for constraint in self._constraints:
440 self.symbols.update(constraint.symbols)
441 self._symbols = tuple(sorted(self._symbols))
442 return self
443
444 @classmethod
445 def fromstring(cls, string):
446 string = string.strip()
447 string = re.sub(r'^\{\s*|\s*\}$', '', string)
448 string = re.sub(r'([^<=>])=([^<=>])', r'\1==\2', string)
449 string = re.sub(r'(\d+|\))\s*([^\W\d_]\w*|\()', r'\1*\2', string)
450 equalities = []
451 inequalities = []
452 for cstr in re.split(r',|;|and|&&|/\\|∧', string, flags=re.I):
453 tree = ast.parse(cstr.strip(), 'eval')
454 if not isinstance(tree, ast.Module) or len(tree.body) != 1:
455 raise SyntaxError('invalid syntax')
456 node = tree.body[0]
457 if not isinstance(node, ast.Expr):
458 raise SyntaxError('invalid syntax')
459 node = node.value
460 if not isinstance(node, ast.Compare):
461 raise SyntaxError('invalid syntax')
462 left = Expression._fromast(node.left)
463 for i in range(len(node.ops)):
464 op = node.ops[i]
465 right = Expression._fromast(node.comparators[i])
466 if isinstance(op, ast.Lt):
467 inequalities.append(right - left - 1)
468 elif isinstance(op, ast.LtE):
469 inequalities.append(right - left)
470 elif isinstance(op, ast.Eq):
471 equalities.append(left - right)
472 elif isinstance(op, ast.GtE):
473 inequalities.append(left - right)
474 elif isinstance(op, ast.Gt):
475 inequalities.append(left - right - 1)
476 else:
477 raise SyntaxError('invalid syntax')
478 left = right
479 return cls(equalities, inequalities)
480
481 @property
482 def equalities(self):
483 return self._equalities
484
485 @property
486 def inequalities(self):
487 return self._inequalities
488
489 @property
490 def constraints(self):
491 return self._constraints
492
493 @property
494 def symbols(self):
495 return self._symbols
496
497 @property
498 def dimension(self):
499 return len(self.symbols)
500
501 def __bool__(self):
502 return not self.is_empty()
503
504 def __contains__(self, value):
505 # is the value in the polyhedron?
506 raise NotImplementedError
507
508 def __eq__(self, other):
509 # works correctly when symbols is not passed
510 # should be equal if values are the same even if symbols are different
511 bset = self._toisl()
512 other = other._toisl()
513 return bool(libisl.isl_basic_set_plain_is_equal(bset, other))
514
515 def isempty(self):
516 bset = self._toisl()
517 return bool(libisl.isl_basic_set_is_empty(bset))
518
519 def isuniverse(self):
520 bset = self._toisl()
521 return bool(libisl.isl_basic_set_is_universe(bset))
522
523 def isdisjoint(self, other):
524 # return true if the polyhedron has no elements in common with other
525 #symbols = self._symbolunion(other)
526 bset = self._toisl()
527 other = other._toisl()
528 return bool(libisl.isl_set_is_disjoint(bset, other))
529
530 def issubset(self, other):
531 # check if self(bset) is a subset of other
532 symbols = self._symbolunion(other)
533 bset = self._toisl(symbols)
534 other = other._toisl(symbols)
535 return bool(libisl.isl_set_is_strict_subset(other, bset))
536
537 def __le__(self, other):
538 return self.issubset(other)
539
540 def __lt__(self, other):
541 symbols = self._symbolunion(other)
542 bset = self._toisl(symbols)
543 other = other._toisl(symbols)
544 return bool(libisl.isl_set_is_strict_subset(other, bset))
545
546 def issuperset(self, other):
547 # test whether every element in other is in the polyhedron
548 raise NotImplementedError
549
550 def __ge__(self, other):
551 return self.issuperset(other)
552
553 def __gt__(self, other):
554 symbols = self._symbolunion(other)
555 bset = self._toisl(symbols)
556 other = other._toisl(symbols)
557 bool(libisl.isl_set_is_strict_subset(other, bset))
558 raise NotImplementedError
559
560 def union(self, *others):
561 # return a new polyhedron with elements from the polyhedron and all
562 # others (convex union)
563 raise NotImplementedError
564
565 def __or__(self, other):
566 return self.union(other)
567
568 def intersection(self, *others):
569 # return a new polyhedron with elements common to the polyhedron and all
570 # others
571 # a poor man's implementation could be:
572 # equalities = list(self.equalities)
573 # inequalities = list(self.inequalities)
574 # for other in others:
575 # equalities.extend(other.equalities)
576 # inequalities.extend(other.inequalities)
577 # return self.__class__(equalities, inequalities)
578 raise NotImplementedError
579
580 def __and__(self, other):
581 return self.intersection(other)
582
583 def difference(self, other):
584 # return a new polyhedron with elements in the polyhedron that are not in the other
585 symbols = self._symbolunion(other)
586 bset = self._toisl(symbols)
587 other = other._toisl(symbols)
588 difference = libisl.isl_set_subtract(bset, other)
589 return difference
590
591 def __sub__(self, other):
592 return self.difference(other)
593
594 def __str__(self):
595 constraints = []
596 for constraint in self.equalities:
597 constraints.append('{} == 0'.format(constraint))
598 for constraint in self.inequalities:
599 constraints.append('{} >= 0'.format(constraint))
600 return '{{{}}}'.format(', '.join(constraints))
601
602 def __repr__(self):
603 if self.isempty():
604 return 'Empty'
605 elif self.isuniverse():
606 return 'Universe'
607 else:
608 equalities = list(self.equalities)
609 inequalities = list(self.inequalities)
610 return '{}(equalities={!r}, inequalities={!r})' \
611 ''.format(self.__class__.__name__, equalities, inequalities)
612
613 def _symbolunion(self, *others):
614 symbols = set(self.symbols)
615 for other in others:
616 symbols.update(other.symbols)
617 return sorted(symbols)
618
619 def _toisl(self, symbols=None):
620 if symbols is None:
621 symbols = self.symbols
622 dimension = len(symbols)
623 space = libisl.isl_space_set_alloc(_main_ctx, 0, dimension)
624 bset = libisl.isl_basic_set_universe(libisl.isl_space_copy(space))
625 ls = libisl.isl_local_space_from_space(space)
626 for equality in self.equalities:
627 ceq = libisl.isl_equality_alloc(libisl.isl_local_space_copy(ls))
628 for symbol, coefficient in equality.coefficients():
629 val = str(coefficient).encode()
630 val = libisl.isl_val_read_from_str(_main_ctx, val)
631 dim = symbols.index(symbol)
632 ceq = libisl.isl_constraint_set_coefficient_val(ceq, libisl.isl_dim_set, dim, val)
633 if equality.constant != 0:
634 val = str(equality.constant).encode()
635 val = libisl.isl_val_read_from_str(_main_ctx, val)
636 ceq = libisl.isl_constraint_set_constant_val(ceq, val)
637 bset = libisl.isl_basic_set_add_constraint(bset, ceq)
638 for inequality in self.inequalities:
639 cin = libisl.isl_inequality_alloc(libisl.isl_local_space_copy(ls))
640 for symbol, coefficient in inequality.coefficients():
641 val = str(coefficient).encode()
642 val = libisl.isl_val_read_from_str(_main_ctx, val)
643 dim = symbols.index(symbol)
644 cin = libisl.isl_constraint_set_coefficient_val(cin, libisl.isl_dim_set, dim, val)
645 if inequality.constant != 0:
646 val = str(inequality.constant).encode()
647 val = libisl.isl_val_read_from_str(_main_ctx, val)
648 cin = libisl.isl_constraint_set_constant_val(cin, val)
649 bset = libisl.isl_basic_set_add_constraint(bset, cin)
650 bset = isl.BasicSet(bset)
651 return bset
652
653 @classmethod
654 def _fromisl(cls, bset, symbols):
655 raise NotImplementedError
656 equalities = ...
657 inequalities = ...
658 return cls(equalities, inequalities)
659 '''takes basic set in isl form and puts back into python version of polyhedron
660 isl example code gives isl form as:
661 "{[i] : exists (a : i = 2a and i >= 10 and i <= 42)}")
662 our printer is giving form as:
663 { [i0, i1] : 2i1 >= -2 - i0 } '''
664
665 Empty = eq(0,1)
666 Universe = Polyhedron()
667
668 if __name__ == '__main__':
669 p1 = Polyhedron('2a + 2b + 1 == 0') # empty
670 print(p1._toisl())
671 p2 = Polyhedron('3x + 2y + 3 == 0') # not empty
672 print(p2._toisl())