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