libpromeki 1.0.0-alpha
PROfessional MEdia toolKIt
 
Loading...
Searching...
No Matches
matrix.h
Go to the documentation of this file.
1
8#pragma once
9
10
11#include <promeki/config.h>
12#if PROMEKI_ENABLE_CORE
13#include <cmath>
14#include <promeki/namespace.h>
15#include <promeki/array.h>
16#include <promeki/error.h>
17
18PROMEKI_NAMESPACE_BEGIN
19
43template <typename T, size_t W, size_t H> class Matrix {
44 public:
46 using RowDataType = Array<T, W>;
47
49 using DataType = Array<RowDataType, H>;
50
55 static Matrix<T, W, H> identity() {
56 Matrix<T, W, H> result;
57 for (size_t i = 0; i < H; ++i) {
58 for (size_t j = 0; j < W; ++j) {
59 result[i][j] = (i == j) ? static_cast<T>(1) : static_cast<T>(0);
60 }
61 }
62 return result;
63 }
64
72 static Matrix<T, W, H> rotationMatrix(T radians, size_t dim, Error *err = nullptr) {
73 static_assert(W == H, "Matrix can only be calculated for square matrices");
74 static_assert(W >= 2, "Matrix size must be at least 2x2 for rotation matrix generation.");
75 if (dim >= W - 1) {
76 if (err != nullptr) *err = Error::InvalidDimension;
77 return Matrix<T, W, H>();
78 }
79
80 Matrix<T, W, H> result;
81 for (size_t i = 0; i < W; ++i) result[i][i] = static_cast<T>(1);
82 result[dim][dim] = std::cos(radians);
83 result[dim][dim + 1] = -std::sin(radians);
84 result[dim + 1][dim] = std::sin(radians);
85 result[dim + 1][dim + 1] = std::cos(radians);
86 return result;
87 }
88
90 Matrix() : d{} {}
91
96 Matrix(const DataType &data) : d(data) {}
97
103 RowDataType &operator[](size_t index) { return d[index]; }
104
110 const RowDataType &operator[](size_t index) const { return d[index]; }
111
113 DataType &data() { return d; }
114
116 const DataType &data() const { return d; }
117
119 size_t width() const { return W; }
120
122 size_t height() const { return H; }
123
125 bool isSquare() const { return W == H; }
126
131 Matrix<T, H, W> transpose() const {
132 Matrix<T, H, W> result;
133 for (size_t i = 0; i < H; ++i) {
134 for (size_t j = 0; j < W; ++j) {
135 result[j][i] = d[i][j];
136 }
137 }
138 return result;
139 }
140
146 Matrix<T, W, H> operator+(const Matrix<T, W, H> &other) const {
147 Matrix<T, W, H> result;
148 for (size_t i = 0; i < H; ++i) {
149 for (size_t j = 0; j < W; ++j) {
150 result[i][j] = d[i][j] + other[i][j];
151 }
152 }
153 return result;
154 }
155
161 Matrix<T, W, H> operator-(const Matrix<T, W, H> &other) const {
162 Matrix<T, W, H> result;
163 for (size_t i = 0; i < H; ++i) {
164 for (size_t j = 0; j < W; ++j) {
165 result[i][j] = d[i][j] - other[i][j];
166 }
167 }
168 return result;
169 }
170
176 Matrix<T, W, H> operator*(const T &scalar) const {
177 Matrix<T, W, H> result;
178 for (size_t i = 0; i < H; ++i) {
179 for (size_t j = 0; j < W; ++j) {
180 result[i][j] = d[i][j] * scalar;
181 }
182 }
183 return result;
184 }
185
192 template <size_t K> Matrix<T, W, K> operator*(const Matrix<T, H, K> &other) const {
193 Matrix<T, W, K> result;
194 for (size_t i = 0; i < W; ++i) {
195 for (size_t j = 0; j < K; ++j) {
196 T sum = 0;
197 for (size_t k = 0; k < H; ++k) {
198 sum += d[i][k] * other[k][j];
199 }
200 result[i][j] = sum;
201 }
202 }
203 return result;
204 }
205
211 Matrix<T, W, H> operator/(const T &scalar) const {
212 Matrix<T, W, H> result;
213 for (size_t i = 0; i < H; ++i) {
214 for (size_t j = 0; j < W; ++j) {
215 result[i][j] = d[i][j] / scalar;
216 }
217 }
218 return result;
219 }
220
226 void LUdecomposition(Matrix<T, W, H> &L, Matrix<T, W, H> &U) const {
227 for (size_t i = 0; i < H; ++i) {
228 for (size_t j = 0; j < W; ++j) {
229 if (j < i) {
230 L[j][i] = 0;
231 } else {
232 L[j][i] = (*this)[j][i];
233 for (size_t k = 0; k < i; ++k) {
234 L[j][i] -= L[j][k] * U[k][i];
235 }
236 }
237 if (i == j) {
238 U[i][i] = 1;
239 }
240 if (i < j) {
241 U[i][j] = 0;
242 } else {
243 U[i][j] = (*this)[i][j] / L[i][i];
244 for (size_t k = 0; k < i; ++k) {
245 U[i][j] -= ((L[i][k] * U[k][j])) / L[i][i];
246 }
247 }
248 }
249 }
250 }
251
260 T determinant() const {
261 static_assert(W == H, "Determinant can only be calculated for square matrices");
262
263 if constexpr (W == 2) {
264 return d[0][0] * d[1][1] - d[0][1] * d[1][0];
265 } else if constexpr (W == 3) {
266 return d[0][0] * (d[1][1] * d[2][2] - d[1][2] * d[2][1]) -
267 d[0][1] * (d[1][0] * d[2][2] - d[1][2] * d[2][0]) +
268 d[0][2] * (d[1][0] * d[2][1] - d[1][1] * d[2][0]);
269 }
270 Matrix<T, W, H> L = Matrix<T, W, H>::identity();
271 Matrix<T, W, H> U = Matrix<T, W, H>::identity();
272
273 LUdecomposition(L, U);
274
275 T det = 1;
276 for (size_t i = 0; i < H; ++i) {
277 det *= L[i][i] * U[i][i];
278 }
279 return det;
280 }
281
292 Matrix<T, W, H> inverse(Error *err = nullptr) const {
293 static_assert(W == H, "Inverse can only be calculated for square matrices");
294
295 T det = determinant();
296 if (det == 0) {
297 if (err != nullptr) *err = Error::SingularMatrix;
298 return Matrix<T, W, H>();
299 }
300
301 Matrix<T, W, H> result;
302 if constexpr (W == 2) {
303 result[0][0] = d[1][1] / det;
304 result[0][1] = -d[0][1] / det;
305 result[1][0] = -d[1][0] / det;
306 result[1][1] = d[0][0] / det;
307 } else if constexpr (W == 3) {
308 result[0][0] = (d[1][1] * d[2][2] - d[1][2] * d[2][1]) / det;
309 result[0][1] = (d[0][2] * d[2][1] - d[0][1] * d[2][2]) / det;
310 result[0][2] = (d[0][1] * d[1][2] - d[0][2] * d[1][1]) / det;
311 result[1][0] = (d[1][2] * d[2][0] - d[1][0] * d[2][2]) / det;
312 result[1][1] = (d[0][0] * d[2][2] - d[0][2] * d[2][0]) / det;
313 result[1][2] = (d[0][2] * d[1][0] - d[0][0] * d[1][2]) / det;
314 result[2][0] = (d[1][0] * d[2][1] - d[1][1] * d[2][0]) / det;
315 result[2][1] = (d[0][1] * d[2][0] - d[0][0] * d[2][1]) / det;
316 result[2][2] = (d[0][0] * d[1][1] - d[0][1] * d[1][0]) / det;
317 } else {
318 // Fallback to Gaussian method
319 result = Matrix<T, W, H>::identity();
320 Matrix<T, W, H> temp = *this;
321
322 for (size_t i = 0; i < W; ++i) {
323 // Find the pivot row
324 size_t pivot = i;
325 for (size_t j = i + 1; j < H; ++j) {
326 if (std::abs(temp[j][i]) > std::abs(temp[pivot][i])) {
327 pivot = j;
328 }
329 }
330
331 // Swap the rows
332 if (pivot != i) {
333 temp[i].swap(temp[pivot]);
334 result[i].swap(result[pivot]);
335 }
336
337 // Divide the pivot row by the pivot element
338 T pivotValue = temp[i][i];
339 if (pivotValue == 0) {
340 if (err != nullptr) *err = Error::SingularMatrix;
341 return Matrix<T, W, H>();
342 }
343
344 for (size_t j = 0; j < W; ++j) {
345 temp[i][j] /= pivotValue;
346 result[i][j] /= pivotValue;
347 }
348
349 // Use the pivot row to eliminate the other rows
350 for (size_t j = 0; j < H; ++j) {
351 if (j != i) {
352 T factor = temp[j][i];
353 for (size_t k = 0; k < W; ++k) {
354 temp[j][k] -= factor * temp[i][k];
355 result[j][k] -= factor * result[i][k];
356 }
357 }
358 }
359 }
360 }
361 return result;
362 }
363
373 template <size_t L> T dot(const Matrix<T, L, 1> &other) const {
374 static_assert(W == 1 || H == 1, "Both matrices must have a single row or a single column");
375 static_assert(L == H || L == W, "The matrices must have the same number of elements");
376
377 size_t size = (W == 1) ? H : W;
378 T result = 0;
379
380 if constexpr (W == 1 && other.width() == 1) { // Both are column vectors
381 for (size_t i = 0; i < size; ++i) {
382 result += (*this)[i][0] * other[i][0];
383 }
384 } else if constexpr (H == 1 && other.height() == 1) { // Both are row vectors
385 for (size_t i = 0; i < size; ++i) {
386 result += (*this)[0][i] * other[0][i];
387 }
388 }
389 return result;
390 }
391
398 template <typename Func> Matrix<T, W, H> apply(Func &&func) const {
399 Matrix<T, W, H> result;
400 for (size_t i = 0; i < H; ++i) {
401 for (size_t j = 0; j < W; ++j) {
402 result[i][j] = func((*this)[i][j]);
403 }
404 }
405 return result;
406 }
407
412 T frobeniusNorm() const {
413 T sum = 0;
414 for (size_t i = 0; i < H; ++i) {
415 for (size_t j = 0; j < W; ++j) {
416 sum += (*this)[i][j] * (*this)[i][j];
417 }
418 }
419 return std::sqrt(sum);
420 }
421
426 T trace() const {
427 static_assert(W == H, "Trace can only be calculated for square matrices");
428 T sum = 0;
429 for (size_t i = 0; i < H; ++i) {
430 sum += (*this)[i][i];
431 }
432 return sum;
433 }
434
439 T sum() const {
440 T total = 0;
441 for (size_t i = 0; i < H; ++i) {
442 for (size_t j = 0; j < W; ++j) {
443 total += (*this)[i][j];
444 }
445 }
446 return total;
447 }
448
454 Matrix<T, W, H> hadamardProduct(const Matrix<T, W, H> &other) const {
455 Matrix<T, W, H> result;
456 for (size_t i = 0; i < H; ++i) {
457 for (size_t j = 0; j < W; ++j) {
458 result[i][j] = (*this)[i][j] * other[i][j];
459 }
460 }
461 return result;
462 }
463
473 template <size_t L> Matrix<T, H, L> outer_product(const Matrix<T, L, 1> &other) const {
474 static_assert(W == 1, "The first matrix must have a single column");
475 Matrix<T, H, L> result;
476 for (size_t i = 0; i < H; ++i) {
477 for (size_t j = 0; j < L; ++j) {
478 result[i][j] = (*this)[i][0] * other[j][0];
479 }
480 }
481 return result;
482 }
483
488 Matrix<T, 1, H> rowSum() const {
489 Matrix<T, 1, H> result;
490 for (size_t i = 0; i < H; ++i) {
491 T total = 0;
492 for (size_t j = 0; j < W; ++j) {
493 total += (*this)[i][j];
494 }
495 result[0][i] = total;
496 }
497 return result;
498 }
499
504 Matrix<T, W, 1> colSum() const {
505 Matrix<T, W, 1> result;
506 for (size_t j = 0; j < W; ++j) {
507 T total = 0;
508 for (size_t i = 0; i < H; ++i) {
509 total += (*this)[i][j];
510 }
511 result[j][0] = total;
512 }
513 return result;
514 }
515
520 Matrix<T, W, 1> diagonal() const {
521 static_assert(W == H, "Diagonal can only be calculated for square matrices");
522
523 Matrix<T, W, 1> result;
524 for (size_t i = 0; i < W; ++i) {
525 result[i][0] = (*this)[i][i];
526 }
527 return result;
528 }
529
534 Matrix<double, 1, H> rowMean() const {
535 Matrix<double, 1, H> result;
536 for (size_t i = 0; i < H; ++i) {
537 double total = 0;
538 for (size_t j = 0; j < W; ++j) {
539 total += (*this)[i][j];
540 }
541 result[0][i] = total / W;
542 }
543 return result;
544 }
545
550 Matrix<double, W, 1> colMean() const {
551 Matrix<double, W, 1> result;
552 for (size_t j = 0; j < W; ++j) {
553 double total = 0;
554 for (size_t i = 0; i < H; ++i) {
555 total += (*this)[i][j];
556 }
557 result[j][0] = total / H;
558 }
559 return result;
560 }
561
566 Matrix<T, H, W> rotateCW() const {
567 Matrix<T, H, W> result;
568 for (size_t i = 0; i < H; ++i) {
569 for (size_t j = 0; j < W; ++j) {
570 result[i][j] = (*this)[H - j - 1][i];
571 }
572 }
573 return result;
574 }
575
580 Matrix<T, H, W> rotateCCW() const {
581 Matrix<T, H, W> result;
582 for (size_t i = 0; i < H; ++i) {
583 for (size_t j = 0; j < W; ++j) {
584 result[i][j] = (*this)[j][W - i - 1];
585 }
586 }
587 return result;
588 }
589
590 private:
591 DataType d;
592};
593
594PROMEKI_NAMESPACE_END
595
596#endif // PROMEKI_ENABLE_CORE