linpy/geometry.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 math
  19 import numbers
  20 import operator
  21
  22 from abc import ABC, abstractproperty, abstractmethod
  23 from collections import OrderedDict, Mapping
  24
  25 from .linexprs import Symbol
  26
  27
  28 __all__ = [
  29     'GeometricObject',
  30     'Point',
  31     'Vector',
  32 ]
  33
  34
  35 class GeometricObject(ABC):
  36     """
  37     GeometricObject is an abstract class to represent objects with a
  38     geometric representation in space. Subclasses of GeometricObject are
  39     Polyhedron, Domain and Point.
  40     """
  41
  42     @abstractproperty
  43     def symbols(self):
  44         """
  45         The tuple of symbols present in the object expression, sorted according
  46         to Symbol.sortkey().
  47         """
  48         pass
  49
  50     @property
  51     def dimension(self):
  52         """
  53         The dimension of the object, i.e. the number of symbols present in it.
  54         """
  55         return len(self.symbols)
  56
  57     @abstractmethod
  58     def aspolyhedron(self):
  59         """
  60         Return a Polyhedron object that approximates the geometric object.
  61         """
  62         pass
  63
  64     def asdomain(self):
  65         """
  66         Return a Domain object that approximates the geometric object.
  67         """
  68         return self.aspolyhedron()
  69
  70
  71 class Coordinates:
  72     """
  73     This class represents coordinate systems.
  74     """
  75
  76     __slots__ = (
  77         '_coordinates',
  78     )
  79
  80     def __new__(cls, coordinates):
  81         """
  82         Create a coordinate system from a dictionary or a sequence that maps the
  83         symbols to their coordinates. Coordinates must be rational numbers.
  84         """
  85         if isinstance(coordinates, Mapping):
  86             coordinates = coordinates.items()
  87         self = object().__new__(cls)
  88         self._coordinates = []
  89         for symbol, coordinate in coordinates:
  90             if not isinstance(symbol, Symbol):
  91                 raise TypeError('symbols must be Symbol instances')
  92             if not isinstance(coordinate, numbers.Real):
  93                 raise TypeError('coordinates must be real numbers')
  94             self._coordinates.append((symbol, coordinate))
  95         self._coordinates.sort(key=lambda item: item[0].sortkey())
  96         self._coordinates = OrderedDict(self._coordinates)
  97         return self
  98
  99     @property
 100     def symbols(self):
 101         """
 102         The tuple of symbols present in the coordinate system, sorted according
 103         to Symbol.sortkey().
 104         """
 105         return tuple(self._coordinates)
 106
 107     @property
 108     def dimension(self):
 109         """
 110         The dimension of the coordinate system, i.e. the number of symbols
 111         present in it.
 112         """
 113         return len(self.symbols)
 114
 115     def coordinate(self, symbol):
 116         """
 117         Return the coordinate value of the given symbol. Raise KeyError if the
 118         symbol is not involved in the coordinate system.
 119         """
 120         if not isinstance(symbol, Symbol):
 121             raise TypeError('symbol must be a Symbol instance')
 122         return self._coordinates[symbol]
 123
 124     __getitem__ = coordinate
 125
 126     def coordinates(self):
 127         """
 128         Iterate over the pairs (symbol, value) of coordinates in the coordinate
 129         system.
 130         """
 131         yield from self._coordinates.items()
 132
 133     def values(self):
 134         """
 135         Iterate over the coordinate values in the coordinate system.
 136         """
 137         yield from self._coordinates.values()
 138
 139     def __bool__(self):
 140         """
 141         Return True if not all coordinates are 0.
 142         """
 143         return any(self._coordinates.values())
 144
 145     def __eq__(self, other):
 146         """
 147         Return True if two coordinate systems are equal.
 148         """
 149         if isinstance(other, self.__class__):
 150             return self._coordinates == other._coordinates
 151         return NotImplemented
 152
 153     def __hash__(self):
 154         return hash(tuple(self.coordinates()))
 155
 156     def __repr__(self):
 157         string = ', '.join(['{!r}: {!r}'.format(symbol, coordinate)
 158             for symbol, coordinate in self.coordinates()])
 159         return '{}({{{}}})'.format(self.__class__.__name__, string)
 160
 161     def _map(self, func):
 162         for symbol, coordinate in self.coordinates():
 163             yield symbol, func(coordinate)
 164
 165     def _iter2(self, other):
 166         if self.symbols != other.symbols:
 167             raise ValueError('arguments must belong to the same space')
 168         coordinates1 = self._coordinates.values()
 169         coordinates2 = other._coordinates.values()
 170         yield from zip(self.symbols, coordinates1, coordinates2)
 171
 172     def _map2(self, other, func):
 173         for symbol, coordinate1, coordinate2 in self._iter2(other):
 174             yield symbol, func(coordinate1, coordinate2)
 175
 176
 177 class Point(Coordinates, GeometricObject):
 178     """
 179     This class represents points in space.
 180
 181     Point instances are hashable and should be treated as immutable.
 182     """
 183
 184     def isorigin(self):
 185         """
 186         Return True if all coordinates are 0.
 187         """
 188         return not bool(self)
 189
 190     def __hash__(self):
 191         return super().__hash__()
 192
 193     def __add__(self, other):
 194         """
 195         Translate the point by a Vector object and return the resulting point.
 196         """
 197         if isinstance(other, Vector):
 198             coordinates = self._map2(other, operator.add)
 199             return Point(coordinates)
 200         return NotImplemented
 201
 202     def __sub__(self, other):
 203         """
 204         If other is a point, substract it from self and return the resulting
 205         vector. If other is a vector, translate the point by the opposite vector
 206         and returns the resulting point.
 207         """
 208         coordinates = []
 209         if isinstance(other, Point):
 210             coordinates = self._map2(other, operator.sub)
 211             return Vector(coordinates)
 212         elif isinstance(other, Vector):
 213             coordinates = self._map2(other, operator.sub)
 214             return Point(coordinates)
 215         return NotImplemented
 216
 217     def aspolyhedron(self):
 218         from .polyhedra import Polyhedron
 219         equalities = []
 220         for symbol, coordinate in self.coordinates():
 221             equalities.append(symbol - coordinate)
 222         return Polyhedron(equalities)
 223
 224
 225 class Vector(Coordinates):
 226     """
 227     This class represents vectors in space.
 228
 229     Vector instances are hashable and should be treated as immutable.
 230     """
 231
 232     def __new__(cls, initial, terminal=None):
 233         """
 234         Create a vector from a dictionary or a sequence that maps the symbols to
 235         their coordinates, or as the displacement between two points.
 236         """
 237         if not isinstance(initial, Point):
 238             initial = Point(initial)
 239         if terminal is None:
 240             coordinates = initial._coordinates
 241         else:
 242             if not isinstance(terminal, Point):
 243                 terminal = Point(terminal)
 244             coordinates = terminal._map2(initial, operator.sub)
 245         return super().__new__(cls, coordinates)
 246
 247     def isnull(self):
 248         """
 249         Return True if all coordinates are 0.
 250         """
 251         return not bool(self)
 252
 253     def __hash__(self):
 254         return super().__hash__()
 255
 256     def __add__(self, other):
 257         """
 258         If other is a point, translate it with the vector self and return the
 259         resulting point. If other is a vector, return the vector self + other.
 260         """
 261         if isinstance(other, (Point, Vector)):
 262             coordinates = self._map2(other, operator.add)
 263             return other.__class__(coordinates)
 264         return NotImplemented
 265
 266     def __sub__(self, other):
 267         """
 268         If other is a point, substract it from the vector self and return the
 269         resulting point. If other is a vector, return the vector self - other.
 270         """
 271         if isinstance(other, (Point, Vector)):
 272             coordinates = self._map2(other, operator.sub)
 273             return other.__class__(coordinates)
 274         return NotImplemented
 275
 276     def __neg__(self):
 277         """
 278         Return the vector -self.
 279         """
 280         coordinates = self._map(operator.neg)
 281         return Vector(coordinates)
 282
 283     def __mul__(self, other):
 284         """
 285         Multiplies a Vector by a scalar value.
 286         """
 287         if isinstance(other, numbers.Real):
 288             coordinates = self._map(lambda coordinate: other * coordinate)
 289             return Vector(coordinates)
 290         return NotImplemented
 291
 292     __rmul__ = __mul__
 293
 294     def __truediv__(self, other):
 295         """
 296         Divide the vector by the specified scalar and returns the result as a
 297         vector.
 298         """
 299         if isinstance(other, numbers.Real):
 300             coordinates = self._map(lambda coordinate: coordinate / other)
 301             return Vector(coordinates)
 302         return NotImplemented
 303
 304     def angle(self, other):
 305         """
 306         Retrieve the angle required to rotate the vector into the vector passed
 307         in argument. The result is an angle in radians, ranging between -pi and
 308         pi.
 309         """
 310         if not isinstance(other, Vector):
 311             raise TypeError('argument must be a Vector instance')
 312         cosinus = self.dot(other) / (self.norm()*other.norm())
 313         return math.acos(cosinus)
 314
 315     def cross(self, other):
 316         """
 317         Compute the cross product of two 3D vectors. If either one of the
 318         vectors is not three-dimensional, a ValueError exception is raised.
 319         """
 320         if not isinstance(other, Vector):
 321             raise TypeError('other must be a Vector instance')
 322         if self.dimension != 3 or other.dimension != 3:
 323             raise ValueError('arguments must be three-dimensional vectors')
 324         if self.symbols != other.symbols:
 325             raise ValueError('arguments must belong to the same space')
 326         x, y, z = self.symbols
 327         coordinates = []
 328         coordinates.append((x, self[y]*other[z] - self[z]*other[y]))
 329         coordinates.append((y, self[z]*other[x] - self[x]*other[z]))
 330         coordinates.append((z, self[x]*other[y] - self[y]*other[x]))
 331         return Vector(coordinates)
 332
 333     def dot(self, other):
 334         """
 335         Compute the dot product of two vectors.
 336         """
 337         if not isinstance(other, Vector):
 338             raise TypeError('argument must be a Vector instance')
 339         result = 0
 340         for symbol, coordinate1, coordinate2 in self._iter2(other):
 341             result += coordinate1 * coordinate2
 342         return result
 343
 344     def __hash__(self):
 345         return super().__hash__()
 346
 347     def norm(self):
 348         """
 349         Return the norm of the vector.
 350         """
 351         return math.sqrt(self.norm2())
 352
 353     def norm2(self):
 354         """
 355         Return the squared norm of the vector.
 356         """
 357         result = 0
 358         for coordinate in self._coordinates.values():
 359             result += coordinate ** 2
 360         return result
 361
 362     def asunit(self):
 363         """
 364         Return the normalized vector, i.e. the vector of same direction but with
 365         norm 1.
 366         """
 367         return self / self.norm()