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