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