DistanceMatrix.h

Go to the documentation of this file.

// Copyright (c) 2002-present, OpenMS Inc. -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Mathias Walzer $
// $Authors: $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/CONCEPT/Colorizer.h>
#include <OpenMS/CONCEPT/Exception.h>
#include <OpenMS/CONCEPT/Macros.h>
#include <OpenMS/CONCEPT/Types.h>
#include <OpenMS/config.h>
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
 
namespace OpenMS
{
 
template<typename Value>
class DistanceMatrix
{
public:
 
  typedef Value value_type;
 
 
  typedef Size SizeType;
  typedef value_type ValueType;
 
  DistanceMatrix(): matrix_(nullptr), init_size_(0), dimensionsize_(0), min_element_(0, 0)
  {
  }
 
  DistanceMatrix(SizeType dimensionsize, Value value = Value()):
      matrix_(new ValueType*[dimensionsize]),
      init_size_(dimensionsize),
      dimensionsize_(dimensionsize),
      min_element_(0, 0)
  {
    matrix_[0] = NULL;
    SizeType i = 1;
    for (i = 1; i < dimensionsize; ++i)
    {
      matrix_[i] = new ValueType[i];
      if (matrix_[i] == NULL)
      {
        SizeType j = i;
        for (i = 1; i < j; i++)
        {
          delete[] matrix_[i];
        }
        delete[] matrix_;
        matrix_ = NULL;
        dimensionsize_ = 0;
        init_size_ = 0;
        throw Exception::OutOfMemory(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
                                     (UInt)((((dimensionsize - 2) * (dimensionsize - 1)) / 2) * sizeof(ValueType)));
      }
    }
    if (matrix_ != NULL)
    {
      for (i = 1; i < dimensionsize; ++i)
      {
        for (SizeType j = 0; j < i; ++j)
        {
          matrix_[i][j] = value;
        }
      }
      min_element_ = std::make_pair(1, 0);
    }
  }
 
  DistanceMatrix(const DistanceMatrix& source):
      matrix_(new ValueType*[source.dimensionsize_]),
      init_size_(source.dimensionsize_),
      dimensionsize_(source.dimensionsize_),
      min_element_(source.min_element_)
  {
    matrix_[0] = NULL;
    SizeType i = 1;
    for (i = 1; i < dimensionsize_; ++i)
    {
      matrix_[i] = new ValueType[i];
      if (matrix_[i] == NULL)
      {
        SizeType j = i;
        for (i = 1; i < j; i++)
        {
          delete[] matrix_[i];
        }
        delete[] matrix_;
        matrix_ = NULL;
        dimensionsize_ = 0;
        init_size_ = 0;
        min_element_ = std::make_pair(0, 0);
        throw Exception::OutOfMemory(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
                                     (UInt)((((dimensionsize_ - 2) * (dimensionsize_ - 1)) / 2) * sizeof(ValueType)));
      }
    }
    if (matrix_ != NULL)
    {
      for (i = 1; i < dimensionsize_; ++i)
      {
        std::copy(source.matrix_[i], source.matrix_[i] + i, matrix_[i]);
      }
    }
  }
 
  ~DistanceMatrix()
  {
    for (SizeType i = 1; i < init_size_; i++)
    {
      delete[] matrix_[i];
    }
    delete[] matrix_;
  }
 
  const ValueType operator()(SizeType i, SizeType j) const
  {
    return getValue(i, j);
  }
 
  ValueType operator()(SizeType i, SizeType j)
  {
    return getValue(i, j);
  }
 
  const ValueType getValue(SizeType i, SizeType j) const
  {
    if (i >= dimensionsize_ || j >= dimensionsize_) { throw Exception::OutOfRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION); }
    // elements on main diagonal are not stored and assumed to be 0
    if (i == j) { return 0; }
    if (i < j) { std::swap(i, j); }
    return (const ValueType)(matrix_[i][j]);
  }
 
