Matrix.h

Go to the documentation of this file.

#ifndef INTARNA_MATRIX_H_
#define INTARNA_MATRIX_H_
 
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <limits>
#include <span>
#include <stdexcept>
#include <utility>
#include <vector>
 
#include "IntaRNA/intarna_config.h"
#if INTARNA_USE_STD_MDSPAN
#include <mdspan>
#else
#include "mdspan/mdspan.hpp"
#endif
 
namespace IntaRNA {
 
namespace matrix_detail {
#if INTARNA_USE_STD_MDSPAN
namespace md = std;
#else
namespace md = MDSPAN_IMPL_STANDARD_NAMESPACE;
#endif
 
inline std::size_t product(std::size_t rows, std::size_t columns) {
    if (columns != 0 && rows > std::numeric_limits<std::size_t>::max() / columns)
        throw std::length_error("matrix dimensions overflow");
    return rows * columns;
}
inline std::size_t triangle(std::size_t n) {
    if (n == std::numeric_limits<std::size_t>::max())
        throw std::length_error("matrix dimensions overflow");
    return n % 2 == 0 ? product(n / 2, n + 1) : product(n, (n + 1) / 2);
}
}
 
template<class T>
class Matrix {
    std::vector<T> values;
    std::size_t rows = 0, columns = 0;
    using Extents = matrix_detail::md::dextents<std::size_t, 2>;
public:
    using value_type = T;
    Matrix() = default;
    Matrix(std::size_t rows, std::size_t columns, const T &value = T{});
    Matrix(const Matrix &) = default;
    Matrix &operator=(const Matrix &) = default;
    Matrix(Matrix &&other) noexcept;
    Matrix &operator=(Matrix &&other) noexcept;
    std::size_t size1() const noexcept;
    std::size_t size2() const noexcept;
    std::size_t storageSize() const noexcept;
    T &operator()(std::size_t i, std::size_t j);
    const T &operator()(std::size_t i, std::size_t j) const;
    void clear();
    void swap(Matrix &other) noexcept;
    void resize(std::size_t newRows, std::size_t newColumns, bool preserve = true);
};
 
template<class T>
inline
Matrix<T>::Matrix(std::size_t rows, std::size_t columns, const T &value)
    : values(matrix_detail::product(rows, columns), value), rows(rows), columns(columns)
{}
 
template<class T>
inline
Matrix<T>::Matrix(Matrix &&other) noexcept : values(std::move(other.values)),
    rows(std::exchange(other.rows, 0)), columns(std::exchange(other.columns, 0))
{}
 
template<class T>
inline
Matrix<T> &Matrix<T>::operator=(Matrix &&other) noexcept
{
    if (this != &other) {
        Matrix moved(std::move(other));
        swap(moved);
    }
    return *this;
}
 
template<class T>
inline
std::size_t Matrix<T>::size1() const noexcept
{ return rows; }
 
template<class T>
inline
std::size_t Matrix<T>::size2() const noexcept
{ return columns; }
 
template<class T>
inline
std::size_t Matrix<T>::storageSize() const noexcept
{ return values.size(); }
 
template<class T>
inline
T &Matrix<T>::operator()(std::size_t i, std::size_t j)
{
    assert(i < rows && j < columns);
    return matrix_detail::md::mdspan<T, Extents>(values.data(), rows, columns)[i, j];
}
 
template<class T>
inline
const T &Matrix<T>::operator()(std::size_t i, std::size_t j) const
{
    assert(i < rows && j < columns);
    return matrix_detail::md::mdspan<const T, Extents>(values.data(), rows, columns)[i, j];
}
 
template<class T>
inline
void Matrix<T>::clear()
{ std::fill(values.begin(), values.end(), T{}); }
 
template<class T>
inline
void Matrix<T>::swap(Matrix &other) noexcept
{
    values.swap(other.values);
    std::swap(rows, other.rows);
    std::swap(columns, other.columns);
}
 
template<class T>
inline
void Matrix<T>::resize(std::size_t newRows, std::size_t newColumns, bool preserve)
{
    if (rows == newRows && columns == newColumns) return;
    if (!preserve) {
        values.resize(matrix_detail::product(newRows, newColumns));
        rows = newRows;
        columns = newColumns;
        return;
    }
    Matrix next(newRows, newColumns);
    for (std::size_t i = 0; i < std::min(rows, newRows); ++i)
        for (std::size_t j = 0; j < std::min(columns, newColumns); ++j)
            next(i, j) = (*this)(i, j);
    swap(next);
}
 
 
template<class T>
class UpperTriangularMatrix {
    std::vector<T> values;
    std::size_t n = 0;
    using Extents = matrix_detail::md::dextents<std::size_t, 1>;
    std::size_t offset(std::size_t i, std::size_t j) const noexcept;
    static std::size_t count(std::size_t rows, std::size_t columns);
public:
    using value_type = T;
    UpperTriangularMatrix() = default;
    UpperTriangularMatrix(std::size_t rows, std::size_t columns);
    UpperTriangularMatrix(const UpperTriangularMatrix &) = default;
    UpperTriangularMatrix &operator=(const UpperTriangularMatrix &) = default;
    UpperTriangularMatrix(UpperTriangularMatrix &&other) noexcept;
    UpperTriangularMatrix &operator=(UpperTriangularMatrix &&other) noexcept;
    std::size_t size1() const noexcept;
    std::size_t size2() const noexcept;
    std::size_t storageSize() const noexcept;
    T &operator()(std::size_t i, std::size_t j);
    const T &operator()(std::size_t i, std::size_t j) const;
    void clear();
    void swap(UpperTriangularMatrix &other) noexcept;
    void resize(std::size_t rows, std::size_t columns, bool preserve = true);
};
 
template<class T>
inline
std::size_t UpperTriangularMatrix<T>::offset(std::size_t i, std::size_t j) const noexcept
{
    const auto remaining = n - i;
    // Constructor has already checked that every triangular size fits.
    const auto tail = remaining % 2 == 0
        ? (remaining / 2) * (remaining + 1) : remaining * ((remaining + 1) / 2);
    return values.size() - tail + j - i;
}
 
template<class T>
inline
std::size_t UpperTriangularMatrix<T>::count(std::size_t rows, std::size_t columns)
{
    if (rows != columns) throw std::invalid_argument("upper-triangular matrix must be square");
    return matrix_detail::triangle(rows);
}
 
template<class T>
inline
UpperTriangularMatrix<T>::UpperTriangularMatrix(std::size_t rows, std::size_t columns)
    : values(count(rows, columns)), n(rows)
{}
 
template<class T>
inline
UpperTriangularMatrix<T>::UpperTriangularMatrix(UpperTriangularMatrix &&other) noexcept
    : values(std::move(other.values)), n(std::exchange(other.n, 0))
{}
 
template<class T>
inline
UpperTriangularMatrix<T> &UpperTriangularMatrix<T>::operator=(UpperTriangularMatrix &&other) noexcept
{
    if (this != &other) {
        UpperTriangularMatrix moved(std::move(other));
        swap(moved);
    }
    return *this;
}
 
template<class T>
inline
std::size_t UpperTriangularMatrix<T>::size1() const noexcept
{ return n; }
 
template<class T>
inline
std::size_t UpperTriangularMatrix<T>::size2() const noexcept
{ return n; }
 
template<class T>
inline
std::size_t UpperTriangularMatrix<T>::storageSize() const noexcept
{ return values.size(); }
 
template<class T>
inline
T &UpperTriangularMatrix<T>::operator()(std::size_t i, std::size_t j)
{
    assert(i <= j && j < n);
    return matrix_detail::md::mdspan<T, Extents>(values.data(), values.size())[offset(i, j)];
}
 
template<class T>
inline
const T &UpperTriangularMatrix<T>::operator()(std::size_t i, std::size_t j) const
{
    assert(i < n && j < n);
    static const T zero{};
    return i > j ? zero : matrix_detail::md::mdspan<const T, Extents>(values.data(), values.size())[offset(i, j)];
}
 
template<class T>
inline
void UpperTriangularMatrix<T>::clear()
{ std::fill(values.begin(), values.end(), T{}); }
 
template<class T>
inline
void UpperTriangularMatrix<T>::swap(UpperTriangularMatrix &other) noexcept
{
    values.swap(other.values);
    std::swap(n, other.n);
}
 
template<class T>
inline
void UpperTriangularMatrix<T>::resize(std::size_t rows, std::size_t columns, bool preserve)
{
    const auto cells = count(rows, columns);
    if (rows == n) return;
    if (!preserve) {
        values.resize(cells);
        n = rows;
        return;
    }
    UpperTriangularMatrix next(rows, columns);
    for (std::size_t i = 0; i < std::min(n, rows); ++i)
        for (std::size_t j = i; j < std::min(n, rows); ++j)
            next(i, j) = (*this)(i, j);
    swap(next);
}
 
 
template<class T>
class UpperBandedMatrix {
    Matrix<T> band;
    std::size_t columns = 0;
    static std::size_t width(std::size_t columns, std::size_t lower, std::size_t upper);
public:
    using value_type = T;
    UpperBandedMatrix() = default;
    UpperBandedMatrix(std::size_t rows, std::size_t columns, std::size_t lower, std::size_t upper);
    UpperBandedMatrix(const UpperBandedMatrix &) = default;
    UpperBandedMatrix &operator=(const UpperBandedMatrix &) = default;
    UpperBandedMatrix(UpperBandedMatrix &&other) noexcept;
    UpperBandedMatrix &operator=(UpperBandedMatrix &&other) noexcept;
    std::size_t size1() const noexcept;
    std::size_t size2() const noexcept;
    std::size_t storageSize() const noexcept;
    T &operator()(std::size_t i, std::size_t j);
    const T &operator()(std::size_t i, std::size_t j) const;
    std::span<T> row(std::size_t i);
    std::span<const T> row(std::size_t i) const;
    void clear();
    void swap(UpperBandedMatrix &other) noexcept;
    void resize(std::size_t rows, std::size_t newColumns, std::size_t lower,
            std::size_t upper, bool preserve = true);
};
 
template<class T>
inline
std::size_t UpperBandedMatrix<T>::width(std::size_t columns, std::size_t lower, std::size_t upper)
{
    if (lower != 0) throw std::invalid_argument("upper-banded matrix requires lower=0");
    return upper >= columns ? columns : upper + 1;
}
 
template<class T>
inline
UpperBandedMatrix<T>::UpperBandedMatrix(std::size_t rows, std::size_t columns, std::size_t lower, std::size_t upper)
    : band(rows, width(columns, lower, upper)), columns(columns)
{}
 
template<class T>
inline
UpperBandedMatrix<T>::UpperBandedMatrix(UpperBandedMatrix &&other) noexcept
    : band(std::move(other.band)), columns(std::exchange(other.columns, 0))
{}
 
template<class T>
inline
UpperBandedMatrix<T> &UpperBandedMatrix<T>::operator=(UpperBandedMatrix &&other) noexcept
{
    if (this != &other) {
        UpperBandedMatrix moved(std::move(other));
        swap(moved);
    }
    return *this;
}
 
template<class T>
inline
std::size_t UpperBandedMatrix<T>::size1() const noexcept
{ return band.size1(); }
 
template<class T>
inline
std::size_t UpperBandedMatrix<T>::size2() const noexcept
{ return columns; }
 
template<class T>
inline
std::size_t UpperBandedMatrix<T>::storageSize() const noexcept
{ return band.storageSize(); }
 
template<class T>
inline
T &UpperBandedMatrix<T>::operator()(std::size_t i, std::size_t j)
{
    assert(i < size1() && j < columns && i <= j && j - i < band.size2());
    return band(i, j - i);
}
 
template<class T>
inline
const T &UpperBandedMatrix<T>::operator()(std::size_t i, std::size_t j) const
{
    assert(i < size1() && j < columns);
    static const T zero{};
    return i > j || j - i >= band.size2() ? zero : band(i, j - i);
}
 
template<class T>
inline
std::span<T> UpperBandedMatrix<T>::row(std::size_t i)
{
    assert(i < size1());
    const auto count = i < columns ? std::min(band.size2(), columns-i) : 0;
    return count == 0 ? std::span<T>{} : std::span<T>{&band(i, 0), count};
}
 
template<class T>
inline
std::span<const T> UpperBandedMatrix<T>::row(std::size_t i) const
{
    assert(i < size1());
    const auto count = i < columns ? std::min(band.size2(), columns-i) : 0;
    return count == 0 ? std::span<const T>{} : std::span<const T>{&band(i, 0), count};
}
 
template<class T>
inline
void UpperBandedMatrix<T>::clear()
{ band.clear(); }
 
template<class T>
inline
void UpperBandedMatrix<T>::swap(UpperBandedMatrix &other) noexcept
{
    band.swap(other.band);
    std::swap(columns, other.columns);
}
 
template<class T>
inline
void UpperBandedMatrix<T>::resize(std::size_t rows, std::size_t newColumns, std::size_t lower,
        std::size_t upper, bool preserve)
{
    UpperBandedMatrix next(rows, newColumns, lower, upper);
    if (preserve) {
        for (std::size_t i = 0; i < std::min(size1(), rows); ++i)
            for (std::size_t d = 0; d < std::min(band.size2(), next.band.size2())
                    && i < std::min(columns, newColumns)
                    && d < std::min(columns, newColumns) - i; ++d)
                next.band(i, d) = band(i, d);
    }
    swap(next);
}
 
 
} // namespace IntaRNA
#endif