SparseMatrix.h

Go to the documentation of this file.

#ifndef SPARSEMATRIX_H
#define SPARSEMATRIX_H
 
#include <iostream>
#include <map>
#include <string>
#include <vector>
#include <stdexcept>
 
template<typename T>
using SparseMatrix = std::map<std::pair<T, T>, double>;
 
template<typename T>
std::vector<T> getDiagonalElements(const SparseMatrix<T>& matrix) {
    std::vector<T> diagonalElements;
    for (const auto& element : matrix) {
        if (element.first.first == element.first.second) {
            diagonalElements.push_back(element.first.first);
        }
    }
    return diagonalElements;
}
 
template<typename T>
SparseMatrix<T> createIdentityMatrix(const std::vector<T>& indices) {
    SparseMatrix<T> identity;
    for (const T& index : indices) {
        identity[{index, index}] = 1.0;
    }
    return identity;
}
 
template<typename T>
double getElement(const SparseMatrix<T>& matrix, const T& row, const T& col) {
    auto it = matrix.find({row, col});
    if (it != matrix.end()) {
        return it->second;
    }
    return 0.0;
}
 
template<typename T>
void setElement(SparseMatrix<T>& matrix, const T& row, const T& col, double value) {
    if (value != 0.0) {
        matrix[{row, col}] = value;
    } else {
        matrix.erase({row, col});
    }
}
 
template<typename T>
SparseMatrix<T> invertMatrix(const SparseMatrix<T>& matrix, const std::vector<T>& indices) {
    SparseMatrix<T> augmentedMatrix = matrix;
    SparseMatrix<T> inverse = createIdentityMatrix(indices);
 
    for (size_t i = 0; i < indices.size(); ++i) {
        T pivot = indices[i];
 
        double pivotValue = getElement(augmentedMatrix, pivot, pivot);
        if (pivotValue == 0.0) {
            throw std::runtime_error("Matrix is singular and cannot be inverted.");
        }
 
        for (size_t k = 0; k < indices.size(); ++k) {
            T col = indices[k];
            setElement(augmentedMatrix, pivot, col, getElement(augmentedMatrix, pivot, col) / pivotValue);
            setElement(inverse, pivot, col, getElement(inverse, pivot, col) / pivotValue);
        }
 
        for (size_t j = 0; j < indices.size(); ++j) {
            if (i == j) continue;
            T row = indices[j];
            double factor = getElement(augmentedMatrix, row, pivot);
            for (size_t k = 0; k < indices.size(); ++k) {
                T col = indices[k];
                setElement(augmentedMatrix, row, col, getElement(augmentedMatrix, row, col) - factor * getElement(augmentedMatrix, pivot, col));
                setElement(inverse, row, col, getElement(inverse, row, col) - factor * getElement(inverse, pivot, col));
            }
        }
    }
 
    return inverse;
}
 
template<typename T>
std::pair<SparseMatrix<T>, std::vector<T>> removeEmptyRowsAndCols(const SparseMatrix<T>& matrix, const std::vector<T>& indices) {
    std::map<T, int> rowCount;
    std::map<T, int> colCount;
 
    // Count non-zero elements per row and column
    for (const auto& [coord, value] : matrix) {
        if (value != 0.0) { // <-- Ignore zero entries explicitly present
            T row = coord.first;
            T col = coord.second;
            rowCount[row]++;
            colCount[col]++;
        }
    }
 
    // Keep only indices where row and column are not empty
    std::vector<T> cleanedIndices;
    for (const T& idx : indices) {
        if (rowCount[idx] > 0 || colCount[idx] > 0) {
            cleanedIndices.push_back(idx);
        }
    }
 
    // Filter matrix to keep only relevant rows/cols
    SparseMatrix<T> cleanedMatrix;
    for (const auto& [coord, value] : matrix) {
        if (value != 0.0) {
            T row = coord.first;
            T col = coord.second;
            if (std::find(cleanedIndices.begin(), cleanedIndices.end(), row) != cleanedIndices.end() &&
                std::find(cleanedIndices.begin(), cleanedIndices.end(), col) != cleanedIndices.end()) {
                cleanedMatrix[coord] = value;
            }
        }
    }
 
    return {cleanedMatrix, cleanedIndices};
}
 
 
template<typename T>
void printMatrix(const SparseMatrix<T>& matrix, const std::vector<T>& indices) {
    for (const auto& i : indices) {
        for (const auto& j : indices) {
            std::cout << getElement(matrix, i, j) << " ";
        }
        std::cout << std::endl;
    }
}
 
template<typename T>
void customPrintMatrix(const SparseMatrix<T>& matrix, const std::vector<T>& indices) {
    // Determine maximum width for any label or value
    size_t maxLabelWidth = 0;
    size_t maxValueWidth = 0;
 
    std::ostringstream oss;
    for (const T& index : indices) {
        oss.str("");
        oss.clear();
        oss << index;
        maxLabelWidth = std::max(maxLabelWidth, oss.str().length());
    }
 
    for (const auto& [coord, value] : matrix) {
        oss.str("");
        oss.clear();
        oss << value;
        maxValueWidth = std::max(maxValueWidth, oss.str().length());
    }
 
    size_t cellWidth = std::max(maxLabelWidth, maxValueWidth) + 2;
 
    // Print top-left empty cell
    std::cout << std::setw(cellWidth) << " ";
 
    // Print column headers
    for (const auto& col : indices) {
        std::cout << std::setw(cellWidth) << col;
    }
    std::cout << "\n";
 
    // Print separator line
    std::cout << std::setw(cellWidth) << " ";
    for (size_t i = 0; i < indices.size(); ++i) {
        std::cout << std::setw(cellWidth) << std::string(cellWidth - 1, '-');
    }
    std::cout << "\n";
 
    // Print rows
    for (const auto& row : indices) {
        std::cout << std::setw(cellWidth) << row;
        for (const auto& col : indices) {
            double val = getElement(matrix, row, col);
            std::cout << std::setw(cellWidth) << val;
        }
        std::cout << "\n";
    }
}
 
template<typename T>
class SparseMatrixWrapper {
public:
    SparseMatrix<T> matrix;
 
    SparseMatrixWrapper() = default;
 
    double getElement(T row, T col) const {
        return ::getElement(matrix, row, col);
    }
 
    void setElement(T row, T col, double value) {
        ::setElement(matrix, row, col, value);
    }
 
    std::vector<T> getDiagonalElements() const {
        return ::getDiagonalElements(matrix);
    }
 
    void print(const std::vector<T>& indices) const {
        ::printMatrix(matrix, indices);
    }
 
    SparseMatrixWrapper<T> invert(const std::vector<T>& indices) const {
        SparseMatrixWrapper<T> result;
        result.matrix = ::invertMatrix(matrix, indices);
        return result;
    }
};
 
 
#endif // __SPARSEMATRIX_H__