Coverage for src/algolib/maths/geometry/geometry.py: 68%
146 statements
« prev ^ index » next coverage.py v7.10.4, created at 2025-08-20 19:37 +0000
1# src/algolib/maths/geometry/geometry.py
2"""
3A professional-grade :math:`N`-dimensional geometry module.
5It implements:
7- **Point**: a location in :math:`\\mathbb{R}^N`.
8- **Vector**: a displacement in :math:`\\mathbb{R}^N`.
9- **Line**: parametric line :math:`P(t) = P_0 + t\\,d`.
10- **Plane**: hyperplane :math:`\\{X:\\; n\\cdot(X-P_0)=0\\}`.
11- **GeometryUtils**: utility routines.
13All classes validate dimensions and numeric inputs, and raise
14:class:`algolib.exceptions.InvalidTypeError` or
15:class:`algolib.exceptions.InvalidValueError` on invalid usage.
16"""
18from __future__ import annotations
20import math
21from typing import Iterable, List, Sequence, Union
23from algolib.exceptions import InvalidTypeError, InvalidValueError
25Number = Union[int, float]
28def _ensure_numbers(xs: Iterable[Number]) -> List[float]:
29 try:
30 out = [float(x) for x in xs]
31 except Exception as e: # noqa: BLE001
32 raise InvalidTypeError("all coordinates/components must be real numbers.") from e
33 return out
36def _same_dim(a_len: int, b_len: int) -> None:
37 if a_len != b_len:
38 raise InvalidValueError("dimensions must match.")
41def _is_zero_vector(v: Sequence[float]) -> bool:
42 return all(c == 0.0 for c in v)
45class Point:
46 r"""
47 Point in :math:`N`-dimensional Euclidean space.
49 Parameters
50 ----------
51 coords : Sequence[Number]
52 Coordinates of the point; length defines the dimension.
54 Notes
55 -----
56 A point has location but no direction or magnitude:
57 :math:`P=(x_1,x_2,\\dots,x_N)`.
58 """
60 def __init__(self, coords: Sequence[Number]):
61 cs = _ensure_numbers(coords)
62 if len(cs) == 0:
63 raise InvalidValueError("point must have at least one coordinate.")
64 self.coords: List[float] = cs
66 def __repr__(self) -> str: # pragma: no cover - trivial
67 return f"Point({self.coords})"
69 def dimension(self) -> int:
70 """Return the dimension (number of coordinates)."""
71 return len(self.coords)
74class Vector:
75 r"""
76 Vector in :math:`N`-dimensional Euclidean space.
78 Parameters
79 ----------
80 comps : Sequence[Number]
81 Components of the vector.
83 Notes
84 -----
85 A vector represents both magnitude and direction.
86 """
88 def __init__(self, comps: Sequence[Number]):
89 cs = _ensure_numbers(comps)
90 if len(cs) == 0:
91 raise InvalidValueError("vector must have at least one component.")
92 self.comps: List[float] = cs
94 def __repr__(self) -> str: # pragma: no cover - trivial
95 return f"Vector({self.comps})"
97 def dimension(self) -> int:
98 """Return the dimension (number of components)."""
99 return len(self.comps)
101 def norm(self) -> float:
102 r"""
103 Return the Euclidean norm :math:`\lVert v \rVert = \sqrt{\sum_i v_i^2}`.
105 Notes
106 -----
107 - Non-finite inputs are handled explicitly:
108 if any component is NaN, the result is NaN;
109 if any component has infinite magnitude, the result is ``inf``.
110 - Uses a scaling strategy for numerical stability to avoid overflow
111 and underflow with mixed-magnitude components.
112 """
113 # Fast paths for non-finite inputs without importing math
114 any_inf = False
115 for c in self.comps:
116 # NaN check: NaN != NaN
117 if c != c:
118 return float("nan")
119 if c == float("inf") or c == float("-inf"):
120 any_inf = True
121 if any_inf:
122 return float("inf")
124 # Stable accumulation with scaling (Smith's method)
125 scale = 0.0
126 sumsq = 1.0
127 nonzero_seen = False
129 for c in self.comps:
130 if c == 0.0:
131 continue
132 nonzero_seen = True
133 abs_c = c if c >= 0.0 else -c
134 if scale < abs_c:
135 if scale == 0.0:
136 # first non-zero term
137 scale = abs_c
138 sumsq = 1.0
139 else:
140 sumsq = 1.0 + sumsq * (scale / abs_c) ** 2
141 scale = abs_c
142 else:
143 # here scale >= abs_c and scale > 0.0, safe to divide
144 sumsq += (abs_c / scale) ** 2
146 if not nonzero_seen:
147 return 0.0
149 # Use built-in power for sqrt without importing math
150 return float(scale * (sumsq ** 0.5))
152 def dot(self, other: "Vector") -> float:
153 r"""
154 Return the dot product with another vector.
156 Formula
157 -------
158 :math:`v\\cdot w = \\sum_i v_i w_i`.
159 """
160 if not isinstance(other, Vector):
161 raise InvalidTypeError("other must be Vector.")
162 _same_dim(self.dimension(), other.dimension())
163 total = 0.0
164 for a, b in zip(self.comps, other.comps):
165 # NaN propagation: if either input is NaN, the dot is NaN
166 if a != a or b != b: # NaN check
167 return float("nan")
169 # Handle indeterminate 0 * inf (or inf * 0) as contributing 0.0
170 if (a == 0.0 and (b == float("inf") or b == float("-inf"))) or \
171 (b == 0.0 and (a == float("inf") or a == float("-inf"))):
172 # mathematically this term is 0; avoid NaN from IEEE 0*inf
173 continue
175 total += a * b
176 return total
179 # convenience ops
180 def __add__(self, other: "Vector") -> "Vector":
181 """Component-wise addition."""
182 if not isinstance(other, Vector):
183 raise InvalidTypeError("other must be Vector.")
184 _same_dim(self.dimension(), other.dimension())
185 return Vector([a + b for a, b in zip(self.comps, other.comps)])
187 def __sub__(self, other: "Vector") -> "Vector":
188 """Component-wise subtraction."""
189 if not isinstance(other, Vector):
190 raise InvalidTypeError("other must be Vector.")
191 _same_dim(self.dimension(), other.dimension())
192 return Vector([a - b for a, b in zip(self.comps, other.comps)])
194 def __mul__(self, k: Number) -> "Vector":
195 """Scalar multiplication ``v * k``.
197 Notes
198 -----
199 - For finite scalars, behaves like standard component-wise multiplication.
200 - For ``k`` being ``±inf``, we avoid the indeterminate form ``0 * inf``
201 producing ``NaN`` by returning ``0.0`` for zero components and
202 ``±inf`` (with proper sign) for non-zero components. This mirrors the
203 typical intent in normalization code where one would divide by a tiny
204 norm rather than multiply by its overflowing reciprocal.
205 - If ``k`` is ``NaN``, propagate ``NaN`` to all components.
206 """
207 if not isinstance(k, (int, float)):
208 raise InvalidTypeError("scalar must be real.")
210 # Handle NaN scalar: propagate NaNs
211 if k != k: # NaN
212 return Vector([float("nan") for _ in self.comps])
214 # Handle infinite scalars: return a signed normalized vector to
215 # avoid generating ±inf components (which later interact with zeros
216 # and produce NaNs in downstream helpers).
217 if k == float("inf") or k == float("-inf"):
218 nrm = self.norm()
219 if nrm == 0.0:
220 # 0 * inf -> keep zero vector (no direction)
221 return Vector([0.0 for _ in self.comps])
222 sign_k = 1.0 if k > 0.0 else -1.0
223 # Avoid forming an infinite scalar; normalize via per-component division
224 out = [ (a / nrm) for a in self.comps ]
225 if sign_k < 0.0:
226 out = [ -v for v in out ]
227 return Vector(out)
229 # Finite scalar: regular multiply
230 return Vector([k * a for a in self.comps])
232 __rmul__ = __mul__ # k * v
235class Line:
236 r"""
237 Parametric line :math:`P(t) = P_0 + t\\,d` in :math:`\\mathbb{R}^N`.
239 Parameters
240 ----------
241 point : Point
242 A point on the line (:math:`P_0`).
243 direction : Vector
244 Direction vector :math:`d` (must be non-zero).
245 """
247 def __init__(self, point: Point, direction: Vector):
248 if not isinstance(point, Point) or not isinstance(direction, Vector):
249 raise InvalidTypeError("point must be Point and direction must be Vector.")
250 _same_dim(point.dimension(), direction.dimension())
251 if _is_zero_vector(direction.comps):
252 raise InvalidValueError("direction vector must be non-zero.")
253 self.point = point
254 self.direction = direction
256 def __repr__(self) -> str: # pragma: no cover - trivial
257 return f"Line(point={self.point}, direction={self.direction})"
259 def point_at(self, t: Number) -> Point:
260 """Return the point :math:`P_0 + t\\,d`."""
261 if not isinstance(t, (int, float)):
262 raise InvalidTypeError("t must be real.")
263 return Point([p + float(t) * d for p, d in zip(self.point.coords, self.direction.comps)])
265 def contains(self, p: Point, tol: float = 1e-12) -> bool:
266 """
267 Return True if ``p`` lies on the line (within tolerance).
269 We check that :math:`p - P_0` is colinear with :math:`d` by comparing ratios;
270 component pairs with :math:`|d_i| \\le \\text{tol}` are skipped.
271 """
272 _same_dim(self.point.dimension(), p.dimension())
273 ratios = []
274 for (pi, p0i), di in zip(zip(p.coords, self.point.coords), self.direction.comps):
275 if abs(di) <= tol:
276 if abs(pi - p0i) > tol:
277 return False # movement along a zero direction -> off-line
278 continue
279 ratios.append((pi - p0i) / di)
280 if not ratios:
281 # direction is almost zero in all components -> treat as aligned if p == P0
282 return all(abs(pi - p0i) <= tol for (pi, p0i) in zip(p.coords, self.point.coords))
283 first = ratios[0]
284 return all(abs(r - first) <= tol for r in ratios[1:])
287class Plane:
288 r"""
289 Hyperplane :math:`\\{X:\\; n\\cdot(X-P_0)=0\\}` in :math:`\\mathbb{R}^N`.
291 Parameters
292 ----------
293 point : Point
294 A reference point :math:`P_0` on the plane.
295 normal : Vector
296 Normal vector :math:`n` (must be non-zero).
297 """
299 def __init__(self, point: Point, normal: Vector):
300 if not isinstance(point, Point) or not isinstance(normal, Vector):
301 raise InvalidTypeError("point must be Point and normal must be Vector.")
302 _same_dim(point.dimension(), normal.dimension())
303 if _is_zero_vector(normal.comps):
304 raise InvalidValueError("normal vector must be non-zero.")
305 self.point = point
306 self.normal = normal
308 def __repr__(self) -> str: # pragma: no cover - trivial
309 return f"Plane(point={self.point}, normal={self.normal})"
311 def signed_distance(self, p: Point) -> float:
312 r"""
313 Return signed distance :math:`\\frac{n\\cdot (p-P_0)}{\\lVert n\\rVert}`.
314 """
315 _same_dim(self.point.dimension(), p.dimension())
316 n2 = self.normal.dot(self.normal)
317 if n2 == 0.0: # guarded at init; defensive
318 raise InvalidValueError("normal vector must be non-zero.")
319 diff = Vector([a - b for a, b in zip(p.coords, self.point.coords)])
320 return self.normal.dot(diff) / math.sqrt(n2)
322 def contains(self, p: Point, tol: float = 1e-12) -> bool:
323 """Return True if ``|n·(p-P0)| <= tol * ||n||``."""
324 return abs(self.signed_distance(p)) <= tol
327class GeometryUtils:
328 """Common geometry utilities."""
330 @staticmethod
331 def distance(p1: Point, p2: Point) -> float:
332 r"""
333 Euclidean distance :math:`\\sqrt{\\sum_i (x_{1i}-x_{2i})^2}`.
334 """
335 _same_dim(p1.dimension(), p2.dimension())
336 return math.sqrt(sum((a - b) ** 2 for a, b in zip(p1.coords, p2.coords)))