Minor improvements in Makefile
[linpy.git] / 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()