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