X-Git-Url: https://scm.cri.mines-paristech.fr/git/linpy.git/blobdiff_plain/50f0fde4dca3f10e535314b1c9d325cd9283735b..4409cffcddd60c95618ad19908fdc4f9097e8caa:/linpy/domains.py diff --git a/linpy/domains.py b/linpy/domains.py index 70e5e4f..f26b2fc 100644 --- a/linpy/domains.py +++ b/linpy/domains.py @@ -15,13 +15,6 @@ # You should have received a copy of the GNU General Public License # along with LinPy. If not, see . -""" -Polyhedral domains - -This module provides classes and functions to deal with polyhedral -domains, i.e. unions of polyhedra. -""" - import ast import functools import re @@ -31,7 +24,7 @@ from fractions import Fraction from . import islhelper from .islhelper import mainctx, libisl -from .linexprs import Expression, Symbol, Rational +from .linexprs import LinExpr, Symbol from .geometry import GeometricObject, Point, Vector @@ -44,7 +37,11 @@ __all__ = [ @functools.total_ordering class Domain(GeometricObject): """ - This class represents polyhedral domains, i.e. unions of polyhedra. + A domain is a union of polyhedra. Unlike polyhedra, domains allow exact + computation of union, subtraction and complementary operations. + + A domain with a unique polyhedron is automatically subclassed as a + Polyhedron instance. """ __slots__ = ( @@ -55,7 +52,28 @@ class Domain(GeometricObject): def __new__(cls, *polyhedra): """ - Create and return a new domain from a string or a list of polyhedra. + Return a domain from a sequence of polyhedra. + + >>> square1 = Polyhedron('0 <= x <= 2, 0 <= y <= 2') + >>> square2 = Polyhedron('1 <= x <= 3, 1 <= y <= 3') + >>> dom = Domain(square1, square2) + >>> dom + Or(And(x <= 2, 0 <= x, y <= 2, 0 <= y), + And(x <= 3, 1 <= x, y <= 3, 1 <= y)) + + It is also possible to build domains from polyhedra using arithmetic + operators Domain.__or__(), Domain.__invert__() or functions Or() and + Not(), using one of the following instructions: + + >>> dom = square1 | square2 + >>> dom = Or(square1, square2) + + Alternatively, a domain can be built from a string: + + >>> dom = Domain('0 <= x <= 2, 0 <= y <= 2; 1 <= x <= 3, 1 <= y <= 3') + + Finally, a domain can be built from a GeometricObject instance, calling + the GeometricObject.asdomain() method. """ from .polyhedra import Polyhedron if len(polyhedra) == 1: @@ -88,35 +106,28 @@ class Domain(GeometricObject): @property def polyhedra(self): """ - The tuple of polyhedra which constitute the domain. + The tuple of polyhedra present in the domain. """ return self._polyhedra @property def symbols(self): """ - The tuple of symbols present in the domain equations. + The tuple of symbols present in the domain equations, sorted according + to Symbol.sortkey(). """ return self._symbols @property def dimension(self): """ - The dimension of the domain, i.e. the number of symbols. + The dimension of the domain, i.e. the number of symbols present in it. """ return self._dimension - def make_disjoint(self): - """ - Return an equivalent domain, whose polyhedra are disjoint. - """ - islset = self._toislset(self.polyhedra, self.symbols) - islset = libisl.isl_set_make_disjoint(mainctx, islset) - return self._fromislset(islset, self.symbols) - def isempty(self): """ - Return True if the domain is empty. + Return True if the domain is empty, that is, equal to Empty. """ islset = self._toislset(self.polyhedra, self.symbols) empty = bool(libisl.isl_set_is_empty(islset)) @@ -131,7 +142,7 @@ class Domain(GeometricObject): def isuniverse(self): """ - Return True if the domain is universal, i.e. with no constraint. + Return True if the domain is universal, that is, equal to Universe. """ islset = self._toislset(self.polyhedra, self.symbols) universe = bool(libisl.isl_set_plain_is_universe(islset)) @@ -149,20 +160,24 @@ class Domain(GeometricObject): def __eq__(self, other): """ - Return True if the two domains are equal. + Return True if two domains are equal. """ - symbols = self._xsymbols([self, other]) - islset1 = self._toislset(self.polyhedra, symbols) - islset2 = other._toislset(other.polyhedra, symbols) - equal = bool(libisl.isl_set_is_equal(islset1, islset2)) - libisl.isl_set_free(islset1) - libisl.isl_set_free(islset2) - return equal + if isinstance(other, Domain): + symbols = self._xsymbols([self, other]) + islset1 = self._toislset(self.polyhedra, symbols) + islset2 = other._toislset(other.polyhedra, symbols) + equal = bool(libisl.isl_set_is_equal(islset1, islset2)) + libisl.isl_set_free(islset1) + libisl.isl_set_free(islset2) + return equal + return NotImplemented def isdisjoint(self, other): """ Return True if two domains have a null intersection. """ + if not isinstance(other, Domain): + raise TypeError('other must be a Domain instance') symbols = self._xsymbols([self, other]) islset1 = self._toislset(self.polyhedra, symbols) islset2 = self._toislset(other.polyhedra, symbols) @@ -175,29 +190,33 @@ class Domain(GeometricObject): """ Report whether another domain contains the domain. """ - symbols = self._xsymbols([self, other]) - islset1 = self._toislset(self.polyhedra, symbols) - islset2 = self._toislset(other.polyhedra, symbols) - equal = bool(libisl.isl_set_is_subset(islset1, islset2)) - libisl.isl_set_free(islset1) - libisl.isl_set_free(islset2) - return equal + return self < other def __le__(self, other): - return self.issubset(other) + if isinstance(other, Domain): + symbols = self._xsymbols([self, other]) + islset1 = self._toislset(self.polyhedra, symbols) + islset2 = self._toislset(other.polyhedra, symbols) + equal = bool(libisl.isl_set_is_subset(islset1, islset2)) + libisl.isl_set_free(islset1) + libisl.isl_set_free(islset2) + return equal + return NotImplemented __le__.__doc__ = issubset.__doc__ def __lt__(self, other): """ Report whether another domain is contained within the domain. """ - symbols = self._xsymbols([self, other]) - islset1 = self._toislset(self.polyhedra, symbols) - islset2 = self._toislset(other.polyhedra, symbols) - equal = bool(libisl.isl_set_is_strict_subset(islset1, islset2)) - libisl.isl_set_free(islset1) - libisl.isl_set_free(islset2) - return equal + if isinstance(other, Domain): + symbols = self._xsymbols([self, other]) + islset1 = self._toislset(self.polyhedra, symbols) + islset2 = self._toislset(other.polyhedra, symbols) + equal = bool(libisl.isl_set_is_strict_subset(islset1, islset2)) + libisl.isl_set_free(islset1) + libisl.isl_set_free(islset2) + return equal + return NotImplemented def complement(self): """ @@ -211,10 +230,18 @@ class Domain(GeometricObject): return self.complement() __invert__.__doc__ = complement.__doc__ + def make_disjoint(self): + """ + Return an equivalent domain, whose polyhedra are disjoint. + """ + islset = self._toislset(self.polyhedra, self.symbols) + islset = libisl.isl_set_make_disjoint(islset) + return self._fromislset(islset, self.symbols) + def coalesce(self): """ Simplify the representation of the domain by trying to combine pairs of - polyhedra into a single polyhedron. + polyhedra into a single polyhedron, and return the resulting domain. """ islset = self._toislset(self.polyhedra, self.symbols) islset = libisl.isl_set_coalesce(islset) @@ -223,7 +250,7 @@ class Domain(GeometricObject): def detect_equalities(self): """ Simplify the representation of the domain by detecting implicit - equalities. + equalities, and return the resulting domain. """ islset = self._toislset(self.polyhedra, self.symbols) islset = libisl.isl_set_detect_equalities(islset) @@ -231,16 +258,14 @@ class Domain(GeometricObject): def remove_redundancies(self): """ - Remove redundant constraints in the domain. + Remove redundant constraints in the domain, and return the resulting + domain. """ islset = self._toislset(self.polyhedra, self.symbols) islset = libisl.isl_set_remove_redundancies(islset) return self._fromislset(islset, self.symbols) def aspolyhedron(self): - """ - Return the polyhedral hull of the domain. - """ from .polyhedra import Polyhedron islset = self._toislset(self.polyhedra, self.symbols) islbset = libisl.isl_set_polyhedral_hull(islset) @@ -251,8 +276,13 @@ class Domain(GeometricObject): def project(self, symbols): """ - Project out the symbols given in arguments. + Project out the sequence of symbols given in arguments, and return the + resulting domain. """ + symbols = list(symbols) + for symbol in symbols: + if not isinstance(symbol, Symbol): + raise TypeError('symbols must be Symbol instances') islset = self._toislset(self.polyhedra, self.symbols) n = 0 for index, symbol in reversed(list(enumerate(self.symbols))): @@ -269,7 +299,8 @@ class Domain(GeometricObject): def sample(self): """ - Return a sample of the domain. + Return a sample of the domain, as an integer instance of Point. If the + domain is empty, a ValueError exception is raised. """ islset = self._toislset(self.polyhedra, self.symbols) islpoint = libisl.isl_set_sample_point(islset) @@ -287,54 +318,62 @@ class Domain(GeometricObject): def intersection(self, *others): """ - Return the intersection of two domains as a new domain. + Return the intersection of two or more domains as a new domain. As an + alternative, function And() can be used. """ - if len(others) == 0: - return self - symbols = self._xsymbols((self,) + others) - islset1 = self._toislset(self.polyhedra, symbols) + result = self for other in others: - islset2 = other._toislset(other.polyhedra, symbols) - islset1 = libisl.isl_set_intersect(islset1, islset2) - return self._fromislset(islset1, symbols) + result &= other + return result def __and__(self, other): - return self.intersection(other) + if isinstance(other, Domain): + symbols = self._xsymbols([self, other]) + islset1 = self._toislset(self.polyhedra, symbols) + islset2 = other._toislset(other.polyhedra, symbols) + islset = libisl.isl_set_intersect(islset1, islset2) + return self._fromislset(islset, symbols) + return NotImplemented __and__.__doc__ = intersection.__doc__ def union(self, *others): """ - Return the union of two domains as a new domain. + Return the union of two or more domains as a new domain. As an + alternative, function Or() can be used. """ - if len(others) == 0: - return self - symbols = self._xsymbols((self,) + others) - islset1 = self._toislset(self.polyhedra, symbols) + result = self for other in others: - islset2 = other._toislset(other.polyhedra, symbols) - islset1 = libisl.isl_set_union(islset1, islset2) - return self._fromislset(islset1, symbols) + result |= other + return result def __or__(self, other): - return self.union(other) + if isinstance(other, Domain): + symbols = self._xsymbols([self, other]) + islset1 = self._toislset(self.polyhedra, symbols) + islset2 = other._toislset(other.polyhedra, symbols) + islset = libisl.isl_set_union(islset1, islset2) + return self._fromislset(islset, symbols) + return NotImplemented __or__.__doc__ = union.__doc__ def __add__(self, other): - return self.union(other) + return self | other __add__.__doc__ = union.__doc__ def difference(self, other): """ Return the difference of two domains as a new domain. """ - symbols = self._xsymbols([self, other]) - islset1 = self._toislset(self.polyhedra, symbols) - islset2 = other._toislset(other.polyhedra, symbols) - islset = libisl.isl_set_subtract(islset1, islset2) - return self._fromislset(islset, symbols) + return self - other def __sub__(self, other): - return self.difference(other) + if isinstance(other, Domain): + symbols = self._xsymbols([self, other]) + islset1 = self._toislset(self.polyhedra, symbols) + islset2 = other._toislset(other.polyhedra, symbols) + islset = libisl.isl_set_subtract(islset1, islset2) + return self._fromislset(islset, symbols) + return NotImplemented __sub__.__doc__ = difference.__doc__ def lexmin(self): @@ -353,11 +392,15 @@ class Domain(GeometricObject): islset = libisl.isl_set_lexmax(islset) return self._fromislset(islset, self.symbols) - _RE_COORDINATE = re.compile(r'\((?P\-?\d+)\)(/(?P\d+))?') + if islhelper.isl_version >= '0.13': + _RE_COORDINATE = re.compile(r'\((?P\-?\d+)\)(/(?P\d+))?') + else: + _RE_COORDINATE = None def vertices(self): """ - Return the vertices of the domain. + Return the vertices of the domain, as a list of rational instances of + Point. """ from .polyhedra import Polyhedron if not self.isbounded(): @@ -368,10 +411,10 @@ class Domain(GeometricObject): vertices = islhelper.isl_vertices_vertices(vertices) points = [] for vertex in vertices: - expr = libisl.isl_vertex_get_expr(vertex) + expression = libisl.isl_vertex_get_expr(vertex) coordinates = [] - if islhelper.isl_version < '0.13': - constraints = islhelper.isl_basic_set_constraints(expr) + if self._RE_COORDINATE is None: + constraints = islhelper.isl_basic_set_constraints(expression) for constraint in constraints: constant = libisl.isl_constraint_get_constant_val(constraint) constant = islhelper.isl_val_to_int(constant) @@ -383,7 +426,7 @@ class Domain(GeometricObject): coordinate = -Fraction(constant, coefficient) coordinates.append((symbol, coordinate)) else: - string = islhelper.isl_multi_aff_to_str(expr) + string = islhelper.isl_multi_aff_to_str(expression) matches = self._RE_COORDINATE.finditer(string) for symbol, match in zip(self.symbols, matches): numerator = int(match.group('num')) @@ -396,7 +439,9 @@ class Domain(GeometricObject): def points(self): """ - Return the points with integer coordinates contained in the domain. + Return the integer points of a bounded domain, as a list of integer + instances of Point. If the domain is not bounded, a ValueError exception + is raised. """ if not self.isbounded(): raise ValueError('domain must be bounded') @@ -414,6 +459,15 @@ class Domain(GeometricObject): points.append(Point(coordinates)) return points + def __contains__(self, point): + """ + Return True if the point if contained within the domain. + """ + for polyhedron in self.polyhedra: + if point in polyhedron: + return True + return False + @classmethod def _polygon_inner_point(cls, points): symbols = points[0].symbols @@ -467,7 +521,9 @@ class Domain(GeometricObject): def faces(self): """ - Return the vertices of the domain, grouped by face. + Return the list of faces of a bounded domain. Each face is represented + by a list of vertices, in the form of rational instances of Point. If + the domain is not bounded, a ValueError exception is raised. """ faces = [] for polyhedron in self.polyhedra: @@ -531,7 +587,11 @@ class Domain(GeometricObject): def plot(self, plot=None, **kwargs): """ - Plot the domain using matplotlib. + Plot a 2D or 3D domain using matplotlib. Draw it to the current plot + object if present, otherwise create a new one. options are keyword + arguments passed to the matplotlib drawing functions, they can be used + to set the drawing color for example. Raise ValueError is the domain is + not 2D or 3D. """ if not self.isbounded(): raise ValueError('domain must be bounded') @@ -540,21 +600,14 @@ class Domain(GeometricObject): elif self.dimension == 3: return self._plot_3d(plot=plot, **kwargs) else: - raise ValueError('polyhedron must be 2 or 3-dimensional') - - def __contains__(self, point): - """ - Return True if point if contained within the domain. - """ - for polyhedron in self.polyhedra: - if point in polyhedron: - return True - return False + raise ValueError('domain must be 2 or 3-dimensional') def subs(self, symbol, expression=None): """ - Subsitute symbol by expression in equations and return the resulting - domain. + Substitute the given symbol by an expression in the domain constraints. + To perform multiple substitutions at once, pass a sequence or a + dictionary of (old, new) pairs to subs. The syntax of this function is + similar to LinExpr.subs(). """ polyhedra = [polyhedron.subs(symbol, expression) for polyhedron in self.polyhedra] @@ -616,10 +669,10 @@ class Domain(GeometricObject): elif isinstance(node, ast.Compare): equalities = [] inequalities = [] - left = Expression._fromast(node.left) + left = LinExpr._fromast(node.left) for i in range(len(node.ops)): op = node.ops[i] - right = Expression._fromast(node.comparators[i]) + right = LinExpr._fromast(node.comparators[i]) if isinstance(op, ast.Lt): inequalities.append(right - left - 1) elif isinstance(op, ast.LtE): @@ -642,25 +695,26 @@ class Domain(GeometricObject): _RE_AND = re.compile(r'\band\b|,|&&|/\\|∧|∩') _RE_OR = re.compile(r'\bor\b|;|\|\||\\/|∨|∪') _RE_NOT = re.compile(r'\bnot\b|!|¬') - _RE_NUM_VAR = Expression._RE_NUM_VAR + _RE_NUM_VAR = LinExpr._RE_NUM_VAR _RE_OPERATORS = re.compile(r'(&|\||~)') @classmethod def fromstring(cls, string): """ - Convert a string into a domain. + Create a domain from a string. Raise SyntaxError if the string is not + properly formatted. """ - # remove curly brackets + # Remove curly brackets. string = cls._RE_BRACES.sub(r'', string) - # replace '=' by '==' + # Replace '=' by '=='. string = cls._RE_EQ.sub(r'\1==\2', string) - # replace 'and', 'or', 'not' + # Replace 'and', 'or', 'not'. string = cls._RE_AND.sub(r' & ', string) string = cls._RE_OR.sub(r' | ', string) string = cls._RE_NOT.sub(r' ~', string) - # add implicit multiplication operators, e.g. '5x' -> '5*x' + # Add implicit multiplication operators, e.g. '5x' -> '5*x'. string = cls._RE_NUM_VAR.sub(r'\1*\2', string) - # add parentheses to force precedence + # Add parentheses to force precedence. tokens = cls._RE_OPERATORS.split(string) for i, token in enumerate(tokens): if i % 2 == 0: @@ -671,23 +725,14 @@ class Domain(GeometricObject): return cls._fromast(tree) def __repr__(self): - """ - Return repr(self). - """ assert len(self.polyhedra) >= 2 strings = [repr(polyhedron) for polyhedron in self.polyhedra] return 'Or({})'.format(', '.join(strings)) - def _repr_latex_(self): - strings = [] - for polyhedron in self.polyhedra: - strings.append('({})'.format(polyhedron._repr_latex_().strip('$'))) - return '${}$'.format(' \\vee '.join(strings)) - @classmethod - def fromsympy(cls, expr): + def fromsympy(cls, expression): """ - Convert a SymPy expression into a domain. + Create a domain from a SymPy expression. """ import sympy from .polyhedra import Lt, Le, Eq, Ne, Ge, Gt @@ -697,16 +742,16 @@ class Domain(GeometricObject): sympy.Eq: Eq, sympy.Ne: Ne, sympy.Ge: Ge, sympy.Gt: Gt, } - if expr.func in funcmap: - args = [Domain.fromsympy(arg) for arg in expr.args] - return funcmap[expr.func](*args) - elif isinstance(expr, sympy.Expr): - return Expression.fromsympy(expr) - raise ValueError('non-domain expression: {!r}'.format(expr)) + if expression.func in funcmap: + args = [Domain.fromsympy(arg) for arg in expression.args] + return funcmap[expression.func](*args) + elif isinstance(expression, sympy.Expr): + return LinExpr.fromsympy(expression) + raise ValueError('non-domain expression: {!r}'.format(expression)) def tosympy(self): """ - Convert a domain into a SymPy expression. + Convert the domain to a SymPy expression. """ import sympy polyhedra = [polyhedron.tosympy() for polyhedron in polyhedra] @@ -714,21 +759,27 @@ class Domain(GeometricObject): def And(*domains): + """ + Create the intersection domain of the domains given in arguments. + """ if len(domains) == 0: from .polyhedra import Universe return Universe else: return domains[0].intersection(*domains[1:]) -And.__doc__ = Domain.intersection.__doc__ def Or(*domains): + """ + Create the union domain of the domains given in arguments. + """ if len(domains) == 0: from .polyhedra import Empty return Empty else: return domains[0].union(*domains[1:]) -Or.__doc__ = Domain.union.__doc__ def Not(domain): + """ + Create the complementary domain of the domain given in argument. + """ return ~domain -Not.__doc__ = Domain.complement.__doc__