DistanceMatrix.h

Go to the documentation of this file.

// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2023.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Mathias Walzer $
// $Authors: $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/CONCEPT/Macros.h>
#include <OpenMS/CONCEPT/Types.h>
#include <OpenMS/CONCEPT/Exception.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);
        }
        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
    {
      OPENMS_PRECONDITION(dimensionsize_ == rhs.dimensionsize_, "DistanceMatrices have different sizes.");
      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)
  {
    typedef typename DistanceMatrix<Value>::SizeType SizeType;
 
    std::ios_base::fmtflags flag_backup = os.setf(std::ios::scientific);
    std::streamsize precision_backup = os.precision();
    //~ os.precision(15);
    os.precision(writtenDigits<double>(0.0)); // #include <OpenMS/CONCEPT/Types.h>
 
    //evtl. color lower triangular matrix o.s.l.t.
    for (SizeType i = 0; i < matrix.dimensionsize(); ++i)
    {
      for (SizeType j = 0; j < matrix.dimensionsize(); ++j)
      {
        os << matrix(i, j) << '\t';
      }
      os << std::endl;
    }
    os.flags(flag_backup);
    os.precision(precision_backup);
    return os;
  }
 
} // namespace OpenMS