From b07d22a8067cac1e3b800117199464f71271b1c3 Mon Sep 17 00:00:00 2001 From: RinRi Date: Tue, 29 Mar 2022 19:18:08 +0600 Subject: [PATCH] initial commit: added SSAD assignment 3 --- .gitignore | 5 ++ CMakeLists.txt | 10 +++ src/matrix.cpp | 174 +++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 189 insertions(+) create mode 100644 .gitignore create mode 100644 CMakeLists.txt create mode 100644 src/matrix.cpp diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..3822669 --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +build/ +CMakeFiles +Makefile +cmake_install.cmake +CMakeCache.txt \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..076eb2c --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,10 @@ +cmake_minimum_required(VERSION 3.22.0) + +project(rinmatrix) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall") + +set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin) + +add_executable(rinmatrix src/matrix.cpp) \ No newline at end of file diff --git a/src/matrix.cpp b/src/matrix.cpp new file mode 100644 index 0000000..15b15a5 --- /dev/null +++ b/src/matrix.cpp @@ -0,0 +1,174 @@ +/* + SSAD Assignment 3 + Implement least square approximation + (matrix multiplication and inverse as well) + + TODO: Refactor the code so that it can be reused (make modular) +*/ + +#include +#include +#include +#include + +class Matrix { + std::vector> m_data; + + public: + Matrix() {} + + // Copy constructor combined with transpose + Matrix(const Matrix& m, bool transpose = 0) { + const std::vector>& data = m.GetData(); + if (transpose) { + if (data.size()) { + int rows = data[0].size(), cols = data.size(); + m_data.resize(rows); + for (auto& i : m_data) i.resize(cols, 0); + + for (int i = 0; i < rows; ++i) + for (int j = 0; j < cols; ++j) m_data[i][j] = data[j][i]; + } + } else { + m_data = data; + } + } + + // Move constructor + Matrix(Matrix&& m) { m_data = std::move(m.GetData()); } + + // Constructor with size + Matrix(size_t rows, size_t cols) { + m_data.resize(rows); + for (auto& i : m_data) i.resize(cols, 0); + } + + // reference-Getter of m_data + std::vector>& GetData() { return m_data; } + + // Getter of m_data + const std::vector>& GetData() const { return m_data; } + + size_t GetRows() const { return m_data.size(); } + + size_t GetCols() const { + if (m_data.size() == 0) return 0; + return m_data[0].size(); + } + + // Matrix multiplication operator + friend Matrix operator*(const Matrix& m1, const Matrix& m2); + + // Inverse matrix using Gaussian elimination + static Matrix Inverse(const Matrix& m) { + int rows = m.GetRows(), cols = m.GetCols(); + if (rows != cols) + throw std::invalid_argument("Inverse is only for square matrices!"); + + Matrix from(m), res(rows, rows); + std::vector>&datares = res.GetData(), + data = from.GetData(); + + for (int i = 0; i < rows; ++i) datares[i][i] = 1; + + for (int i = 0; i < rows; ++i) { + if (data[i][i] == 0.0) + throw std::invalid_argument("No inverse for this matrix!"); + + for (int j = 0; j < rows; ++j) { + if (i != j) { + double ratio = data[j][i] / data[i][i]; + + for (int k = 0; k < rows; ++k) + data[j][k] -= ratio * data[i][k]; + + for (int k = 0; k < rows; ++k) + datares[j][k] -= ratio * datares[i][k]; + } + } + } + + for (int i = 0; i < rows; ++i) + for (int j = 0; j < rows; ++j) datares[i][j] /= data[i][i]; + + return res; + } + + void print(std::ostream& out) const { + for (const auto& i : m_data) { + for (const auto& j : i) out << j << ' '; + out << '\n'; + } + } +}; + +Matrix operator*(const Matrix& m1, const Matrix& m2) { + size_t rows = m1.GetRows(), cols = m2.GetCols(), mid = m1.GetCols(); + if (mid != m2.GetRows()) throw std::invalid_argument("Couldn't multiply!"); + + Matrix res(rows, cols); + const std::vector>&data1 = m1.GetData(), + &data2 = m2.GetData(); + std::vector>& datares = res.GetData(); + + for (size_t i = 0; i < rows; ++i) { + for (size_t j = 0; j < cols; ++j) { + datares[i][j] = 0; + for (size_t k = 0; k < mid; ++k) { + datares[i][j] += data1[i][k] * data2[k][j]; + } + } + } + + return res; +} + +inline std::ostream& operator<<(std::ostream& out, const Matrix& m) { + m.print(out); + return out; +} + +int main() { + std::ifstream inputtxt("input.txt"); + std::ofstream outputtxt("output.txt"); + +#ifdef DEBUG + std::istream& in = inputtxt; + std::ostream& out = std::cout; +#else + std::istream& in = inputtxt; + std::ostream& out = outputtxt; +#endif + out.precision(2); + + int n, m; + in >> n >> m; + Matrix A(m, n + 1), Y(m, 1); + + // input + std::vector>&a = A.GetData(), &y = Y.GetData(); + for (int i = 0; i < m; ++i) { + for (int j = 0; j < n + 1; ++j) { + if (j == 0) a[i][j] = 1; + if (j < n) + in >> a[i][j + 1]; + else + in >> y[i][0]; + } + } + + out << "A:\n" << std::fixed << A << '\n'; + out << "b:\n" << std::fixed << Y << '\n'; + + Matrix A_T(A, true), // A transpose + A_TA(std::move(A_T * A)), // A transpose multiplied by A + A_TA_inv( + std::move(Matrix::Inverse(A_TA))), // inverse of the previous one + A_TA_invA_T(std::move(A_TA_inv * A_T)), // multiplied by A transpose + ans(std::move(A_TA_invA_T * Y)); // answer x + + out << "A_T*A:\n" << std::fixed << A_TA << '\n'; + out << "(A_T*A)_-1:\n" << std::fixed << A_TA_inv << '\n'; + out << "(A_T*A)_-1*A_T:\n" << std::fixed << A_TA_invA_T << '\n'; + out << "x:\n" << std::fixed << ans << '\n'; +}