UpperTriangularMatrix.h

Go to the documentation of this file.

#ifndef INCLUDE_MOLASSEMBLER_TEMPLE_CONSTEXPR_TRIANGULAR_MATRIX_H
#define INCLUDE_MOLASSEMBLER_TEMPLE_CONSTEXPR_TRIANGULAR_MATRIX_H

#include "Molassembler/Temple/constexpr/FloatingPointComparison.h"

namespace Scine {
namespace Molassembler {
namespace Temple {
namespace UpperTriangularMatrixImpl {

template<typename T, size_t size>
using ArrayType = std::array<T, size>;

namespace index_conversion {
  template<size_t N, typename UnsignedType>
  PURITY_STRONG constexpr unsigned toSingleIndex(const UnsignedType i, const UnsignedType j) {
    return (
      N * (N - 1) / 2
      - (N - i) * ((N - i) - 1) / 2
      + j - i - 1
    );
  }

  template<size_t N, typename UnsignedType>
  PURITY_STRONG constexpr std::pair<UnsignedType, UnsignedType> toDoubleIndex(const UnsignedType k) {
    UnsignedType i {
      N - 2
      - static_cast<unsigned>(
        Math::floor(
          Math::sqrt(
            0.0
            + 4 * N * (N - 1)
            - 8 * k
            - 7
          ) / 2.0
          - 0.5
        )
      )
    };

    return {
      i,
      (
        k + i + 1
        - N * (N - 1) / 2
        + (N - i) * ((N - i) - 1) / 2
      )
    };
  }
} // namespace index_conversion

PURITY_STRONG constexpr bool isValidMatrixSize(const size_t dataSize) {
  auto root = Math::sqrt(1.0 + 8.0 * dataSize);

  return(
    Floating::isCloseRelative(
      root,
      static_cast<double>(Math::floor(root)),
      1e-6
    )
  );
}

PURITY_STRONG constexpr size_t getMatrixSize(const size_t dataSize) {
  return (
    1 + Math::sqrt(1.0 + 8.0 * dataSize)
  ) / 2;
}

} // namespace UpperTriangularMatrixImpl

template<typename ValueType, size_t dataSize>
class UpperTriangularMatrix {
public:
  using DataType = UpperTriangularMatrixImpl::ArrayType<ValueType, dataSize>;

  static constexpr unsigned N = UpperTriangularMatrixImpl::getMatrixSize(dataSize);

  constexpr UpperTriangularMatrix() : data_ {} {}

  template<
    template<typename, size_t> class ArrayType
  > constexpr explicit UpperTriangularMatrix(
    const ArrayType<ValueType, dataSize>& data
  ) : data_(data)
  {
    static_assert(
      UpperTriangularMatrixImpl::isValidMatrixSize(dataSize) && dataSize > 0,
      "Passed data size is not consistent with an upper triangular matrix!"
    );
  }


  constexpr ValueType& at(const unsigned i, const unsigned j) {
    if(i >= j || i >= N || j >= N) {
      throw "Index out of bounds!";
    }

    return data_.at(
      UpperTriangularMatrixImpl::index_conversion::toSingleIndex<N>(i, j)
    );
  }

  PURITY_WEAK constexpr const ValueType& at(
    const unsigned i,
    const unsigned j
  ) const {
    if(i >= j || i >= N || j >= N) {
      throw "Index out of bounds!";
    }

    return data_.at(
      UpperTriangularMatrixImpl::index_conversion::toSingleIndex<N>(i, j)
    );
  }

  constexpr const DataType& getData() const {
    return data_;
  }

private:
/* State */
  DataType data_;
};

template<
  template<typename, size_t> class ArrayType,
  typename ValueType,
  size_t size
> constexpr UpperTriangularMatrix<ValueType, size> makeUpperTriangularMatrix(
  const ArrayType<ValueType, size>& data
) {
  return UpperTriangularMatrix<ValueType, size>(data);
}

} // namespace Temple
} // namespace Molassembler
} // namespace Scine

#endif