Coverage for src/algolib/maths/algebra/matrix_dense.py: 91%
147 statements
« prev ^ index » next coverage.py v7.10.4, created at 2025-08-20 19:37 +0000
1# src/algolib/maths/algebra/matrix_dense.py
2from __future__ import annotations
4from dataclasses import dataclass
5from typing import Iterable, List, Sequence, Tuple, Union, overload
7from algolib.exceptions import InvalidTypeError, InvalidValueError
9Number = Union[int, float]
12@dataclass(frozen=True)
13class MatrixDense:
14 """
15 A simple dense matrix (row-major) without external dependencies.
17 Designed for correctness and clarity (teaching/learning oriented).
18 Not optimized for large-scale numerical workloads.
20 Parameters
21 ----------
22 rows : Sequence[Sequence[Number]]
23 Row-major data. Must be rectangular (all rows same length).
25 Raises
26 ------
27 InvalidTypeError
28 If rows are not sequences of numbers.
29 InvalidValueError
30 If matrix is empty or not rectangular.
32 Examples
33 --------
34 >>> A = MatrixDense([[1, 2], [3, 4]])
35 >>> B = MatrixDense.identity(2)
36 >>> C = (A * B).rows == A.rows
37 True
38 """
40 rows: Tuple[Tuple[Number, ...], ...] # immutable, rectangular
42 # ---------- Constructors ----------
44 def __init__(self, rows: Sequence[Sequence[Number]]):
45 # Validate type & rectangular
46 if not isinstance(rows, Sequence) or len(rows) == 0:
47 raise InvalidValueError("rows must be a non-empty sequence of sequences.")
48 n_cols = None
49 norm_rows: List[Tuple[Number, ...]] = []
50 for r in rows:
51 if not isinstance(r, Sequence) or len(r) == 0:
52 raise InvalidValueError("each row must be a non-empty sequence.")
53 if n_cols is None:
54 n_cols = len(r)
55 elif len(r) != n_cols:
56 raise InvalidValueError("matrix must be rectangular (same number of columns per row).")
57 # check elements are numbers
58 for x in r:
59 if not isinstance(x, (int, float)):
60 raise InvalidTypeError("matrix elements must be int or float.")
61 norm_rows.append(tuple(x for x in r))
62 object.__setattr__(self, "rows", tuple(norm_rows))
64 @staticmethod
65 def from_rows(rows: Sequence[Sequence[Number]]) -> "MatrixDense":
66 """Construct from row-major data (alias of the constructor)."""
67 return MatrixDense(rows)
69 @staticmethod
70 def zeros(n_rows: int, n_cols: int) -> "MatrixDense":
71 """Return an :math:`(n_{\\mathrm{rows}} × n_{\\mathrm{cols}})` zero matrix."""
72 if n_rows <= 0 or n_cols <= 0:
73 raise InvalidValueError("matrix dimensions must be positive.")
74 return MatrixDense([[0.0] * n_cols for _ in range(n_rows)])
76 @staticmethod
77 def identity(n: int) -> "MatrixDense":
78 """Return the :math:`n \\times n` identity matrix."""
79 if n <= 0:
80 raise InvalidValueError("size must be positive.")
81 return MatrixDense([[1.0 if i == j else 0.0 for j in range(n)] for i in range(n)])
83 eye = identity
84 """Alias of :func:`identity`."""
86 # ---------- Basic properties ----------
88 @property
89 def shape(self) -> Tuple[int, int]:
90 """(n_rows, n_cols)"""
91 return (len(self.rows), len(self.rows[0]))
93 def copy(self) -> "MatrixDense":
94 """Shallow copy (rows are tuples; safe)."""
95 return MatrixDense([list(r) for r in self.rows])
97 # ---------- Equality (tolerant for floats) ----------
99 def equals(self, other: "MatrixDense", tol: float = 1e-12) -> bool:
100 """Element-wise comparison with tolerance for floats."""
101 self._check_same_shape(other)
102 for ra, rb in zip(self.rows, other.rows):
103 for a, b in zip(ra, rb):
104 if abs(a - b) > tol:
105 return False
106 return True
108 # ---------- Arithmetic ----------
110 def _check_same_shape(self, other: "MatrixDense") -> None:
111 if self.shape != other.shape:
112 raise InvalidValueError(f"shape mismatch: {self.shape} vs {other.shape}")
114 def __add__(self, other: "MatrixDense") -> "MatrixDense":
115 """Matrix addition."""
116 self._check_same_shape(other)
117 return MatrixDense([[a + b for a, b in zip(ra, rb)] for ra, rb in zip(self.rows, other.rows)])
119 def __sub__(self, other: "MatrixDense") -> "MatrixDense":
120 """Matrix subtraction."""
121 self._check_same_shape(other)
122 return MatrixDense([[a - b for a, b in zip(ra, rb)] for ra, rb in zip(self.rows, other.rows)])
124 @overload
125 def __mul__(self, other: Number) -> "MatrixDense":
126 ...
128 @overload
129 def __mul__(self, other: "MatrixDense") -> "MatrixDense":
130 ...
132 def __mul__(self, other: Union[Number, "MatrixDense"]) -> "MatrixDense":
133 """
134 Scalar or matrix multiplication.
136 - If `other` is a number: scale this matrix.
137 - If `other` is a MatrixDense: matrix-matrix product.
138 """
139 if isinstance(other, (int, float)):
140 return MatrixDense([[a * other for a in r] for r in self.rows])
141 if isinstance(other, MatrixDense):
142 a_rows, a_cols = self.shape
143 b_rows, b_cols = other.shape
144 if a_cols != b_rows:
145 raise InvalidValueError(f"incompatible shapes for matmul: {self.shape} * {other.shape}")
146 out = []
147 # naive triple loop (clear & correct)
148 for i in range(a_rows):
149 row = []
150 for j in range(b_cols):
151 s = 0.0
152 for k in range(a_cols):
153 s += self.rows[i][k] * other.rows[k][j]
154 row.append(s)
155 out.append(row)
156 return MatrixDense(out)
157 raise InvalidTypeError("unsupported operand for *: MatrixDense and non-(number|MatrixDense)")
159 def __rmul__(self, other: Number) -> "MatrixDense":
160 """Right scalar multiplication."""
161 return self.__mul__(other)
163 def matvec(self, v: Sequence[Number]) -> List[float]:
164 """
165 Matrix-vector product (Ax).
167 Parameters
168 ----------
169 v : Sequence[Number]
170 Vector of length equal to number of columns.
172 Returns
173 -------
174 list[float]
175 Resulting vector.
176 """
177 n_rows, n_cols = self.shape
178 if not isinstance(v, Sequence) or len(v) != n_cols:
179 raise InvalidValueError(f"vector length must be {n_cols}")
180 for x in v:
181 if not isinstance(x, (int, float)):
182 raise InvalidTypeError("vector elements must be numbers.")
183 out: List[float] = []
184 for i in range(n_rows):
185 s = 0.0
186 for j in range(n_cols):
187 s += self.rows[i][j] * v[j]
188 out.append(s)
189 return out
191 # ---------- Linear algebra helpers (small n) ----------
193 def T(self) -> "MatrixDense":
194 """Transpose."""
195 n_rows, n_cols = self.shape
196 return MatrixDense([[self.rows[i][j] for i in range(n_rows)] for j in range(n_cols)])
198 def det(self) -> float:
199 """
200 Determinant via Gaussian elimination with partial pivoting (O(n^3)).
202 Raises
203 ------
204 InvalidValueError
205 If matrix is not square.
206 """
207 n, m = self.shape
208 if n != m:
209 raise InvalidValueError("determinant is defined for square matrices.")
210 # Work on a mutable copy
211 a = [list(r) for r in self.rows]
212 sign = 1.0
213 det_val = 1.0
214 for i in range(n):
215 # pivot
216 pivot_row = max(range(i, n), key=lambda r: abs(a[r][i]))
217 if abs(a[pivot_row][i]) < 1e-15:
218 return 0.0
219 if pivot_row != i:
220 a[i], a[pivot_row] = a[pivot_row], a[i]
221 sign *= -1.0
222 pivot = a[i][i]
223 det_val *= pivot
224 # eliminate
225 for r in range(i + 1, n):
226 factor = a[r][i] / pivot
227 if factor == 0:
228 continue
229 for c in range(i, n):
230 a[r][c] -= factor * a[i][c]
231 return float(sign * det_val)
233 def inv(self) -> "MatrixDense":
234 """
235 Inverse via Gauss-Jordan elimination (O(n^3)).
237 Raises
238 ------
239 InvalidValueError
240 If matrix is not square or singular.
241 """
242 n, m = self.shape
243 if n != m:
244 raise InvalidValueError("inverse is defined for square matrices.")
245 # augmented [A | I]
246 A = [list(row) + [1.0 if i == j else 0.0 for j in range(n)] for i, row in enumerate(self.rows)]
248 # Gauss-Jordan with partial pivoting
249 for col in range(n):
250 # pivot
251 pivot_row = max(range(col, n), key=lambda r: abs(A[r][col]))
252 if abs(A[pivot_row][col]) < 1e-15:
253 raise InvalidValueError("matrix is singular; cannot invert.")
254 if pivot_row != col:
255 A[col], A[pivot_row] = A[pivot_row], A[col]
256 pivot = A[col][col]
257 # normalize pivot row
258 inv_p = 1.0 / pivot
259 for j in range(2 * n):
260 A[col][j] *= inv_p
261 # eliminate other rows
262 for r in range(n):
263 if r == col:
264 continue
265 factor = A[r][col]
266 if factor == 0:
267 continue
268 for j in range(2 * n):
269 A[r][j] -= factor * A[col][j]
271 inv_rows = [row[n:] for row in A]
272 return MatrixDense(inv_rows)