#pragma once #include "global.h" #include #include namespace pcpo { /// Matrix. Row and column are indexed beginning at 0 template class Matrix { protected: std::vector> vectors; int width; int height; public: //MARK: - Initializers /// Creates a matrix with dimensions height x width and initalizes its values to `value` /// @param width Width of the matrix /// @param height Height of the matrix /// @param value Initial value for each element Matrix(int height, int width, T value) { this->width = width; this->height = height; this->vectors.reserve(width); for (int i = 0; i < height; i++) { std::vector vector(width,value); vectors.push_back(vector); } }; /// Creates a matrix with dimensions height x width and initalizes its values to 0 /// @param width Width of the matrix /// @param height Height of the matrix Matrix(int height, int width) : Matrix(height, width, 0) {}; Matrix() = default; Matrix(Matrix const& matrix) = default; /// Creates an identity matrix with dimension eye x eye /// @param eye dimension of the matrix Matrix(int eye){ this->width = eye; this->height = eye; this->vectors.reserve(width); for (int i = 0; i < height; i++) { std::vector vector(width,0); vector[i] = 1; vectors.push_back(vector); } }; /// Creates a matrix from a 2D vector /// @param vectors 2D vector containing columns with rows Matrix(std::vector> const &vectors) { assert(all_of(vectors.begin(), vectors.end(), [&vectors](std::vector vec){ return vec.size() == vectors.front().size(); })); this->width = vectors.empty() ? 0 : vectors.front().size(); this->height = vectors.size(); this->vectors = vectors; }; /// Creates a vector from a std::vector /// @param vector the vector Matrix(std::vector const &vector) { std::vector> vectors = {vector}; this->width = vector.size(); this->height = vector.empty() ? 0 : 1; this->vectors = vectors; }; Matrix(std::vector const &values, int rows, int columns) { assert(int(values.size()) == rows * columns); std::vector> result; result.reserve(rows); for (int row = 0; row < rows; row++) { std::vector rowVector; rowVector.reserve(columns); for (int column = 0; column < columns; column++) { rowVector.push_back(values[row * rows + column]); } result.push_back(rowVector); } this->vectors = result; this->width = columns; this->height = rows; }; // MARK: - Properties /// The height of the matrix (number of rows) int getHeight() const { return height; }; /// The width of the matrix (number of columns) int getWidth() const { return width; }; /// Returns true when the matrix is empty bool empty() const { return getWidth() == 0 && getHeight() == 0; }; // MARK: - Matrix operations /// Transposes the matrix Matrix transpose() const { Matrix result = Matrix(width, height); for (int i = 0; i < height; i++) { for (int j = 0; j < width; j++) { result.value(j,i) = value(i,j); } } return result; }; /// Transforms the matrix to reduced row echelon form Matrix echelonForm() const { Matrix result = Matrix(*this); int pivot = 0; for (int row = 0; row < height; row++) { if (pivot >= width) { return result; } int i = row; while (result.value(i,pivot) == 0) { if (++i >= height) { i = row; if (++pivot >= width) { return result; } } } result.swap_rows(i, row); result.divide_row(row, result.value(row,pivot)); for (i = 0; i < height; i++) { if (i != row) { result.add_multiple_row(i, row, -result.value(i,pivot)); } } } return result; }; /// The rank of the matrix int getRank() const { Matrix e = echelonForm(); int rank = 0; for (int row = 0; row < height; row++) { for (int column = 0; column < width; column++) { if ((e.value(row,column) == 0) && (column == width - 1)) { return rank; } else if (e.value(row,column) != 0) { break; } } rank++; } return rank; } /// Basis of the linear span of the column vectors static Matrix span(Matrix const& matrix, bool transposed = false) { std::vector> rows; // if matrix is already transposed don't do it again Matrix t = transposed ? matrix : matrix.transpose(); Matrix te = t.echelonForm(); int rank = te.getRank(); rows.reserve(rank); for (int row = 0; row < rank; row++) { rows.push_back(te.row(row)); } return Matrix(rows).transpose(); } /// Computes the null space for the column vectors static Matrix null(Matrix const& matrix); /// Converts the matrix to a 1D Vector by stacking the column vectors std::vector toVector() const { std::vector result; result.reserve(getWidth() * getHeight()); for (int column = 0; column < getWidth(); column++) { for (int row = 0; row < getHeight(); row++) { result.push_back(value(row,column)); } } return result; } /// Converts a 1D Vector to a Matrix with given dimensions /// @param rows number of rows /// @param columns number of columns Matrix reshape(int rows, int columns) const { assert(rows > 0 && columns > 0); // FIXME: Performance Matrix t = transpose(); return Matrix(t.vectors.front(), rows, columns).transpose(); }; std::vector> reshapeColumns(int height, int width) const { std::vector> result; for (int c = 0; c < getWidth(); c++) { result.push_back(Matrix(column(c), height, width).transpose()); } return result; } /// Returns the value at row i and column j /// @param row /// @param column T& value(int row, int column) { assert(row < getHeight() && column < getWidth()); return vectors[row][column]; }; /// Returns the value at row i and column j /// @param row /// @param column T const& value(int row, int column) const { assert(row < getHeight() && column < getWidth()); return vectors[row][column]; }; /// Returns a vector with the elements of the row at index i. The returned row can be modified. /// @param i Index of the row to return. std::vector& row(int i) { assert(i < getHeight()); return vectors[i]; }; std::vector row(int i) const { assert(i < getHeight()); return vectors[i]; }; /// Returns the column at index i. The returned column cannot be modified. /// @param i Index of the column to return std::vector column(int i) const { assert(i < getWidth()); std::vector row; row.reserve(width); for (std::vector const& x : vectors) { row.push_back(x[i]); } return row; } void setColumn(std::vector const& vector, int column) { assert(int(vector.size()) == height); for (int row = 0; row < height; row++) { value(row,column) = vector[row]; } } // MARK: - Operators T& operator()(int i, int j) { return value(i,j); }; Matrix& operator *=(Matrix const& rhs) { assert(width == rhs.height); Matrix result = Matrix(height,rhs.width); for (int i = 0; i < height; i++) { for (int k = 0; k < width; k++) { for (int j = 0; j < rhs.width; j++) { result.value(i,j) += value(i,k) * rhs.value(k,j); } } } this->vectors = result.vectors; this->width = result.width; this->height = result.height; return *this; }; Matrix& operator *=(T scalar) { for (int i = 0; i < height; i++) { for (int j = 0; j < width; j++) { value(i,j) *= scalar; } } return *this; }; Matrix& operator +=(Matrix const& rhs) { assert(rhs.width == width && rhs.height == height); for (int i=0;i& operator +=(T scalar) { for (int i = 0; i < height; i++) { for (int j = 0; j < width; j++) { value(i,j) += scalar; } } return *this; }; Matrix& operator -=(Matrix const& rhs) { assert(rhs.width == width && rhs.height == height); for (int i = 0; i < height; i++) { for (int j= 0; j < width; j++) { value(i,j) -= rhs.value(i,j); } } return *this; }; Matrix& operator -=(int scalar) { for (int i = 0; i < height; i++) { for (int j = 0; j < width; j++) { value(i,j) -= scalar; } } return *this; }; bool operator==(Matrix const& rhs) const { return rhs.vectors == rhs.vectors && width == rhs.width && height == rhs.height; }; void print() const { dbgs(4) << *this; } protected: // MARK: - Echelon helpers /// Swaps two rows /// @param a index of the first row /// @param b index of the second row void swap_rows(int a, int b) { vectors[a].swap(vectors[b]); }; /// Divides a row by a constant /// @param row index of the row to divide /// @param quotient quotient to divide the row by void divide_row(int row, T quotient) { for (int column = 0; column < width; column++) { value(row,column) /= quotient; } }; /// Adds a multiple of row b to a row a /// @param a Row to add a multiple of b to /// @param b Row to be added to row a /// @param factor Factor to multiply row b with when adding it to row a void add_multiple_row(int a, int b, T factor) { for (int column = 0; column < width; column++) { value(a,column) += value(b,column) * factor; } }; // MARK: - Utils /// Greates common divisor. /// @param lhs /// @param rhs static int ggT(int lhs, int rhs) { int h; if (lhs == 0) { return abs(rhs); } if (rhs == 0) { return abs(lhs); } do { h = lhs % rhs; lhs = rhs; rhs = h; } while (rhs != 0); return abs(lhs); }; /// Least common multiple. /// @param lhs /// @param rhs static int kgV(int lhs, int rhs) { return (lhs * rhs) / ggT(lhs, rhs); }; }; template inline Matrix operator*(Matrix lhs, Matrix const& rhs) { return lhs *= rhs; }; template inline Matrix operator*(Matrix lhs, T scalar) { return lhs *= scalar; }; template inline Matrix operator+(Matrix lhs, Matrix const& rhs) { return lhs += rhs; }; template inline Matrix operator+(Matrix lhs, T scalar) { return lhs += scalar; }; template inline Matrix operator-(Matrix lhs, Matrix const& rhs) { return lhs -= rhs; }; template inline Matrix operator-(Matrix lhs, T scalar) { return lhs -= scalar; }; template llvm::raw_ostream& operator<<(llvm::raw_ostream& os, Matrix const& matrix) { for (int row = 0; row < matrix.getHeight(); row++) { os << "[ "; for (int column = 0; column < matrix.getWidth(); column++) { if constexpr (std::is_floating_point_v) { if (column == matrix.getWidth() - 1) { os << llvm::format("%g", matrix.value(row,column)); } else { os << llvm::format("%-6g", matrix.value(row,column)); } } else { if (column == matrix.getWidth() - 1) { os << llvm::format("%d", matrix.value(row,column)); } else { os << llvm::format("%-6d", matrix.value(row,column)); } } } os << " ]\n"; } return os; }; }