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