Fix LinExpr.subs() implementation
[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 = 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 a point from another and returns the
196 resulting vector. If other is a vector, translate the point by the
197 opposite vector 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()