X-Git-Url: https://scm.cri.ensmp.fr/git/linpy.git/blobdiff_plain/66e41ccd173874b3309e4c24be64c7b6b2ac6298..1d494bb187b70135df721c13306d7f26fdf33f50:/pypol/linear.py diff --git a/pypol/linear.py b/pypol/linear.py deleted file mode 100644 index b40415f..0000000 --- a/pypol/linear.py +++ /dev/null @@ -1,784 +0,0 @@ -import ast -import functools -import numbers -import re - -from fractions import Fraction, gcd - -from . import isl -from .isl import libisl - - -__all__ = [ - 'Expression', 'Constant', 'Symbol', 'symbols', - 'eq', 'le', 'lt', 'ge', 'gt', - 'Polyhedron', - 'Empty', 'Universe' -] - - -def _polymorphic_method(func): - @functools.wraps(func) - def wrapper(a, b): - if isinstance(b, Expression): - return func(a, b) - if isinstance(b, numbers.Rational): - b = Constant(b) - return func(a, b) - return NotImplemented - return wrapper - -def _polymorphic_operator(func): - # A polymorphic operator should call a polymorphic method, hence we just - # have to test the left operand. - @functools.wraps(func) - def wrapper(a, b): - if isinstance(a, numbers.Rational): - a = Constant(a) - return func(a, b) - elif isinstance(a, Expression): - return func(a, b) - raise TypeError('arguments must be linear expressions') - return wrapper - - -_main_ctx = isl.Context() - - -class Expression: - """ - This class implements linear expressions. - """ - - __slots__ = ( - '_coefficients', - '_constant', - '_symbols', - '_dimension', - ) - - def __new__(cls, coefficients=None, constant=0): - if isinstance(coefficients, str): - if constant: - raise TypeError('too many arguments') - return cls.fromstring(coefficients) - if isinstance(coefficients, dict): - coefficients = coefficients.items() - if coefficients is None: - return Constant(constant) - coefficients = [(symbol, coefficient) - for symbol, coefficient in coefficients if coefficient != 0] - if len(coefficients) == 0: - return Constant(constant) - elif len(coefficients) == 1 and constant == 0: - symbol, coefficient = coefficients[0] - if coefficient == 1: - return Symbol(symbol) - self = object().__new__(cls) - self._coefficients = {} - for symbol, coefficient in coefficients: - if isinstance(symbol, Symbol): - symbol = symbol.name - elif not isinstance(symbol, str): - raise TypeError('symbols must be strings or Symbol instances') - if isinstance(coefficient, Constant): - coefficient = coefficient.constant - if not isinstance(coefficient, numbers.Rational): - raise TypeError('coefficients must be rational numbers or Constant instances') - self._coefficients[symbol] = coefficient - if isinstance(constant, Constant): - constant = constant.constant - if not isinstance(constant, numbers.Rational): - raise TypeError('constant must be a rational number or a Constant instance') - self._constant = constant - self._symbols = tuple(sorted(self._coefficients)) - self._dimension = len(self._symbols) - return self - - @classmethod - def _fromast(cls, node): - if isinstance(node, ast.Module) and len(node.body) == 1: - return cls._fromast(node.body[0]) - elif isinstance(node, ast.Expr): - return cls._fromast(node.value) - elif isinstance(node, ast.Name): - return Symbol(node.id) - elif isinstance(node, ast.Num): - return Constant(node.n) - elif isinstance(node, ast.UnaryOp) and isinstance(node.op, ast.USub): - return -cls._fromast(node.operand) - elif isinstance(node, ast.BinOp): - left = cls._fromast(node.left) - right = cls._fromast(node.right) - if isinstance(node.op, ast.Add): - return left + right - elif isinstance(node.op, ast.Sub): - return left - right - elif isinstance(node.op, ast.Mult): - return left * right - elif isinstance(node.op, ast.Div): - return left / right - raise SyntaxError('invalid syntax') - - @classmethod - def fromstring(cls, string): - string = re.sub(r'(\d+|\))\s*([^\W\d_]\w*|\()', r'\1*\2', string) - tree = ast.parse(string, 'eval') - return cls._fromast(tree) - - @property - def symbols(self): - return self._symbols - - @property - def dimension(self): - return self._dimension - - def coefficient(self, symbol): - if isinstance(symbol, Symbol): - symbol = str(symbol) - elif not isinstance(symbol, str): - raise TypeError('symbol must be a string or a Symbol instance') - try: - return self._coefficients[symbol] - except KeyError: - return 0 - - __getitem__ = coefficient - - def coefficients(self): - for symbol in self.symbols: - yield symbol, self.coefficient(symbol) - - @property - def constant(self): - return self._constant - - def isconstant(self): - return False - - def values(self): - for symbol in self.symbols: - yield self.coefficient(symbol) - yield self.constant - - def issymbol(self): - return False - - def __bool__(self): - return True - - def __pos__(self): - return self - - def __neg__(self): - return self * -1 - - @_polymorphic_method - def __add__(self, other): - coefficients = dict(self.coefficients()) - for symbol, coefficient in other.coefficients(): - if symbol in coefficients: - coefficients[symbol] += coefficient - else: - coefficients[symbol] = coefficient - constant = self.constant + other.constant - return Expression(coefficients, constant) - - __radd__ = __add__ - - @_polymorphic_method - def __sub__(self, other): - coefficients = dict(self.coefficients()) - for symbol, coefficient in other.coefficients(): - if symbol in coefficients: - coefficients[symbol] -= coefficient - else: - coefficients[symbol] = -coefficient - constant = self.constant - other.constant - return Expression(coefficients, constant) - - def __rsub__(self, other): - return -(self - other) - - @_polymorphic_method - def __mul__(self, other): - if other.isconstant(): - coefficients = dict(self.coefficients()) - for symbol in coefficients: - coefficients[symbol] *= other.constant - constant = self.constant * other.constant - return Expression(coefficients, constant) - if isinstance(other, Expression) and not self.isconstant(): - raise ValueError('non-linear expression: ' - '{} * {}'.format(self._parenstr(), other._parenstr())) - return NotImplemented - - __rmul__ = __mul__ - - @_polymorphic_method - def __truediv__(self, other): - if other.isconstant(): - coefficients = dict(self.coefficients()) - for symbol in coefficients: - coefficients[symbol] = \ - Fraction(coefficients[symbol], other.constant) - constant = Fraction(self.constant, other.constant) - return Expression(coefficients, constant) - if isinstance(other, Expression): - raise ValueError('non-linear expression: ' - '{} / {}'.format(self._parenstr(), other._parenstr())) - return NotImplemented - - def __rtruediv__(self, other): - if isinstance(other, self): - if self.isconstant(): - constant = Fraction(other, self.constant) - return Expression(constant=constant) - else: - raise ValueError('non-linear expression: ' - '{} / {}'.format(other._parenstr(), self._parenstr())) - return NotImplemented - - def __str__(self): - string = '' - i = 0 - for symbol in self.symbols: - coefficient = self.coefficient(symbol) - if coefficient == 1: - if i == 0: - string += symbol - else: - string += ' + {}'.format(symbol) - elif coefficient == -1: - if i == 0: - string += '-{}'.format(symbol) - else: - string += ' - {}'.format(symbol) - else: - if i == 0: - string += '{}*{}'.format(coefficient, symbol) - elif coefficient > 0: - string += ' + {}*{}'.format(coefficient, symbol) - else: - assert coefficient < 0 - coefficient *= -1 - string += ' - {}*{}'.format(coefficient, symbol) - i += 1 - constant = self.constant - if constant != 0 and i == 0: - string += '{}'.format(constant) - elif constant > 0: - string += ' + {}'.format(constant) - elif constant < 0: - constant *= -1 - string += ' - {}'.format(constant) - if string == '': - string = '0' - return string - - def _parenstr(self, always=False): - string = str(self) - if not always and (self.isconstant() or self.issymbol()): - return string - else: - return '({})'.format(string) - - def __repr__(self): - return '{}({!r})'.format(self.__class__.__name__, str(self)) - - @_polymorphic_method - def __eq__(self, other): - # "normal" equality - # see http://docs.sympy.org/dev/tutorial/gotchas.html#equals-signs - return isinstance(other, Expression) and \ - self._coefficients == other._coefficients and \ - self.constant == other.constant - - def __hash__(self): - return hash((tuple(sorted(self._coefficients.items())), self._constant)) - - def _toint(self): - lcm = functools.reduce(lambda a, b: a*b // gcd(a, b), - [value.denominator for value in self.values()]) - return self * lcm - - @_polymorphic_method - def _eq(self, other): - return Polyhedron(equalities=[(self - other)._toint()]) - - @_polymorphic_method - def __le__(self, other): - return Polyhedron(inequalities=[(other - self)._toint()]) - - @_polymorphic_method - def __lt__(self, other): - return Polyhedron(inequalities=[(other - self)._toint() - 1]) - - @_polymorphic_method - def __ge__(self, other): - return Polyhedron(inequalities=[(self - other)._toint()]) - - @_polymorphic_method - def __gt__(self, other): - return Polyhedron(inequalities=[(self - other)._toint() - 1]) - - @classmethod - def fromsympy(cls, expr): - import sympy - coefficients = {} - constant = 0 - for symbol, coefficient in expr.as_coefficients_dict().items(): - coefficient = Fraction(coefficient.p, coefficient.q) - if symbol == sympy.S.One: - constant = coefficient - elif isinstance(symbol, sympy.Symbol): - symbol = symbol.name - coefficients[symbol] = coefficient - else: - raise ValueError('non-linear expression: {!r}'.format(expr)) - return cls(coefficients, constant) - - def tosympy(self): - import sympy - expr = 0 - for symbol, coefficient in self.coefficients(): - term = coefficient * sympy.Symbol(symbol) - expr += term - expr += self.constant - return expr - - -class Constant(Expression): - - def __new__(cls, numerator=0, denominator=None): - self = object().__new__(cls) - if denominator is None: - if isinstance(numerator, numbers.Rational): - self._constant = numerator - elif isinstance(numerator, Constant): - self._constant = numerator.constant - else: - raise TypeError('constant must be a rational number or a Constant instance') - else: - self._constant = Fraction(numerator, denominator) - self._coefficients = {} - self._symbols = () - self._dimension = 0 - return self - - def isconstant(self): - return True - - def __bool__(self): - return bool(self.constant) - - def __repr__(self): - if self.constant.denominator == 1: - return '{}({!r})'.format(self.__class__.__name__, self.constant) - else: - return '{}({!r}, {!r})'.format(self.__class__.__name__, - self.constant.numerator, self.constant.denominator) - - @classmethod - def fromsympy(cls, expr): - import sympy - if isinstance(expr, sympy.Rational): - return cls(expr.p, expr.q) - elif isinstance(expr, numbers.Rational): - return cls(expr) - else: - raise TypeError('expr must be a sympy.Rational instance') - - -class Symbol(Expression): - - __slots__ = Expression.__slots__ + ( - '_name', - ) - - def __new__(cls, name): - if isinstance(name, Symbol): - name = name.name - elif not isinstance(name, str): - raise TypeError('name must be a string or a Symbol instance') - self = object().__new__(cls) - self._coefficients = {name: 1} - self._constant = 0 - self._symbols = tuple(name) - self._name = name - self._dimension = 1 - return self - - @property - def name(self): - return self._name - - def issymbol(self): - return True - - def __repr__(self): - return '{}({!r})'.format(self.__class__.__name__, self._name) - - @classmethod - def fromsympy(cls, expr): - import sympy - if isinstance(expr, sympy.Symbol): - return cls(expr.name) - else: - raise TypeError('expr must be a sympy.Symbol instance') - - -def symbols(names): - if isinstance(names, str): - names = names.replace(',', ' ').split() - return (Symbol(name) for name in names) - - -@_polymorphic_operator -def eq(a, b): - return a.__eq__(b) - -@_polymorphic_operator -def le(a, b): - return a.__le__(b) - -@_polymorphic_operator -def lt(a, b): - return a.__lt__(b) - -@_polymorphic_operator -def ge(a, b): - return a.__ge__(b) - -@_polymorphic_operator -def gt(a, b): - return a.__gt__(b) - - -class Polyhedron: - """ - This class implements polyhedrons. - """ - - __slots__ = ( - '_equalities', - '_inequalities', - '_constraints', - '_symbols', - ) - - def __new__(cls, equalities=None, inequalities=None): - if isinstance(equalities, str): - if inequalities is not None: - raise TypeError('too many arguments') - return cls.fromstring(equalities) - self = super().__new__(cls) - self._equalities = [] - if equalities is not None: - for constraint in equalities: - for value in constraint.values(): - if value.denominator != 1: - raise TypeError('non-integer constraint: ' - '{} == 0'.format(constraint)) - self._equalities.append(constraint) - self._equalities = tuple(self._equalities) - self._inequalities = [] - if inequalities is not None: - for constraint in inequalities: - for value in constraint.values(): - if value.denominator != 1: - raise TypeError('non-integer constraint: ' - '{} <= 0'.format(constraint)) - self._inequalities.append(constraint) - self._inequalities = tuple(self._inequalities) - self._constraints = self._equalities + self._inequalities - self._symbols = set() - for constraint in self._constraints: - self.symbols.update(constraint.symbols) - self._symbols = tuple(sorted(self._symbols)) - return self - - @classmethod - def _fromast(cls, node): - if isinstance(node, ast.Module) and len(node.body) == 1: - return cls._fromast(node.body[0]) - elif isinstance(node, ast.Expr): - return cls._fromast(node.value) - elif isinstance(node, ast.BinOp) and isinstance(node.op, ast.BitAnd): - equalities1, inequalities1 = cls._fromast(node.left) - equalities2, inequalities2 = cls._fromast(node.right) - equalities = equalities1 + equalities2 - inequalities = inequalities1 + inequalities2 - return equalities, inequalities - elif isinstance(node, ast.Compare): - equalities = [] - inequalities = [] - left = Expression._fromast(node.left) - for i in range(len(node.ops)): - op = node.ops[i] - right = Expression._fromast(node.comparators[i]) - if isinstance(op, ast.Lt): - inequalities.append(right - left - 1) - elif isinstance(op, ast.LtE): - inequalities.append(right - left) - elif isinstance(op, ast.Eq): - equalities.append(left - right) - elif isinstance(op, ast.GtE): - inequalities.append(left - right) - elif isinstance(op, ast.Gt): - inequalities.append(left - right - 1) - else: - break - left = right - else: - return equalities, inequalities - raise SyntaxError('invalid syntax') - - @classmethod - def fromstring(cls, string): - string = string.strip() - string = re.sub(r'^\{\s*|\s*\}$', '', string) - string = re.sub(r'([^<=>])=([^<=>])', r'\1==\2', string) - string = re.sub(r'(\d+|\))\s*([^\W\d_]\w*|\()', r'\1*\2', string) - tokens = re.split(r',|;|and|&&|/\\|∧', string, flags=re.I) - tokens = ['({})'.format(token) for token in tokens] - string = ' & '.join(tokens) - tree = ast.parse(string, 'eval') - equalities, inequalities = cls._fromast(tree) - return cls(equalities, inequalities) - - @property - def equalities(self): - return self._equalities - - @property - def inequalities(self): - return self._inequalities - - @property - def constraints(self): - return self._constraints - - @property - def symbols(self): - return self._symbols - - @property - def dimension(self): - return len(self.symbols) - - def __bool__(self): - return not self.is_empty() - - def __contains__(self, value): - # is the value in the polyhedron? - raise NotImplementedError - - def __eq__(self, other): - # works correctly when symbols is not passed - # should be equal if values are the same even if symbols are different - bset = self._toisl() - other = other._toisl() - return bool(libisl.isl_basic_set_plain_is_equal(bset, other)) - - def isempty(self): - bset = self._toisl() - return bool(libisl.isl_basic_set_is_empty(bset)) - - def isuniverse(self): - bset = self._toisl() - return bool(libisl.isl_basic_set_is_universe(bset)) - - def isdisjoint(self, other): - # return true if the polyhedron has no elements in common with other - #symbols = self._symbolunion(other) - bset = self._toisl() - other = other._toisl() - return bool(libisl.isl_set_is_disjoint(bset, other)) - - def issubset(self, other): - # check if self(bset) is a subset of other - symbols = self._symbolunion(other) - bset = self._toisl(symbols) - other = other._toisl(symbols) - return bool(libisl.isl_set_is_strict_subset(other, bset)) - - def __le__(self, other): - return self.issubset(other) - - def __lt__(self, other): - symbols = self._symbolunion(other) - bset = self._toisl(symbols) - other = other._toisl(symbols) - return bool(libisl.isl_set_is_strict_subset(other, bset)) - - def issuperset(self, other): - # test whether every element in other is in the polyhedron - raise NotImplementedError - - def __ge__(self, other): - return self.issuperset(other) - - def __gt__(self, other): - symbols = self._symbolunion(other) - bset = self._toisl(symbols) - other = other._toisl(symbols) - bool(libisl.isl_set_is_strict_subset(other, bset)) - raise NotImplementedError - - def union(self, *others): - # return a new polyhedron with elements from the polyhedron and all - # others (convex union) - raise NotImplementedError - - def __or__(self, other): - return self.union(other) - - def intersection(self, *others): - # return a new polyhedron with elements common to the polyhedron and all - # others - # a poor man's implementation could be: - # equalities = list(self.equalities) - # inequalities = list(self.inequalities) - # for other in others: - # equalities.extend(other.equalities) - # inequalities.extend(other.inequalities) - # return self.__class__(equalities, inequalities) - raise NotImplementedError - - def __and__(self, other): - return self.intersection(other) - - def difference(self, other): - # return a new polyhedron with elements in the polyhedron that are not in the other - symbols = self._symbolunion(other) - bset = self._toisl(symbols) - other = other._toisl(symbols) - difference = libisl.isl_set_subtract(bset, other) - return difference - - def __sub__(self, other): - return self.difference(other) - - def __str__(self): - constraints = [] - for constraint in self.equalities: - constraints.append('{} == 0'.format(constraint)) - for constraint in self.inequalities: - constraints.append('{} >= 0'.format(constraint)) - return '{}'.format(', '.join(constraints)) - - def __repr__(self): - if self.isempty(): - return 'Empty' - elif self.isuniverse(): - return 'Universe' - else: - return '{}({!r})'.format(self.__class__.__name__, str(self)) - - @classmethod - def _fromsympy(cls, expr): - import sympy - equalities = [] - inequalities = [] - if expr.func == sympy.And: - for arg in expr.args: - arg_eqs, arg_ins = cls._fromsympy(arg) - equalities.extend(arg_eqs) - inequalities.extend(arg_ins) - elif expr.func == sympy.Eq: - expr = Expression.fromsympy(expr.args[0] - expr.args[1]) - equalities.append(expr) - else: - if expr.func == sympy.Lt: - expr = Expression.fromsympy(expr.args[1] - expr.args[0] - 1) - elif expr.func == sympy.Le: - expr = Expression.fromsympy(expr.args[1] - expr.args[0]) - elif expr.func == sympy.Ge: - expr = Expression.fromsympy(expr.args[0] - expr.args[1]) - elif expr.func == sympy.Gt: - expr = Expression.fromsympy(expr.args[0] - expr.args[1] - 1) - else: - raise ValueError('non-polyhedral expression: {!r}'.format(expr)) - inequalities.append(expr) - return equalities, inequalities - - @classmethod - def fromsympy(cls, expr): - import sympy - equalities, inequalities = cls._fromsympy(expr) - return cls(equalities, inequalities) - - def tosympy(self): - import sympy - constraints = [] - for equality in self.equalities: - constraints.append(sympy.Eq(equality.tosympy(), 0)) - for inequality in self.inequalities: - constraints.append(sympy.Ge(inequality.tosympy(), 0)) - return sympy.And(*constraints) - - def _symbolunion(self, *others): - symbols = set(self.symbols) - for other in others: - symbols.update(other.symbols) - return sorted(symbols) - - def _toisl(self, symbols=None): - if symbols is None: - symbols = self.symbols - dimension = len(symbols) - space = libisl.isl_space_set_alloc(_main_ctx, 0, dimension) - bset = libisl.isl_basic_set_universe(libisl.isl_space_copy(space)) - ls = libisl.isl_local_space_from_space(space) - for equality in self.equalities: - ceq = libisl.isl_equality_alloc(libisl.isl_local_space_copy(ls)) - for symbol, coefficient in equality.coefficients(): - val = str(coefficient).encode() - val = libisl.isl_val_read_from_str(_main_ctx, val) - dim = symbols.index(symbol) - ceq = libisl.isl_constraint_set_coefficient_val(ceq, libisl.isl_dim_set, dim, val) - if equality.constant != 0: - val = str(equality.constant).encode() - val = libisl.isl_val_read_from_str(_main_ctx, val) - ceq = libisl.isl_constraint_set_constant_val(ceq, val) - bset = libisl.isl_basic_set_add_constraint(bset, ceq) - for inequality in self.inequalities: - cin = libisl.isl_inequality_alloc(libisl.isl_local_space_copy(ls)) - for symbol, coefficient in inequality.coefficients(): - val = str(coefficient).encode() - val = libisl.isl_val_read_from_str(_main_ctx, val) - dim = symbols.index(symbol) - cin = libisl.isl_constraint_set_coefficient_val(cin, libisl.isl_dim_set, dim, val) - if inequality.constant != 0: - val = str(inequality.constant).encode() - val = libisl.isl_val_read_from_str(_main_ctx, val) - cin = libisl.isl_constraint_set_constant_val(cin, val) - bset = libisl.isl_basic_set_add_constraint(bset, cin) - bset = isl.BasicSet(bset) - return bset - - @classmethod - def _fromisl(cls, bset, symbols): - raise NotImplementedError - equalities = ... - inequalities = ... - return cls(equalities, inequalities) - '''takes basic set in isl form and puts back into python version of polyhedron - isl example code gives isl form as: - "{[i] : exists (a : i = 2a and i >= 10 and i <= 42)}") - our printer is giving form as: - { [i0, i1] : 2i1 >= -2 - i0 } ''' - -Empty = eq(0,1) - -Universe = Polyhedron() - - -if __name__ == '__main__': - #p = Polyhedron('2a + 2b + 1 == 0') # empty - p = Polyhedron('3x + 2y + 3 == 0, y == 0') # not empty - ip = p._toisl() - print(ip) - print(ip.constraints())