From 0995f79a8db531aa3a2ca951af65f3f0fa4943e5 Mon Sep 17 00:00:00 2001 From: Vivien Maisonneuve Date: Fri, 11 Jul 2014 14:10:50 +0200 Subject: [PATCH] New module coordinates.py --- pypol/__init__.py | 2 + pypol/coordinates.py | 244 ++++++++++++++++++++++++++++++++ pypol/domains.py | 15 +- pypol/tests/test_coordinates.py | 87 ++++++++++++ 4 files changed, 341 insertions(+), 7 deletions(-) create mode 100644 pypol/coordinates.py create mode 100644 pypol/tests/test_coordinates.py diff --git a/pypol/__init__.py b/pypol/__init__.py index e39bce7..be1440b 100644 --- a/pypol/__init__.py +++ b/pypol/__init__.py @@ -3,12 +3,14 @@ A polyhedral library based on ISL. """ from .linexprs import Expression, Symbol, Dummy, symbols, Rational +from .coordinates import Point, Vector from .polyhedra import Polyhedron, Eq, Ne, Le, Lt, Ge, Gt, Ne, Empty, Universe from .domains import Domain, And, Or, Not __all__ = [ 'Expression', 'Symbol', 'Dummy', 'symbols', 'Rational', + 'Point', 'Vector', 'Polyhedron', 'Eq', 'Ne', 'Le', 'Lt', 'Ge', 'Gt', 'Empty', 'Universe', 'Domain', 'And', 'Or', 'Not', ] diff --git a/pypol/coordinates.py b/pypol/coordinates.py new file mode 100644 index 0000000..b5c23c9 --- /dev/null +++ b/pypol/coordinates.py @@ -0,0 +1,244 @@ +import math +import numbers +import operator + +from collections import OrderedDict + +from .linexprs import Symbol + + +__all__ = [ + 'Point', + 'Vector', +] + + +def _map(obj, func): + for symbol, coordinate in obj.coordinates(): + yield symbol, func(coordinate) + +def _iter2(obj1, obj2): + if obj1.symbols != obj2.symbols: + raise ValueError('arguments must belong to the same space') + coordinates1 = obj1._coordinates.values() + coordinates2 = obj2._coordinates.values() + yield from zip(obj1.symbols, coordinates1, coordinates2) + +def _map2(obj1, obj2, func): + for symbol, coordinate1, coordinate2 in _iter2(obj1, obj2): + yield symbol, func(coordinate1, coordinate2) + + +class Point: + """ + This class represents points in space. + """ + + def __new__(cls, coordinates=None): + if isinstance(coordinates, dict): + coordinates = coordinates.items() + self = object().__new__(cls) + self._coordinates = OrderedDict() + for symbol, coordinate in sorted(coordinates, + key=lambda item: item[0].sortkey()): + if not isinstance(symbol, Symbol): + raise TypeError('symbols must be Symbol instances') + if not isinstance(coordinate, numbers.Real): + raise TypeError('coordinates must be real numbers') + self._coordinates[symbol] = coordinate + return self + + @property + def symbols(self): + return tuple(self._coordinates) + + @property + def dimension(self): + return len(self.symbols) + + def coordinates(self): + yield from self._coordinates.items() + + def coordinate(self, symbol): + if not isinstance(symbol, Symbol): + raise TypeError('symbol must be a Symbol instance') + return self._coordinates[symbol] + + __getitem__ = coordinate + + def isorigin(self): + return not bool(self) + + def __bool__(self): + return any(self._coordinates.values()) + + def __add__(self, other): + if not isinstance(other, Vector): + return NotImplemented + coordinates = _map2(self, other, operator.add) + return Point(coordinates) + + def __sub__(self, other): + coordinates = [] + if isinstance(other, Point): + coordinates = _map2(self, other, operator.sub) + return Vector(coordinates) + elif isinstance(other, Vector): + coordinates = _map2(self, other, operator.sub) + return Point(coordinates) + else: + return NotImplemented + + def __eq__(self, other): + return isinstance(other, Point) and \ + self._coordinates == other._coordinates + + def __hash__(self): + return hash(tuple(self.coordinates())) + + def __repr__(self): + string = ', '.join(['{!r}: {!r}'.format(symbol, coordinate) + for symbol, coordinate in self.coordinates()]) + return '{}({{{}}})'.format(self.__class__.__name__, string) + + +class Vector: + """ + This class represents displacements in space. + """ + + __slots__ = ( + '_coordinates', + ) + + def __new__(cls, initial, terminal=None): + self = object().__new__(cls) + if not isinstance(initial, Point): + initial = Point(initial) + if terminal is None: + self._coordinates = initial._coordinates + elif not isinstance(terminal, Point): + terminal = Point(terminal) + self._coordinates = _map2(terminal, initial, operator.sub) + return self + + @property + def symbols(self): + return tuple(self._coordinates) + + @property + def dimension(self): + return len(self.symbols) + + def coordinates(self): + yield from self._coordinates.items() + + def coordinate(self, symbol): + if not isinstance(symbol, Symbol): + raise TypeError('symbol must be a Symbol instance') + return self._coordinates[symbol] + + __getitem__ = coordinate + + def isnull(self): + return not bool(self) + + def __bool__(self): + return any(self._coordinates.values()) + + def __add__(self, other): + if isinstance(other, (Point, Vector)): + coordinates = _map2(self, other, operator.add) + return other.__class__(coordinates) + return NotImplemented + + def angle(self, other): + """ + Retrieve the angle required to rotate the vector into the vector passed + in argument. The result is an angle in radians, ranging between -pi and + pi. + """ + if not isinstance(other, Vector): + raise TypeError('argument must be a Vector instance') + cosinus = self.dot(other) / (self.norm() * other.norm()) + return math.acos(cosinus) + + def cross(self, other): + """ + Calculate the cross product of two Vector3D structures. + """ + if not isinstance(other, Vector): + raise TypeError('other must be a Vector instance') + if self.dimension != 3 or other.dimension != 3: + raise ValueError('arguments must be three-dimensional vectors') + if self.symbols != other.symbols: + raise ValueError('arguments must belong to the same space') + x, y, z = self.symbols + coordinates = [] + coordinates.append((x, self[y]*other[z] - self[z]*other[y])) + coordinates.append((y, self[z]*other[x] - self[x]*other[z])) + coordinates.append((z, self[x]*other[y] - self[y]*other[x])) + return Vector(coordinates) + + def __truediv__(self, other): + """ + Divide the vector by the specified scalar and returns the result as a + vector. + """ + if not isinstance(other, numbers.Real): + return NotImplemented + coordinates = _map(self, lambda coordinate: coordinate / other) + return Vector(coordinates) + + def dot(self, other): + """ + Calculate the dot product of two vectors. + """ + if not isinstance(other, Vector): + raise TypeError('argument must be a Vector instance') + result = 0 + for symbol, coordinate1, coordinate2 in _iter2(self, other): + result += coordinate1 * coordinate2 + return result + + def __eq__(self, other): + return isinstance(other, Vector) and \ + self._coordinates == other._coordinates + + def __hash__(self): + return hash(tuple(self.coordinates())) + + def __mul__(self, other): + if not isinstance(other, numbers.Real): + return NotImplemented + coordinates = _map(self, lambda coordinate: other * coordinate) + return Vector(coordinates) + + __rmul__ = __mul__ + + def __neg__(self): + coordinates = _map(self, operator.neg) + return Vector(coordinates) + + def norm(self): + return math.sqrt(self.norm2()) + + def norm2(self): + result = 0 + for coordinate in self._coordinates.values(): + result += coordinate ** 2 + return result + + def asunit(self): + return self / self.norm() + + def __sub__(self, other): + if isinstance(other, (Point, Vector)): + coordinates = _map2(self, other, operator.sub) + return other.__class__(coordinates) + return NotImplemented + + def __repr__(self): + string = ', '.join(['{!r}: {!r}'.format(symbol, coordinate) + for symbol, coordinate in self.coordinates()]) + return '{}({{{}}})'.format(self.__class__.__name__, string) diff --git a/pypol/domains.py b/pypol/domains.py index 45040c8..3ab7d44 100644 --- a/pypol/domains.py +++ b/pypol/domains.py @@ -6,6 +6,7 @@ from fractions import Fraction from . import islhelper from .islhelper import mainctx, libisl, isl_set_basic_sets +from .coordinates import Point from .linexprs import Expression, Symbol @@ -268,7 +269,7 @@ class Domain: points = [] for vertex in vertices: expr = libisl.isl_vertex_get_expr(vertex) - point = {} + coordinates = [] if islhelper.isl_version < '0.13': constraints = islhelper.isl_basic_set_constraints(expr) for constraint in constraints: @@ -280,7 +281,7 @@ class Domain: coefficient = islhelper.isl_val_to_int(coefficient) if coefficient != 0: coordinate = -Fraction(constant, coefficient) - point[symbol]= coordinate + coordinates.append((symbol, coordinate)) else: # horrible hack, find a cleaner solution string = islhelper.isl_multi_aff_to_str(expr) @@ -290,8 +291,8 @@ class Domain: denominator = match.group('den') denominator = 1 if denominator is None else int(denominator) coordinate = Fraction(numerator, denominator) - point[symbol] = coordinate - points.append(point) + coordinates.append((symbol, coordinate)) + points.append(Point(coordinates)) return points def points(self): @@ -302,13 +303,13 @@ class Domain: islpoints = islhelper.isl_set_points(islset) points = [] for islpoint in islpoints: - point = {} + coordinates = {} for index, symbol in enumerate(self.symbols): coordinate = libisl.isl_point_get_coordinate_val(islpoint, libisl.isl_dim_set, index) coordinate = islhelper.isl_val_to_int(coordinate) - point[symbol] = coordinate - points.append(point) + coordinates[symbol] = coordinate + points.append(Point(coordinates)) return points def subs(self, symbol, expression=None): diff --git a/pypol/tests/test_coordinates.py b/pypol/tests/test_coordinates.py new file mode 100644 index 0000000..2e14032 --- /dev/null +++ b/pypol/tests/test_coordinates.py @@ -0,0 +1,87 @@ +import math +import unittest + +from ..linexprs import Symbol +from ..coordinates import * + + +class TestPoint(unittest.TestCase): + + def setUp(self): + self.x = Symbol('x') + self.y = Symbol('y') + self.z = Symbol('z') + self.pt1 = Point({self.x: 10, self.y: 5, self.z: 1}) + self.pt2 = Point({self.x: 15, self.y: 40, self.z: 60}) + self.vec1 = Vector({self.x: 20, self.y: 30, self.z: 40}) + + def test_add(self): + self.assertEqual(self.pt1 + self.vec1, Point({self.x: 30, self.y: 35, self.z: 41})) + with self.assertRaises(TypeError): + self.pt1 + self.pt2 + + def test_eq(self): + self.assertEqual(self.pt1, self.pt1) + self.assertNotEqual(self.pt1, self.pt2) + self.assertNotEqual(self.pt1, self.vec1) + + def test_sub(self): + self.assertEqual(self.pt1 - self.pt2, Vector({self.x: -5, self.y: -35, self.z: -59})) + self.assertEqual(self.pt1 - self.vec1, Point({self.x: -10, self.y: -25, self.z: -39})) + + +class TestVector(unittest.TestCase): + + def setUp(self): + self.x = Symbol('x') + self.y = Symbol('y') + self.z = Symbol('z') + self.pt1 = Point({self.x: 10, self.y: 5, self.z: 1}) + self.pt2 = Point({self.x: 15, self.y: 40, self.z: 60}) + self.vec1 = Vector({self.x: 20, self.y: 30, self.z: 40}) + self.vec2 = Vector({self.x: 45, self.y: 70, self.z: 80}) + + def test_add(self): + self.assertEqual(self.vec1 + self.pt1, Point({self.x: 30, self.y: 35, self.z: 41})) + self.assertEqual(self.vec1 + self.vec2, Vector({self.x: 65, self.y: 100, self.z: 120})) + + def test_angle(self): + self.assertEqual(math.degrees(self.vec1.angle(self.vec1)), 0) + self.assertAlmostEqual(math.degrees(self.vec1.angle(self.vec2)), 4.15129, places=5) + self.assertAlmostEqual(math.degrees(self.vec2.angle(self.vec1)), 4.15129, places=5) + + def test_cross(self): + self.assertEqual(self.vec1.cross(self.vec2), Vector({self.x: -400, self.y: 200, self.z: 50})) + + def test_div(self): + self.assertEqual(self.vec1 / 10, Vector({self.x: 2, self.y: 3, self.z: 4})) + + def test_dot(self): + self.assertEqual(self.vec1.dot(self.vec2), 6200) + + def test_eq(self): + self.assertEqual(self.vec1, self.vec1) + self.assertNotEqual(self.vec1, self.vec2) + + def test_mul(self): + self.assertEqual(75 * self.vec1, Vector({self.x: 1500, self.y: 2250, self.z: 3000})) + self.assertEqual(self.vec1 * 75, Vector({self.x: 1500, self.y: 2250, self.z: 3000})) + + def test_neg(self): + self.assertEqual(-self.vec1, Vector({self.x: -20, self.y: -30, self.z: -40})) + + def test_norm(self): + self.assertAlmostEqual(self.vec1.norm(), 53.85165, places=5) + + def test_norm2(self): + self.assertEqual(self.vec1.norm2(), 2900) + + def test_asunit(self): + unit = self.vec1.asunit() + self.assertAlmostEqual(unit[self.x], 0.37139, 5) + self.assertAlmostEqual(unit[self.y], 0.55709, 5) + self.assertAlmostEqual(unit[self.z], 0.74278, 5) + + def test_sub(self): + self.assertEqual(self.vec1 - self.pt1, Point({self.x: 10, self.y: 25, self.z: 39})) + self.assertEqual(self.vec1 - self.vec2, Vector({self.x: -25, self.y: -40, self.z: -40})) -- 2.20.1