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