From: Vivien Maisonneuve <v.maisonneuve@gmail.com>
Date: Fri, 11 Jul 2014 12:10:50 +0000 (+0200)
Subject: New module coordinates.py
X-Git-Tag: 1.0~133
X-Git-Url: https://svn.cri.ensmp.fr/git/linpy.git/commitdiff_plain/0995f79a8db531aa3a2ca951af65f3f0fa4943e5

New module 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}))