Matrix.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1

#ifndef BPP_NUMERIC_MATRIX_MATRIX_H
#define BPP_NUMERIC_MATRIX_MATRIX_H

#include <iostream>
#include <vector>

#include "../../Clonable.h"
#include "../NumConstants.h"
#include "../NumTools.h"
#include "../VectorExceptions.h"

namespace bpp
{
template<class Scalar>
class Matrix :
  public Clonable
{
public:
  Matrix() {}
  virtual ~Matrix() {}

public:
  virtual const Scalar& operator()(size_t i, size_t j) const = 0;

  virtual Scalar& operator()(size_t i, size_t j) = 0;

  virtual bool equals(const Matrix& m, double threshold = NumConstants::TINY())
  {
    if (m.getNumberOfRows() != getNumberOfRows() || m.getNumberOfColumns() != getNumberOfColumns())
      return false;
    for (size_t i = 0; i < getNumberOfRows(); i++)
    {
      for (size_t j = 0; j < getNumberOfColumns(); j++)
      {
        if (NumTools::abs<double>(static_cast<double>(operator()(i, j)) - static_cast<double>(m(i, j))) > threshold) return false;
      }
    }
    return true;
  }
  virtual size_t getNumberOfRows() const = 0;
  virtual size_t getNumberOfColumns() const = 0;
  virtual std::vector<Scalar> row(size_t i) const = 0;
  virtual std::vector<Scalar> col(size_t j) const = 0;
  virtual void resize(size_t nRows, size_t nCols) = 0;
};


template<class Scalar>
class RowMatrix :
  public Matrix<Scalar>
{
protected:
  std::vector< std::vector<Scalar>> m_;

public:
  RowMatrix() : m_() {}

  RowMatrix(size_t nRow, size_t nCol) : m_(nRow)
  {
    for (size_t i = 0; i < nRow; i++)
    {
      m_[i].resize(nCol);
    }
  }

  RowMatrix(const Matrix<Scalar>& m) : m_(m.getNumberOfRows())
  {
    size_t nr = m.getNumberOfRows();
    size_t nc = m.getNumberOfColumns();
    for (size_t i = 0; i < nr; i++)
    {
      m_[i].resize(nc);
      for (size_t j = 0; j < nc; j++)
      {
        m_[i][j] = m(i, j);
      }
    }
  }

  RowMatrix& operator=(const Matrix<Scalar>& m)
  {
    size_t nr = m.getNumberOfRows();
    m_.resize(nr);
    size_t nc = m.getNumberOfColumns();
    for (size_t i = 0; i < nr; i++)
    {
      m_[i].resize(nc);
      for (size_t j = 0; j < nc; j++)
      {
        m_[i][j] = m(i, j);
      }
    }
    return *this;
  }

//  virtual ~RowMatrix() {}

public:
  RowMatrix* clone() const { return new RowMatrix(*this); }

  const Scalar& operator()(size_t i, size_t j) const { return m_[i][j]; }

  Scalar& operator()(size_t i, size_t j) { return m_[i][j]; }

  size_t getNumberOfRows() const { return m_.size(); }

  size_t getNumberOfColumns() const { return m_.size() == 0 ? 0 : m_[0].size(); }

  std::vector<Scalar> row(size_t i) const
  {
    std::vector<Scalar> r(getNumberOfColumns());
    for (size_t j = 0; j < getNumberOfColumns(); j++) { r[j] = operator()(i, j); }
    return r;
  }

  const std::vector<Scalar>& getRow(size_t i) const
  {
    return m_[i];
  }

  std::vector<Scalar>& getRow(size_t i)
  {
    return m_[i];
  }

  std::vector<Scalar> col(size_t j) const
  {
    std::vector<Scalar> c(getNumberOfRows());
    for (size_t i = 0; i < getNumberOfRows(); i++) { c[i] = operator()(i, j); }
    return c;
  }

  void resize(size_t nRows, size_t nCols)
  {
    m_.resize(nRows);
    for (size_t i = 0; i < nRows; i++)
    {
      m_[i].resize(nCols);
    }
  }

  void addRow(const std::vector<Scalar>& newRow)
  {
    if (getNumberOfColumns() != 0 && newRow.size() != getNumberOfColumns())
      throw DimensionException("RowMatrix::addRow: invalid row dimension", newRow.size(), getNumberOfColumns());
    m_.push_back(newRow);
  }
};

template<class Scalar>
class ColMatrix :
  public Matrix<Scalar>
{
private:
  std::vector< std::vector<Scalar>> m_;

public:
  ColMatrix() : m_() {}

  ColMatrix(size_t nRow, size_t nCol) : m_(nCol)
  {
    for (size_t i = 0; i < nCol; i++)
    {
      m_[i].resize(nRow);
    }
  }

  ColMatrix(const Matrix<Scalar>& m) : m_(m.getNumberOfColumns())
  {
    size_t nr = m.getNumberOfRows();
    size_t nc = m.getNumberOfColumns();
    for (size_t i = 0; i < nc; i++)
    {
      m_[i].resize(nr);
      for (size_t j = 0; j < nr; j++)
      {
        m_[i][j] = m(j, i);
      }
    }
  }

  ColMatrix& operator=(const Matrix<Scalar>& m)
  {
    size_t nc = m.getNumberOfColumns();
    m_.resize(nc);
    size_t nr = m.getNumberOfRows();
    for (size_t i = 0; i < nc; i++)
    {
      m_[i].resize(nr);
      for (size_t j = 0; j < nr; j++)
      {
        m_[i][j] = m(j, i);
      }
    }
    return *this;
  }

  virtual ~ColMatrix() {}

public:
  ColMatrix* clone() const { return new ColMatrix(*this); }

  const Scalar& operator()(size_t i, size_t j) const { return m_[j][i]; }

  Scalar& operator()(size_t i, size_t j) { return m_[j][i]; }

  size_t getNumberOfColumns() const { return m_.size(); }

  size_t getNumberOfRows() const { return m_.size() == 0 ? 0 : m_[0].size(); }

  std::vector<Scalar> row(size_t i) const
  {
    std::vector<Scalar> r(getNumberOfColumns());
    for (size_t j = 0; j < getNumberOfColumns(); j++) { r[j] = operator()(i, j); }
    return r;
  }

  const std::vector<Scalar>& getCol(size_t i) const
  {
    return m_[i];
  }

  std::vector<Scalar>& getCol(size_t i)
  {
    return m_[i];
  }

  std::vector<Scalar> col(size_t j) const
  {
    std::vector<Scalar> c(getNumberOfRows());
    for (size_t i = 0; i < getNumberOfRows(); i++) { c[i] = operator()(i, j); }
    return c;
  }

  void resize(size_t nRows, size_t nCols)
  {
    m_.resize(nCols);
    for (size_t i = 0; i < nCols; i++)
    {
      m_[i].resize(nRows);
    }
  }

  void addCol(const std::vector<Scalar>& newCol)
  {
    if (getNumberOfRows() != 0 && newCol.size() != getNumberOfRows())
      throw DimensionException("ColMatrix::addCol: invalid column dimension", newCol.size(), getNumberOfRows());
    m_.push_back(newCol);
  }
};

template<class Scalar>
class LinearMatrix :
  public Matrix<Scalar>
{
private:
  std::vector<Scalar> m_;
  size_t rows_;
  size_t cols_;

public:
  LinearMatrix() : m_(),
    rows_(0),
    cols_(0) { resize_(0, 0); }

  LinearMatrix(size_t nRow, size_t nCol) : m_(),
    rows_(nRow),
    cols_(nCol) { resize_(nRow, nCol); }

  LinearMatrix(const Matrix<Scalar>& m) : m_(m.getNumberOfRows() * m.getNumberOfColumns()), rows_(m.getNumberOfRows()), cols_(m.getNumberOfColumns())
  {
    for (size_t i = 0; i < rows_; i++)
    {
      for (size_t j = 0; j < cols_; j++)
      {
        m_[i * cols_ + j] = m(i, j);
      }
    }
  }

  LinearMatrix& operator=(const Matrix<Scalar>& m)
  {
    rows_ = m.getNumberOfRows();
    cols_ = m.getNumberOfColumns();
    m_.resize(rows_ * cols_);
    for (size_t i = 0; i < rows_; i++)
    {
      for (size_t j = 0; j < cols_; j++)
      {
        m_[i * cols_ + j] = m(i, j);
      }
    }

    return *this;
  }

  virtual ~LinearMatrix() {}

public:
  LinearMatrix* clone() const { return new LinearMatrix(*this); }

  const Scalar& operator()(size_t i, size_t j) const { return m_[i * cols_ + j]; }

  Scalar& operator()(size_t i, size_t j) { return m_[i * cols_ + j]; }

  size_t getNumberOfRows() const { return rows_; }

  size_t getNumberOfColumns() const { return cols_; }

  std::vector<Scalar> row(size_t i) const
  {
    std::vector<Scalar> r(getNumberOfColumns());
    for (size_t j = 0; j < getNumberOfColumns(); j++)
    {
      r[j] = operator()(i, j);
    }
    return r;
  }

  std::vector<Scalar> col(size_t j) const
  {
    std::vector<Scalar> c(getNumberOfRows());
    for (size_t i = 0; i < getNumberOfRows(); i++)
    {
      c[i] = operator()(i, j);
    }
    return c;
  }

  void resize(size_t nRows, size_t nCols)
  {
    resize(nRows, nCols, true);
  }

  void resize(size_t nRows, size_t nCols, bool keepValues)
  {
    LinearMatrix<Scalar> tmpM;
    if (keepValues)
      tmpM = *this;
    resize_(nRows, nCols);
    if (keepValues)
    {
      for (size_t i = 0; i < nRows; i++)
      {
        for (size_t j = 0; j < nCols; j++)
        {
          if (i < tmpM.getNumberOfRows() && j < tmpM.getNumberOfColumns())
          {
            operator()(i, j) = tmpM(i, j);
          }
          else
          {
            operator()(i, j) = 0;
          }
        }
      }
    }
  }

private:
  void resize_(size_t nRows, size_t nCols)
  {
    m_.resize(nRows * nCols);
    rows_ = nRows;
    cols_ = nCols;
  }
};

template<class Scalar>
bool operator==(const Matrix<Scalar>& m1, const Matrix<Scalar>& m2)
{
  if (m1.getNumberOfRows() != m2.getNumberOfRows() || m1.getNumberOfColumns() != m2.getNumberOfColumns())
    return false;
  for (size_t i = 0; i < m1.getNumberOfRows(); i++)
  {
    for (size_t j = 0; j < m1.getNumberOfColumns(); j++)
    {
      if (m1(i, j) != m2(i, j))
        return false;
    }
  }
  return true;
}
} // end of namespace bpp.
#endif // BPP_NUMERIC_MATRIX_MATRIX_H