rinmatrix/src/matrix.cpp

/*
    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 <fstream>
#include <iostream>
#include <stdexcept>
#include <vector>

class Matrix {
    std::vector<std::vector<double>> m_data;

   public:
    Matrix() {}

    // Copy constructor combined with transpose
    Matrix(const Matrix& m, bool transpose = 0) {
        const std::vector<std::vector<double>>& 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<std::vector<double>>& GetData() { return m_data; }

    // Getter of m_data
    const std::vector<std::vector<double>>& 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<std::vector<double>>&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<std::vector<double>>&data1 = m1.GetData(),
          &data2 = m2.GetData();
    std::vector<std::vector<double>>& 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<std::vector<double>>&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';
}