spectree.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/processing/spectree/spectree.cpp
 * \date 11/12/2023
 * \author Olivier Langella
 * \brief Matthieu David's SpecTree structure
 *
 * C++ implementation of algorithm already described in :
 * 1. David, M., Fertin, G., Rogniaux, H. & Tessier, D. SpecOMS: A Full Open
 * Modification Search Method Performing All-to-All Spectra Comparisons within
 * Minutes. J. Proteome Res. 16, 3030–3038 (2017).
 *
 * https://www.theses.fr/2019NANT4092
 */
 
 
/*
 * SpecTree
 * Copyright (C) 2023  Olivier Langella
 * <olivier.langella@universite-paris-saclay.fr>
 *
 * This program is free software: you can redistribute ipetide to spectrum
 * alignmentt and/or modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
 
#include "spectree.h"
#include <QDebug>
#include <QObject>
#include "pappsomspp/core/exception/exceptioninterrupted.h"
#include "pappsomspp/core/exception/exceptionoutofrange.h"
 
 
namespace pappso
{
namespace spectree
{
 
SpecTree::SpecTree(const BucketClustering &bucket_clustering)
{
 
 
  std::vector<Bucket> bins = bucket_clustering.asSortedList();
  SpecTreeNode node;
  node.value = bins[0].getCartList().at(0);
  node.count = 1;
  qDebug() << node.value;
  // Specific processing for the first bucket : the first node and all its
  // children are added in SpecTrees and the transversal accessor is updated
  // accordingly
 
 
  std::vector<std::size_t> spectrumLastNodeIndex;
  spectrumLastNodeIndex.resize(bucket_clustering.getItemCartCount(),
                               index_not_defined);
  m_spectrumFirstNodeIndex.resize(bucket_clustering.getItemCartCount(),
                                  index_not_defined);
 
  m_spectrumFirstNodeIndex[node.value] = m_nodeList.size() - 1;
 
 
  m_nodeList.push_back(node);
  manageSideAccess(spectrumLastNodeIndex);
  for(std::size_t cpt = 1; cpt < bins[0].size(); ++cpt)
    {
      node.parentIndex = m_nodeList.size() - 1;
      node.nextIndex   = index_not_defined;
      node.value       = bins[0].getCartList().at(cpt);
      node.count       = 1;
      m_nodeList.push_back(node);
      manageSideAccess(spectrumLastNodeIndex);
    }
 
 
  qDebug();
  // For the rest of the buckets, we increment the node counter in prefix
  // common parts without creating new nodes. Once out of the prefix common
  // part, new nodes are created similarily as they were for the first bucket.
  for(std::size_t bins_index = 1; bins_index < bins.size(); ++bins_index)
    {
      qDebug() << "bins_index=" << bins_index;
      std::size_t pos        = 0;
      bool common            = false;
      std::size_t parentNode = index_not_defined;
 
      std::size_t previous_bin_size = bins[bins_index - 1].size();
 
      // detection of a common first element
      if(pos < previous_bin_size)
        {
          if(bins[bins_index].getCartList().at(pos) ==
             bins[bins_index - 1].getCartList().at(pos))
            {
              common = true;
            }
        }
      // insertion procedure for a common prefix part insertion : we go
      // through existing nodes and increment their counter value by one
      while(common)
        {
          qDebug() << "common pos=" << pos;
          std::size_t spectrum_index = bins[bins_index].getCartList().at(pos);
          qDebug() << "common spectrum_index=" << spectrum_index;
          std::size_t tempNodeIndex = spectrumLastNodeIndex[spectrum_index];
          qDebug() << "common tempNodeIndex=" << tempNodeIndex;
          m_nodeList[tempNodeIndex].count++;
          // checking on the prefix common part persistence conditions ; if
          // not valid, we are not in a common prefix part anymore for the
          // currently inserted bucket
          if((previous_bin_size - 1 == pos) ||
             (bins[bins_index].size() - 1 == pos))
            {
              qDebug() << "last one";
              common = false;
            }
          else
            {
              qDebug() << "bins[bins_index - 1].getSpectrumIndex(pos + 1)="
                       << bins[bins_index - 1].getCartList().at(pos + 1);
              qDebug() << "bins[bins_index].getSpectrumIndex(pos + 1)="
                       << bins[bins_index].getCartList().at(pos + 1);
              if(bins[bins_index - 1].getCartList().at(pos + 1) !=
                 bins[bins_index].getCartList().at(pos + 1))
                {
                  common = false;
                }
            }
          parentNode = tempNodeIndex;
          ++pos;
 
          qDebug() << "common=" << common;
        }
      qDebug();
      // insertion procedure for all the remaining nodes not in the prefix
      // common part : we create new nodes with a counter value of 1 and
      // update the transversal accessor and last inserted node list
      while(pos < bins[bins_index].size())
        {
          std::size_t spectrum_index = bins[bins_index].getCartList().at(pos);
 
          node.parentIndex = parentNode;
          node.nextIndex   = index_not_defined;
          node.value       = spectrum_index;
          node.count       = 1;
          m_nodeList.push_back(node);
          manageSideAccess(spectrumLastNodeIndex);
          parentNode = m_nodeList.size() - 1;
 
          ++pos;
        }
    }
}
 
SpecTree::~SpecTree()
{
}
 
void
SpecTree::addNewNode(const SpecTree::SpecTreeNode &node)
{
  m_nodeList.push_back(node);
}
 
void
SpecTree::manageSideAccess(std::vector<std::size_t> &spectrumLastNodeIndex)
{
  auto spectrum_index = m_nodeList.back().value;
  auto node_position  = m_nodeList.size() - 1;
  if(m_spectrumFirstNodeIndex[spectrum_index] == index_not_defined)
    {
      // first occurence of this spectrum
      m_spectrumFirstNodeIndex[spectrum_index] = node_position;
      spectrumLastNodeIndex[spectrum_index]    = node_position;
    }
  else
    {
      if(spectrumLastNodeIndex[spectrum_index] == node_position)
        {
        }
      else
        {
          m_nodeList[spectrumLastNodeIndex[spectrum_index]].nextIndex =
            node_position;
          spectrumLastNodeIndex[spectrum_index] = node_position;
        }
    }
}
 
QString
SpecTree::toString() const
{
  QString node_str("node|spectrum|parent|next|count|\n");
  std::size_t i = 0;
  for(auto node : m_nodeList)
    {
      int parent = -1;
      if(node.parentIndex < index_not_defined)
        parent = node.parentIndex;
      int next = -1;
      if(node.nextIndex < index_not_defined)
        next = node.nextIndex;
      node_str.append(QString("node_%1|%2|%3|%4|%5end\n")
                        .arg(i)
                        .arg(node.value)
                        .arg(parent)
                        .arg(next)
                        .arg(node.count));
 
      i++;
    }
  return node_str;
}
 
const std::vector<std::size_t> &
SpecTree::getSpectrumFirstNodeIndex() const
{
  return m_spectrumFirstNodeIndex;
}
 
void
SpecTree::xtract(UiMonitorInterface &monitor,
                 SpecXtractInterface &spec_xtract,
                 std::size_t minimum_count,
                 std::size_t cart_id_range_max,
                 std::size_t cart_id_range_min,
                 std::size_t target_cart_id_max,
                 std::size_t target_cart_id_min) const
{
 
  if(cart_id_range_max < cart_id_range_min)
    {
      std::swap(cart_id_range_max, cart_id_range_min);
    }
 
  if(cart_id_range_max < m_spectrumFirstNodeIndex.size())
    {
 
      monitor.setStatus(QObject::tr("starting SpecXtract operation"));
 
      MapSimilarityCount map_count;
      map_count.map_id_count.resize(cart_id_range_max);
 
      // initialisation of the structure to store the compared pairs
      // int[][] results = new int[m_count][3];
      // int count       = 0;
 
      // iteration on the transversal accessor starting from the deepest spetrum
      // and climbing up to limit
      std::size_t spectra1 = cart_id_range_max;
      monitor.setTotalSteps(cart_id_range_max - cart_id_range_min);
      while(true)
        {
 
          if(monitor.shouldIstop())
            {
              throw pappso::ExceptionInterrupted(
                QObject::tr("SpecXtract job stopped"));
            }
          qDebug() << "spectra1=" << spectra1;
          spec_xtract.beginItemCartExtraction(spectra1);
          monitor.count();
          // System.out.println(start);
 
          map_count.keys.clear();
          map_count.aboveThreshold.clear();
 
          extractSpectrumSimilarityCount(map_count,
                                         minimum_count,
                                         spectra1,
                                         target_cart_id_max,
                                         target_cart_id_min);
 
          qDebug() << "spectra1=" << spectra1
                   << " map_count.aboveThreshold.size()="
                   << map_count.aboveThreshold.size()
                   << " map_count.keys.size()=" << map_count.keys.size();
 
          for(std::size_t spectra2 : map_count.aboveThreshold)
            {
 
              // std::size_t spectra2 = *it;
              qDebug() << "spectra2=" << spectra2;
 
              // std::size_t total_extracted_count =
              // map_count.map_id_count.at(spectra2);
              spec_xtract.reportSimilarity(
                spectra1, spectra2, map_count.map_id_count.at(spectra2).count);
            }
 
          spec_xtract.endItemCartExtraction(spectra1);
          // monitor.count();
          if(spectra1 == cart_id_range_min)
            break;
          spectra1--;
        }
    }
  else
    {
      throw pappso::ExceptionOutOfRange(
        QObject::tr(
          "item cart index %1 out of range (item cart id max size=%2)")
          .arg(cart_id_range_max)
          .arg(m_spectrumFirstNodeIndex.size()));
    }
  qDebug();
}
 
void
SpecTree::walkBackInBranchFromNode(const SpecTree::SpecTreeNode &start_node,
                                   MapSimilarityCount &map_count,
                                   std::size_t minimum_count,
                                   std::size_t target_cart_id_max,
                                   std::size_t target_cart_id_min) const
{
 
  std::size_t parent_id = start_node.parentIndex;
 
  qDebug() << " start node spectrum=" << start_node.value
           << " parent_id=" << parent_id << " init_count=" << start_node.count
           << " minimum_count=" << minimum_count
           << " target_cart_id_min=" << target_cart_id_min;
 
  while(parent_id != index_not_defined)
    {
      // next_node :
      const SpecTree::SpecTreeNode &current_node = m_nodeList[parent_id];
 
      if(current_node.value < target_cart_id_min)
        break;
      if(current_node.value <= target_cart_id_max)
        {
          auto &map_id_count = map_count.map_id_count[current_node.value];
          if(map_id_count.lastWitness == start_node.value)
            {
              map_id_count.count += start_node.count;
              // if at one time, a similarity exceeds the threshold value for
              // the first time, it is written into the proper stack
              if((map_id_count.count >= minimum_count) &&
                 (map_id_count.aboveThreshold == false))
                {
                  map_id_count.aboveThreshold = true;
                  map_count.aboveThreshold.push_back(current_node.value);
                }
            }
          else
            {
              map_id_count.lastWitness    = start_node.value;
              map_id_count.count          = start_node.count;
              map_id_count.aboveThreshold = false;
              // if at one time, a similarity exceeds the threshold value for
              // the first time, it is written into the proper stack
              if(map_id_count.count >= minimum_count)
                {
                  map_id_count.aboveThreshold = true;
                  map_count.aboveThreshold.push_back(current_node.value);
                }
              map_count.keys.push_back(current_node.value);
            }
        }
      parent_id = current_node.parentIndex;
    }
}
 
std::size_t
SpecTree::walkBackInBranchFromNodeToTarget(
  const SpecTree::SpecTreeNode &start_node,
  std::size_t spectrum_index_target) const
{
  std::size_t parent_id = start_node.parentIndex;
 
  while(parent_id != index_not_defined)
    {
      // next_node :
      const SpecTree::SpecTreeNode &current_node = m_nodeList[parent_id];
 
      if(current_node.value < spectrum_index_target)
        return 0;
      if(current_node.value == spectrum_index_target)
        return start_node.count;
 
      parent_id = current_node.parentIndex;
    }
  return 0;
}
 
 
/**
 * Extraction of the similarities above a given threshold between a given
 * spectrum and all other spectra higher located in SpecTrees.
 * @param tree_indices SpecTrees the data has to be extracted from
 * @param spectrum_index The spectrum with ated in SpecTrees.
 * @param tree_indices SpecTrees the data has to be extracted from
 * @param spectrum_index The spectrum with which we want to extract similarities
 * @param minimum_count The threshold to use for the extraction
 * @since 0.1
 */