  ValueType getValue(SizeType i, SizeType j)
  {
    if (i >= dimensionsize_ || j >= dimensionsize_) { throw Exception::OutOfRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION); }
    // elements on main diagonal are not stored and assumed to be 0
    if (i == j) { return 0; }
    if (i < j) { std::swap(i, j); }
    return matrix_[i][j];
  }
 
  void setValue(SizeType i, SizeType j, ValueType value)
  {
    if (i >= dimensionsize_ || j >= dimensionsize_) { throw Exception::OutOfRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION); }
    // elements on main diagonal are not stored and assumed to be 0
    if (i != j)
    {
      if (i < j) { std::swap(i, j); }
      // Fast path whenever (i,j) is NOT the current minimum cell: writing into any
      // other cell can only introduce a new (smaller) minimum, never invalidate the
      // existing one. Using '||' (i.e. "not the min cell") instead of '&&' is required:
      // with '&&' a cell that shares exactly one index with the min coordinate fell into
      // the else-branch, which never refreshes min_element_ when the new value is smaller,
      // leaving a stale cached minimum (issue #9488, DATA-24).
      if (i != min_element_.first || j != min_element_.second)
      {
        matrix_[i][j] = value;
        if (value < matrix_[min_element_.first][min_element_.second]) // keep min_element_ up-to-date
        {
          min_element_ = std::make_pair(i, j);
        }
      }
      else
      {
        if (value <= matrix_[min_element_.first][min_element_.second]) { matrix_[i][j] = value; }
        else
        {
          matrix_[i][j] = value;
          updateMinElement();
        }
      }
    }
  }
 
  void setValueQuick(SizeType i, SizeType j, ValueType value)
  {
    if (i >= dimensionsize_ || j >= dimensionsize_) { throw Exception::OutOfRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION); }
    // elements on main diagonal are not stored and assumed to be 0
    if (i != j)
    {
      if (i < j) { std::swap(i, j); }
      matrix_[i][j] = value;
    }
  }
 
  void clear()
  {
    for (SizeType i = 1; i < init_size_; i++)
    {
      delete[] matrix_[i];
    }
    delete[] matrix_;
    matrix_ = nullptr;
    min_element_ = std::make_pair(0, 0);
    dimensionsize_ = 0;
    init_size_ = 0;
  }
 
  void resize(SizeType dimensionsize, Value value = Value())
  {
    for (SizeType j = 1; j < init_size_; j++)
    {
      delete[] matrix_[j];
    }
    delete[] matrix_;
    dimensionsize_ = dimensionsize;
    init_size_ = dimensionsize;
    min_element_ = std::make_pair(0, 0);
    matrix_ = new ValueType*[dimensionsize_];
    for (SizeType j = 1; j < dimensionsize_; ++j)
    {
      matrix_[j] = new ValueType[j];
      if (matrix_[j] == nullptr)
      {
        for (SizeType k = 1; k < j; ++k)
        {
          delete[] matrix_[k];
        }
        delete[] matrix_;
        matrix_ = nullptr;
        dimensionsize_ = 0;
        init_size_ = 0;
        throw Exception::OutOfMemory(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
                                     (UInt)((((dimensionsize_ - 2) * (dimensionsize_ - 1)) / 2) * sizeof(Value)));
      }
    }
    if (matrix_ != nullptr)
    {
      for (SizeType j = 0; j < dimensionsize; ++j)
      {
        for (SizeType k = 0; k < j; ++k)
        {
          matrix_[j][k] = value;
        }
      }
      min_element_ = std::make_pair(1, 0);
    }
  }
 
  void reduce(SizeType j)
  {
    if (j >= dimensionsize_) { throw Exception::OutOfRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION); }
    // delete row j and therefor overwrite with row j+1 and iterate like this to last row
    SizeType i = j + 1;
    while (i < dimensionsize_ && matrix_[i] != nullptr)
    {
      // left out in the copy is each rows jth element, pointer working here as iterators just fine
      std::copy(matrix_[i] + j + 1, matrix_[i] + i, std::copy(matrix_[i], matrix_[i] + j, matrix_[i - 1]));
      ++i;
    }
    // last row is freed and the pointer set to NULL (outer array's size is not changed)
    delete[] matrix_[i - 1];
    matrix_[i - 1] = nullptr;
    --dimensionsize_;
  }
 
  SizeType dimensionsize() const
  {
    return dimensionsize_;
  }
 
  void updateMinElement()
  {
    min_element_ = std::make_pair(1, 0);
    // error if dimensionsize_<1, return if dimensionsize_ == 1, else
    if (dimensionsize_ < 1) { throw Exception::OutOfRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION); }
    if (dimensionsize_ != 1) // else matrix has one element: (1,0)
    {
      ValueType* row_min_;
      for (SizeType r = 2; r < dimensionsize_ && matrix_[r] != nullptr; ++r)
      {
        row_min_ = std::min_element(matrix_[r], matrix_[r] + r);
        if (*row_min_ < matrix_[min_element_.first][min_element_.second]) { min_element_ = std::make_pair(r, row_min_ - matrix_[r]); }
      }
    }
  }
 
  bool operator==(DistanceMatrix<ValueType> const& rhs) const
  {
    if (dimensionsize_ != rhs.dimensionsize_) { return false; }
    for (Size i = 1; i < rhs.dimensionsize(); ++i)
    {
      for (Size j = 0; j < i; ++j)
      {
        if (matrix_[i][j] != rhs.matrix_[i][j]) { return false; }
      }
    }
    return true;
  }
 
  std::pair<SizeType, SizeType> getMinElementCoordinates() const
  {
    if (dimensionsize_ == 0) { throw Exception::OutOfRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION); }
    return min_element_;
  }
 
protected:
  ValueType** matrix_;
  SizeType init_size_; // actual size of outer array
  SizeType dimensionsize_; // number of virtual elements: ((dimensionsize-1)*(dimensionsize))/2
  std::pair<SizeType, SizeType> min_element_;
 
private:
  DistanceMatrix& operator=(const DistanceMatrix& rhs)
  {
    matrix_ = rhs.matrix_;
    init_size_ = rhs.init_size_;
    dimensionsize_ = rhs.dimensionsize_;
    min_element_ = rhs.min_element_;
 
    return *this;
  }
 
}; // class DistanceMatrix
 
template<typename Value>
std::ostream& operator<<(std::ostream& os, const DistanceMatrix<Value>& matrix)
{
  using SizeType = typename DistanceMatrix<Value>::SizeType;
 
  // we need to print a square matrix. So we set the width
  //std::ios_base::fmtflags flag_backup = os.setf(std::ios::scientific); // 'scientific' messes with the width; don't do it
  std::streamsize precision_backup = os.precision(6); // we could go with `writtenDigits<Value>(Value())`, but it becomes unreadable...
  auto width_backup = os.width(8);
 
  for (SizeType i = 0; i < matrix.dimensionsize(); ++i)
  {
    for (SizeType j = 0; j < matrix.dimensionsize(); ++j)
    {
      if (i == j)
      { // color the diagonal in red (conditional, see Colorizer)
        os << red(matrix(i, j)) << '\t';
      }
      else 
      {
        os << matrix(i, j) << '\t';
      }
    }
    os << '\n';
  }
  //os.flags(flag_backup);
  os.precision(precision_backup);
  os.width(width_backup);
  return os;
}
 
} // namespace OpenMS