pypol/polyhedra.py

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