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