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 expression.
 120     """
 121         equalities = [equality.subs(symbol, expression)
 122             for equality in self.equalities]
 123         inequalities = [inequality.subs(symbol, expression)
 124             for inequality in self.inequalities]
 125         return Polyhedron(equalities, inequalities)
 126
 127     def _asinequalities(self):
 128         inequalities = list(self.equalities)
 129         inequalities.extend([-expression for expression in self.equalities])
 130         inequalities.extend(self.inequalities)
 131         return inequalities
 132
 133     def widen(self, other):
 134         if not isinstance(other, Polyhedron):
 135             raise ValueError('argument must be a Polyhedron instance')
 136         inequalities1 = self._asinequalities()
 137         inequalities2 = other._asinequalities()
 138         inequalities = []
 139         for inequality1 in inequalities1:
 140             if other <= Polyhedron(inequalities=[inequality1]):
 141                 inequalities.append(inequality1)
 142         for inequality2 in inequalities2:
 143             for i in range(len(inequalities1)):
 144                 inequalities3 = inequalities1[:i] + inequalities[i + 1:]
 145                 inequalities3.append(inequality2)
 146                 polyhedron3 = Polyhedron(inequalities=inequalities3)
 147                 if self == polyhedron3:
 148                     inequalities.append(inequality2)
 149                     break
 150         return Polyhedron(inequalities=inequalities)
 151
 152     @classmethod
 153     def _fromislbasicset(cls, islbset, symbols):
 154         islconstraints = islhelper.isl_basic_set_constraints(islbset)
 155         equalities = []
 156         inequalities = []
 157         for islconstraint in islconstraints:
 158             constant = libisl.isl_constraint_get_constant_val(islconstraint)
 159             constant = islhelper.isl_val_to_int(constant)
 160             coefficients = {}
 161             for index, symbol in enumerate(symbols):
 162                 coefficient = libisl.isl_constraint_get_coefficient_val(islconstraint,
 163                     libisl.isl_dim_set, index)
 164                 coefficient = islhelper.isl_val_to_int(coefficient)
 165                 if coefficient != 0:
 166                     coefficients[symbol] = coefficient
 167             expression = Expression(coefficients, constant)
 168             if libisl.isl_constraint_is_equality(islconstraint):
 169                 equalities.append(expression)
 170             else:
 171                 inequalities.append(expression)
 172         libisl.isl_basic_set_free(islbset)
 173         self = object().__new__(Polyhedron)
 174         self._equalities = tuple(equalities)
 175         self._inequalities = tuple(inequalities)
 176         self._constraints = tuple(equalities + inequalities)
 177         self._symbols = cls._xsymbols(self._constraints)
 178         self._dimension = len(self._symbols)
 179         return self
 180
 181     @classmethod
 182     def _toislbasicset(cls, equalities, inequalities, symbols):
 183         dimension = len(symbols)
 184         indices = {symbol: index for index, symbol in enumerate(symbols)}
 185         islsp = libisl.isl_space_set_alloc(mainctx, 0, dimension)
 186         islbset = libisl.isl_basic_set_universe(libisl.isl_space_copy(islsp))
 187         islls = libisl.isl_local_space_from_space(islsp)
 188         for equality in equalities:
 189             isleq = libisl.isl_equality_alloc(libisl.isl_local_space_copy(islls))
 190             for symbol, coefficient in equality.coefficients():
 191                 islval = str(coefficient).encode()
 192                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 193                 index = indices[symbol]
 194                 isleq = libisl.isl_constraint_set_coefficient_val(isleq,
 195                     libisl.isl_dim_set, index, islval)
 196             if equality.constant != 0:
 197                 islval = str(equality.constant).encode()
 198                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 199                 isleq = libisl.isl_constraint_set_constant_val(isleq, islval)
 200             islbset = libisl.isl_basic_set_add_constraint(islbset, isleq)
 201         for inequality in inequalities:
 202             islin = libisl.isl_inequality_alloc(libisl.isl_local_space_copy(islls))
 203             for symbol, coefficient in inequality.coefficients():
 204                 islval = str(coefficient).encode()
 205                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 206                 index = indices[symbol]
 207                 islin = libisl.isl_constraint_set_coefficient_val(islin,
 208                     libisl.isl_dim_set, index, islval)
 209             if inequality.constant != 0:
 210                 islval = str(inequality.constant).encode()
 211                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 212                 islin = libisl.isl_constraint_set_constant_val(islin, islval)
 213             islbset = libisl.isl_basic_set_add_constraint(islbset, islin)
 214         return islbset
 215
 216     @classmethod
 217     def fromstring(cls, string):
 218         domain = Domain.fromstring(string)
 219         if not isinstance(domain, Polyhedron):
 220             raise ValueError('non-polyhedral expression: {!r}'.format(string))
 221         return domain
 222
 223     def __repr__(self):
 224         strings = []
 225         for equality in self.equalities:
 226             strings.append('Eq({}, 0)'.format(equality))
 227         for inequality in self.inequalities:
 228             strings.append('Ge({}, 0)'.format(inequality))
 229         if len(strings) == 1:
 230             return strings[0]
 231         else:
 232             return 'And({})'.format(', '.join(strings))
 233
 234
 235     def _repr_latex_(self):
 236         strings = []
 237         for equality in self.equalities:
 238             strings.append('{} = 0'.format(equality._repr_latex_().strip('$')))
 239         for inequality in self.inequalities:
 240             strings.append('{} \\ge 0'.format(inequality._repr_latex_().strip('$')))
 241         return '$${}$$'.format(' \\wedge '.join(strings))
 242
 243     @classmethod
 244     def fromsympy(cls, expr):
 245     """
 246     Convert a sympy object to an expression.
 247     """
 248         domain = Domain.fromsympy(expr)
 249         if not isinstance(domain, Polyhedron):
 250             raise ValueError('non-polyhedral expression: {!r}'.format(expr))
 251         return domain
 252
 253     def tosympy(self):
 254     """
 255     Return an expression as a sympy object.
 256     """
 257         import sympy
 258         constraints = []
 259         for equality in self.equalities:
 260             constraints.append(sympy.Eq(equality.tosympy(), 0))
 261         for inequality in self.inequalities:
 262             constraints.append(sympy.Ge(inequality.tosympy(), 0))
 263         return sympy.And(*constraints)
 264
 265 class EmptyType(Polyhedron):
 266
 267     __slots__ = Polyhedron.__slots__
 268
 269     def __new__(cls):
 270         self = object().__new__(cls)
 271         self._equalities = (Rational(1),)
 272         self._inequalities = ()
 273         self._constraints = self._equalities
 274         self._symbols = ()
 275         self._dimension = 0
 276         return self
 277
 278     def widen(self, other):
 279         if not isinstance(other, Polyhedron):
 280             raise ValueError('argument must be a Polyhedron instance')
 281         return other
 282
 283     def __repr__(self):
 284         return 'Empty'
 285
 286     def _repr_latex_(self):
 287         return '$$\\emptyset$$'
 288
 289 Empty = EmptyType()
 290
 291
 292 class UniverseType(Polyhedron):
 293
 294     __slots__ = Polyhedron.__slots__
 295
 296     def __new__(cls):
 297         self = object().__new__(cls)
 298         self._equalities = ()
 299         self._inequalities = ()
 300         self._constraints = ()
 301         self._symbols = ()
 302         self._dimension = ()
 303         return self
 304
 305     def __repr__(self):
 306         return 'Universe'
 307
 308     def _repr_latex_(self):
 309         return '$$\\Omega$$'
 310
 311 Universe = UniverseType()
 312
 313
 314 def _polymorphic(func):
 315     @functools.wraps(func)
 316     def wrapper(left, right):
 317         if not isinstance(left, Expression):
 318             if isinstance(left, numbers.Rational):
 319                 left = Rational(left)
 320             else:
 321                 raise TypeError('left must be a a rational number '
 322                     'or a linear expression')
 323         if not isinstance(right, Expression):
 324             if isinstance(right, numbers.Rational):
 325                 right = Rational(right)
 326             else:
 327                 raise TypeError('right must be a a rational number '
 328                     'or a linear expression')
 329         return func(left, right)
 330     return wrapper
 331
 332 @_polymorphic
 333 def Lt(left, right):
 334     """
 335     Assert first set is less than the second set.
 336     """
 337     return Polyhedron([], [right - left - 1])
 338
 339 @_polymorphic
 340 def Le(left, right):
 341     """
 342     Assert first set is less than or equal to the second set.
 343     """
 344     return Polyhedron([], [right - left])
 345
 346 @_polymorphic
 347 def Eq(left, right):
 348     """
 349     Assert first set is equal to the second set.
 350     """
 351     return Polyhedron([left - right], [])
 352
 353 @_polymorphic
 354 def Ne(left, right):
 355     """
 356     Assert first set is not equal to the second set.
 357     """
 358     return ~Eq(left, right)
 359
 360 @_polymorphic
 361 def Gt(left, right):
 362     """
 363     Assert first set is greater than the second set.
 364     """
 365     return Polyhedron([], [left - right - 1])
 366
 367 @_polymorphic
 368 def Ge(left, right):
 369     """
 370     Assert first set is greater than or equal to the second set.
 371     """
 372     return Polyhedron([], [left - right])
 373