X-Git-Url: https://scm.cri.ensmp.fr/git/linpy.git/blobdiff_plain/82e100a9d5e7db532fda649849dc784148e55069..25ce908cffca380f930182a77c1e5a4491042a1c:/pypol/polyhedra.py~ diff --git a/pypol/polyhedra.py~ b/pypol/polyhedra.py~ new file mode 100644 index 0000000..ccb1a8c --- /dev/null +++ b/pypol/polyhedra.py~ @@ -0,0 +1,358 @@ +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])