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 numbers
  20
  21 from . import islhelper
  22
  23 from .domains import Domain
  24 from .geometry import GeometricObject, Point
  25 from .islhelper import libisl, mainctx
  26 from .linexprs import LinExpr, Rational
  27
  28
  29 __all__ = [
  30     'Empty',
  31     'Eq',
  32     'Ge',
  33     'Gt',
  34     'Le',
  35     'Lt',
  36     'Ne',
  37     'Polyhedron',
  38     'Universe',
  39 ]
  40
  41
  42 class Polyhedron(Domain):
  43     """
  44     A convex polyhedron (or simply "polyhedron") is the space defined by a
  45     system of linear equalities and inequalities. This space can be unbounded.
  46     A Z-polyhedron (simply called "polyhedron" in LinPy) is the set of integer
  47     points in a convex polyhedron.
  48     """
  49
  50     __slots__ = (
  51         '_equalities',
  52         '_inequalities',
  53         '_symbols',
  54         '_dimension',
  55     )
  56
  57     def __new__(cls, equalities=None, inequalities=None):
  58         """
  59         Return a polyhedron from two sequences of linear expressions:
  60         equalities is a list of expressions equal to 0, and inequalities is a
  61         list of expressions greater or equal to 0. For example, the polyhedron
  62         0 <= x <= 2, 0 <= y <= 2 can be constructed with:
  63
  64         >>> x, y = symbols('x y')
  65         >>> square1 = Polyhedron([], [x, 2 - x, y, 2 - y])
  66         >>> square1
  67         And(0 <= x, x <= 2, 0 <= y, y <= 2)
  68
  69         It may be easier to use comparison operators LinExpr.__lt__(),
  70         LinExpr.__le__(), LinExpr.__ge__(), LinExpr.__gt__(), or
  71         functions Lt(), Le(), Eq(), Ge() and Gt(), using one of the following
  72         instructions:
  73
  74         >>> x, y = symbols('x y')
  75         >>> square1 = (0 <= x) & (x <= 2) & (0 <= y) & (y <= 2)
  76         >>> square1 = Le(0, x, 2) & Le(0, y, 2)
  77
  78         It is also possible to build a polyhedron from a string.
  79
  80         >>> square1 = Polyhedron('0 <= x <= 2, 0 <= y <= 2')
  81
  82         Finally, a polyhedron can be constructed from a GeometricObject
  83         instance, calling the GeometricObject.aspolyedron() method. This way,
  84         it is possible to compute the polyhedral hull of a Domain instance,
  85         i.e., the convex hull of two polyhedra:
  86
  87         >>> square1 = Polyhedron('0 <= x <= 2, 0 <= y <= 2')
  88         >>> square2 = Polyhedron('1 <= x <= 3, 1 <= y <= 3')
  89         >>> Polyhedron(square1 | square2)
  90         And(0 <= x, 0 <= y, x <= y + 2, y <= x + 2, x <= 3, y <= 3)
  91         """
  92         if isinstance(equalities, str):
  93             if inequalities is not None:
  94                 raise TypeError('too many arguments')
  95             return cls.fromstring(equalities)
  96         elif isinstance(equalities, GeometricObject):
  97             if inequalities is not None:
  98                 raise TypeError('too many arguments')
  99             return equalities.aspolyhedron()
 100         sc_equalities = []
 101         if equalities is not None:
 102             for equality in equalities:
 103                 if isinstance(equality, LinExpr):
 104                     sc_equalities.append(equality.scaleint())
 105                 elif isinstance(equality, numbers.Rational):
 106                     sc_equalities.append(Rational(equality).scaleint())
 107                 else:
 108                     raise TypeError('equalities must be linear expressions '
 109                                     'or rational numbers')
 110         sc_inequalities = []
 111         if inequalities is not None:
 112             for inequality in inequalities:
 113                 if isinstance(inequality, LinExpr):
 114                     sc_inequalities.append(inequality.scaleint())
 115                 elif isinstance(inequality, numbers.Rational):
 116                     sc_inequalities.append(Rational(inequality).scaleint())
 117                 else:
 118                     raise TypeError('inequalities must be linear expressions '
 119                                     'or rational numbers')
 120         symbols = cls._xsymbols(sc_equalities + sc_inequalities)
 121         islbset = cls._toislbasicset(sc_equalities, sc_inequalities, symbols)
 122         return cls._fromislbasicset(islbset, symbols)
 123
 124     @property
 125     def equalities(self):
 126         """
 127         The tuple of equalities. This is a list of LinExpr instances that are
 128         equal to 0 in the polyhedron.
 129         """
 130         return self._equalities
 131
 132     @property
 133     def inequalities(self):
 134         """
 135         The tuple of inequalities. This is a list of LinExpr instances that are
 136         greater or equal to 0 in the polyhedron.
 137         """
 138         return self._inequalities
 139
 140     @property
 141     def constraints(self):
 142         """
 143         The tuple of constraints, i.e., equalities and inequalities. This is
 144         semantically equivalent to: equalities + inequalities.
 145         """
 146         return self._equalities + self._inequalities
 147
 148     @property
 149     def polyhedra(self):
 150         return self,
 151
 152     def make_disjoint(self):
 153         return self
 154
 155     def isuniverse(self):
 156         islbset = self._toislbasicset(self.equalities, self.inequalities,
 157                                       self.symbols)
 158         universe = bool(libisl.isl_basic_set_is_universe(islbset))
 159         libisl.isl_basic_set_free(islbset)
 160         return universe
 161
 162     def aspolyhedron(self):
 163         return self
 164
 165     def convex_union(self, *others):
 166         """
 167         Return the convex union of two or more polyhedra.
 168         """
 169         for other in others:
 170             if not isinstance(other, Polyhedron):
 171                 raise TypeError('arguments must be Polyhedron instances')
 172         return Polyhedron(self.union(*others))
 173
 174     def __contains__(self, point):
 175         if not isinstance(point, Point):
 176             raise TypeError('point must be a Point instance')
 177         if self.symbols != point.symbols:
 178             raise ValueError('arguments must belong to the same space')
 179         for equality in self.equalities:
 180             if equality.subs(point.coordinates()) != 0:
 181                 return False
 182         for inequality in self.inequalities:
 183             if inequality.subs(point.coordinates()) < 0:
 184                 return False
 185         return True
 186
 187     def subs(self, symbol, expression=None):
 188         equalities = [equality.subs(symbol, expression)
 189                       for equality in self.equalities]
 190         inequalities = [inequality.subs(symbol, expression)
 191                         for inequality in self.inequalities]
 192         return Polyhedron(equalities, inequalities)
 193
 194     def asinequalities(self):
 195         """
 196         Express the polyhedron using inequalities, given as a list of
 197         expressions greater or equal to 0.
 198         """
 199         inequalities = list(self.equalities)
 200         inequalities.extend([-expression for expression in self.equalities])
 201         inequalities.extend(self.inequalities)
 202         return inequalities
 203
 204     def widen(self, other):
 205         """
 206         Compute the standard widening of two polyhedra, à la Halbwachs.
 207
 208         In its current implementation, this method is slow and should not be
 209         used on large polyhedra.
 210         """
 211         if not isinstance(other, Polyhedron):
 212             raise TypeError('argument must be a Polyhedron instance')
 213         inequalities1 = self.asinequalities()
 214         inequalities2 = other.asinequalities()
 215         inequalities = []
 216         for inequality1 in inequalities1:
 217             if other <= Polyhedron(inequalities=[inequality1]):
 218                 inequalities.append(inequality1)
 219         for inequality2 in inequalities2:
 220             for i in range(len(inequalities1)):
 221                 inequalities3 = inequalities1[:i] + inequalities[i + 1:]
 222                 inequalities3.append(inequality2)
 223                 polyhedron3 = Polyhedron(inequalities=inequalities3)
 224                 if self == polyhedron3:
 225                     inequalities.append(inequality2)
 226                     break
 227         return Polyhedron(inequalities=inequalities)
 228
 229     @classmethod
 230     def _fromislbasicset(cls, islbset, symbols):
 231         if bool(libisl.isl_basic_set_is_empty(islbset)):
 232             return Empty
 233         if bool(libisl.isl_basic_set_is_universe(islbset)):
 234             return Universe
 235         islconstraints = islhelper.isl_basic_set_constraints(islbset)
 236         equalities = []
 237         inequalities = []
 238         for islconstraint in islconstraints:
 239             constant = libisl.isl_constraint_get_constant_val(islconstraint)
 240             constant = islhelper.isl_val_to_int(constant)
 241             coefficients = {}
 242             for index, symbol in enumerate(symbols):
 243                 coefficient = libisl.isl_constraint_get_coefficient_val(
 244                     islconstraint, libisl.isl_dim_set, index)
 245                 coefficient = islhelper.isl_val_to_int(coefficient)
 246                 if coefficient != 0:
 247                     coefficients[symbol] = coefficient
 248             expression = LinExpr(coefficients, constant)
 249             if libisl.isl_constraint_is_equality(islconstraint):
 250                 equalities.append(expression)
 251             else:
 252                 inequalities.append(expression)
 253         libisl.isl_basic_set_free(islbset)
 254         self = object().__new__(Polyhedron)
 255         self._equalities = tuple(equalities)
 256         self._inequalities = tuple(inequalities)
 257         self._symbols = cls._xsymbols(self.constraints)
 258         self._dimension = len(self._symbols)
 259         return self
 260
 261     @classmethod
 262     def _toislbasicset(cls, equalities, inequalities, symbols):
 263         dimension = len(symbols)
 264         indices = {symbol: index for index, symbol in enumerate(symbols)}
 265         islsp = libisl.isl_space_set_alloc(mainctx, 0, dimension)
 266         islbset = libisl.isl_basic_set_universe(libisl.isl_space_copy(islsp))
 267         islls = libisl.isl_local_space_from_space(islsp)
 268         for equality in equalities:
 269             isleq = libisl.isl_equality_alloc(
 270                 libisl.isl_local_space_copy(islls))
 271             for symbol, coefficient in equality.coefficients():
 272                 islval = str(coefficient).encode()
 273                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 274                 index = indices[symbol]
 275                 isleq = libisl.isl_constraint_set_coefficient_val(
 276                     isleq, libisl.isl_dim_set, index, islval)
 277             if equality.constant != 0:
 278                 islval = str(equality.constant).encode()
 279                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 280                 isleq = libisl.isl_constraint_set_constant_val(isleq, islval)
 281             islbset = libisl.isl_basic_set_add_constraint(islbset, isleq)
 282         for inequality in inequalities:
 283             islin = libisl.isl_inequality_alloc(
 284                 libisl.isl_local_space_copy(islls))
 285             for symbol, coefficient in inequality.coefficients():
 286                 islval = str(coefficient).encode()
 287                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 288                 index = indices[symbol]
 289                 islin = libisl.isl_constraint_set_coefficient_val(
 290                     islin, libisl.isl_dim_set, index, islval)
 291             if inequality.constant != 0:
 292                 islval = str(inequality.constant).encode()
 293                 islval = libisl.isl_val_read_from_str(mainctx, islval)
 294                 islin = libisl.isl_constraint_set_constant_val(islin, islval)
 295             islbset = libisl.isl_basic_set_add_constraint(islbset, islin)
 296         return islbset
 297
 298     @classmethod
 299     def fromstring(cls, string):
 300         domain = Domain.fromstring(string)
 301         if not isinstance(domain, Polyhedron):
 302             raise ValueError('non-polyhedral expression: {!r}'.format(string))
 303         return domain
 304
 305     def __repr__(self):
 306         strings = []
 307         for equality in self.equalities:
 308             left, right, swap = 0, 0, False
 309             for i, (symbol, coefficient) in enumerate(equality.coefficients()):
 310                 if coefficient > 0:
 311                     left += coefficient * symbol
 312                 else:
 313                     right -= coefficient * symbol
 314                     if i == 0:
 315                         swap = True
 316             if equality.constant > 0:
 317                 left += equality.constant
 318             else:
 319                 right -= equality.constant
 320             if swap:
 321                 left, right = right, left
 322             strings.append('{} == {}'.format(left, right))
 323         for inequality in self.inequalities:
 324             left, right = 0, 0
 325             for symbol, coefficient in inequality.coefficients():
 326                 if coefficient < 0:
 327                     left -= coefficient * symbol
 328                 else:
 329                     right += coefficient * symbol
 330             if inequality.constant < 0:
 331                 left -= inequality.constant
 332             else:
 333                 right += inequality.constant
 334             strings.append('{} <= {}'.format(left, right))
 335         if len(strings) == 1:
 336             return strings[0]
 337         else:
 338             return 'And({})'.format(', '.join(strings))
 339
 340     @classmethod
 341     def fromsympy(cls, expression):
 342         domain = Domain.fromsympy(expression)
 343         if not isinstance(domain, Polyhedron):
 344             raise ValueError('non-polyhedral expression: {!r}'.format(
 345                 expression))
 346         return domain
 347
 348     def tosympy(self):
 349         import sympy
 350         constraints = []
 351         for equality in self.equalities:
 352             constraints.append(sympy.Eq(equality.tosympy(), 0))
 353         for inequality in self.inequalities:
 354             constraints.append(sympy.Ge(inequality.tosympy(), 0))
 355         return sympy.And(*constraints)
 356
 357
 358 class EmptyType(Polyhedron):
 359     """
 360     The empty polyhedron, whose set of constraints is not satisfiable.
 361     """
 362
 363     def __new__(cls):
 364         self = object().__new__(cls)
 365         self._equalities = (Rational(1),)
 366         self._inequalities = ()
 367         self._symbols = ()
 368         self._dimension = 0
 369         return self
 370
 371     def widen(self, other):
 372         if not isinstance(other, Polyhedron):
 373             raise ValueError('argument must be a Polyhedron instance')
 374         return other
 375
 376     def __repr__(self):
 377         return 'Empty'
 378
 379 Empty = EmptyType()
 380
 381
 382 class UniverseType(Polyhedron):
 383     """
 384     The universe polyhedron, whose set of constraints is always satisfiable,
 385     i.e. is empty.
 386     """
 387
 388     def __new__(cls):
 389         self = object().__new__(cls)
 390         self._equalities = ()
 391         self._inequalities = ()
 392         self._symbols = ()
 393         self._dimension = ()
 394         return self
 395
 396     def __repr__(self):
 397         return 'Universe'
 398
 399 Universe = UniverseType()
 400
 401
 402 def _pseudoconstructor(func):
 403     @functools.wraps(func)
 404     def wrapper(expression1, expression2, *expressions):
 405         expressions = (expression1, expression2) + expressions
 406         for expression in expressions:
 407             if not isinstance(expression, LinExpr):
 408                 if isinstance(expression, numbers.Rational):
 409                     expression = Rational(expression)
 410                 else:
 411                     raise TypeError('arguments must be rational numbers '
 412                                     'or linear expressions')
 413         return func(*expressions)
 414     return wrapper
 415
 416
 417 @_pseudoconstructor
 418 def Lt(*expressions):
 419     """
 420     Create the polyhedron with constraints expr1 < expr2 < expr3 ...
 421     """
 422     inequalities = []
 423     for left, right in zip(expressions, expressions[1:]):
 424         inequalities.append(right - left - 1)
 425     return Polyhedron([], inequalities)
 426
 427
 428 @_pseudoconstructor
 429 def Le(*expressions):
 430     """
 431     Create the polyhedron with constraints expr1 <= expr2 <= expr3 ...
 432     """
 433     inequalities = []
 434     for left, right in zip(expressions, expressions[1:]):
 435         inequalities.append(right - left)
 436     return Polyhedron([], inequalities)
 437
 438
 439 @_pseudoconstructor
 440 def Eq(*expressions):
 441     """
 442     Create the polyhedron with constraints expr1 == expr2 == expr3 ...
 443     """
 444     equalities = []
 445     for left, right in zip(expressions, expressions[1:]):
 446         equalities.append(left - right)
 447     return Polyhedron(equalities, [])
 448
 449
 450 @_pseudoconstructor
 451 def Ne(*expressions):
 452     """
 453     Create the domain such that expr1 != expr2 != expr3 ... The result is a
 454     Domain object, not a Polyhedron.
 455     """
 456     domain = Universe
 457     for left, right in zip(expressions, expressions[1:]):
 458         domain &= ~Eq(left, right)
 459     return domain
 460
 461
 462 @_pseudoconstructor
 463 def Ge(*expressions):
 464     """
 465     Create the polyhedron with constraints expr1 >= expr2 >= expr3 ...
 466     """
 467     inequalities = []
 468     for left, right in zip(expressions, expressions[1:]):
 469         inequalities.append(left - right)
 470     return Polyhedron([], inequalities)
 471
 472
 473 @_pseudoconstructor
 474 def Gt(*expressions):
 475     """
 476     Create the polyhedron with constraints expr1 > expr2 > expr3 ...
 477     """
 478     inequalities = []
 479     for left, right in zip(expressions, expressions[1:]):
 480         inequalities.append(left - right - 1)
 481     return Polyhedron([], inequalities)