libpromeki main
PROfessional MEdia toolKIt
 
Loading...
Searching...
No Matches
matrix.h
Go to the documentation of this file.
1
8#pragma once
9
10#include <cmath>
11#include <promeki/core/namespace.h>
12#include <promeki/core/array.h>
13#include <promeki/core/error.h>
14
15PROMEKI_NAMESPACE_BEGIN
16
35template <typename T, size_t W, size_t H>
36class Matrix {
37 public:
39 using RowDataType = Array<T, W>;
40
42 using DataType = Array<RowDataType, H>;
43
48 static Matrix<T, W, H> identity() {
49 Matrix<T, W, H> result;
50 for (size_t i = 0; i < H; ++i) {
51 for (size_t j = 0; j < W; ++j) {
52 result[i][j] = (i == j) ? static_cast<T>(1) : static_cast<T>(0);
53 }
54 }
55 return result;
56 }
57
65 static Matrix<T, W, H> rotationMatrix(T radians, size_t dim, Error *err = nullptr) {
66 static_assert(W == H, "Matrix can only be calculated for square matrices");
67 static_assert(W >= 2, "Matrix size must be at least 2x2 for rotation matrix generation.");
68 if(dim >= W - 1) {
69 if(err != nullptr) *err = Error::InvalidDimension;
70 return Matrix<T, W, H>();
71 }
72
73 Matrix<T, W, H> result;
74 for (size_t i = 0; i < W; ++i) result[i][i] = static_cast<T>(1);
75 result[dim][dim] = std::cos(radians);
76 result[dim][dim + 1] = -std::sin(radians);
77 result[dim + 1][dim] = std::sin(radians);
78 result[dim + 1][dim + 1] = std::cos(radians);
79 return result;
80 }
81
83 Matrix() : d{} {}
84
89 Matrix(const DataType &data) : d(data) {}
90
96 RowDataType &operator[](size_t index) {
97 return d[index];
98 }
99
105 const RowDataType &operator[](size_t index) const {
106 return d[index];
107 }
108
110 DataType& data() {
111 return d;
112 }
113
115 const DataType& data() const {
116 return d;
117 }
118
120 size_t width() const {
121 return W;
122 }
123
125 size_t height() const {
126 return H;
127 }
128
130 bool isSquare() const {
131 return W == H;
132 }
133
138 Matrix<T, H, W> transpose() const {
139 Matrix<T, H, W> result;
140 for (size_t i = 0; i < H; ++i) {
141 for (size_t j = 0; j < W; ++j) {
142 result[j][i] = d[i][j];
143 }
144 }
145 return result;
146 }
147
153 Matrix<T, W, H> operator+(const Matrix<T, W, H>& other) const {
154 Matrix<T, W, H> result;
155 for (size_t i = 0; i < H; ++i) {
156 for (size_t j = 0; j < W; ++j) {
157 result[i][j] = d[i][j] + other[i][j];
158 }
159 }
160 return result;
161 }
162
168 Matrix<T, W, H> operator-(const Matrix<T, W, H>& other) const {
169 Matrix<T, W, H> result;
170 for (size_t i = 0; i < H; ++i) {
171 for (size_t j = 0; j < W; ++j) {
172 result[i][j] = d[i][j] - other[i][j];
173 }
174 }
175 return result;
176 }
177
183 Matrix<T, W, H> operator*(const T& scalar) const {
184 Matrix<T, W, H> result;
185 for (size_t i = 0; i < H; ++i) {
186 for (size_t j = 0; j < W; ++j) {
187 result[i][j] = d[i][j] * scalar;
188 }
189 }
190 return result;
191 }
192
199 template<size_t K> Matrix<T, W, K> operator*(const Matrix<T, H, K>& other) const {
200 Matrix<T, W, K> result;
201 for (size_t i = 0; i < W; ++i) {
202 for (size_t j = 0; j < K; ++j) {
203 T sum = 0;
204 for (size_t k = 0; k < H; ++k) {
205 sum += d[i][k] * other[k][j];
206 }
207 result[i][j] = sum;
208 }
209 }
210 return result;
211 }
212
218 Matrix<T, W, H> operator/(const T& scalar) const {
219 Matrix<T, W, H> result;
220 for (size_t i = 0; i < H; ++i) {
221 for (size_t j = 0; j < W; ++j) {
222 result[i][j] = d[i][j] / scalar;
223 }
224 }
225 return result;
226 }
227
233 void LUdecomposition(Matrix<T, W, H>& L, Matrix<T, W, H>& U) const {
234 for (size_t i = 0; i < H; ++i) {
235 for (size_t j = 0; j < W; ++j) {
236 if (j < i) {
237 L[j][i] = 0;
238 } else {
239 L[j][i] = (*this)[j][i];
240 for (size_t k = 0; k < i; ++k) {
241 L[j][i] -= L[j][k] * U[k][i];
242 }
243 }
244 if (i == j) {
245 U[i][i] = 1;
246 }
247 if (i < j) {
248 U[i][j] = 0;
249 } else {
250 U[i][j] = (*this)[i][j] / L[i][i];
251 for (size_t k = 0; k < i; ++k) {
252 U[i][j] -= ((L[i][k] * U[k][j])) / L[i][i];
253 }
254 }
255 }
256 }
257 }
258
267 T determinant() const {
268 static_assert(W == H, "Determinant can only be calculated for square matrices");
269
270 if constexpr (W == 2) {
271 return d[0][0] * d[1][1] - d[0][1] * d[1][0];
272 } else if constexpr (W == 3) {
273 return d[0][0] * (d[1][1] * d[2][2] - d[1][2] * d[2][1]) -
274 d[0][1] * (d[1][0] * d[2][2] - d[1][2] * d[2][0]) +
275 d[0][2] * (d[1][0] * d[2][1] - d[1][1] * d[2][0]);
276 }
277 Matrix<T, W, H> L = Matrix<T, W, H>::identity();
278 Matrix<T, W, H> U = Matrix<T, W, H>::identity();
279
280 LUdecomposition(L, U);
281
282 T det = 1;
283 for (size_t i = 0; i < H; ++i) {
284 det *= L[i][i] * U[i][i];
285 }
286 return det;
287 }
288
299 Matrix<T, W, H> inverse(Error *err = nullptr) const {
300 static_assert(W == H, "Inverse can only be calculated for square matrices");
301
302 T det = determinant();
303 if(det == 0) {
304 if(err != nullptr) *err = Error::SingularMatrix;
305 return Matrix<T, W, H>();
306 }
307
308 Matrix<T, W, H> result;
309 if constexpr (W == 2) {
310 result[0][0] = d[1][1] / det;
311 result[0][1] = -d[0][1] / det;
312 result[1][0] = -d[1][0] / det;
313 result[1][1] = d[0][0] / det;
314 } else if constexpr (W == 3) {
315 result[0][0] = (d[1][1] * d[2][2] - d[1][2] * d[2][1]) / det;
316 result[0][1] = (d[0][2] * d[2][1] - d[0][1] * d[2][2]) / det;
317 result[0][2] = (d[0][1] * d[1][2] - d[0][2] * d[1][1]) / det;
318 result[1][0] = (d[1][2] * d[2][0] - d[1][0] * d[2][2]) / det;
319 result[1][1] = (d[0][0] * d[2][2] - d[0][2] * d[2][0]) / det;
320 result[1][2] = (d[0][2] * d[1][0] - d[0][0] * d[1][2]) / det;
321 result[2][0] = (d[1][0] * d[2][1] - d[1][1] * d[2][0]) / det;
322 result[2][1] = (d[0][1] * d[2][0] - d[0][0] * d[2][1]) / det;
323 result[2][2] = (d[0][0] * d[1][1] - d[0][1] * d[1][0]) / det;
324 } else {
325 // Fallback to Gaussian method
326 result = Matrix<T, W, H>::identity();
327 Matrix<T, W, H> temp = *this;
328
329 for (size_t i = 0; i < W; ++i) {
330 // Find the pivot row
331 size_t pivot = i;
332 for (size_t j = i + 1; j < H; ++j) {
333 if (std::abs(temp[j][i]) > std::abs(temp[pivot][i])) {
334 pivot = j;
335 }
336 }
337
338 // Swap the rows
339 if (pivot != i) {
340 temp[i].swap(temp[pivot]);
341 result[i].swap(result[pivot]);
342 }
343
344 // Divide the pivot row by the pivot element
345 T pivotValue = temp[i][i];
346 if(pivotValue == 0) {
347 if(err != nullptr) *err = Error::SingularMatrix;
348 return Matrix<T, W, H>();
349 }
350
351 for (size_t j = 0; j < W; ++j) {
352 temp[i][j] /= pivotValue;
353 result[i][j] /= pivotValue;
354 }
355
356 // Use the pivot row to eliminate the other rows
357 for (size_t j = 0; j < H; ++j) {
358 if (j != i) {
359 T factor = temp[j][i];
360 for (size_t k = 0; k < W; ++k) {
361 temp[j][k] -= factor * temp[i][k];
362 result[j][k] -= factor * result[i][k];
363 }
364 }
365 }
366 }
367 }
368 return result;
369 }
370
380 template<size_t L> T dot(const Matrix<T, L, 1>& other) const {
381 static_assert(W == 1 || H == 1, "Both matrices must have a single row or a single column");
382 static_assert(L == H || L == W, "The matrices must have the same number of elements");
383
384 size_t size = (W == 1) ? H : W;
385 T result = 0;
386
387 if constexpr (W == 1 && other.width() == 1) { // Both are column vectors
388 for (size_t i = 0; i < size; ++i) {
389 result += (*this)[i][0] * other[i][0];
390 }
391 } else if constexpr (H == 1 && other.height() == 1) { // Both are row vectors
392 for (size_t i = 0; i < size; ++i) {
393 result += (*this)[0][i] * other[0][i];
394 }
395 }
396 return result;
397 }
398
405 template<typename Func> Matrix<T, W, H> apply(Func&& func) const {
406 Matrix<T, W, H> result;
407 for (size_t i = 0; i < H; ++i) {
408 for (size_t j = 0; j < W; ++j) {
409 result[i][j] = func((*this)[i][j]);
410 }
411 }
412 return result;
413 }
414
419 T frobeniusNorm() const {
420 T sum = 0;
421 for (size_t i = 0; i < H; ++i) {
422 for (size_t j = 0; j < W; ++j) {
423 sum += (*this)[i][j] * (*this)[i][j];
424 }
425 }
426 return std::sqrt(sum);
427 }
428
433 T trace() const {
434 static_assert(W == H, "Trace can only be calculated for square matrices");
435 T sum = 0;
436 for (size_t i = 0; i < H; ++i) {
437 sum += (*this)[i][i];
438 }
439 return sum;
440 }
441
446 T sum() const {
447 T total = 0;
448 for (size_t i = 0; i < H; ++i) {
449 for (size_t j = 0; j < W; ++j) {
450 total += (*this)[i][j];
451 }
452 }
453 return total;
454 }
455
461 Matrix<T, W, H> hadamardProduct(const Matrix<T, W, H>& other) const {
462 Matrix<T, W, H> result;
463 for (size_t i = 0; i < H; ++i) {
464 for (size_t j = 0; j < W; ++j) {
465 result[i][j] = (*this)[i][j] * other[i][j];
466 }
467 }
468 return result;
469 }
470
480 template<size_t L> Matrix<T, H, L> outer_product(const Matrix<T, L, 1>& other) const {
481 static_assert(W == 1, "The first matrix must have a single column");
482 Matrix<T, H, L> result;
483 for (size_t i = 0; i < H; ++i) {
484 for (size_t j = 0; j < L; ++j) {
485 result[i][j] = (*this)[i][0] * other[j][0];
486 }
487 }
488 return result;
489 }
490
495 Matrix<T, 1, H> rowSum() const {
496 Matrix<T, 1, H> result;
497 for (size_t i = 0; i < H; ++i) {
498 T total = 0;
499 for (size_t j = 0; j < W; ++j) {
500 total += (*this)[i][j];
501 }
502 result[0][i] = total;
503 }
504 return result;
505 }
506
511 Matrix<T, W, 1> colSum() const {
512 Matrix<T, W, 1> result;
513 for (size_t j = 0; j < W; ++j) {
514 T total = 0;
515 for (size_t i = 0; i < H; ++i) {
516 total += (*this)[i][j];
517 }
518 result[j][0] = total;
519 }
520 return result;
521 }
522
527 Matrix<T, W, 1> diagonal() const {
528 static_assert(W == H, "Diagonal can only be calculated for square matrices");
529
530 Matrix<T, W, 1> result;
531 for (size_t i = 0; i < W; ++i) {
532 result[i][0] = (*this)[i][i];
533 }
534 return result;
535 }
536
541 Matrix<double, 1, H> rowMean() const {
542 Matrix<double, 1, H> result;
543 for (size_t i = 0; i < H; ++i) {
544 double total = 0;
545 for (size_t j = 0; j < W; ++j) {
546 total += (*this)[i][j];
547 }
548 result[0][i] = total / W;
549 }
550 return result;
551 }
552
557 Matrix<double, W, 1> colMean() const {
558 Matrix<double, W, 1> result;
559 for (size_t j = 0; j < W; ++j) {
560 double total = 0;
561 for (size_t i = 0; i < H; ++i) {
562 total += (*this)[i][j];
563 }
564 result[j][0] = total / H;
565 }
566 return result;
567 }
568
573 Matrix<T, H, W> rotateCW() const {
574 Matrix<T, H, W> result;
575 for (size_t i = 0; i < H; ++i) {
576 for (size_t j = 0; j < W; ++j) {
577 result[i][j] = (*this)[H - j - 1][i];
578 }
579 }
580 return result;
581 }
582
587 Matrix<T, H, W> rotateCCW() const {
588 Matrix<T, H, W> result;
589 for (size_t i = 0; i < H; ++i) {
590 for (size_t j = 0; j < W; ++j) {
591 result[i][j] = (*this)[j][W - i - 1];
592 }
593 }
594 return result;
595 }
596
597 private:
598 DataType d;
599};
600
601PROMEKI_NAMESPACE_END
602
Error
Lightweight error code wrapper for the promeki library.
Definition error.h:39
Error::InvalidDimension
@ InvalidDimension
Invalid dimension value.
Definition error.h:87
Error::SingularMatrix
@ SingularMatrix
Matrix is singular and cannot be inverted.
Definition error.h:54
List
Dynamic array container wrapping std::vector.
Definition list.h:40
List::swap
void swap(List< T > &other) noexcept
Swaps the list data with another list of the same type.
Definition list.h:530
Matrix
Generic fixed-size matrix with compile-time dimensions.
Definition matrix.h:36
Matrix::determinant
T determinant() const
Computes the determinant of a square matrix.
Definition matrix.h:267
Matrix::identity
static Matrix< T, W, H > identity()
Creates an identity matrix.
Definition matrix.h:48
Matrix::operator/
Matrix< T, W, H > operator/(const T &scalar) const
Scalar division.
Definition matrix.h:218
Matrix::rotateCCW
Matrix< T, H, W > rotateCCW() const
Rotates the matrix 90 degrees counter-clockwise.
Definition matrix.h:587
Matrix::apply
Matrix< T, W, H > apply(Func &&func) const
Applies a function to each element of the matrix.
Definition matrix.h:405
Matrix::colMean
Matrix< double, W, 1 > colMean() const
Computes the mean of each column.
Definition matrix.h:557
Matrix::rowMean
Matrix< double, 1, H > rowMean() const
Computes the mean of each row.
Definition matrix.h:541
Matrix::data
DataType & data()
Returns a mutable reference to the underlying data.
Definition matrix.h:110
Matrix::operator[]
RowDataType & operator[](size_t index)
Accesses a row by index.
Definition matrix.h:96
Matrix::LUdecomposition
void LUdecomposition(Matrix< T, W, H > &L, Matrix< T, W, H > &U) const
Computes the LU decomposition of this matrix.
Definition matrix.h:233
Matrix::height
size_t height() const
Returns the number of rows.
Definition matrix.h:125
Matrix::dot
T dot(const Matrix< T, L, 1 > &other) const
Computes the dot product between two matrices treated as vectors.
Definition matrix.h:380
Matrix::diagonal
Matrix< T, W, 1 > diagonal() const
Extracts the diagonal elements of a square matrix.
Definition matrix.h:527
Matrix::data
const DataType & data() const
Returns a const reference to the underlying data.
Definition matrix.h:115
Matrix::transpose
Matrix< T, H, W > transpose() const
Returns the transpose of this matrix.
Definition matrix.h:138
Matrix::outer_product
Matrix< T, H, L > outer_product(const Matrix< T, L, 1 > &other) const
Computes the outer product of two column vectors.
Definition matrix.h:480
Matrix::Matrix
Matrix(const DataType &data)
Constructs a matrix from existing data.
Definition matrix.h:89
Matrix::operator*
Matrix< T, W, K > operator*(const Matrix< T, H, K > &other) const
Matrix multiplication.
Definition matrix.h:199
Matrix::colSum
Matrix< T, W, 1 > colSum() const
Computes the sum of each column.
Definition matrix.h:511
Matrix::operator*
Matrix< T, W, H > operator*(const T &scalar) const
Scalar multiplication.
Definition matrix.h:183
Matrix::sum
T sum() const
Computes the sum of all elements in the matrix.
Definition matrix.h:446
Matrix::width
size_t width() const
Returns the number of columns.
Definition matrix.h:120
Matrix::isSquare
bool isSquare() const
Returns true if the matrix is square (W == H).
Definition matrix.h:130
Matrix::operator+
Matrix< T, W, H > operator+(const Matrix< T, W, H > &other) const
Element-wise matrix addition.
Definition matrix.h:153
Matrix::trace
T trace() const
Computes the trace (sum of diagonal elements) of a square matrix.
Definition matrix.h:433
Matrix::DataType
Array< RowDataType, H > DataType
Type alias for the entire matrix storage (array of rows).
Definition matrix.h:42
Matrix::hadamardProduct
Matrix< T, W, H > hadamardProduct(const Matrix< T, W, H > &other) const
Computes the Hadamard (element-wise) product of two matrices.
Definition matrix.h:461
Matrix::operator[]
const RowDataType & operator[](size_t index) const
Accesses a row by index (const).
Definition matrix.h:105
Matrix::inverse
Matrix< T, W, H > inverse(Error *err=nullptr) const
Computes the inverse of a square matrix.
Definition matrix.h:299
Matrix::rowSum
Matrix< T, 1, H > rowSum() const
Computes the sum of each row.
Definition matrix.h:495
Matrix::frobeniusNorm
T frobeniusNorm() const
Computes the Frobenius norm of the matrix.
Definition matrix.h:419
Matrix::operator-
Matrix< T, W, H > operator-(const Matrix< T, W, H > &other) const
Element-wise matrix subtraction.
Definition matrix.h:168
Matrix::rotationMatrix
static Matrix< T, W, H > rotationMatrix(T radians, size_t dim, Error *err=nullptr)
Creates a rotation matrix for the given angle and dimension pair.
Definition matrix.h:65
Matrix::Matrix
Matrix()
Constructs a zero-initialized matrix.
Definition matrix.h:83
Matrix::rotateCW
Matrix< T, H, W > rotateCW() const
Rotates the matrix 90 degrees clockwise.
Definition matrix.h:573
PROMEKI_NAMESPACE_BEGIN
#define PROMEKI_NAMESPACE_BEGIN
Starts a promeki namespace block.
Definition namespace.h:14
PROMEKI_NAMESPACE_END
#define PROMEKI_NAMESPACE_END
Ends a promeki namespace block.
Definition namespace.h:19