import functools import math import numbers from . import islhelper from .islhelper import mainctx, libisl from .geometry import GeometricObject, Point from .linexprs import Expression, Rational from .domains import Domain __all__ = [ 'Polyhedron', 'Lt', 'Le', 'Eq', 'Ne', 'Ge', 'Gt', 'Empty', 'Universe', ] class Polyhedron(Domain): __slots__ = ( '_equalities', '_inequalities', '_constraints', '_symbols', '_dimension', ) 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) elif isinstance(equalities, GeometricObject): if inequalities is not None: raise TypeError('too many arguments') return equalities.aspolyhedron() if equalities is None: equalities = [] else: for i, equality in enumerate(equalities): if not isinstance(equality, Expression): raise TypeError('equalities must be linear expressions') equalities[i] = equality.scaleint() if inequalities is None: inequalities = [] else: for i, inequality in enumerate(inequalities): if not isinstance(inequality, Expression): raise TypeError('inequalities must be linear expressions') inequalities[i] = inequality.scaleint() symbols = cls._xsymbols(equalities + inequalities) islbset = cls._toislbasicset(equalities, inequalities, symbols) return cls._fromislbasicset(islbset, symbols) @property def equalities(self): return self._equalities @property def inequalities(self): return self._inequalities @property def constraints(self): return self._constraints @property def polyhedra(self): return self, def disjoint(self): """ Return this set as disjoint. """ return self def isuniverse(self): """ Return true if this set is the Universe set. """ islbset = self._toislbasicset(self.equalities, self.inequalities, self.symbols) universe = bool(libisl.isl_basic_set_is_universe(islbset)) libisl.isl_basic_set_free(islbset) return universe def aspolyhedron(self): """ Return polyhedral hull of this set. """ return self def __contains__(self, point): if not isinstance(point, Point): raise TypeError('point must be a Point instance') if self.symbols != point.symbols: raise ValueError('arguments must belong to the same space') for equality in self.equalities: if equality.subs(point.coordinates()) != 0: return False for inequality in self.inequalities: if inequality.subs(point.coordinates()) < 0: return False return True def subs(self, symbol, expression=None): equalities = [equality.subs(symbol, expression) for equality in self.equalities] inequalities = [inequality.subs(symbol, expression) for inequality in self.inequalities] return Polyhedron(equalities, inequalities) def _asinequalities(self): inequalities = list(self.equalities) inequalities.extend([-expression for expression in self.equalities]) inequalities.extend(self.inequalities) return inequalities def widen(self, other): if not isinstance(other, Polyhedron): raise ValueError('argument must be a Polyhedron instance') inequalities1 = self._asinequalities() inequalities2 = other._asinequalities() inequalities = [] for inequality1 in inequalities1: if other <= Polyhedron(inequalities=[inequality1]): inequalities.append(inequality1) for inequality2 in inequalities2: for i in range(len(inequalities1)): inequalities3 = inequalities1[:i] + inequalities[i + 1:] inequalities3.append(inequality2) polyhedron3 = Polyhedron(inequalities=inequalities3) if self == polyhedron3: inequalities.append(inequality2) break return Polyhedron(inequalities=inequalities) @classmethod def _fromislbasicset(cls, islbset, symbols): if libisl.isl_basic_set_is_empty(islbset): return Empty if libisl.isl_basic_set_is_universe(islbset): return Universe islconstraints = islhelper.isl_basic_set_constraints(islbset) equalities = [] inequalities = [] for islconstraint in islconstraints: constant = libisl.isl_constraint_get_constant_val(islconstraint) constant = islhelper.isl_val_to_int(constant) coefficients = {} for index, symbol in enumerate(symbols): coefficient = libisl.isl_constraint_get_coefficient_val(islconstraint, libisl.isl_dim_set, index) coefficient = islhelper.isl_val_to_int(coefficient) if coefficient != 0: coefficients[symbol] = coefficient expression = Expression(coefficients, constant) if libisl.isl_constraint_is_equality(islconstraint): equalities.append(expression) else: inequalities.append(expression) libisl.isl_basic_set_free(islbset) self = object().__new__(Polyhedron) self._equalities = tuple(equalities) self._inequalities = tuple(inequalities) self._constraints = tuple(equalities + inequalities) self._symbols = cls._xsymbols(self._constraints) self._dimension = len(self._symbols) return self @classmethod def _toislbasicset(cls, equalities, inequalities, symbols): dimension = len(symbols) indices = {symbol: index for index, symbol in enumerate(symbols)} islsp = libisl.isl_space_set_alloc(mainctx, 0, dimension) islbset = libisl.isl_basic_set_universe(libisl.isl_space_copy(islsp)) islls = libisl.isl_local_space_from_space(islsp) for equality in equalities: isleq = libisl.isl_equality_alloc(libisl.isl_local_space_copy(islls)) for symbol, coefficient in equality.coefficients(): islval = str(coefficient).encode() islval = libisl.isl_val_read_from_str(mainctx, islval) index = indices[symbol] isleq = libisl.isl_constraint_set_coefficient_val(isleq, libisl.isl_dim_set, index, islval) if equality.constant != 0: islval = str(equality.constant).encode() islval = libisl.isl_val_read_from_str(mainctx, islval) isleq = libisl.isl_constraint_set_constant_val(isleq, islval) islbset = libisl.isl_basic_set_add_constraint(islbset, isleq) for inequality in inequalities: islin = libisl.isl_inequality_alloc(libisl.isl_local_space_copy(islls)) for symbol, coefficient in inequality.coefficients(): islval = str(coefficient).encode() islval = libisl.isl_val_read_from_str(mainctx, islval) index = indices[symbol] islin = libisl.isl_constraint_set_coefficient_val(islin, libisl.isl_dim_set, index, islval) if inequality.constant != 0: islval = str(inequality.constant).encode() islval = libisl.isl_val_read_from_str(mainctx, islval) islin = libisl.isl_constraint_set_constant_val(islin, islval) islbset = libisl.isl_basic_set_add_constraint(islbset, islin) return islbset @classmethod def fromstring(cls, string): domain = Domain.fromstring(string) if not isinstance(domain, Polyhedron): raise ValueError('non-polyhedral expression: {!r}'.format(string)) return domain def __repr__(self): strings = [] for equality in self.equalities: strings.append('Eq({}, 0)'.format(equality)) for inequality in self.inequalities: strings.append('Ge({}, 0)'.format(inequality)) if len(strings) == 1: return strings[0] else: return 'And({})'.format(', '.join(strings)) def _repr_latex_(self): strings = [] for equality in self.equalities: strings.append('{} = 0'.format(equality._repr_latex_().strip('$'))) for inequality in self.inequalities: strings.append('{} \\ge 0'.format(inequality._repr_latex_().strip('$'))) return '$${}$$'.format(' \\wedge '.join(strings)) @classmethod def fromsympy(cls, expr): domain = Domain.fromsympy(expr) if not isinstance(domain, Polyhedron): raise ValueError('non-polyhedral expression: {!r}'.format(expr)) return domain 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) class EmptyType(Polyhedron): __slots__ = Polyhedron.__slots__ def __new__(cls): self = object().__new__(cls) self._equalities = (Rational(1),) self._inequalities = () self._constraints = self._equalities self._symbols = () self._dimension = 0 return self def widen(self, other): if not isinstance(other, Polyhedron): raise ValueError('argument must be a Polyhedron instance') return other def __repr__(self): return 'Empty' def _repr_latex_(self): return '$$\\emptyset$$' Empty = EmptyType() class UniverseType(Polyhedron): __slots__ = Polyhedron.__slots__ def __new__(cls): self = object().__new__(cls) self._equalities = () self._inequalities = () self._constraints = () self._symbols = () self._dimension = () return self def __repr__(self): return 'Universe' def _repr_latex_(self): return '$$\\Omega$$' Universe = UniverseType() def _polymorphic(func): @functools.wraps(func) def wrapper(left, right): if not isinstance(left, Expression): if isinstance(left, numbers.Rational): left = Rational(left) else: raise TypeError('left must be a a rational number ' 'or a linear expression') if not isinstance(right, Expression): if isinstance(right, numbers.Rational): right = Rational(right) else: raise TypeError('right must be a a rational number ' 'or a linear expression') return func(left, right) return wrapper @_polymorphic def Lt(left, right): """ Return true if the first set is less than the second. """ return Polyhedron([], [right - left - 1]) @_polymorphic def Le(left, right): """ Return true the first set is less than or equal to the second. """ return Polyhedron([], [right - left]) @_polymorphic def Eq(left, right): """ Return true if the sets are equal. """ return Polyhedron([left - right], []) @_polymorphic def Ne(left, right): """ Return true if the sets are NOT equal. """ return ~Eq(left, right) @_polymorphic def Gt(left, right): """ Return true if the first set is greater than the second set. """ return Polyhedron([], [left - right - 1]) @_polymorphic def Ge(left, right): """ Return true if the first set is greater than or equal the second set. """ return Polyhedron([], [left - right])