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         islset = self._toislset(self.polyhedra, self.symbols)
 282         n = 0
 283         for index, symbol in reversed(list(enumerate(self.symbols))):
 284             if symbol in symbols:
 285                 n += 1
 286             elif n > 0:
 287                 islset = libisl.isl_set_project_out(islset,
 288                     libisl.isl_dim_set, index + 1, n)
 289                 n = 0
 290         if n > 0:
 291             islset = libisl.isl_set_project_out(islset, libisl.isl_dim_set, 0, n)
 292         symbols = [symbol for symbol in self.symbols if symbol not in symbols]
 293         return Domain._fromislset(islset, symbols)
 294
 295     def sample(self):
 296         """
 297         Return a sample of the domain, as an integer instance of Point. If the
 298         domain is empty, a ValueError exception is raised.
 299         """
 300         islset = self._toislset(self.polyhedra, self.symbols)
 301         islpoint = libisl.isl_set_sample_point(islset)
 302         if bool(libisl.isl_point_is_void(islpoint)):
 303             libisl.isl_point_free(islpoint)
 304             raise ValueError('domain must be non-empty')
 305         point = {}
 306         for index, symbol in enumerate(self.symbols):
 307             coordinate = libisl.isl_point_get_coordinate_val(islpoint,
 308                 libisl.isl_dim_set, index)
 309             coordinate = islhelper.isl_val_to_int(coordinate)
 310             point[symbol] = coordinate
 311         libisl.isl_point_free(islpoint)
 312         return point
 313
 314     def intersection(self, *others):
 315         """
 316         Return the intersection of two or more domains as a new domain. As an
 317         alternative, function And() can be used.
 318         """
 319         if len(others) == 0:
 320             return self
 321         symbols = self._xsymbols((self,) + others)
 322         islset1 = self._toislset(self.polyhedra, symbols)
 323         for other in others:
 324             islset2 = other._toislset(other.polyhedra, symbols)
 325             islset1 = libisl.isl_set_intersect(islset1, islset2)
 326         return self._fromislset(islset1, symbols)
 327
 328     def __and__(self, other):
 329         return self.intersection(other)
 330     __and__.__doc__ = intersection.__doc__
 331
 332     def union(self, *others):
 333         """
 334         Return the union of two or more domains as a new domain. As an
 335         alternative, function Or() can be used.
 336         """
 337         if len(others) == 0:
 338             return self
 339         symbols = self._xsymbols((self,) + others)
 340         islset1 = self._toislset(self.polyhedra, symbols)
 341         for other in others:
 342             islset2 = other._toislset(other.polyhedra, symbols)
 343             islset1 = libisl.isl_set_union(islset1, islset2)
 344         return self._fromislset(islset1, symbols)
 345
 346     def __or__(self, other):
 347         return self.union(other)
 348     __or__.__doc__ = union.__doc__
 349
 350     def __add__(self, other):
 351         return self.union(other)
 352     __add__.__doc__ = union.__doc__
 353
 354     def difference(self, other):
 355         """
 356         Return the difference of two domains as a new domain.
 357         """
 358         symbols = self._xsymbols([self, other])
 359         islset1 = self._toislset(self.polyhedra, symbols)
 360         islset2 = other._toislset(other.polyhedra, symbols)
 361         islset = libisl.isl_set_subtract(islset1, islset2)
 362         return self._fromislset(islset, symbols)
 363
 364     def __sub__(self, other):
 365         return self.difference(other)
 366     __sub__.__doc__ = difference.__doc__
 367
 368     def lexmin(self):
 369         """
 370         Return the lexicographic minimum of the elements in the domain.
 371         """
 372         islset = self._toislset(self.polyhedra, self.symbols)
 373         islset = libisl.isl_set_lexmin(islset)
 374         return self._fromislset(islset, self.symbols)
 375
 376     def lexmax(self):
 377         """
 378         Return the lexicographic maximum of the elements in the domain.
 379         """
 380         islset = self._toislset(self.polyhedra, self.symbols)
 381         islset = libisl.isl_set_lexmax(islset)
 382         return self._fromislset(islset, self.symbols)
 383
 384     _RE_COORDINATE = re.compile(r'\((?P<num>\-?\d+)\)(/(?P<den>\d+))?')
 385
 386     def vertices(self):
 387         """
 388         Return the vertices of the domain, as a list of rational instances of
 389         Point.
 390         """
 391         from .polyhedra import Polyhedron
 392         if not self.isbounded():
 393             raise ValueError('domain must be bounded')
 394         islbset = self._toislbasicset(self.equalities, self.inequalities,
 395             self.symbols)
 396         vertices = libisl.isl_basic_set_compute_vertices(islbset);
 397         vertices = islhelper.isl_vertices_vertices(vertices)
 398         points = []
 399         for vertex in vertices:
 400             expr = libisl.isl_vertex_get_expr(vertex)
 401             coordinates = []
 402             if islhelper.isl_version < '0.13':
 403                 constraints = islhelper.isl_basic_set_constraints(expr)
 404                 for constraint in constraints:
 405                     constant = libisl.isl_constraint_get_constant_val(constraint)
 406                     constant = islhelper.isl_val_to_int(constant)
 407                     for index, symbol in enumerate(self.symbols):
 408                         coefficient = libisl.isl_constraint_get_coefficient_val(constraint,
 409                             libisl.isl_dim_set, index)
 410                         coefficient = islhelper.isl_val_to_int(coefficient)
 411                         if coefficient != 0:
 412                             coordinate = -Fraction(constant, coefficient)
 413                             coordinates.append((symbol, coordinate))
 414             else:
 415                 string = islhelper.isl_multi_aff_to_str(expr)
 416                 matches = self._RE_COORDINATE.finditer(string)
 417                 for symbol, match in zip(self.symbols, matches):
 418                     numerator = int(match.group('num'))
 419                     denominator = match.group('den')
 420                     denominator = 1 if denominator is None else int(denominator)
 421                     coordinate = Fraction(numerator, denominator)
 422                     coordinates.append((symbol, coordinate))
 423             points.append(Point(coordinates))
 424         return points
 425
 426     def points(self):
 427         """
 428         Return the integer points of a bounded domain, as a list of integer
 429         instances of Point. If the domain is not bounded, a ValueError exception
 430         is raised.
 431         """
 432         if not self.isbounded():
 433             raise ValueError('domain must be bounded')
 434         from .polyhedra import Universe, Eq
 435         islset = self._toislset(self.polyhedra, self.symbols)
 436         islpoints = islhelper.isl_set_points(islset)
 437         points = []
 438         for islpoint in islpoints:
 439             coordinates = {}
 440             for index, symbol in enumerate(self.symbols):
 441                 coordinate = libisl.isl_point_get_coordinate_val(islpoint,
 442                     libisl.isl_dim_set, index)
 443                 coordinate = islhelper.isl_val_to_int(coordinate)
 444                 coordinates[symbol] = coordinate
 445             points.append(Point(coordinates))
 446         return points
 447
 448     def __contains__(self, point):
 449         """
 450         Return True if the point if contained within the domain.
 451         """
 452         for polyhedron in self.polyhedra:
 453             if point in polyhedron:
 454                 return True
 455         return False
 456
 457     @classmethod
 458     def _polygon_inner_point(cls, points):
 459         symbols = points[0].symbols
 460         coordinates = {symbol: 0 for symbol in symbols}
 461         for point in points:
 462             for symbol, coordinate in point.coordinates():
 463                 coordinates[symbol] += coordinate
 464         for symbol in symbols:
 465             coordinates[symbol] /= len(points)
 466         return Point(coordinates)
 467
 468     @classmethod
 469     def _sort_polygon_2d(cls, points):
 470         if len(points) <= 3:
 471             return points
 472         o = cls._polygon_inner_point(points)
 473         angles = {}
 474         for m in points:
 475             om = Vector(o, m)
 476             dx, dy = (coordinate for symbol, coordinate in om.coordinates())
 477             angle = math.atan2(dy, dx)
 478             angles[m] = angle
 479         return sorted(points, key=angles.get)
 480
 481     @classmethod
 482     def _sort_polygon_3d(cls, points):
 483         if len(points) <= 3:
 484             return points
 485         o = cls._polygon_inner_point(points)
 486         a = points[0]
 487         oa = Vector(o, a)
 488         norm_oa = oa.norm()
 489         for b in points[1:]:
 490             ob = Vector(o, b)
 491             u = oa.cross(ob)
 492             if not u.isnull():
 493                 u = u.asunit()
 494                 break
 495         else:
 496             raise ValueError('degenerate polygon')
 497         angles = {a: 0.}
 498         for m in points[1:]:
 499             om = Vector(o, m)
 500             normprod = norm_oa * om.norm()
 501             cosinus = max(oa.dot(om) / normprod, -1.)
 502             sinus = u.dot(oa.cross(om)) / normprod
 503             angle = math.acos(cosinus)
 504             angle = math.copysign(angle, sinus)
 505             angles[m] = angle
 506         return sorted(points, key=angles.get)
 507
 508     def faces(self):
 509         """
 510         Return the list of faces of a bounded domain. Each face is represented
 511         by a list of vertices, in the form of rational instances of Point. If
 512         the domain is not bounded, a ValueError exception is raised.
 513         """
 514         faces = []
 515         for polyhedron in self.polyhedra:
 516             vertices = polyhedron.vertices()
 517             for constraint in polyhedron.constraints:
 518                 face = []
 519                 for vertex in vertices:
 520                     if constraint.subs(vertex.coordinates()) == 0:
 521                         face.append(vertex)
 522                 if len(face) >= 3:
 523                     faces.append(face)
 524         return faces
 525
 526     def _plot_2d(self, plot=None, **kwargs):
 527         import matplotlib.pyplot as plt
 528         from matplotlib.patches import Polygon
 529         if plot is None:
 530             fig = plt.figure()
 531             plot = fig.add_subplot(1, 1, 1)
 532         xmin, xmax = plot.get_xlim()
 533         ymin, ymax = plot.get_ylim()
 534         for polyhedron in self.polyhedra:
 535             vertices = polyhedron._sort_polygon_2d(polyhedron.vertices())
 536             xys = [tuple(vertex.values()) for vertex in vertices]
 537             xs, ys = zip(*xys)
 538             xmin, xmax = min(xmin, float(min(xs))), max(xmax, float(max(xs)))
 539             ymin, ymax = min(ymin, float(min(ys))), max(ymax, float(max(ys)))
 540             plot.add_patch(Polygon(xys, closed=True, **kwargs))
 541         plot.set_xlim(xmin, xmax)
 542         plot.set_ylim(ymin, ymax)
 543         return plot
 544
 545     def _plot_3d(self, plot=None, **kwargs):
 546         import matplotlib.pyplot as plt
 547         from mpl_toolkits.mplot3d import Axes3D
 548         from mpl_toolkits.mplot3d.art3d import Poly3DCollection
 549         if plot is None:
 550             fig = plt.figure()
 551             axes = Axes3D(fig)
 552         else:
 553             axes = plot
 554         xmin, xmax = axes.get_xlim()
 555         ymin, ymax = axes.get_ylim()
 556         zmin, zmax = axes.get_zlim()
 557         poly_xyzs = []
 558         for vertices in self.faces():
 559             vertices = self._sort_polygon_3d(vertices)
 560             vertices.append(vertices[0])
 561             face_xyzs = [tuple(vertex.values()) for vertex in vertices]
 562             xs, ys, zs = zip(*face_xyzs)
 563             xmin, xmax = min(xmin, float(min(xs))), max(xmax, float(max(xs)))
 564             ymin, ymax = min(ymin, float(min(ys))), max(ymax, float(max(ys)))
 565             zmin, zmax = min(zmin, float(min(zs))), max(zmax, float(max(zs)))
 566             poly_xyzs.append(face_xyzs)
 567         collection = Poly3DCollection(poly_xyzs, **kwargs)
 568         axes.add_collection3d(collection)
 569         axes.set_xlim(xmin, xmax)
 570         axes.set_ylim(ymin, ymax)
 571         axes.set_zlim(zmin, zmax)
 572         return axes
 573
 574     def plot(self, plot=None, **kwargs):
 575         """
 576         Plot a 2D or 3D domain using matplotlib. Draw it to the current plot
 577         object if present, otherwise create a new one. options are keyword
 578         arguments passed to the matplotlib drawing functions, they can be used
 579         to set the drawing color for example. Raise ValueError is the domain is
 580         not 2D or 3D.
 581         """
 582         if not self.isbounded():
 583             raise ValueError('domain must be bounded')
 584         elif self.dimension == 2:
 585             return self._plot_2d(plot=plot, **kwargs)
 586         elif self.dimension == 3:
 587             return self._plot_3d(plot=plot, **kwargs)
 588         else:
 589             raise ValueError('polyhedron must be 2 or 3-dimensional')
 590
 591     def subs(self, symbol, expression=None):
 592         """
 593         Substitute the given symbol by an expression in the domain constraints.
 594         To perform multiple substitutions at once, pass a sequence or a
 595         dictionary of (old, new) pairs to subs. The syntax of this function is
 596         similar to LinExpr.subs().
 597         """
 598         polyhedra = [polyhedron.subs(symbol, expression)
 599             for polyhedron in self.polyhedra]
 600         return Domain(*polyhedra)
 601
 602     @classmethod
 603     def _fromislset(cls, islset, symbols):
 604         from .polyhedra import Polyhedron
 605         islset = libisl.isl_set_remove_divs(islset)
 606         islbsets = islhelper.isl_set_basic_sets(islset)
 607         libisl.isl_set_free(islset)
 608         polyhedra = []
 609         for islbset in islbsets:
 610             polyhedron = Polyhedron._fromislbasicset(islbset, symbols)
 611             polyhedra.append(polyhedron)
 612         if len(polyhedra) == 0:
 613             from .polyhedra import Empty
 614             return Empty
 615         elif len(polyhedra) == 1:
 616             return polyhedra[0]
 617         else:
 618             self = object().__new__(Domain)
 619             self._polyhedra = tuple(polyhedra)
 620             self._symbols = cls._xsymbols(polyhedra)
 621             self._dimension = len(self._symbols)
 622             return self
 623
 624     @classmethod
 625     def _toislset(cls, polyhedra, symbols):
 626         polyhedron = polyhedra[0]
 627         islbset = polyhedron._toislbasicset(polyhedron.equalities,
 628             polyhedron.inequalities, symbols)
 629         islset1 = libisl.isl_set_from_basic_set(islbset)
 630         for polyhedron in polyhedra[1:]:
 631             islbset = polyhedron._toislbasicset(polyhedron.equalities,
 632                 polyhedron.inequalities, symbols)
 633             islset2 = libisl.isl_set_from_basic_set(islbset)
 634             islset1 = libisl.isl_set_union(islset1, islset2)
 635         return islset1
 636
 637     @classmethod
 638     def _fromast(cls, node):
 639         from .polyhedra import Polyhedron
 640         if isinstance(node, ast.Module) and len(node.body) == 1:
 641             return cls._fromast(node.body[0])
 642         elif isinstance(node, ast.Expr):
 643             return cls._fromast(node.value)
 644         elif isinstance(node, ast.UnaryOp):
 645             domain = cls._fromast(node.operand)
 646             if isinstance(node.operand, ast.invert):
 647                 return Not(domain)
 648         elif isinstance(node, ast.BinOp):
 649             domain1 = cls._fromast(node.left)
 650             domain2 = cls._fromast(node.right)
 651             if isinstance(node.op, ast.BitAnd):
 652                 return And(domain1, domain2)
 653             elif isinstance(node.op, ast.BitOr):
 654                 return Or(domain1, domain2)
 655         elif isinstance(node, ast.Compare):
 656             equalities = []
 657             inequalities = []
 658             left = LinExpr._fromast(node.left)
 659             for i in range(len(node.ops)):
 660                 op = node.ops[i]
 661                 right = LinExpr._fromast(node.comparators[i])
 662                 if isinstance(op, ast.Lt):
 663                     inequalities.append(right - left - 1)
 664                 elif isinstance(op, ast.LtE):
 665                     inequalities.append(right - left)
 666                 elif isinstance(op, ast.Eq):
 667                     equalities.append(left - right)
 668                 elif isinstance(op, ast.GtE):
 669                     inequalities.append(left - right)
 670                 elif isinstance(op, ast.Gt):
 671                     inequalities.append(left - right - 1)
 672                 else:
 673                     break
 674                 left = right
 675             else:
 676                 return Polyhedron(equalities, inequalities)
 677         raise SyntaxError('invalid syntax')
 678
 679     _RE_BRACES = re.compile(r'^\{\s*|\s*\}$')
 680     _RE_EQ = re.compile(r'([^<=>])=([^<=>])')
 681     _RE_AND = re.compile(r'\band\b|,|&&|/\\|∧|∩')
 682     _RE_OR = re.compile(r'\bor\b|;|\|\||\\/|∨|∪')
 683     _RE_NOT = re.compile(r'\bnot\b|!|¬')
 684     _RE_NUM_VAR = LinExpr._RE_NUM_VAR
 685     _RE_OPERATORS = re.compile(r'(&|\||~)')
 686
 687     @classmethod
 688     def fromstring(cls, string):
 689         """
 690         Create a domain from a string. Raise SyntaxError if the string is not
 691         properly formatted.
 692         """
 693         # remove curly brackets
 694         string = cls._RE_BRACES.sub(r'', string)
 695         # replace '=' by '=='
 696         string = cls._RE_EQ.sub(r'\1==\2', string)
 697         # replace 'and', 'or', 'not'
 698         string = cls._RE_AND.sub(r' & ', string)
 699         string = cls._RE_OR.sub(r' | ', string)
 700         string = cls._RE_NOT.sub(r' ~', string)
 701         # add implicit multiplication operators, e.g. '5x' -> '5*x'
 702         string = cls._RE_NUM_VAR.sub(r'\1*\2', string)
 703         # add parentheses to force precedence
 704         tokens = cls._RE_OPERATORS.split(string)
 705         for i, token in enumerate(tokens):
 706             if i % 2 == 0:
 707                 token = '({})'.format(token)
 708                 tokens[i] = token
 709         string = ''.join(tokens)
 710         tree = ast.parse(string, 'eval')
 711         return cls._fromast(tree)
 712
 713     def __repr__(self):
 714         assert len(self.polyhedra) >= 2
 715         strings = [repr(polyhedron) for polyhedron in self.polyhedra]
 716         return 'Or({})'.format(', '.join(strings))
 717
 718     def _repr_latex_(self):
 719         strings = []
 720         for polyhedron in self.polyhedra:
 721             strings.append('({})'.format(polyhedron._repr_latex_().strip('$')))
 722         return '${}$'.format(' \\vee '.join(strings))
 723
 724     @classmethod
 725     def fromsympy(cls, expr):
 726         """
 727         Create a domain from a sympy expression.
 728         """
 729         import sympy
 730         from .polyhedra import Lt, Le, Eq, Ne, Ge, Gt
 731         funcmap = {
 732             sympy.And: And, sympy.Or: Or, sympy.Not: Not,
 733             sympy.Lt: Lt, sympy.Le: Le,
 734             sympy.Eq: Eq, sympy.Ne: Ne,
 735             sympy.Ge: Ge, sympy.Gt: Gt,
 736         }
 737         if expr.func in funcmap:
 738             args = [Domain.fromsympy(arg) for arg in expr.args]
 739             return funcmap[expr.func](*args)
 740         elif isinstance(expr, sympy.Expr):
 741             return LinExpr.fromsympy(expr)
 742         raise ValueError('non-domain expression: {!r}'.format(expr))
 743
 744     def tosympy(self):
 745         """
 746         Convert the domain to a sympy expression.
 747         """
 748         import sympy
 749         polyhedra = [polyhedron.tosympy() for polyhedron in polyhedra]
 750         return sympy.Or(*polyhedra)
 751
 752
 753 def And(*domains):
 754     """
 755     Create the intersection domain of the domains given in arguments.
 756     """
 757     if len(domains) == 0:
 758         from .polyhedra import Universe
 759         return Universe
 760     else:
 761         return domains[0].intersection(*domains[1:])
 762 And.__doc__ = Domain.intersection.__doc__
 763
 764 def Or(*domains):
 765     """
 766     Create the union domain of the domains given in arguments.
 767     """
 768     if len(domains) == 0:
 769         from .polyhedra import Empty
 770         return Empty
 771     else:
 772         return domains[0].union(*domains[1:])
 773 Or.__doc__ = Domain.union.__doc__
 774
 775 def Not(domain):
 776     """
 777     Create the complementary domain of the domain given in argument.
 778     """
 779     return ~domain
 780 Not.__doc__ = Domain.complement.__doc__