Merge pull request #1 from Cynopolis/Add-Matrix

Added a matrix library and unit tests for it
This commit is contained in:
Quinn
2024-12-14 19:54:47 -05:00
committed by GitHub
16 changed files with 1383 additions and 0 deletions

Binary file not shown.

1
.gitignore vendored Normal file
View File

@@ -0,0 +1 @@
build/

43
.vscode/launch.json vendored Normal file
View File

@@ -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"
}
]
}

77
.vscode/settings.json vendored Normal file
View File

@@ -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
}

16
.vscode/tasks.json vendored Normal file
View File

@@ -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"
}
]
}

40
CMakeLists.txt Normal file
View File

@@ -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
.
)

443
Matrix.cpp Normal file
View File

@@ -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 <algorithm>
#include <cmath>
#include <cstdlib>
#include <type_traits>
template <uint8_t rows, uint8_t columns>
Matrix<rows, columns>::Matrix(float value) {
this->Fill(value);
}
template <uint8_t rows, uint8_t columns>
Matrix<rows, columns>::Matrix(const std::array<float, rows * columns> &array) {
this->setMatrixToArray(array);
}
template <uint8_t rows, uint8_t columns>
template <typename... Args>
Matrix<rows, columns>::Matrix(Args... args) {
constexpr uint16_t arraySize{static_cast<uint16_t>(rows) *
static_cast<uint16_t>(columns)};
std::initializer_list<float> initList{static_cast<float>(args)...};
// choose whichever buffer size is smaller for the copy length
uint32_t minSize =
std::min(arraySize, static_cast<uint16_t>(initList.size()));
memcpy(this->matrix.begin(), initList.begin(), minSize * sizeof(float));
}
template <uint8_t rows, uint8_t columns>
Matrix<rows, columns>::Matrix(const Matrix<rows, columns> &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 <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::setMatrixToArray(
const std::array<float, rows * columns> &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<uint16_t>(row_idx) * static_cast<uint16_t>(columns) +
static_cast<uint16_t>(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 <uint8_t rows, uint8_t columns>
Matrix<rows, columns> &
Matrix<rows, columns>::Add(const Matrix<rows, columns> &other,
Matrix<rows, columns> &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 <uint8_t rows, uint8_t columns>
Matrix<rows, columns> &
Matrix<rows, columns>::Sub(const Matrix<rows, columns> &other,
Matrix<rows, columns> &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 <uint8_t rows, uint8_t columns>
template <uint8_t other_columns>
Matrix<rows, columns> &
Matrix<rows, columns>::Mult(const Matrix<columns, other_columns> &other,
Matrix<rows, other_columns> &result) const {
// allocate some buffers for all of our dot products
Matrix<1, columns> this_row;
Matrix<rows, 1> 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<rows, columns>::dotProduct(this_row, other_column_t);
}
}
return result;
}
template <uint8_t rows, uint8_t columns>
Matrix<rows, columns> &
Matrix<rows, columns>::Mult(float scalar, Matrix<rows, columns> &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 <uint8_t rows, uint8_t columns>
Matrix<rows, columns> &
Matrix<rows, columns>::Invert(Matrix<rows, columns> &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<rows, columns> 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 <uint8_t rows, uint8_t columns>
Matrix<columns, rows> &
Matrix<rows, columns>::Transpose(Matrix<columns, rows> &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 <uint8_t rows, uint8_t columns>
float Matrix<rows, columns>::Det() const {
static_assert(rows == columns,
"You can't take the determinant of a non-square matrix.");
Matrix<rows - 1, columns - 1> 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 <uint8_t rows, uint8_t columns>
Matrix<rows, columns> &
Matrix<rows, columns>::ElementMultiply(const Matrix<rows, columns> &other,
Matrix<rows, columns> &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 <uint8_t rows, uint8_t columns>
Matrix<rows, columns> &
Matrix<rows, columns>::ElementDivide(const Matrix<rows, columns> &other,
Matrix<rows, columns> &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 <uint8_t rows, uint8_t columns>
float Matrix<rows, columns>::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 <uint8_t rows, uint8_t columns>
Matrix<1, columns> &
Matrix<rows, columns>::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 <uint8_t rows, uint8_t columns>
Matrix<rows, 1> &
Matrix<rows, columns>::GetColumn(uint8_t column_index,
Matrix<rows, 1> &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 <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::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 <uint8_t rows, uint8_t columns>
std::array<float, columns> &Matrix<rows, columns>::
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<std::array<float, columns> *>(
&(this->matrix[row_index * columns]));
}
template <uint8_t rows, uint8_t columns>
Matrix<rows, columns> &Matrix<rows, columns>::
operator=(const Matrix<rows, columns> &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 <uint8_t rows, uint8_t columns>
Matrix<rows, columns> Matrix<rows, columns>::
operator+(const Matrix<rows, columns> &other) const {
Matrix<rows, columns> buffer{};
this->Add(other, buffer);
return buffer;
}
template <uint8_t rows, uint8_t columns>
Matrix<rows, columns> Matrix<rows, columns>::
operator-(const Matrix<rows, columns> &other) const {
Matrix<rows, columns> buffer{};
this->Sub(other, buffer);
return buffer;
}
template <uint8_t rows, uint8_t columns>
Matrix<rows, columns> Matrix<rows, columns>::
operator*(const Matrix<rows, columns> &other) const {
Matrix<rows, columns> buffer{};
this->Mult(other, buffer);
return buffer;
}
template <uint8_t rows, uint8_t columns>
Matrix<rows, columns> Matrix<rows, columns>::operator*(float scalar) const {
Matrix<rows, columns> buffer{};
this->Mult(scalar, buffer);
return buffer;
}
template <uint8_t rows, uint8_t columns>
template <uint8_t vector_size>
float Matrix<rows, columns>::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 <uint8_t rows, uint8_t columns>
template <uint8_t vector_size>
float Matrix<rows, columns>::dotProduct(const Matrix<vector_size, 1> &vec1,
const Matrix<vector_size, 1> &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 <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::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 <uint8_t rows, uint8_t columns>
Matrix<rows, columns> &
Matrix<rows, columns>::MatrixOfMinors(Matrix<rows, columns> &result) const {
Matrix<rows - 1, columns - 1> 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 <uint8_t rows, uint8_t columns>
Matrix<rows - 1, columns - 1> &
Matrix<rows, columns>::MinorMatrix(Matrix<rows - 1, columns - 1> &result,
uint8_t row_idx, uint8_t column_idx) const {
std::array<float, (rows - 1) * (columns - 1)> 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<rows - 1, columns - 1>{subArray};
return result;
}
template <uint8_t rows, uint8_t columns>
Matrix<rows, columns> &
Matrix<rows, columns>::adjugate(Matrix<rows, columns> &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 <uint8_t rows, uint8_t columns>
Matrix<rows, columns> &
Matrix<rows, columns>::Normalize(Matrix<rows, columns> &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_

206
Matrix.hpp Normal file
View File

@@ -0,0 +1,206 @@
#ifndef MATRIX_H_
#define MATRIX_H_
#include <array>
#include <cstdint>
// 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 <uint8_t rows, uint8_t columns> 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<float, rows * columns> &array);
/**
* @brief Initialize a matrix as a copy of another matrix
*/
Matrix(const Matrix<rows, columns> &other);
/**
* @brief Initialize a matrix directly with any number of arguments
*/
template <typename... Args> 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<rows, columns> &Add(const Matrix<rows, columns> &other,
Matrix<rows, columns> &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<rows, columns> &Sub(const Matrix<rows, columns> &other,
Matrix<rows, columns> &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 <uint8_t other_columns>
Matrix<rows, columns> &Mult(const Matrix<columns, other_columns> &other,
Matrix<rows, other_columns> &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<rows, columns> &Mult(float scalar,
Matrix<rows, columns> &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<rows, columns> &ElementMultiply(const Matrix<rows, columns> &other,
Matrix<rows, columns> &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<rows, columns> &ElementDivide(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const;
Matrix<rows - 1, columns - 1> &
MinorMatrix(Matrix<rows - 1, columns - 1> &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<rows, columns> &MatrixOfMinors(Matrix<rows, columns> &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<rows, columns> &Invert(Matrix<rows, columns> &result) const;
/**
* @brief Transpose this matrix
* @param result A buffer to store the result into
*/
Matrix<columns, rows> &Transpose(Matrix<columns, rows> &result) const;
/**
* @brief reduce the matrix so the sum of its elements equal 1
* @param result a buffer to store the result into
*/
Matrix<rows, columns> &Normalize(Matrix<rows, columns> &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<rows, 1> &GetColumn(uint8_t column_index,
Matrix<rows, 1> &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<float, columns> &operator[](uint8_t row_index);
/**
* @brief Copy the contents of other into this matrix
*/
Matrix<rows, columns> &operator=(const Matrix<rows, columns> &other);
/**
* @brief Return a new matrix that is the sum of this matrix and other matrix
*/
Matrix<rows, columns> operator+(const Matrix<rows, columns> &other) const;
Matrix<rows, columns> operator-(const Matrix<rows, columns> &other) const;
Matrix<rows, columns> operator*(const Matrix<rows, columns> &other) const;
Matrix<rows, columns> operator*(float scalar) const;
private:
/**
* @brief take the dot product of the two vectors
*/
template <uint8_t vector_size>
static float dotProduct(const Matrix<1, vector_size> &vec1,
const Matrix<1, vector_size> &vec2);
template <uint8_t vector_size>
static float dotProduct(const Matrix<vector_size, 1> &vec1,
const Matrix<vector_size, 1> &vec2);
Matrix<rows, columns> &adjugate(Matrix<rows, columns> &result) const;
void setMatrixToArray(const std::array<float, rows * columns> &array);
std::array<float, rows * columns> matrix;
};
#include "Matrix.cpp"
#endif // MATRIX_H_

View File

@@ -2,6 +2,7 @@
#include <cstdint> #include <cstdint>
#include <cmath> #include <cmath>
#include <type_traits>
template <typename Type> template <typename Type>
class V3D{ class V3D{

1
unit-tests/.gitignore vendored Normal file
View File

@@ -0,0 +1 @@
matrix-test-timings-temp.txt

21
unit-tests/CMakeLists.txt Normal file
View File

@@ -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
)

View File

@@ -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

405
unit-tests/matrix-tests.cpp Normal file
View File

@@ -0,0 +1,405 @@
// include the unit test framework first
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
// include the module you're going to test next
#include "Matrix.hpp"
// any other libraries
#include <array>
#include <cmath>
#include <iostream>
TEST_CASE("Elementary Matrix Operations", "Matrix") {
std::array<float, 4> 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<float, 50 * 50> arr1{};
for (uint16_t i{0}; i < 50 * 50; i++) {
arr1[i] = i;
}
std::array<float, 50 * 50> 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);
}
}
}

View File

@@ -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

View File

@@ -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.")