X-Git-Url: https://scm.cri.ensmp.fr/git/linpy.git/blobdiff_plain/2a1055d4f4754fa33c53d3f155cc050e4911a2a3..7b93cea1daf2889e9ee10ca9c22a1b5124404937:/pypol/polyhedra.py diff --git a/pypol/polyhedra.py b/pypol/polyhedra.py deleted file mode 100644 index a5048d1..0000000 --- a/pypol/polyhedra.py +++ /dev/null @@ -1,417 +0,0 @@ -import functools -import math -import numbers - -from . import islhelper - -from .islhelper import mainctx, libisl -from .geometry import GeometricObject, Point, Vector -from .linexprs import Expression, Symbol, Rational -from .domains import Domain - - -__all__ = [ - 'Polyhedron', - 'Lt', 'Le', 'Eq', 'Ne', 'Ge', 'Gt', - 'Empty', 'Universe', -] - - -class Polyhedron(Domain): - - __slots__ = ( - '_equalities', - '_inequalities', - '_constraints', - '_symbols', - '_dimension', - ) - - def __new__(cls, equalities=None, inequalities=None): - if isinstance(equalities, str): - if inequalities is not None: - raise TypeError('too many arguments') - return cls.fromstring(equalities) - elif isinstance(equalities, GeometricObject): - if inequalities is not None: - raise TypeError('too many arguments') - return equalities.aspolyhedron() - if equalities is None: - equalities = [] - else: - for i, equality in enumerate(equalities): - if not isinstance(equality, Expression): - raise TypeError('equalities must be linear expressions') - equalities[i] = equality.scaleint() - if inequalities is None: - inequalities = [] - else: - for i, inequality in enumerate(inequalities): - if not isinstance(inequality, Expression): - raise TypeError('inequalities must be linear expressions') - inequalities[i] = inequality.scaleint() - symbols = cls._xsymbols(equalities + inequalities) - islbset = cls._toislbasicset(equalities, inequalities, symbols) - return cls._fromislbasicset(islbset, symbols) - - @property - def equalities(self): - return self._equalities - - @property - def inequalities(self): - return self._inequalities - - @property - def constraints(self): - return self._constraints - - @property - def polyhedra(self): - return self, - - def disjoint(self): - """ - Return this set as disjoint. - """ - return self - - def isuniverse(self): - """ - Return true if this set is the Universe set. - """ - islbset = self._toislbasicset(self.equalities, self.inequalities, - self.symbols) - universe = bool(libisl.isl_basic_set_is_universe(islbset)) - libisl.isl_basic_set_free(islbset) - return universe - - def aspolyhedron(self): - """ - Return polyhedral hull of this set. - """ - return self - - def __contains__(self, point): - if not isinstance(point, Point): - raise TypeError('point must be a Point instance') - if self.symbols != point.symbols: - raise ValueError('arguments must belong to the same space') - for equality in self.equalities: - if equality.subs(point.coordinates()) != 0: - return False - for inequality in self.inequalities: - if inequality.subs(point.coordinates()) < 0: - return False - return True - - def subs(self, symbol, expression=None): - equalities = [equality.subs(symbol, expression) - for equality in self.equalities] - inequalities = [inequality.subs(symbol, expression) - for inequality in self.inequalities] - return Polyhedron(equalities, inequalities) - - @classmethod - def _fromislbasicset(cls, islbset, symbols): - islconstraints = islhelper.isl_basic_set_constraints(islbset) - equalities = [] - inequalities = [] - for islconstraint in islconstraints: - constant = libisl.isl_constraint_get_constant_val(islconstraint) - constant = islhelper.isl_val_to_int(constant) - coefficients = {} - for index, symbol in enumerate(symbols): - coefficient = libisl.isl_constraint_get_coefficient_val(islconstraint, - libisl.isl_dim_set, index) - coefficient = islhelper.isl_val_to_int(coefficient) - if coefficient != 0: - coefficients[symbol] = coefficient - expression = Expression(coefficients, constant) - if libisl.isl_constraint_is_equality(islconstraint): - equalities.append(expression) - else: - inequalities.append(expression) - libisl.isl_basic_set_free(islbset) - self = object().__new__(Polyhedron) - self._equalities = tuple(equalities) - self._inequalities = tuple(inequalities) - self._constraints = tuple(equalities + inequalities) - self._symbols = cls._xsymbols(self._constraints) - self._dimension = len(self._symbols) - return self - - @classmethod - def _toislbasicset(cls, equalities, inequalities, symbols): - dimension = len(symbols) - indices = {symbol: index for index, symbol in enumerate(symbols)} - islsp = libisl.isl_space_set_alloc(mainctx, 0, dimension) - islbset = libisl.isl_basic_set_universe(libisl.isl_space_copy(islsp)) - islls = libisl.isl_local_space_from_space(islsp) - for equality in equalities: - isleq = libisl.isl_equality_alloc(libisl.isl_local_space_copy(islls)) - for symbol, coefficient in equality.coefficients(): - islval = str(coefficient).encode() - islval = libisl.isl_val_read_from_str(mainctx, islval) - index = indices[symbol] - isleq = libisl.isl_constraint_set_coefficient_val(isleq, - libisl.isl_dim_set, index, islval) - if equality.constant != 0: - islval = str(equality.constant).encode() - islval = libisl.isl_val_read_from_str(mainctx, islval) - isleq = libisl.isl_constraint_set_constant_val(isleq, islval) - islbset = libisl.isl_basic_set_add_constraint(islbset, isleq) - for inequality in inequalities: - islin = libisl.isl_inequality_alloc(libisl.isl_local_space_copy(islls)) - for symbol, coefficient in inequality.coefficients(): - islval = str(coefficient).encode() - islval = libisl.isl_val_read_from_str(mainctx, islval) - index = indices[symbol] - islin = libisl.isl_constraint_set_coefficient_val(islin, - libisl.isl_dim_set, index, islval) - if inequality.constant != 0: - islval = str(inequality.constant).encode() - islval = libisl.isl_val_read_from_str(mainctx, islval) - islin = libisl.isl_constraint_set_constant_val(islin, islval) - islbset = libisl.isl_basic_set_add_constraint(islbset, islin) - return islbset - - @classmethod - def fromstring(cls, string): - domain = Domain.fromstring(string) - if not isinstance(domain, Polyhedron): - raise ValueError('non-polyhedral expression: {!r}'.format(string)) - return domain - - def __repr__(self): - if self.isempty(): - return 'Empty' - elif self.isuniverse(): - return 'Universe' - else: - strings = [] - for equality in self.equalities: - strings.append('Eq({}, 0)'.format(equality)) - for inequality in self.inequalities: - strings.append('Ge({}, 0)'.format(inequality)) - if len(strings) == 1: - return strings[0] - else: - return 'And({})'.format(', '.join(strings)) - - def _repr_latex_(self): - if self.isempty(): - return '$\\emptyset$' - elif self.isuniverse(): - return '$\\Omega$' - else: - strings = [] - for equality in self.equalities: - strings.append('{} = 0'.format(equality._repr_latex_().strip('$'))) - for inequality in self.inequalities: - strings.append('{} \\ge 0'.format(inequality._repr_latex_().strip('$'))) - return '${}$'.format(' \\wedge '.join(strings)) - - @classmethod - def fromsympy(cls, expr): - domain = Domain.fromsympy(expr) - if not isinstance(domain, Polyhedron): - raise ValueError('non-polyhedral expression: {!r}'.format(expr)) - return domain - - def tosympy(self): - import sympy - constraints = [] - for equality in self.equalities: - constraints.append(sympy.Eq(equality.tosympy(), 0)) - for inequality in self.inequalities: - constraints.append(sympy.Ge(inequality.tosympy(), 0)) - return sympy.And(*constraints) - - @classmethod - def _polygon_inner_point(cls, points): - symbols = points[0].symbols - coordinates = {symbol: 0 for symbol in symbols} - for point in points: - for symbol, coordinate in point.coordinates(): - coordinates[symbol] += coordinate - for symbol in symbols: - coordinates[symbol] /= len(points) - return Point(coordinates) - - @classmethod - def _sort_polygon_2d(cls, points): - if len(points) <= 3: - return points - o = cls._polygon_inner_point(points) - angles = {} - for m in points: - om = Vector(o, m) - dx, dy = (coordinate for symbol, coordinate in om.coordinates()) - angle = math.atan2(dy, dx) - angles[m] = angle - return sorted(points, key=angles.get) - - @classmethod - def _sort_polygon_3d(cls, points): - if len(points) <= 3: - return points - o = cls._polygon_inner_point(points) - a = points[0] - oa = Vector(o, a) - norm_oa = oa.norm() - for b in points[1:]: - ob = Vector(o, b) - u = oa.cross(ob) - if not u.isnull(): - u = u.asunit() - break - else: - raise ValueError('degenerate polygon') - angles = {a: 0.} - for m in points[1:]: - om = Vector(o, m) - normprod = norm_oa * om.norm() - cosinus = max(oa.dot(om) / normprod, -1.) - sinus = u.dot(oa.cross(om)) / normprod - angle = math.acos(cosinus) - angle = math.copysign(angle, sinus) - angles[m] = angle - return sorted(points, key=angles.get) - - def faces(self): - vertices = self.vertices() - faces = [] - for constraint in self.constraints: - face = [] - for vertex in vertices: - if constraint.subs(vertex.coordinates()) == 0: - face.append(vertex) - faces.append(face) - return faces - - def _plot_2d(self, plot=None, **kwargs): - from matplotlib import pylab - import matplotlib.pyplot as plt - from matplotlib.axes import Axes - from matplotlib.patches import Polygon - vertices = self._sort_polygon_2d(self.vertices()) - xys = [tuple(vertex.values()) for vertex in vertices] - if plot is None: - fig = plt.figure() - plot = fig.add_subplot(1, 1, 1) - xs, ys = zip(*xys) - plot.set_xlim(float(min(xs)), float(max(xs))) - plot.set_ylim(float(min(ys)), float(max(ys))) - plot.add_patch(Polygon(xys, closed=True, **kwargs)) - return plot - - def _plot_3d(self, plot=None, **kwargs): - import matplotlib.pyplot as plt - from mpl_toolkits.mplot3d import Axes3D - from mpl_toolkits.mplot3d.art3d import Poly3DCollection - if plot is None: - fig = plt.figure() - axes = Axes3D(fig) - xmin, xmax = float('inf'), float('-inf') - ymin, ymax = float('inf'), float('-inf') - zmin, zmax = float('inf'), float('-inf') - else: - axes = plot - poly_xyzs = [] - for vertices in self.faces(): - if len(vertices) == 0: - continue - vertices = Polyhedron._sort_polygon_3d(vertices) - vertices.append(vertices[0]) - face_xyzs = [tuple(vertex.values()) for vertex in vertices] - if plot is None: - xs, ys, zs = zip(*face_xyzs) - xmin, xmax = min(xmin, float(min(xs))), max(xmax, float(max(xs))) - ymin, ymax = min(ymin, float(min(ys))), max(ymax, float(max(ys))) - zmin, zmax = min(zmin, float(min(zs))), max(zmax, float(max(zs))) - poly_xyzs.append(face_xyzs) - collection = Poly3DCollection(poly_xyzs, **kwargs) - axes.add_collection3d(collection) - if plot is None: - axes.set_xlim(xmin, xmax) - axes.set_ylim(ymin, ymax) - axes.set_zlim(zmin, zmax) - return axes - - def plot(self, plot=None, **kwargs): - """ - Display 3D plot of set. - """ - if self.dimension == 2: - return self._plot_2d(plot=plot, **kwargs) - elif self.dimension == 3: - return self._plot_3d(plot=plot, **kwargs) - else: - raise ValueError('polyhedron must be 2 or 3-dimensional') - - -def _polymorphic(func): - @functools.wraps(func) - def wrapper(left, right): - if not isinstance(left, Expression): - if isinstance(left, numbers.Rational): - left = Rational(left) - else: - raise TypeError('left must be a a rational number ' - 'or a linear expression') - if not isinstance(right, Expression): - if isinstance(right, numbers.Rational): - right = Rational(right) - else: - raise TypeError('right must be a a rational number ' - 'or a linear expression') - return func(left, right) - return wrapper - -@_polymorphic -def Lt(left, right): - """ - Return true if the first set is less than the second. - """ - return Polyhedron([], [right - left - 1]) - -@_polymorphic -def Le(left, right): - """ - Return true the first set is less than or equal to the second. - """ - return Polyhedron([], [right - left]) - -@_polymorphic -def Eq(left, right): - """ - Return true if the sets are equal. - """ - return Polyhedron([left - right], []) - -@_polymorphic -def Ne(left, right): - """ - Return true if the sets are NOT equal. - """ - return ~Eq(left, right) - -@_polymorphic -def Gt(left, right): - """ - Return true if the first set is greater than the second set. - """ - return Polyhedron([], [left - right - 1]) - -@_polymorphic -def Ge(left, right): - """ - Return true if the first set is greater than or equal the second set. - """ - return Polyhedron([], [left - right]) - - -Empty = Eq(1, 0) - -Universe = Polyhedron([])