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