void
SpecTree::extractSpectrumSimilarityCount(MapSimilarityCount &map_count,
                                         std::size_t minimum_count,
                                         std::size_t spectrum_index,
                                         std::size_t target_cart_id_max,
                                         std::size_t target_cart_id_min) const
{
 
  qDebug() << " spectrum_index=" << spectrum_index;
  // sOne;
  // we position the current node at the first occurence of the spectrum in
  // the transversal accessor and initialise variablethresholds accordingly
  std::size_t node_index = m_spectrumFirstNodeIndex[spectrum_index];
  while(node_index != index_not_defined)
    {
      qDebug() << "node_index=" << node_index;
      const SpecTreeNode &current_node = m_nodeList[node_index];
      // go back in the tree branch from the current_node
      walkBackInBranchFromNode(current_node,
                               map_count,
                               minimum_count,
                               target_cart_id_max,
                               target_cart_id_min);
      node_index = current_node.nextIndex;
    }
  qDebug();
}
 
std::size_t
SpecTree::extractSpectrumPairSimilarityCount(std::size_t spectrum_a_index,
                                             std::size_t spectrum_b_index) const
{
  if(spectrum_a_index > spectrum_b_index)
    std::swap(spectrum_a_index, spectrum_b_index);
 
  std::size_t count_target = 0;
  if(spectrum_b_index < m_spectrumFirstNodeIndex.size())
    {
      std::size_t node_index = m_spectrumFirstNodeIndex[spectrum_b_index];
      while(node_index != index_not_defined)
        {
          qDebug() << "node_index=" << node_index;
          const SpecTreeNode &current_node = m_nodeList[node_index];
          // go back in the tree branch from the current_node
          count_target +=
            walkBackInBranchFromNodeToTarget(current_node, spectrum_a_index);
          node_index = current_node.nextIndex;
        }
    }
  else
    {
      throw pappso::ExceptionOutOfRange(
        QObject::tr(
          "spectrum_index %1 out of range (spectrum_index max size=%2)")
          .arg(spectrum_b_index)
          .arg(m_spectrumFirstNodeIndex.size()));
    }
  return count_target;
}
} // namespace spectree
} // namespace pappso