diff --git a/.cache/clangd/index/Matrix.hpp.B9C1D36F0ACA1714.idx b/.cache/clangd/index/Matrix.hpp.B9C1D36F0ACA1714.idx new file mode 100644 index 0000000..34b906c Binary files /dev/null and b/.cache/clangd/index/Matrix.hpp.B9C1D36F0ACA1714.idx differ diff --git a/.cache/clangd/index/matrix-tests.cpp.88D512D6040A12A3.idx b/.cache/clangd/index/matrix-tests.cpp.88D512D6040A12A3.idx new file mode 100644 index 0000000..f7707e0 Binary files /dev/null and b/.cache/clangd/index/matrix-tests.cpp.88D512D6040A12A3.idx differ diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d163863 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +build/ \ No newline at end of file diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 0000000..d860e76 --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,43 @@ +{ + // Use IntelliSense to learn about possible attributes. + // Hover to view descriptions of existing attributes. + // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 + "version": "0.2.0", + "configurations": [ + { + "name": "Debug Matrix Unit Tests", + "type": "cppdbg", + "request": "launch", + "program": "${workspaceFolder}/build/unit-tests/matrix-tests", + "args": [], + "stopAtEntry": false, + "cwd": "${workspaceFolder}", + "environment": [], + "externalConsole": false, + "MIMode": "gdb", + "miDebuggerPath": "/usr/bin/gdb", // Adjust to your debugger path + "setupCommands": [ + { + "description": "Enable pretty-printing for gdb", + "text": "-enable-pretty-printing", + "ignoreFailures": true + } + ], + "preLaunchTask": "build_tests", // Task to compile unit tests + "internalConsoleOptions": "openOnSessionStart" + }, + { + "name": "Run Matrix Unit Tests", + "type": "cpp", + "request": "launch", + "program": "${workspaceFolder}/build/unit-tests/matrix-tests", + "args": [], + "stopAtEntry": false, + "cwd": "${workspaceFolder}", + "environment": [], + "externalConsole": false, + "preLaunchTask": "build_tests", // Compile unit tests before running + "internalConsoleOptions": "openOnSessionStart" + } + ] +} \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..71e5c54 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,77 @@ +{ + "C_Cpp.intelliSenseEngine": "default", + "clangd.arguments": [ + "--include-directory=build/unit-tests" + ], + "C_Cpp.default.intelliSenseMode": "linux-gcc-x64", + "files.associations": { + "*.h": "cpp", + "array": "cpp", + "atomic": "cpp", + "bit": "cpp", + "*.tcc": "cpp", + "bitset": "cpp", + "cctype": "cpp", + "chrono": "cpp", + "clocale": "cpp", + "cmath": "cpp", + "compare": "cpp", + "concepts": "cpp", + "condition_variable": "cpp", + "cstdarg": "cpp", + "cstddef": "cpp", + "cstdint": "cpp", + "cstdio": "cpp", + "cstdlib": "cpp", + "cstring": "cpp", + "ctime": "cpp", + "cwchar": "cpp", + "cwctype": "cpp", + "deque": "cpp", + "list": "cpp", + "map": "cpp", + "set": "cpp", + "string": "cpp", + "unordered_map": "cpp", + "vector": "cpp", + "exception": "cpp", + "algorithm": "cpp", + "functional": "cpp", + "iterator": "cpp", + "memory": "cpp", + "memory_resource": "cpp", + "numeric": "cpp", + "optional": "cpp", + "random": "cpp", + "ratio": "cpp", + "regex": "cpp", + "string_view": "cpp", + "system_error": "cpp", + "tuple": "cpp", + "type_traits": "cpp", + "utility": "cpp", + "fstream": "cpp", + "future": "cpp", + "initializer_list": "cpp", + "iomanip": "cpp", + "iosfwd": "cpp", + "iostream": "cpp", + "istream": "cpp", + "limits": "cpp", + "mutex": "cpp", + "new": "cpp", + "numbers": "cpp", + "ostream": "cpp", + "semaphore": "cpp", + "sstream": "cpp", + "stdexcept": "cpp", + "stop_token": "cpp", + "streambuf": "cpp", + "thread": "cpp", + "typeinfo": "cpp", + "variant": "cpp", + "shared_mutex": "cpp" + }, + "clangd.enable": false, + "C_Cpp.dimInactiveRegions": false +} \ No newline at end of file diff --git a/.vscode/tasks.json b/.vscode/tasks.json new file mode 100644 index 0000000..45f53aa --- /dev/null +++ b/.vscode/tasks.json @@ -0,0 +1,16 @@ +{ + "version": "2.0.0", + "tasks": [ + { + "label": "build_tests", + "type": "shell", + "command": "cd build && ninja matrix-tests", + "group": { + "kind": "build", + "isDefault": true + }, + "problemMatcher": ["$gcc"], + "detail": "Generated task to build unit test executable" + } + ] +} \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..2fb853f --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,40 @@ +cmake_minimum_required(VERSION 3.6) +project(Vector3D) + +add_subdirectory(unit-tests) + +set(CMAKE_CXX_STANDARD 11) + +add_compile_options(-fdiagnostics-color=always) + +# Vector3d +add_library(Vector3D + STATIC + Vector3D.hpp +) + +set_target_properties(Vector3D + PROPERTIES + LINKER_LANGUAGE CXX +) + +target_include_directories(Vector3D PUBLIC + include +) + +# Matrix +add_library(Matrix + STATIC + Matrix.hpp + Matrix.cpp +) + +set_target_properties(Matrix + PROPERTIES + LINKER_LANGUAGE CXX +) + +target_include_directories(Matrix + PUBLIC + . +) \ No newline at end of file diff --git a/Matrix.cpp b/Matrix.cpp new file mode 100644 index 0000000..11954b7 --- /dev/null +++ b/Matrix.cpp @@ -0,0 +1,443 @@ +#ifdef MATRIX_H_ // since the .cpp file has to be included by the .hpp file this + // will evaluate to true +#include "Matrix.hpp" + +#include +#include +#include +#include + +template +Matrix::Matrix(float value) { + this->Fill(value); +} + +template +Matrix::Matrix(const std::array &array) { + this->setMatrixToArray(array); +} + +template +template +Matrix::Matrix(Args... args) { + constexpr uint16_t arraySize{static_cast(rows) * + static_cast(columns)}; + + std::initializer_list initList{static_cast(args)...}; + // choose whichever buffer size is smaller for the copy length + uint32_t minSize = + std::min(arraySize, static_cast(initList.size())); + memcpy(this->matrix.begin(), initList.begin(), minSize * sizeof(float)); +} + +template +Matrix::Matrix(const Matrix &other) { + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + this->matrix[row_idx * columns + column_idx] = + other.Get(row_idx, column_idx); + } + } +} + +template +void Matrix::setMatrixToArray( + const std::array &array) { + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + uint16_t array_idx = + static_cast(row_idx) * static_cast(columns) + + static_cast(column_idx); + if (array_idx < array.size()) { + this->matrix[row_idx * columns + column_idx] = array[array_idx]; + } else { + this->matrix[row_idx * columns + column_idx] = 0; + } + } + } +} + +template +Matrix & +Matrix::Add(const Matrix &other, + Matrix &result) const { + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + result[row_idx][column_idx] = + this->Get(row_idx, column_idx) + other.Get(row_idx, column_idx); + } + } + return result; +} + +template +Matrix & +Matrix::Sub(const Matrix &other, + Matrix &result) const { + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + result[row_idx][column_idx] = + this->Get(row_idx, column_idx) - other.Get(row_idx, column_idx); + } + } + + return result; +} + +template +template +Matrix & +Matrix::Mult(const Matrix &other, + Matrix &result) const { + // allocate some buffers for all of our dot products + Matrix<1, columns> this_row; + Matrix other_column; + Matrix<1, rows> other_column_t; + + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + // get our row + this->GetRow(row_idx, this_row); + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + // get the other matrix'ss column + other.GetColumn(column_idx, other_column); + // transpose the other matrix's column + other_column.Transpose(other_column_t); + + // the result's index is equal to the dot product of these two vectors + result[row_idx][column_idx] = + Matrix::dotProduct(this_row, other_column_t); + } + } + + return result; +} + +template +Matrix & +Matrix::Mult(float scalar, Matrix &result) const { + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + result[row_idx][column_idx] = this->Get(row_idx, column_idx) * scalar; + } + } + + return result; +} + +template +Matrix & +Matrix::Invert(Matrix &result) const { + // since all matrix sizes have to be statically specified at compile time we + // can do this + static_assert(rows == columns, + "Your matrix isn't square and can't be inverted"); + + // unfortunately we can't calculate this at compile time so we'll just reurn + // zeros + float determinant{this->Det()}; + if (determinant == 0) { + // you can't invert a matrix with a negative determinant + result.Fill(0); + return result; + } + + // TODO: This algorithm is really inneficient because of the matrix of minors. + // We should make a different algorithm how to calculate the inverse: + // https://www.mathsisfun.com/algebra/matrix-inverse-minors-cofactors-adjugate.html + + // calculate the matrix of minors + Matrix minors{}; + this->MatrixOfMinors(minors); + + // now adjugate the matrix and save it in our output + minors.adjugate(result); + + // scale the result by 1/determinant and we have our answer + result = result * (1 / determinant); + // result.Mult(1 / determinant, result); + + return result; +} + +template +Matrix & +Matrix::Transpose(Matrix &result) const { + for (uint8_t column_idx{0}; column_idx < rows; column_idx++) { + for (uint8_t row_idx{0}; row_idx < columns; row_idx++) { + result[row_idx][column_idx] = this->Get(column_idx, row_idx); + } + } + + return result; +} + +// explicitly define the determinant for a 2x2 matrix because it is definitely +// the fastest way to calculate a 2x2 matrix determinant +template <> float Matrix<0, 0>::Det() const { return 1e+6; } +template <> float Matrix<1, 1>::Det() const { return this->matrix[0]; } +template <> float Matrix<2, 2>::Det() const { + return this->matrix[0] * this->matrix[3] - this->matrix[1] * this->matrix[2]; +} + +template +float Matrix::Det() const { + static_assert(rows == columns, + "You can't take the determinant of a non-square matrix."); + + Matrix MinorMatrix{}; + float determinant{0}; + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + // for odd indices the sign is negative + float sign = (column_idx % 2 == 0) ? 1 : -1; + determinant += sign * this->matrix[column_idx] * + this->MinorMatrix(MinorMatrix, 0, column_idx).Det(); + } + + return determinant; +} + +template +Matrix & +Matrix::ElementMultiply(const Matrix &other, + Matrix &result) const { + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + result[row_idx][column_idx] = + this->Get(row_idx, column_idx) * other.Get(row_idx, column_idx); + } + } + + return result; +} + +template +Matrix & +Matrix::ElementDivide(const Matrix &other, + Matrix &result) const { + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + result[row_idx][column_idx] = + this->Get(row_idx, column_idx) / other.Get(row_idx, column_idx); + } + } + + return result; +} + +template +float Matrix::Get(uint8_t row_index, + uint8_t column_index) const { + if (row_index > rows - 1 || column_index > columns - 1) { + return 1e+10; // TODO: We should throw something here instead of failing + // quietly + } + return this->matrix[row_index * columns + column_index]; +} + +template +Matrix<1, columns> & +Matrix::GetRow(uint8_t row_index, + Matrix<1, columns> &row) const { + memcpy(&(row[0]), this->matrix.begin() + row_index * columns, + columns * sizeof(float)); + + return row; +} + +template +Matrix & +Matrix::GetColumn(uint8_t column_index, + Matrix &column) const { + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + column[row_idx][0] = this->Get(row_idx, column_index); + } + + return column; +} + +template +void Matrix::ToString(std::string &stringBuffer) const { + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + stringBuffer += "|"; + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + stringBuffer += + std::to_string(this->matrix[row_idx * columns + column_idx]); + if (column_idx != columns - 1) { + stringBuffer += "\t"; + } + } + stringBuffer += "|\n"; + } +} + +template +std::array &Matrix:: +operator[](uint8_t row_index) { + if (row_index > rows - 1) { + // TODO: We should throw something here instead of failing quietly. + row_index = 0; + } + // cursed reinterpret_cast that will help us fake having a nested array when + // we really don't + return *reinterpret_cast *>( + &(this->matrix[row_index * columns])); +} + +template +Matrix &Matrix:: +operator=(const Matrix &other) { + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + this->matrix[row_idx * columns + column_idx] = + other.Get(row_idx, column_idx); + } + } + // return a reference to ourselves so you can chain together these functions + return *this; +} + +template +Matrix Matrix:: +operator+(const Matrix &other) const { + Matrix buffer{}; + this->Add(other, buffer); + return buffer; +} + +template +Matrix Matrix:: +operator-(const Matrix &other) const { + Matrix buffer{}; + this->Sub(other, buffer); + return buffer; +} + +template +Matrix Matrix:: +operator*(const Matrix &other) const { + Matrix buffer{}; + this->Mult(other, buffer); + return buffer; +} + +template +Matrix Matrix::operator*(float scalar) const { + Matrix buffer{}; + this->Mult(scalar, buffer); + return buffer; +} + +template +template +float Matrix::dotProduct(const Matrix<1, vector_size> &vec1, + const Matrix<1, vector_size> &vec2) { + float sum{0}; + for (uint8_t i{0}; i < vector_size; i++) { + sum += vec1.Get(0, i) * vec2.Get(0, i); + } + + return sum; +} + +template +template +float Matrix::dotProduct(const Matrix &vec1, + const Matrix &vec2) { + float sum{0}; + for (uint8_t i{0}; i < vector_size; i++) { + sum += vec1.Get(i, 0) * vec2.Get(i, 0); + } + + return sum; +} + +template +void Matrix::Fill(float value) { + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + this->matrix[row_idx * columns + column_idx] = value; + } + } +} + +template +Matrix & +Matrix::MatrixOfMinors(Matrix &result) const { + Matrix MinorMatrix{}; + + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + this->MinorMatrix(MinorMatrix, row_idx, column_idx); + result[row_idx][column_idx] = MinorMatrix.Det(); + } + } + + return result; +} + +template +Matrix & +Matrix::MinorMatrix(Matrix &result, + uint8_t row_idx, uint8_t column_idx) const { + std::array subArray{}; + uint16_t array_idx{0}; + for (uint8_t row_iter{0}; row_iter < rows; row_iter++) { + if (row_iter == row_idx) { + continue; + } + for (uint8_t column_iter{0}; column_iter < columns; column_iter++) { + if (column_iter == column_idx) { + continue; + } + subArray[array_idx] = this->Get(row_iter, column_iter); + array_idx++; + } + } + + result = Matrix{subArray}; + return result; +} + +template +Matrix & +Matrix::adjugate(Matrix &result) const { + for (uint8_t row_iter{0}; row_iter < rows; row_iter++) { + for (uint8_t column_iter{0}; column_iter < columns; column_iter++) { + float sign = ((row_iter + 1) % 2) == 0 ? -1 : 1; + sign *= ((column_iter + 1) % 2) == 0 ? -1 : 1; + result[column_iter][row_iter] = this->Get(row_iter, column_iter) * sign; + } + } + + return result; +} + +template +Matrix & +Matrix::Normalize(Matrix &result) const { + float sum{0}; + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + float val{this->Get(row_idx, column_idx)}; + sum += val * val; + } + } + + if (sum == 0) { + // this wouldn't do anything anyways + result.Fill(1e+6); + return result; + } + + sum = sqrt(sum); + + for (uint8_t row_idx{0}; row_idx < rows; row_idx++) { + for (uint8_t column_idx{0}; column_idx < columns; column_idx++) { + result[row_idx][column_idx] = this->Get(row_idx, column_idx) / sum; + } + } + + return result; +} + +#endif // MATRIX_H_ \ No newline at end of file diff --git a/Matrix.hpp b/Matrix.hpp new file mode 100644 index 0000000..4b5479b --- /dev/null +++ b/Matrix.hpp @@ -0,0 +1,206 @@ +#ifndef MATRIX_H_ +#define MATRIX_H_ + +#include +#include + +// TODO: Add a function to calculate eigenvalues/vectors +// TODO: Add a function to compute RREF +// TODO: Add a function for SVD decomposition +// TODO: Add a function for LQ decomposition + +template class Matrix { +public: + /** + * @brief create a matrix but leave all of its values unitialized + */ + Matrix() = default; + + /** + * @brief Create a matrix but fill all of its entries with one value + */ + Matrix(float value); + + /** + * @brief Initialize a matrix with an array + */ + Matrix(const std::array &array); + + /** + * @brief Initialize a matrix as a copy of another matrix + */ + Matrix(const Matrix &other); + + /** + * @brief Initialize a matrix directly with any number of arguments + */ + template Matrix(Args... args); + /** + * @brief Set all elements in this to value + */ + void Fill(float value); + + /** + * @brief Element-wise matrix addition + * @param other the other matrix to add to this one + * @param result A buffer to store the result into + * @note there is no problem if result == this + */ + Matrix &Add(const Matrix &other, + Matrix &result) const; + + /** + * @brief Element-wise subtract matrix + * @param other the other matrix to subtract from this one + * @param result A buffer to store the result into + * @note there is no problem if result == this + */ + Matrix &Sub(const Matrix &other, + Matrix &result) const; + + /** + * @brief Matrix multiply the two matrices + * @param other the other matrix to multiply into this one + * @param result A buffer to store the result into + */ + template + Matrix &Mult(const Matrix &other, + Matrix &result) const; + + /** + * @brief Multiply the matrix by a scalar + * @param scalar the the scalar to multiply by + * @param result A buffer to store the result into + * @note there is no problem if result == this + */ + Matrix &Mult(float scalar, + Matrix &result) const; + + /** + * @brief Element-wise multiply the two matrices + * @param other the other matrix to multiply into this one + * @param result A buffer to store the result into + * @note there is no problem if result == this + */ + Matrix &ElementMultiply(const Matrix &other, + Matrix &result) const; + + /** + * @brief Element-wise divide the two matrices + * @param other the other matrix to multiply into this one + * @param result A buffer to store the result into + * @note there is no problem if result == this + */ + Matrix &ElementDivide(const Matrix &other, + Matrix &result) const; + + Matrix & + MinorMatrix(Matrix &result, uint8_t row_idx, + uint8_t column_idx) const; + + /** + * @return Get the determinant of the matrix + * @note for right now only 2x2 and 3x3 matrices are supported + */ + float Det() const; + + Matrix &MatrixOfMinors(Matrix &result) const; + + /** + * @brief Invert this matrix + * @param result A buffer to store the result into + * @warning this is super slow! Only call it if you absolutely have to!!! + */ + Matrix &Invert(Matrix &result) const; + + /** + * @brief Transpose this matrix + * @param result A buffer to store the result into + */ + Matrix &Transpose(Matrix &result) const; + + /** + * @brief reduce the matrix so the sum of its elements equal 1 + * @param result a buffer to store the result into + */ + Matrix &Normalize(Matrix &result) const; + + /** + * @brief Get a row from the matrix + * @param row_index the row index to get + * @param row a buffer to write the row into + */ + Matrix<1, columns> &GetRow(uint8_t row_index, Matrix<1, columns> &row) const; + + /** + * @brief Get a row from the matrix + * @param column_index the row index to get + * @param column a buffer to write the row into + */ + Matrix &GetColumn(uint8_t column_index, + Matrix &column) const; + + /** + * @brief Get the number of rows in this matrix + */ + constexpr uint8_t GetRowSize() { return rows; } + + /** + * @brief Get the number of columns in this matrix + */ + constexpr uint8_t GetColumnSize() { return columns; } + + void ToString(std::string &stringBuffer) const; + + /** + * @brief Get an element from the matrix + * @param row the row index of the element + * @param column the column index of the element + * @return The value of the element you want to get + */ + float Get(uint8_t row_index, uint8_t column_index) const; + + /** + * @brief get the specified row of the matrix returned as a reference to the + * internal array + */ + std::array &operator[](uint8_t row_index); + + /** + * @brief Copy the contents of other into this matrix + */ + Matrix &operator=(const Matrix &other); + + /** + * @brief Return a new matrix that is the sum of this matrix and other matrix + */ + Matrix operator+(const Matrix &other) const; + + Matrix operator-(const Matrix &other) const; + + Matrix operator*(const Matrix &other) const; + + Matrix operator*(float scalar) const; + +private: + /** + * @brief take the dot product of the two vectors + */ + template + static float dotProduct(const Matrix<1, vector_size> &vec1, + const Matrix<1, vector_size> &vec2); + + template + static float dotProduct(const Matrix &vec1, + const Matrix &vec2); + + Matrix &adjugate(Matrix &result) const; + + void setMatrixToArray(const std::array &array); + + std::array matrix; +}; + +#include "Matrix.cpp" + +#endif // MATRIX_H_ \ No newline at end of file diff --git a/Vector3D.h b/Vector3D.hpp similarity index 98% rename from Vector3D.h rename to Vector3D.hpp index cb60300..f625f16 100644 --- a/Vector3D.h +++ b/Vector3D.hpp @@ -2,6 +2,7 @@ #include #include +#include template class V3D{ diff --git a/unit-tests/.gitignore b/unit-tests/.gitignore new file mode 100644 index 0000000..746a662 --- /dev/null +++ b/unit-tests/.gitignore @@ -0,0 +1 @@ +matrix-test-timings-temp.txt \ No newline at end of file diff --git a/unit-tests/CMakeLists.txt b/unit-tests/CMakeLists.txt new file mode 100644 index 0000000..83f3807 --- /dev/null +++ b/unit-tests/CMakeLists.txt @@ -0,0 +1,21 @@ +cmake_minimum_required (VERSION 3.11) + +project ("test_driver") + +include(FetchContent) + +FetchContent_Declare( + Catch2 + GIT_REPOSITORY https://github.com/catchorg/Catch2.git + GIT_TAG v3.0.1 # or a later release +) + +FetchContent_MakeAvailable(Catch2) + +add_executable(matrix-tests matrix-tests.cpp) + +target_link_libraries(matrix-tests + PRIVATE + Matrix + Catch2::Catch2WithMain +) \ No newline at end of file diff --git a/unit-tests/matrix-test-timings.txt b/unit-tests/matrix-test-timings.txt new file mode 100644 index 0000000..1e2dff2 --- /dev/null +++ b/unit-tests/matrix-test-timings.txt @@ -0,0 +1,14 @@ +Addition: 0.419 s +Subtraction: 0.421 s +Multiplication: 3.297 s +Scalar Multiplication: 0.329 s +Element Multiply: 0.306 s +Element Divide: 0.302 s +Minor Matrix: 0.331 s +Determinant: 0.177 s +Matrix of Minors: 0.766 s +Invert: 0.183 s +Transpose: 0.215 s +Normalize: 0.315 s +GET ROW: 0.008 s +GET COLUMN: 0.43 s diff --git a/unit-tests/matrix-tests.cpp b/unit-tests/matrix-tests.cpp new file mode 100644 index 0000000..09a3da9 --- /dev/null +++ b/unit-tests/matrix-tests.cpp @@ -0,0 +1,405 @@ +// include the unit test framework first +#include +#include + +// include the module you're going to test next +#include "Matrix.hpp" + +// any other libraries +#include +#include +#include + +TEST_CASE("Elementary Matrix Operations", "Matrix") { + std::array arr2{5, 6, 7, 8}; + Matrix<2, 2> mat1{1, 2, 3, 4}; + Matrix<2, 2> mat2{arr2}; + Matrix<2, 2> mat3{}; + + SECTION("Initialization") { + // array initialization + REQUIRE(mat1.Get(0, 0) == 1); + REQUIRE(mat1.Get(0, 1) == 2); + REQUIRE(mat1.Get(1, 0) == 3); + REQUIRE(mat1.Get(1, 1) == 4); + + // empty initialization + REQUIRE(mat3.Get(0, 0) == 0); + REQUIRE(mat3.Get(0, 1) == 0); + REQUIRE(mat3.Get(1, 0) == 0); + REQUIRE(mat3.Get(1, 1) == 0); + + // template pack initialization + REQUIRE(mat2.Get(0, 0) == 5); + REQUIRE(mat2.Get(0, 1) == 6); + REQUIRE(mat2.Get(1, 0) == 7); + REQUIRE(mat2.Get(1, 1) == 8); + + // large matrix + Matrix<255, 255> mat6{}; + mat6.Fill(4); + for (uint8_t row{0}; row < 255; row++) { + for (uint8_t column{0}; column < 255; column++) { + REQUIRE(mat6.Get(row, column) == 4); + } + } + } + + SECTION("Fill") { + mat1.Fill(0); + REQUIRE(mat1.Get(0, 0) == 0); + REQUIRE(mat1.Get(0, 1) == 0); + REQUIRE(mat1.Get(1, 0) == 0); + REQUIRE(mat1.Get(1, 1) == 0); + + mat2.Fill(100000); + REQUIRE(mat2.Get(0, 0) == 100000); + REQUIRE(mat2.Get(0, 1) == 100000); + REQUIRE(mat2.Get(1, 0) == 100000); + REQUIRE(mat2.Get(1, 1) == 100000); + + mat3.Fill(-20); + REQUIRE(mat3.Get(0, 0) == -20); + REQUIRE(mat3.Get(0, 1) == -20); + REQUIRE(mat3.Get(1, 0) == -20); + REQUIRE(mat3.Get(1, 1) == -20); + } + + SECTION("Addition") { + std::string strBuf1 = ""; + mat1.ToString(strBuf1); + std::cout << "Matrix 1:\n" << strBuf1 << std::endl; + + mat1.Add(mat2, mat3); + + REQUIRE(mat3.Get(0, 0) == 6); + REQUIRE(mat3.Get(0, 1) == 8); + REQUIRE(mat3.Get(1, 0) == 10); + REQUIRE(mat3.Get(1, 1) == 12); + + // try out addition with overloaded operators + mat3.Fill(0); + mat3 = mat1 + mat2; + REQUIRE(mat3.Get(0, 0) == 6); + REQUIRE(mat3.Get(0, 1) == 8); + REQUIRE(mat3.Get(1, 0) == 10); + REQUIRE(mat3.Get(1, 1) == 12); + } + + SECTION("Subtraction") { + mat1.Sub(mat2, mat3); + + REQUIRE(mat3.Get(0, 0) == -4); + REQUIRE(mat3.Get(0, 1) == -4); + REQUIRE(mat3.Get(1, 0) == -4); + REQUIRE(mat3.Get(1, 1) == -4); + + // try out subtraction with operators + mat3.Fill(0); + mat3 = mat1 - mat2; + REQUIRE(mat3.Get(0, 0) == -4); + REQUIRE(mat3.Get(0, 1) == -4); + REQUIRE(mat3.Get(1, 0) == -4); + REQUIRE(mat3.Get(1, 1) == -4); + } + + SECTION("Multiplication") { + mat1.Mult(mat2, mat3); + + REQUIRE(mat3.Get(0, 0) == 19); + REQUIRE(mat3.Get(0, 1) == 22); + REQUIRE(mat3.Get(1, 0) == 43); + REQUIRE(mat3.Get(1, 1) == 50); + + // try out multiplication with operators + mat3.Fill(0); + mat3 = mat1 * mat2; + REQUIRE(mat3.Get(0, 0) == 19); + REQUIRE(mat3.Get(0, 1) == 22); + REQUIRE(mat3.Get(1, 0) == 43); + REQUIRE(mat3.Get(1, 1) == 50); + } + + SECTION("Scalar Multiplication") { + mat1.Mult(2, mat3); + + REQUIRE(mat3.Get(0, 0) == 2); + REQUIRE(mat3.Get(0, 1) == 4); + REQUIRE(mat3.Get(1, 0) == 6); + REQUIRE(mat3.Get(1, 1) == 8); + } + + SECTION("Element Multiply") { + mat1.ElementMultiply(mat2, mat3); + + REQUIRE(mat3.Get(0, 0) == 5); + REQUIRE(mat3.Get(0, 1) == 12); + REQUIRE(mat3.Get(1, 0) == 21); + REQUIRE(mat3.Get(1, 1) == 32); + } + + SECTION("Element Divide") { + mat1.ElementDivide(mat2, mat3); + + REQUIRE_THAT(mat3.Get(0, 0), Catch::Matchers::WithinRel(0.2f, 1e-6f)); + REQUIRE_THAT(mat3.Get(0, 1), Catch::Matchers::WithinRel(0.3333333f, 1e-6f)); + REQUIRE_THAT(mat3.Get(1, 0), Catch::Matchers::WithinRel(0.4285714f, 1e-6f)); + REQUIRE_THAT(mat3.Get(1, 1), Catch::Matchers::WithinRel(0.5f, 1e-6f)); + } + + SECTION("Minor Matrix") { + // what about matrices of 0,0 or 1,1? + // minor matrix for 2x2 matrix + Matrix<1, 1> minorMat1{}; + mat1.MinorMatrix(minorMat1, 0, 0); + REQUIRE(minorMat1.Get(0, 0) == 4); + mat1.MinorMatrix(minorMat1, 0, 1); + REQUIRE(minorMat1.Get(0, 0) == 3); + mat1.MinorMatrix(minorMat1, 1, 0); + REQUIRE(minorMat1.Get(0, 0) == 2); + mat1.MinorMatrix(minorMat1, 1, 1); + REQUIRE(minorMat1.Get(0, 0) == 1); + + // minor matrix for 3x3 matrix + Matrix<3, 3> mat4{1, 2, 3, 4, 5, 6, 7, 8, 9}; + Matrix<2, 2> minorMat4{}; + + mat4.MinorMatrix(minorMat4, 0, 0); + REQUIRE(minorMat4.Get(0, 0) == 5); + REQUIRE(minorMat4.Get(0, 1) == 6); + REQUIRE(minorMat4.Get(1, 0) == 8); + REQUIRE(minorMat4.Get(1, 1) == 9); + + mat4.MinorMatrix(minorMat4, 1, 1); + REQUIRE(minorMat4.Get(0, 0) == 1); + REQUIRE(minorMat4.Get(0, 1) == 3); + REQUIRE(minorMat4.Get(1, 0) == 7); + REQUIRE(minorMat4.Get(1, 1) == 9); + + mat4.MinorMatrix(minorMat4, 2, 2); + REQUIRE(minorMat4.Get(0, 0) == 1); + REQUIRE(minorMat4.Get(0, 1) == 2); + REQUIRE(minorMat4.Get(1, 0) == 4); + REQUIRE(minorMat4.Get(1, 1) == 5); + } + + SECTION("Determinant") { + float det1 = mat1.Det(); + + REQUIRE_THAT(det1, Catch::Matchers::WithinRel(-2.0F, 1e-6f)); + + Matrix<3, 3> mat4{1, 2, 3, 4, 5, 6, 7, 8, 9}; + + float det4 = mat4.Det(); + + REQUIRE_THAT(det4, Catch::Matchers::WithinRel(0.0F, 1e-6f)); + + Matrix<3, 3> mat5{1, 0, 0, 0, 2, 0, 0, 0, 3}; + + float det5 = mat5.Det(); + REQUIRE_THAT(det5, Catch::Matchers::WithinRel(6.0F, 1e-6f)); + } + + SECTION("Matrix of Minors") { + mat1.MatrixOfMinors(mat3); + REQUIRE_THAT(mat3.Get(0, 0), Catch::Matchers::WithinRel(4.0F, 1e-6f)); + REQUIRE_THAT(mat3.Get(0, 1), Catch::Matchers::WithinRel(3.0F, 1e-6f)); + REQUIRE_THAT(mat3.Get(1, 0), Catch::Matchers::WithinRel(2.0F, 1e-6f)); + REQUIRE_THAT(mat3.Get(1, 1), Catch::Matchers::WithinRel(1.0F, 1e-6f)); + + Matrix<3, 3> mat4{1, 2, 3, 4, 5, 6, 7, 8, 9}; + Matrix<3, 3> mat5{0}; + mat4.MatrixOfMinors(mat5); + REQUIRE_THAT(mat5.Get(0, 0), Catch::Matchers::WithinRel(-3.0F, 1e-6f)); + REQUIRE_THAT(mat5.Get(0, 1), Catch::Matchers::WithinRel(-6.0F, 1e-6f)); + REQUIRE_THAT(mat5.Get(0, 2), Catch::Matchers::WithinRel(-3.0F, 1e-6f)); + REQUIRE_THAT(mat5.Get(1, 0), Catch::Matchers::WithinRel(-6.0F, 1e-6f)); + REQUIRE_THAT(mat5.Get(1, 1), Catch::Matchers::WithinRel(-12.0F, 1e-6f)); + REQUIRE_THAT(mat5.Get(1, 2), Catch::Matchers::WithinRel(-6.0F, 1e-6f)); + REQUIRE_THAT(mat5.Get(2, 0), Catch::Matchers::WithinRel(-3.0F, 1e-6f)); + REQUIRE_THAT(mat5.Get(2, 1), Catch::Matchers::WithinRel(-6.0F, 1e-6f)); + REQUIRE_THAT(mat5.Get(2, 2), Catch::Matchers::WithinRel(-3.0F, 1e-6f)); + } + + SECTION("Invert") { + mat1.Invert(mat3); + REQUIRE_THAT(mat3.Get(0, 0), Catch::Matchers::WithinRel(-2.0F, 1e-6f)); + REQUIRE_THAT(mat3.Get(0, 1), Catch::Matchers::WithinRel(1.0F, 1e-6f)); + REQUIRE_THAT(mat3.Get(1, 0), Catch::Matchers::WithinRel(1.5F, 1e-6f)); + REQUIRE_THAT(mat3.Get(1, 1), Catch::Matchers::WithinRel(-0.5F, 1e-6f)); + }; + + SECTION("Transpose") { + // transpose a square matrix + mat1.Transpose(mat3); + + REQUIRE(mat3.Get(0, 0) == 1); + REQUIRE(mat3.Get(0, 1) == 3); + REQUIRE(mat3.Get(1, 0) == 2); + REQUIRE(mat3.Get(1, 1) == 4); + + // transpose a non-square matrix + Matrix<2, 3> mat4{1, 2, 3, 4, 5, 6}; + Matrix<3, 2> mat5{}; + + mat4.Transpose(mat5); + + REQUIRE(mat5.Get(0, 0) == 1); + REQUIRE(mat5.Get(0, 1) == 4); + REQUIRE(mat5.Get(1, 0) == 2); + REQUIRE(mat5.Get(1, 1) == 5); + REQUIRE(mat5.Get(2, 0) == 3); + REQUIRE(mat5.Get(2, 1) == 6); + } + + SECTION("Normalize") { + mat1.Normalize(mat3); + + float sqrt_30{sqrt(30)}; + + REQUIRE(mat3.Get(0, 0) == 1 / sqrt_30); + REQUIRE(mat3.Get(0, 1) == 2 / sqrt_30); + REQUIRE(mat3.Get(1, 0) == 3 / sqrt_30); + REQUIRE(mat3.Get(1, 1) == 4 / sqrt_30); + + Matrix<2, 1> mat4{-0.878877044, 2.92092276}; + Matrix<2, 1> mat5{}; + mat4.Normalize(mat5); + + REQUIRE_THAT(mat5.Get(0, 0), + Catch::Matchers::WithinRel(-0.288129855179f, 1e-6f)); + REQUIRE_THAT(mat5.Get(1, 0), + Catch::Matchers::WithinRel(0.957591346325f, 1e-6f)); + } + + SECTION("GET ROW") { + Matrix<1, 2> mat1Rows{}; + mat1.GetRow(0, mat1Rows); + REQUIRE(mat1Rows.Get(0, 0) == 1); + REQUIRE(mat1Rows.Get(0, 1) == 2); + + mat1.GetRow(1, mat1Rows); + REQUIRE(mat1Rows.Get(0, 0) == 3); + REQUIRE(mat1Rows.Get(0, 1) == 4); + } + + SECTION("GET COLUMN") { + Matrix<2, 1> mat1Columns{}; + mat1.GetColumn(0, mat1Columns); + REQUIRE(mat1Columns.Get(0, 0) == 1); + REQUIRE(mat1Columns.Get(1, 0) == 3); + + mat1.GetColumn(1, mat1Columns); + REQUIRE(mat1Columns.Get(0, 0) == 2); + REQUIRE(mat1Columns.Get(1, 0) == 4); + } +} + +// basically re-run all of the previous tests with huge matrices and time the +// results. +TEST_CASE("Timing Tests", "Matrix") { + std::array arr1{}; + for (uint16_t i{0}; i < 50 * 50; i++) { + arr1[i] = i; + } + std::array arr2{5, 6, 7, 8}; + for (uint16_t i{50 * 50}; i < 2 * 50 * 50; i++) { + arr2[i] = i; + } + Matrix<50, 50> mat1{arr1}; + Matrix<50, 50> mat2{arr2}; + Matrix<50, 50> mat3{}; + + // A smaller matrix to use for really badly optimized operations + Matrix<4, 4> mat4{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}; + Matrix<4, 4> mat5{}; + + SECTION("Addition") { + for (uint32_t i{0}; i < 10000; i++) { + mat3 = mat1 + mat2; + } + } + + SECTION("Subtraction") { + for (uint32_t i{0}; i < 10000; i++) { + mat3 = mat1 - mat2; + } + } + + SECTION("Multiplication") { + for (uint32_t i{0}; i < 1000; i++) { + mat3 = mat1 * mat2; + } + } + + SECTION("Scalar Multiplication") { + for (uint32_t i{0}; i < 10000; i++) { + mat3 = mat1 * 3; + } + } + + SECTION("Element Multiply") { + for (uint32_t i{0}; i < 10000; i++) { + mat1.ElementMultiply(mat2, mat3); + } + } + + SECTION("Element Divide") { + for (uint32_t i{0}; i < 10000; i++) { + mat1.ElementDivide(mat2, mat3); + } + } + + SECTION("Minor Matrix") { + // what about matrices of 0,0 or 1,1? + // minor matrix for 2x2 matrix + Matrix<49, 49> minorMat1{}; + for (uint32_t i{0}; i < 10000; i++) { + mat1.MinorMatrix(minorMat1, 0, 0); + } + } + + SECTION("Determinant") { + for (uint32_t i{0}; i < 100000; i++) { + float det1 = mat4.Det(); + } + } + + SECTION("Matrix of Minors") { + for (uint32_t i{0}; i < 100000; i++) { + mat4.MatrixOfMinors(mat5); + } + } + + SECTION("Invert") { + for (uint32_t i{0}; i < 100000; i++) { + mat4.Invert(mat5); + } + }; + + SECTION("Transpose") { + for (uint32_t i{0}; i < 10000; i++) { + mat1.Transpose(mat3); + } + } + + SECTION("Normalize") { + for (uint32_t i{0}; i < 10000; i++) { + mat1.Normalize(mat3); + } + } + + SECTION("GET ROW") { + Matrix<1, 50> mat1Rows{}; + for (uint32_t i{0}; i < 1000000; i++) { + mat1.GetRow(0, mat1Rows); + } + } + + SECTION("GET COLUMN") { + Matrix<50, 1> mat1Columns{}; + for (uint32_t i{0}; i < 1000000; i++) { + mat1.GetColumn(0, mat1Columns); + } + } +} \ No newline at end of file diff --git a/unit-tests/run-matrix-timing-test.sh b/unit-tests/run-matrix-timing-test.sh new file mode 100755 index 0000000..f6d6638 --- /dev/null +++ b/unit-tests/run-matrix-timing-test.sh @@ -0,0 +1,7 @@ +# be in the root folder of this project when you run this +cd build/ +ninja matrix-tests +echo "Running tests. This will take a while." +./unit-tests/matrix-tests -n "Timing Tests" -d yes > ../unit-tests/matrix-test-timings-temp.txt +cd ../unit-tests/ +python3 test-timing-post-process.py \ No newline at end of file diff --git a/unit-tests/test-timing-post-process.py b/unit-tests/test-timing-post-process.py new file mode 100644 index 0000000..3b56360 --- /dev/null +++ b/unit-tests/test-timing-post-process.py @@ -0,0 +1,108 @@ +class Timing: + def __init__(self, time: float, test_name: str): + self.time = time + self.test_name: str = test_name + self.difference = -1 + + def __eq__(self, other: "Timing"): + return self.test_name.lower() == other.test_name.lower() + + def __sub__(self, other: "Timing"): + return self.time - other.time + + def to_string(self): + return f"{self.test_name}: {self.time} s" + + def to_string_w_diff(self): + diff = self.difference + if diff == -1: + diff = 0 + return f"{self.test_name}: {round(self.time,3)} s, Difference: {round(diff,3)}" + +def create_timing_from_test_line(line: str) -> Timing: + time_end_idx: int = line.find(" ") + if time_end_idx == -1: + return None + try: + time: float = float(line[0:time_end_idx]) + except: + print("Couldn't convert: " + line[0:time_end_idx] + " to a float") + return None + + test_name = line[time_end_idx+4:-1] + return Timing(time, test_name) + +def parse_test_file(file_path: str) -> list[Timing]: + timings: list[Timing] = [] + with open(file_path, 'r') as file: + previous_line = "" + for line in file: + if line.find("Timing Tests") != -1: + timing = create_timing_from_test_line(previous_line) + if timing is not None: + timings.append(timing) + previous_line = line[:] # deep copy line + return timings + +def parse_timing_file(file_path: str) -> list[Timing]: + timings: list[Timing] = [] + with open(file_path, 'r') as file: + for line in file: + seperator_idx = line.find(":") + if seperator_idx == -1: + continue + test_name = line[:seperator_idx] + try: + time = float(line[seperator_idx+2:-2]) + print(time) + except: + print("Couldn't convert: " + line[seperator_idx:-2] + " to a float") + continue + timings.append(Timing(time, test_name)) + + + return timings + +def save_timings(timings: list[Timing], file_path: str): + with open(file_path, 'w') as file: + for timing in timings: + file.write(f"{timing.to_string()}\n") + +parse_file_path = "matrix-test-timings-temp.txt" +save_file_path = "matrix-test-timings.txt" + +# get the new timings +new_timings = parse_test_file(parse_file_path) + +# get the old timings +old_timings = parse_timing_file(save_file_path) + +difference_increased = "" +# calculate the timing difference +for new_timing in new_timings: + for old_timing in old_timings: + if new_timing == old_timing: + new_timing.difference = new_timing - old_timing + if abs(new_timing.difference) >= 0.03: + difference_increased += f"{new_timing.test_name}, " + +def save_option(): + # save the new timings + while True: + option = input("Save Results? (y/n)") + if option[0].lower() == 'y': + save_timings(new_timings, save_file_path) + print("Saved.") + break + elif option[0].lower() == 'n': + break + +# print the new timing results along with the difference +for timing in new_timings: + print(timing.to_string_w_diff()) + +if len(difference_increased) > 0: + print("You've made major timing changes for:" + difference_increased) + save_option() +else: + print("No times have changed outside the margin of error.")