linpy/domains.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 ast
  19 import functools
  20 import re
  21 import math
  22
  23 from fractions import Fraction
  24
  25 from . import islhelper
  26 from .islhelper import mainctx, libisl
  27 from .linexprs import LinExpr, Symbol, Rational
  28 from .geometry import GeometricObject, Point, Vector
  29
  30
  31 __all__ = [
  32     'Domain',
  33     'And', 'Or', 'Not',
  34 ]
  35
  36
  37 @functools.total_ordering
  38 class Domain(GeometricObject):
  39     """
  40     A domain is a union of polyhedra. Unlike polyhedra, domains allow exact
  41     computation of union and complementary operations.
  42
  43     A domain with a unique polyhedron is automatically subclassed as a
  44     Polyhedron instance.
  45     """
  46
  47     __slots__ = (
  48         '_polyhedra',
  49         '_symbols',
  50         '_dimension',
  51     )
  52
  53     def __new__(cls, *polyhedra):
  54         """
  55         Return a domain from a sequence of polyhedra.
  56
  57         >>> square = Polyhedron('0 <= x <= 2, 0 <= y <= 2')
  58         >>> square2 = Polyhedron('2 <= x <= 4, 2 <= y <= 4')
  59         >>> dom = Domain([square, square2])
  60
  61         It is also possible to build domains from polyhedra using arithmetic
  62         operators Domain.__and__(), Domain.__or__() or functions And() and Or(),
  63         using one of the following instructions:
  64
  65         >>> square = Polyhedron('0 <= x <= 2, 0 <= y <= 2')
  66         >>> square2 = Polyhedron('2 <= x <= 4, 2 <= y <= 4')
  67         >>> dom = square | square2
  68         >>> dom = Or(square, square2)
  69
  70         Alternatively, a domain can be built from a string:
  71
  72         >>> dom = Domain('0 <= x <= 2, 0 <= y <= 2; 2 <= x <= 4, 2 <= y <= 4')
  73
  74         Finally, a domain can be built from a GeometricObject instance, calling
  75         the GeometricObject.asdomain() method.
  76         """
  77         from .polyhedra import Polyhedron
  78         if len(polyhedra) == 1:
  79             argument = polyhedra[0]
  80             if isinstance(argument, str):
  81                 return cls.fromstring(argument)
  82             elif isinstance(argument, GeometricObject):
  83                 return argument.aspolyhedron()
  84             else:
  85                 raise TypeError('argument must be a string '
  86                     'or a GeometricObject instance')
  87         else:
  88             for polyhedron in polyhedra:
  89                 if not isinstance(polyhedron, Polyhedron):
  90                     raise TypeError('arguments must be Polyhedron instances')
  91             symbols = cls._xsymbols(polyhedra)
  92             islset = cls._toislset(polyhedra, symbols)
  93             return cls._fromislset(islset, symbols)
  94
  95     @classmethod
  96     def _xsymbols(cls, iterator):
  97         """
  98         Return the ordered tuple of symbols present in iterator.
  99         """
 100         symbols = set()
 101         for item in iterator:
 102             symbols.update(item.symbols)
 103         return tuple(sorted(symbols, key=Symbol.sortkey))
 104
 105     @property
 106     def polyhedra(self):
 107         """
 108         The tuple of polyhedra present in the domain.
 109         """
 110         return self._polyhedra
 111
 112     @property
 113     def symbols(self):
 114         """
 115         The tuple of symbols present in the domain equations, sorted according
 116         to Symbol.sortkey().
 117         """
 118         return self._symbols
 119
 120     @property
 121     def dimension(self):
 122         """
 123         The dimension of the domain, i.e. the number of symbols present in it.
 124         """
 125         return self._dimension
 126
 127     def isempty(self):
 128         """
 129         Return True if the domain is empty, that is, equal to Empty.
 130         """
 131         islset = self._toislset(self.polyhedra, self.symbols)
 132         empty = bool(libisl.isl_set_is_empty(islset))
 133         libisl.isl_set_free(islset)
 134         return empty
 135
 136     def __bool__(self):
 137         """
 138         Return True if the domain is non-empty.
 139         """
 140         return not self.isempty()
 141
 142     def isuniverse(self):
 143         """
 144         Return True if the domain is universal, that is, equal to Universe.
 145         """
 146         islset = self._toislset(self.polyhedra, self.symbols)
 147         universe = bool(libisl.isl_set_plain_is_universe(islset))
 148         libisl.isl_set_free(islset)
 149         return universe
 150
 151     def isbounded(self):
 152         """
 153         Return True if the domain is bounded.
 154         """
 155         islset = self._toislset(self.polyhedra, self.symbols)
 156         bounded = bool(libisl.isl_set_is_bounded(islset))
 157         libisl.isl_set_free(islset)
 158         return bounded
 159
 160     def __eq__(self, other):
 161         """
 162         Return True if two domains are equal.
 163         """
 164         if isinstance(other, Domain):
 165             symbols = self._xsymbols([self, other])
 166             islset1 = self._toislset(self.polyhedra, symbols)
 167             islset2 = other._toislset(other.polyhedra, symbols)
 168             equal = bool(libisl.isl_set_is_equal(islset1, islset2))
 169             libisl.isl_set_free(islset1)
 170             libisl.isl_set_free(islset2)
 171             return equal
 172         return NotImplemented
 173
 174     def isdisjoint(self, other):
 175         """
 176         Return True if two domains have a null intersection.
 177         """
 178         if not isinstance(other, Domain):
 179             raise TypeError('other must be a Domain instance')
 180         symbols = self._xsymbols([self, other])
 181         islset1 = self._toislset(self.polyhedra, symbols)
 182         islset2 = self._toislset(other.polyhedra, symbols)
 183         equal = bool(libisl.isl_set_is_disjoint(islset1, islset2))
 184         libisl.isl_set_free(islset1)
 185         libisl.isl_set_free(islset2)
 186         return equal
 187
 188     def issubset(self, other):
 189         """
 190         Report whether another domain contains the domain.
 191         """
 192         return self < other
 193
 194     def __le__(self, other):
 195         if isinstance(other, Domain):
 196             symbols = self._xsymbols([self, other])
 197             islset1 = self._toislset(self.polyhedra, symbols)
 198             islset2 = self._toislset(other.polyhedra, symbols)
 199             equal = bool(libisl.isl_set_is_subset(islset1, islset2))
 200             libisl.isl_set_free(islset1)
 201             libisl.isl_set_free(islset2)
 202             return equal
 203         return NotImplemented
 204     __le__.__doc__ = issubset.__doc__
 205
 206     def __lt__(self, other):
 207         """
 208         Report whether another domain is contained within the domain.
 209         """
 210         if isinstance(other, Domain):
 211             symbols = self._xsymbols([self, other])
 212             islset1 = self._toislset(self.polyhedra, symbols)
 213             islset2 = self._toislset(other.polyhedra, symbols)
 214             equal = bool(libisl.isl_set_is_strict_subset(islset1, islset2))
 215             libisl.isl_set_free(islset1)
 216             libisl.isl_set_free(islset2)
 217             return equal
 218         return NotImplemented
 219
 220     def complement(self):
 221         """
 222         Return the complementary domain of the domain.
 223         """
 224         islset = self._toislset(self.polyhedra, self.symbols)
 225         islset = libisl.isl_set_complement(islset)
 226         return self._fromislset(islset, self.symbols)
 227
 228     def __invert__(self):
 229         return self.complement()
 230     __invert__.__doc__ = complement.__doc__
 231
 232     def make_disjoint(self):
 233         """
 234         Return an equivalent domain, whose polyhedra are disjoint.
 235         """
 236         islset = self._toislset(self.polyhedra, self.symbols)
 237         islset = libisl.isl_set_make_disjoint(mainctx, islset)
 238         return self._fromislset(islset, self.symbols)
 239
 240     def coalesce(self):
 241         """
 242         Simplify the representation of the domain by trying to combine pairs of
 243         polyhedra into a single polyhedron, and return the resulting domain.
 244         """
 245         islset = self._toislset(self.polyhedra, self.symbols)
 246         islset = libisl.isl_set_coalesce(islset)
 247         return self._fromislset(islset, self.symbols)
 248
 249     def detect_equalities(self):
 250         """
 251         Simplify the representation of the domain by detecting implicit
 252         equalities, and return the resulting domain.
 253         """
 254         islset = self._toislset(self.polyhedra, self.symbols)
 255         islset = libisl.isl_set_detect_equalities(islset)
 256         return self._fromislset(islset, self.symbols)
 257
 258     def remove_redundancies(self):
 259         """
 260         Remove redundant constraints in the domain, and return the resulting
 261         domain.
 262         """
 263         islset = self._toislset(self.polyhedra, self.symbols)
 264         islset = libisl.isl_set_remove_redundancies(islset)
 265         return self._fromislset(islset, self.symbols)
 266
 267     def aspolyhedron(self):
 268         from .polyhedra import Polyhedron
 269         islset = self._toislset(self.polyhedra, self.symbols)
 270         islbset = libisl.isl_set_polyhedral_hull(islset)
 271         return Polyhedron._fromislbasicset(islbset, self.symbols)
 272
 273     def asdomain(self):
 274         return self
 275
 276     def project(self, symbols):
 277         """
 278         Project out the sequence of symbols given in arguments, and return the
 279         resulting domain.
 280         """
 281         symbols = list(symbols)
 282         for symbol in symbols:
 283             if not isinstance(symbol, Symbol):
 284                 raise TypeError('symbols must be Symbol instances')
 285         islset = self._toislset(self.polyhedra, self.symbols)
 286         n = 0
 287         for index, symbol in reversed(list(enumerate(self.symbols))):
 288             if symbol in symbols:
 289                 n += 1
 290             elif n > 0:
 291                 islset = libisl.isl_set_project_out(islset,
 292                     libisl.isl_dim_set, index + 1, n)
 293                 n = 0
 294         if n > 0:
 295             islset = libisl.isl_set_project_out(islset, libisl.isl_dim_set, 0, n)
 296         symbols = [symbol for symbol in self.symbols if symbol not in symbols]
 297         return Domain._fromislset(islset, symbols)
 298
 299     def sample(self):
 300         """
 301         Return a sample of the domain, as an integer instance of Point. If the
 302         domain is empty, a ValueError exception is raised.
 303         """
 304         islset = self._toislset(self.polyhedra, self.symbols)
 305         islpoint = libisl.isl_set_sample_point(islset)
 306         if bool(libisl.isl_point_is_void(islpoint)):
 307             libisl.isl_point_free(islpoint)
 308             raise ValueError('domain must be non-empty')
 309         point = {}
 310         for index, symbol in enumerate(self.symbols):
 311             coordinate = libisl.isl_point_get_coordinate_val(islpoint,
 312                 libisl.isl_dim_set, index)
 313             coordinate = islhelper.isl_val_to_int(coordinate)
 314             point[symbol] = coordinate
 315         libisl.isl_point_free(islpoint)
 316         return point
 317
 318     def intersection(self, *others):
 319         """
 320         Return the intersection of two or more domains as a new domain. As an
 321         alternative, function And() can be used.
 322         """
 323         result = self
 324         for other in others:
 325             result &= other
 326         return result
 327
 328     def __and__(self, other):
 329         if isinstance(other, Domain):
 330             symbols = self._xsymbols([self, other])
 331             islset1 = self._toislset(self.polyhedra, symbols)
 332             islset2 = other._toislset(other.polyhedra, symbols)
 333             islset = libisl.isl_set_intersect(islset1, islset2)
 334             return self._fromislset(islset, symbols)
 335         return NotImplemented
 336     __and__.__doc__ = intersection.__doc__
 337
 338     def union(self, *others):
 339         """
 340         Return the union of two or more domains as a new domain. As an
 341         alternative, function Or() can be used.
 342         """
 343         if len(others) == 0:
 344             return self
 345         symbols = self._xsymbols((self,) + others)
 346         islset1 = self._toislset(self.polyhedra, symbols)
 347         for other in others:
 348             islset2 = other._toislset(other.polyhedra, symbols)
 349             islset1 = libisl.isl_set_union(islset1, islset2)
 350         return self._fromislset(islset1, symbols)
 351
 352     def __or__(self, other):
 353         return self.union(other)
 354     __or__.__doc__ = union.__doc__
 355
 356     def __add__(self, other):
 357         return self.union(other)
 358     __add__.__doc__ = union.__doc__
 359
 360     def difference(self, other):
 361         """
 362         Return the difference of two domains as a new domain.
 363         """
 364         symbols = self._xsymbols([self, other])
 365         islset1 = self._toislset(self.polyhedra, symbols)
 366         islset2 = other._toislset(other.polyhedra, symbols)
 367         islset = libisl.isl_set_subtract(islset1, islset2)
 368         return self._fromislset(islset, symbols)
 369
 370     def __sub__(self, other):
 371         return self.difference(other)
 372     __sub__.__doc__ = difference.__doc__
 373
 374     def lexmin(self):
 375         """
 376         Return the lexicographic minimum of the elements in the domain.
 377         """
 378         islset = self._toislset(self.polyhedra, self.symbols)
 379         islset = libisl.isl_set_lexmin(islset)
 380         return self._fromislset(islset, self.symbols)
 381
 382     def lexmax(self):
 383         """
 384         Return the lexicographic maximum of the elements in the domain.
 385         """
 386         islset = self._toislset(self.polyhedra, self.symbols)
 387         islset = libisl.isl_set_lexmax(islset)
 388         return self._fromislset(islset, self.symbols)
 389
 390     _RE_COORDINATE = re.compile(r'\((?P<num>\-?\d+)\)(/(?P<den>\d+))?')
 391
 392     def vertices(self):
 393         """
 394         Return the vertices of the domain, as a list of rational instances of
 395         Point.
 396         """
 397         from .polyhedra import Polyhedron
 398         if not self.isbounded():
 399             raise ValueError('domain must be bounded')
 400         islbset = self._toislbasicset(self.equalities, self.inequalities,
 401             self.symbols)
 402         vertices = libisl.isl_basic_set_compute_vertices(islbset);
 403         vertices = islhelper.isl_vertices_vertices(vertices)
 404         points = []
 405         for vertex in vertices:
 406             expr = libisl.isl_vertex_get_expr(vertex)
 407             coordinates = []
 408             if islhelper.isl_version < '0.13':
 409                 constraints = islhelper.isl_basic_set_constraints(expr)
 410                 for constraint in constraints:
 411                     constant = libisl.isl_constraint_get_constant_val(constraint)
 412                     constant = islhelper.isl_val_to_int(constant)
 413                     for index, symbol in enumerate(self.symbols):
 414                         coefficient = libisl.isl_constraint_get_coefficient_val(constraint,
 415                             libisl.isl_dim_set, index)
 416                         coefficient = islhelper.isl_val_to_int(coefficient)
 417                         if coefficient != 0:
 418                             coordinate = -Fraction(constant, coefficient)
 419                             coordinates.append((symbol, coordinate))
 420             else:
 421                 string = islhelper.isl_multi_aff_to_str(expr)
 422                 matches = self._RE_COORDINATE.finditer(string)
 423                 for symbol, match in zip(self.symbols, matches):
 424                     numerator = int(match.group('num'))
 425                     denominator = match.group('den')
 426                     denominator = 1 if denominator is None else int(denominator)
 427                     coordinate = Fraction(numerator, denominator)
 428                     coordinates.append((symbol, coordinate))
 429             points.append(Point(coordinates))
 430         return points
 431
 432     def points(self):
 433         """
 434         Return the integer points of a bounded domain, as a list of integer
 435         instances of Point. If the domain is not bounded, a ValueError exception
 436         is raised.
 437         """
 438         if not self.isbounded():
 439             raise ValueError('domain must be bounded')
 440         from .polyhedra import Universe, Eq
 441         islset = self._toislset(self.polyhedra, self.symbols)
 442         islpoints = islhelper.isl_set_points(islset)
 443         points = []
 444         for islpoint in islpoints:
 445             coordinates = {}
 446             for index, symbol in enumerate(self.symbols):
 447                 coordinate = libisl.isl_point_get_coordinate_val(islpoint,
 448                     libisl.isl_dim_set, index)
 449                 coordinate = islhelper.isl_val_to_int(coordinate)
 450                 coordinates[symbol] = coordinate
 451             points.append(Point(coordinates))
 452         return points
 453
 454     def __contains__(self, point):
 455         """
 456         Return True if the point if contained within the domain.
 457         """
 458         for polyhedron in self.polyhedra:
 459             if point in polyhedron:
 460                 return True
 461         return False
 462
 463     @classmethod
 464     def _polygon_inner_point(cls, points):
 465         symbols = points[0].symbols
 466         coordinates = {symbol: 0 for symbol in symbols}
 467         for point in points:
 468             for symbol, coordinate in point.coordinates():
 469                 coordinates[symbol] += coordinate
 470         for symbol in symbols:
 471             coordinates[symbol] /= len(points)
 472         return Point(coordinates)
 473
 474     @classmethod
 475     def _sort_polygon_2d(cls, points):
 476         if len(points) <= 3:
 477             return points
 478         o = cls._polygon_inner_point(points)
 479         angles = {}
 480         for m in points:
 481             om = Vector(o, m)
 482             dx, dy = (coordinate for symbol, coordinate in om.coordinates())
 483             angle = math.atan2(dy, dx)
 484             angles[m] = angle
 485         return sorted(points, key=angles.get)
 486
 487     @classmethod
 488     def _sort_polygon_3d(cls, points):
 489         if len(points) <= 3:
 490             return points
 491         o = cls._polygon_inner_point(points)
 492         a = points[0]
 493         oa = Vector(o, a)
 494         norm_oa = oa.norm()
 495         for b in points[1:]:
 496             ob = Vector(o, b)
 497             u = oa.cross(ob)
 498             if not u.isnull():
 499                 u = u.asunit()
 500                 break
 501         else:
 502             raise ValueError('degenerate polygon')
 503         angles = {a: 0.}
 504         for m in points[1:]:
 505             om = Vector(o, m)
 506             normprod = norm_oa * om.norm()
 507             cosinus = max(oa.dot(om) / normprod, -1.)
 508             sinus = u.dot(oa.cross(om)) / normprod
 509             angle = math.acos(cosinus)
 510             angle = math.copysign(angle, sinus)
 511             angles[m] = angle
 512         return sorted(points, key=angles.get)
 513
 514     def faces(self):
 515         """
 516         Return the list of faces of a bounded domain. Each face is represented
 517         by a list of vertices, in the form of rational instances of Point. If
 518         the domain is not bounded, a ValueError exception is raised.
 519         """
 520         faces = []
 521         for polyhedron in self.polyhedra:
 522             vertices = polyhedron.vertices()
 523             for constraint in polyhedron.constraints:
 524                 face = []
 525                 for vertex in vertices:
 526                     if constraint.subs(vertex.coordinates()) == 0:
 527                         face.append(vertex)
 528                 if len(face) >= 3:
 529                     faces.append(face)
 530         return faces
 531
 532     def _plot_2d(self, plot=None, **kwargs):
 533         import matplotlib.pyplot as plt
 534         from matplotlib.patches import Polygon
 535         if plot is None:
 536             fig = plt.figure()
 537             plot = fig.add_subplot(1, 1, 1)
 538         xmin, xmax = plot.get_xlim()
 539         ymin, ymax = plot.get_ylim()
 540         for polyhedron in self.polyhedra:
 541             vertices = polyhedron._sort_polygon_2d(polyhedron.vertices())
 542             xys = [tuple(vertex.values()) for vertex in vertices]
 543             xs, ys = zip(*xys)
 544             xmin, xmax = min(xmin, float(min(xs))), max(xmax, float(max(xs)))
 545             ymin, ymax = min(ymin, float(min(ys))), max(ymax, float(max(ys)))
 546             plot.add_patch(Polygon(xys, closed=True, **kwargs))
 547         plot.set_xlim(xmin, xmax)
 548         plot.set_ylim(ymin, ymax)
 549         return plot
 550
 551     def _plot_3d(self, plot=None, **kwargs):
 552         import matplotlib.pyplot as plt
 553         from mpl_toolkits.mplot3d import Axes3D
 554         from mpl_toolkits.mplot3d.art3d import Poly3DCollection
 555         if plot is None:
 556             fig = plt.figure()
 557             axes = Axes3D(fig)
 558         else:
 559             axes = plot
 560         xmin, xmax = axes.get_xlim()
 561         ymin, ymax = axes.get_ylim()
 562         zmin, zmax = axes.get_zlim()
 563         poly_xyzs = []
 564         for vertices in self.faces():
 565             vertices = self._sort_polygon_3d(vertices)
 566             vertices.append(vertices[0])
 567             face_xyzs = [tuple(vertex.values()) for vertex in vertices]
 568             xs, ys, zs = zip(*face_xyzs)
 569             xmin, xmax = min(xmin, float(min(xs))), max(xmax, float(max(xs)))
 570             ymin, ymax = min(ymin, float(min(ys))), max(ymax, float(max(ys)))
 571             zmin, zmax = min(zmin, float(min(zs))), max(zmax, float(max(zs)))
 572             poly_xyzs.append(face_xyzs)
 573         collection = Poly3DCollection(poly_xyzs, **kwargs)
 574         axes.add_collection3d(collection)
 575         axes.set_xlim(xmin, xmax)
 576         axes.set_ylim(ymin, ymax)
 577         axes.set_zlim(zmin, zmax)
 578         return axes
 579
 580     def plot(self, plot=None, **kwargs):
 581         """
 582         Plot a 2D or 3D domain using matplotlib. Draw it to the current plot
 583         object if present, otherwise create a new one. options are keyword
 584         arguments passed to the matplotlib drawing functions, they can be used
 585         to set the drawing color for example. Raise ValueError is the domain is
 586         not 2D or 3D.
 587         """
 588         if not self.isbounded():
 589             raise ValueError('domain must be bounded')
 590         elif self.dimension == 2:
 591             return self._plot_2d(plot=plot, **kwargs)
 592         elif self.dimension == 3:
 593             return self._plot_3d(plot=plot, **kwargs)
 594         else:
 595             raise ValueError('polyhedron must be 2 or 3-dimensional')
 596
 597     def subs(self, symbol, expression=None):
 598         """
 599         Substitute the given symbol by an expression in the domain constraints.
 600         To perform multiple substitutions at once, pass a sequence or a
 601         dictionary of (old, new) pairs to subs. The syntax of this function is
 602         similar to LinExpr.subs().
 603         """
 604         polyhedra = [polyhedron.subs(symbol, expression)
 605             for polyhedron in self.polyhedra]
 606         return Domain(*polyhedra)
 607
 608     @classmethod
 609     def _fromislset(cls, islset, symbols):
 610         from .polyhedra import Polyhedron
 611         islset = libisl.isl_set_remove_divs(islset)
 612         islbsets = islhelper.isl_set_basic_sets(islset)
 613         libisl.isl_set_free(islset)
 614         polyhedra = []
 615         for islbset in islbsets:
 616             polyhedron = Polyhedron._fromislbasicset(islbset, symbols)
 617             polyhedra.append(polyhedron)
 618         if len(polyhedra) == 0:
 619             from .polyhedra import Empty
 620             return Empty
 621         elif len(polyhedra) == 1:
 622             return polyhedra[0]
 623         else:
 624             self = object().__new__(Domain)
 625             self._polyhedra = tuple(polyhedra)
 626             self._symbols = cls._xsymbols(polyhedra)
 627             self._dimension = len(self._symbols)
 628             return self
 629
 630     @classmethod
 631     def _toislset(cls, polyhedra, symbols):
 632         polyhedron = polyhedra[0]
 633         islbset = polyhedron._toislbasicset(polyhedron.equalities,
 634             polyhedron.inequalities, symbols)
 635         islset1 = libisl.isl_set_from_basic_set(islbset)
 636         for polyhedron in polyhedra[1:]:
 637             islbset = polyhedron._toislbasicset(polyhedron.equalities,
 638                 polyhedron.inequalities, symbols)
 639             islset2 = libisl.isl_set_from_basic_set(islbset)
 640             islset1 = libisl.isl_set_union(islset1, islset2)
 641         return islset1
 642
 643     @classmethod
 644     def _fromast(cls, node):
 645         from .polyhedra import Polyhedron
 646         if isinstance(node, ast.Module) and len(node.body) == 1:
 647             return cls._fromast(node.body[0])
 648         elif isinstance(node, ast.Expr):
 649             return cls._fromast(node.value)
 650         elif isinstance(node, ast.UnaryOp):
 651             domain = cls._fromast(node.operand)
 652             if isinstance(node.operand, ast.invert):
 653                 return Not(domain)
 654         elif isinstance(node, ast.BinOp):
 655             domain1 = cls._fromast(node.left)
 656             domain2 = cls._fromast(node.right)
 657             if isinstance(node.op, ast.BitAnd):
 658                 return And(domain1, domain2)
 659             elif isinstance(node.op, ast.BitOr):
 660                 return Or(domain1, domain2)
 661         elif isinstance(node, ast.Compare):
 662             equalities = []
 663             inequalities = []
 664             left = LinExpr._fromast(node.left)
 665             for i in range(len(node.ops)):
 666                 op = node.ops[i]
 667                 right = LinExpr._fromast(node.comparators[i])
 668                 if isinstance(op, ast.Lt):
 669                     inequalities.append(right - left - 1)
 670                 elif isinstance(op, ast.LtE):
 671                     inequalities.append(right - left)
 672                 elif isinstance(op, ast.Eq):
 673                     equalities.append(left - right)
 674                 elif isinstance(op, ast.GtE):
 675                     inequalities.append(left - right)
 676                 elif isinstance(op, ast.Gt):
 677                     inequalities.append(left - right - 1)
 678                 else:
 679                     break
 680                 left = right
 681             else:
 682                 return Polyhedron(equalities, inequalities)
 683         raise SyntaxError('invalid syntax')
 684
 685     _RE_BRACES = re.compile(r'^\{\s*|\s*\}$')
 686     _RE_EQ = re.compile(r'([^<=>])=([^<=>])')
 687     _RE_AND = re.compile(r'\band\b|,|&&|/\\|∧|∩')
 688     _RE_OR = re.compile(r'\bor\b|;|\|\||\\/|∨|∪')
 689     _RE_NOT = re.compile(r'\bnot\b|!|¬')
 690     _RE_NUM_VAR = LinExpr._RE_NUM_VAR
 691     _RE_OPERATORS = re.compile(r'(&|\||~)')
 692
 693     @classmethod
 694     def fromstring(cls, string):
 695         """
 696         Create a domain from a string. Raise SyntaxError if the string is not
 697         properly formatted.
 698         """
 699         # remove curly brackets
 700         string = cls._RE_BRACES.sub(r'', string)
 701         # replace '=' by '=='
 702         string = cls._RE_EQ.sub(r'\1==\2', string)
 703         # replace 'and', 'or', 'not'
 704         string = cls._RE_AND.sub(r' & ', string)
 705         string = cls._RE_OR.sub(r' | ', string)
 706         string = cls._RE_NOT.sub(r' ~', string)
 707         # add implicit multiplication operators, e.g. '5x' -> '5*x'
 708         string = cls._RE_NUM_VAR.sub(r'\1*\2', string)
 709         # add parentheses to force precedence
 710         tokens = cls._RE_OPERATORS.split(string)
 711         for i, token in enumerate(tokens):
 712             if i % 2 == 0:
 713                 token = '({})'.format(token)
 714                 tokens[i] = token
 715         string = ''.join(tokens)
 716         tree = ast.parse(string, 'eval')
 717         return cls._fromast(tree)
 718
 719     def __repr__(self):
 720         assert len(self.polyhedra) >= 2
 721         strings = [repr(polyhedron) for polyhedron in self.polyhedra]
 722         return 'Or({})'.format(', '.join(strings))
 723
 724     def _repr_latex_(self):
 725         strings = []
 726         for polyhedron in self.polyhedra:
 727             strings.append('({})'.format(polyhedron._repr_latex_().strip('$')))
 728         return '${}$'.format(' \\vee '.join(strings))
 729
 730     @classmethod
 731     def fromsympy(cls, expr):
 732         """
 733         Create a domain from a sympy expression.
 734         """
 735         import sympy
 736         from .polyhedra import Lt, Le, Eq, Ne, Ge, Gt
 737         funcmap = {
 738             sympy.And: And, sympy.Or: Or, sympy.Not: Not,
 739             sympy.Lt: Lt, sympy.Le: Le,
 740             sympy.Eq: Eq, sympy.Ne: Ne,
 741             sympy.Ge: Ge, sympy.Gt: Gt,
 742         }
 743         if expr.func in funcmap:
 744             args = [Domain.fromsympy(arg) for arg in expr.args]
 745             return funcmap[expr.func](*args)
 746         elif isinstance(expr, sympy.Expr):
 747             return LinExpr.fromsympy(expr)
 748         raise ValueError('non-domain expression: {!r}'.format(expr))
 749
 750     def tosympy(self):
 751         """
 752         Convert the domain to a sympy expression.
 753         """
 754         import sympy
 755         polyhedra = [polyhedron.tosympy() for polyhedron in polyhedra]
 756         return sympy.Or(*polyhedra)
 757
 758
 759 def And(*domains):
 760     """
 761     Create the intersection domain of the domains given in arguments.
 762     """
 763     if len(domains) == 0:
 764         from .polyhedra import Universe
 765         return Universe
 766     else:
 767         return domains[0].intersection(*domains[1:])
 768 And.__doc__ = Domain.intersection.__doc__
 769
 770 def Or(*domains):
 771     """
 772     Create the union domain of the domains given in arguments.
 773     """
 774     if len(domains) == 0:
 775         from .polyhedra import Empty
 776         return Empty
 777     else:
 778         return domains[0].union(*domains[1:])
 779 Or.__doc__ = Domain.union.__doc__
 780
 781 def Not(domain):
 782     """
 783     Create the complementary domain of the domain given in argument.
 784     """
 785     return ~domain
 786 Not.__doc__ = Domain.complement.__doc__