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 isinstance(other, Vector):
 189             coordinates = self._map2(other, operator.add)
 190             return Point(coordinates)
 191         return NotImplemented
 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         return NotImplemented
 207
 208     def __eq__(self, other):
 209         """
 210         Test whether two points are equal.
 211         """
 212         if isinstance(other, Point):
 213             return self._coordinates == other._coordinates
 214         return NotImplemented
 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 displacement 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 isinstance(other, numbers.Real):
 287             coordinates = self._map(lambda coordinate: other * coordinate)
 288             return Vector(coordinates)
 289         return NotImplemented
 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 isinstance(other, numbers.Real):
 299             coordinates = self._map(lambda coordinate: coordinate / other)
 300             return Vector(coordinates)
 301         return NotImplemented
 302
 303     def __eq__(self, other):
 304         """
 305         Test whether two vectors are equal.
 306         """
 307         if isinstance(other, Vector):
 308             return self._coordinates == other._coordinates
 309         return NotImplemented
 310
 311     def angle(self, other):
 312         """
 313         Retrieve the angle required to rotate the vector into the vector passed
 314         in argument. The result is an angle in radians, ranging between -pi and
 315         pi.
 316         """
 317         if not isinstance(other, Vector):
 318             raise TypeError('argument must be a Vector instance')
 319         cosinus = self.dot(other) / (self.norm()*other.norm())
 320         return math.acos(cosinus)
 321
 322     def cross(self, other):
 323         """
 324         Compute the cross product of two 3D vectors. If either one of the
 325         vectors is not three-dimensional, a ValueError exception is raised.
 326         """
 327         if not isinstance(other, Vector):
 328             raise TypeError('other must be a Vector instance')
 329         if self.dimension != 3 or other.dimension != 3:
 330             raise ValueError('arguments must be three-dimensional vectors')
 331         if self.symbols != other.symbols:
 332             raise ValueError('arguments must belong to the same space')
 333         x, y, z = self.symbols
 334         coordinates = []
 335         coordinates.append((x, self[y]*other[z] - self[z]*other[y]))
 336         coordinates.append((y, self[z]*other[x] - self[x]*other[z]))
 337         coordinates.append((z, self[x]*other[y] - self[y]*other[x]))
 338         return Vector(coordinates)
 339
 340     def dot(self, other):
 341         """
 342         Compute the dot product of two vectors.
 343         """
 344         if not isinstance(other, Vector):
 345             raise TypeError('argument must be a Vector instance')
 346         result = 0
 347         for symbol, coordinate1, coordinate2 in self._iter2(other):
 348             result += coordinate1 * coordinate2
 349         return result
 350
 351     def __hash__(self):
 352         return super().__hash__()
 353
 354     def norm(self):
 355         """
 356         Return the norm of the vector.
 357         """
 358         return math.sqrt(self.norm2())
 359
 360     def norm2(self):
 361         """
 362         Return the squared norm of the vector.
 363         """
 364         result = 0
 365         for coordinate in self._coordinates.values():
 366             result += coordinate ** 2
 367         return result
 368
 369     def asunit(self):
 370         """
 371         Return the normalized vector, i.e. the vector of same direction but with
 372         norm 1.
 373         """
 374         return self / self.norm()