linpy/domains.py

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