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