linpy/domains.py

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