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