pypol/polyhedra.py

   1 import functools
   2 import math
   3 import numbers
   4
   5 from . import islhelper
   6
   7 from .islhelper import mainctx, libisl
   8 from .geometry import GeometricObject, Point
   9 from .linexprs import Expression, Rational
  10 from .domains import Domain
  11
  12
  13 __all__ = [
  14     'Polyhedron',
  15     'Lt', 'Le', 'Eq', 'Ne', 'Ge', 'Gt',
  16     'Empty', 'Universe',
  17 ]
  18
  19
  20 class Polyhedron(Domain):
  21
  22     __slots__ = (
  23         '_equalities',
  24         '_inequalities',
  25         '_constraints',
  26         '_symbols',
  27         '_dimension',
  28     )
  29
  30     def __new__(cls, equalities=None, inequalities=None):
  31         if isinstance(equalities, str):
  32             if inequalities is not None:
  33                 raise TypeError('too many arguments')
  34             return cls.fromstring(equalities)
  35         elif isinstance(equalities, GeometricObject):
  36             if inequalities is not None:
  37                 raise TypeError('too many arguments')
  38             return equalities.aspolyhedron()
  39         if equalities is None:
  40             equalities = []
  41         else:
  42             for i, equality in enumerate(equalities):
  43                 if not isinstance(equality, Expression):
  44                     raise TypeError('equalities must be linear expressions')
  45                 equalities[i] = equality.scaleint()
  46         if inequalities is None:
  47             inequalities = []
  48         else:
  49             for i, inequality in enumerate(inequalities):
  50                 if not isinstance(inequality, Expression):
  51                     raise TypeError('inequalities must be linear expressions')
  52                 inequalities[i] = inequality.scaleint()
  53         symbols = cls._xsymbols(equalities + inequalities)
  54         islbset = cls._toislbasicset(equalities, inequalities, symbols)
  55         return cls._fromislbasicset(islbset, symbols)
  56
  57     @property
  58     def equalities(self):
  59         """
  60         Return a list of the equalities in a set.
  61         """
  62         return self._equalities
  63
  64     @property
  65     def inequalities(self):
  66         """
  67         Return a list of the inequalities in a set.
  68         """
  69         return self._inequalities
  70
  71     @property
  72     def constraints(self):
  73         """
  74         Return ta list of the constraints of a set.
  75         """
  76         return self._constraints
  77
  78     @property
  79     def polyhedra(self):
  80         return self,
  81
  82     def disjoint(self):
  83         """
  84         Return a set as disjoint.
  85         """
  86         return self
  87
  88     def isuniverse(self):
  89         """
  90         Return true if a set is the Universe set.
  91         """
  92         islbset = self._toislbasicset(self.equalities, self.inequalities,
  93             self.symbols)
  94         universe = bool(libisl.isl_basic_set_is_universe(islbset))
  95         libisl.isl_basic_set_free(islbset)
  96         return universe
  97
  98     def aspolyhedron(self):
  99         """
 100         Return polyhedral hull of a set.
 101         """
 102         return self
 103
 104     def __contains__(self, point):
 105         if not isinstance(point, Point):
 106             raise TypeError('point must be a Point instance')
 107         if self.symbols != point.symbols:
 108             raise ValueError('arguments must belong to the same space')
 109         for equality in self.equalities:
 110             if equality.subs(point.coordinates()) != 0:
 111                 return False
 112         for inequality in self.inequalities:
 113             if inequality.subs(point.coordinates()) < 0:
 114                 return False
 115         return True
 116
 117     def subs(self, symbol, expression=None):
 118         """
 119         Subsitute the given value into an expression and return the resulting
 120         expression.
 121         """
 122         equalities = [equality.subs(symbol, expression)
 123             for equality in self.equalities]
 124         inequalities = [inequality.subs(symbol, expression)
 125             for inequality in self.inequalities]
 126         return Polyhedron(equalities, inequalities)
 127
 128     def _asinequalities(self):
 129         inequalities = list(self.equalities)
 130         inequalities.extend([-expression for expression in self.equalities])
 131         inequalities.extend(self.inequalities)
 132         return inequalities
 133
 134     def widen(self, other):
 135         if not isinstance(other, Polyhedron):
 136             raise ValueError('argument must be a Polyhedron instance')
 137         inequalities1 = self._asinequalities()
 138         inequalities2 = other._asinequalities()
 139         inequalities = []
 140         for inequality1 in inequalities1:
 141             if other <= Polyhedron(inequalities=[inequality1]):
 142                 inequalities.append(inequality1)
 143         for inequality2 in inequalities2:
 144             for i in range(len(inequalities1)):
 145                 inequalities3 = inequalities1[:i] + inequalities[i + 1:]
 146                 inequalities3.append(inequality2)
 147                 polyhedron3 = Polyhedron(inequalities=inequalities3)
 148                 if self == polyhedron3:
 149                     inequalities.append(inequality2)
 150                     break
 151         return Polyhedron(inequalities=inequalities)
 152
 153     @classmethod
 154     def _fromislbasicset(cls, islbset, symbols):
 155         islconstraints = islhelper.isl_basic_set_constraints(islbset)
 156         equalities = []
 157         inequalities = []
 158         for islconstraint in islconstraints:
 159             constant = libisl.isl_constraint_get_constant_val(islconstraint)
 160             constant = islhelper.isl_val_to_int(constant)
 161             coefficients = {}
 162             for index, symbol in enumerate(symbols):
 163                 coefficient = libisl.isl_constraint_get_coefficient_val(islconstraint,
 164                     libisl.isl_dim_set, index)
 165                 coefficient = islhelper.isl_val_to_int(coefficient)
 166                 if coefficient != 0:
 167                     coefficients[symbol] = coefficient
 168             expression = Expression(coefficients, constant)
 169             if libisl.isl_constraint_is_equality(islconstraint):
 170                 equalities.append(expression)
 171             else:
 172                 inequalities.append(expression)
 173         libisl.isl_basic_set_free(islbset)
 174         self = object().__new__(Polyhedron)
 175         self._equalities = tuple(equalities)
 176         self._inequalities = tuple(inequalities)
 177         self._constraints = tuple(equalities + inequalities)
 178         self._symbols = cls._xsymbols(self._constraints)
 179         self._dimension = len(self._symbols)
 180         return self
 181
 182     @classmethod
 183     def _toislbasicset(cls, equalities, inequalities, symbols):
 184         dimension = len(symbols)
 185         indices = {symbol: index for index, symbol in enumerate(symbols)}
 186         islsp = libisl.isl_space_set_alloc(mainctx, 0, dimension)
 187         islbset = libisl.isl_basic_set_universe(libisl.isl_space_copy(islsp))
 188         islls = libisl.isl_local_space_from_space(islsp)
 189         for equality in equalities:
 190             isleq = libisl.isl_equality_alloc(libisl.isl_local_space_copy(islls))
 191             for symbol, coefficient in equality.coefficients():
 192                 islval = str(coefficient).encode()
 193                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 194                 index = indices[symbol]
 195                 isleq = libisl.isl_constraint_set_coefficient_val(isleq,
 196                     libisl.isl_dim_set, index, islval)
 197             if equality.constant != 0:
 198                 islval = str(equality.constant).encode()
 199                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 200                 isleq = libisl.isl_constraint_set_constant_val(isleq, islval)
 201             islbset = libisl.isl_basic_set_add_constraint(islbset, isleq)
 202         for inequality in inequalities:
 203             islin = libisl.isl_inequality_alloc(libisl.isl_local_space_copy(islls))
 204             for symbol, coefficient in inequality.coefficients():
 205                 islval = str(coefficient).encode()
 206                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 207                 index = indices[symbol]
 208                 islin = libisl.isl_constraint_set_coefficient_val(islin,
 209                     libisl.isl_dim_set, index, islval)
 210             if inequality.constant != 0:
 211                 islval = str(inequality.constant).encode()
 212                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 213                 islin = libisl.isl_constraint_set_constant_val(islin, islval)
 214             islbset = libisl.isl_basic_set_add_constraint(islbset, islin)
 215         return islbset
 216
 217     @classmethod
 218     def fromstring(cls, string):
 219         domain = Domain.fromstring(string)
 220         if not isinstance(domain, Polyhedron):
 221             raise ValueError('non-polyhedral expression: {!r}'.format(string))
 222         return domain
 223
 224     def __repr__(self):
 225         strings = []
 226         for equality in self.equalities:
 227             strings.append('Eq({}, 0)'.format(equality))
 228         for inequality in self.inequalities:
 229             strings.append('Ge({}, 0)'.format(inequality))
 230         if len(strings) == 1:
 231             return strings[0]
 232         else:
 233             return 'And({})'.format(', '.join(strings))
 234
 235
 236     def _repr_latex_(self):
 237         strings = []
 238         for equality in self.equalities:
 239             strings.append('{} = 0'.format(equality._repr_latex_().strip('$')))
 240         for inequality in self.inequalities:
 241             strings.append('{} \\ge 0'.format(inequality._repr_latex_().strip('$')))
 242         return '$${}$$'.format(' \\wedge '.join(strings))
 243
 244     @classmethod
 245     def fromsympy(cls, expr):
 246         """
 247         Convert a sympy object to an expression.
 248         """
 249         domain = Domain.fromsympy(expr)
 250         if not isinstance(domain, Polyhedron):
 251             raise ValueError('non-polyhedral expression: {!r}'.format(expr))
 252         return domain
 253
 254     def tosympy(self):
 255         """
 256         Return an expression as a sympy object.
 257         """
 258         import sympy
 259         constraints = []
 260         for equality in self.equalities:
 261             constraints.append(sympy.Eq(equality.tosympy(), 0))
 262         for inequality in self.inequalities:
 263             constraints.append(sympy.Ge(inequality.tosympy(), 0))
 264         return sympy.And(*constraints)
 265
 266 class EmptyType(Polyhedron):
 267
 268     __slots__ = Polyhedron.__slots__
 269
 270     def __new__(cls):
 271         self = object().__new__(cls)
 272         self._equalities = (Rational(1),)
 273         self._inequalities = ()
 274         self._constraints = self._equalities
 275         self._symbols = ()
 276         self._dimension = 0
 277         return self
 278
 279     def widen(self, other):
 280         if not isinstance(other, Polyhedron):
 281             raise ValueError('argument must be a Polyhedron instance')
 282         return other
 283
 284     def __repr__(self):
 285         return 'Empty'
 286
 287     def _repr_latex_(self):
 288         return '$$\\emptyset$$'
 289
 290 Empty = EmptyType()
 291
 292
 293 class UniverseType(Polyhedron):
 294
 295     __slots__ = Polyhedron.__slots__
 296
 297     def __new__(cls):
 298         self = object().__new__(cls)
 299         self._equalities = ()
 300         self._inequalities = ()
 301         self._constraints = ()
 302         self._symbols = ()
 303         self._dimension = ()
 304         return self
 305
 306     def __repr__(self):
 307         return 'Universe'
 308
 309     def _repr_latex_(self):
 310         return '$$\\Omega$$'
 311
 312 Universe = UniverseType()
 313
 314
 315 def _polymorphic(func):
 316     @functools.wraps(func)
 317     def wrapper(left, right):
 318         if not isinstance(left, Expression):
 319             if isinstance(left, numbers.Rational):
 320                 left = Rational(left)
 321             else:
 322                 raise TypeError('left must be a a rational number '
 323                     'or a linear expression')
 324         if not isinstance(right, Expression):
 325             if isinstance(right, numbers.Rational):
 326                 right = Rational(right)
 327             else:
 328                 raise TypeError('right must be a a rational number '
 329                     'or a linear expression')
 330         return func(left, right)
 331     return wrapper
 332
 333 @_polymorphic
 334 def Lt(left, right):
 335     """
 336     Assert first set is less than the second set.
 337     """
 338     return Polyhedron([], [right - left - 1])
 339
 340 @_polymorphic
 341 def Le(left, right):
 342     """
 343     Assert first set is less than or equal to the second set.
 344     """
 345     return Polyhedron([], [right - left])
 346
 347 @_polymorphic
 348 def Eq(left, right):
 349     """
 350     Assert first set is equal to the second set.
 351     """
 352     return Polyhedron([left - right], [])
 353
 354 @_polymorphic
 355 def Ne(left, right):
 356     """
 357     Assert first set is not equal to the second set.
 358     """
 359     return ~Eq(left, right)
 360
 361 @_polymorphic
 362 def Gt(left, right):
 363     """
 364     Assert first set is greater than the second set.
 365     """
 366     return Polyhedron([], [left - right - 1])
 367
 368 @_polymorphic
 369 def Ge(left, right):
 370     """
 371     Assert first set is greater than or equal to the second set.
 372     """
 373     return Polyhedron([], [left - right])
 374