msrundatasettree.cpp

Go to the documentation of this file.

// GPL 3+
// Filippo Rusconi
 
#include <map>
#include <limits>
#include <iostream>
#include <iomanip>
 
#include "msrundatasettree.h"
 
#include "../pappsoexception.h"
#include "../exception/exceptionnotpossible.h"
 
#include <QVariant>
 
 
namespace pappso
{
 
 
MsRunDataSetTree::MsRunDataSetTree(MsRunIdCstSPtr ms_run_id_csp) : mcsp_msRunId(ms_run_id_csp)
{
}
 
 
MsRunDataSetTree::~MsRunDataSetTree()
{
  // qDebug();
 
  for(auto &&node : m_rootNodes)
    {
      // Each node is responsible for freeing its children nodes!
 
      delete node;
    }
 
  m_rootNodes.clear();
 
  // Beware not to delete the node member of the map, as we have already
  // destroyed them above!
  //
  // for(auto iterator = m_indexNodeMap.begin(); iterator !=
  // m_indexNodeMap.end();
  //++iterator)
  //{
  // delete(iterator->second);
  //}
 
  // qDebug();
}
 
 
MsRunDataSetTreeNode *
MsRunDataSetTree::addMassSpectrum(QualifiedMassSpectrumCstSPtr mass_spectrum_csp)
{
  qDebug();
 
  if(mass_spectrum_csp == nullptr)
    qFatal("Cannot be nullptr");
 
  if(mass_spectrum_csp.get() == nullptr)
    qFatal("Cannot be nullptr");
 
  // We need to get the precursor spectrum index, in case this spectrum is a
  // fragmentation index.
 
  MsRunDataSetTreeNode *new_node_p = nullptr;
 
  std::size_t precursor_spectrum_index = mass_spectrum_csp->getPrecursorSpectrumIndex();
 
  qDebug() << "The precursor_spectrum_index:" << precursor_spectrum_index;
 
  if(precursor_spectrum_index == std::numeric_limits<std::size_t>::max())
    {
      // This spectrum is a full scan spectrum, not a fragmentation spectrum.
      // Create a new node with no parent and push it back to the root nodes
      // vector.
 
      new_node_p = new MsRunDataSetTreeNode(mass_spectrum_csp, nullptr);
 
      // Since there is no parent in this overload, it is assumed that the node
      // to be populated with the new node is the root node.
 
      m_rootNodes.push_back(new_node_p);
 
      // true: with_data
      // qDebug().noquote() << "Pushed back to the roots node vector node:"
      //<< new_node_p->toString(true);
    }
  else
    {
      // This spectrum is a fragmentation spectrum.
 
      // Sanity check
 
      if(mass_spectrum_csp->getMsLevel() <= 1)
        {
          throw ExceptionNotPossible(
            "msrundatasettree.cpp -- ERROR the MS level needs to be > 1 in a "
            "fragmentation spectrum.");
        }
 
      // Get the node that contains the precursor ion mass spectrum.
      MsRunDataSetTreeNode *parent_node_p = findNode(precursor_spectrum_index);
 
      if(parent_node_p == nullptr)
        {
          throw ExceptionNotPossible(QString("%1 @ %2 - %3\n")
                                       .arg(__FILE__)
                                       .arg(__LINE__)
                                       .arg("msrundatasettree.cpp -- ERROR could not find "
                                            "a tree node matching the index."));
        }
 
      // qDebug() << "Fragmentation spectrum"
      //<< "Found parent node:" << parent_node_p
      //<< "for precursor index:" << precursor_spectrum_index;
 
      // At this point, create a new node with the right parent.
 
      new_node_p = new MsRunDataSetTreeNode(mass_spectrum_csp, parent_node_p);
 
      parent_node_p->m_children.push_back(new_node_p);
    }
 
  // And now document that addition in the node index map.
  m_indexNodeMap.insert(std::pair<std::size_t, MsRunDataSetTreeNode *>(
    mass_spectrum_csp->getMassSpectrumId().getSpectrumIndex(), new_node_p));
 
  // We also want to document the new node relating to the
  // retention time.
 
  documentNodeInDtRtMap(mass_spectrum_csp->getRtInMinutes(), new_node_p, DataKind::rt);
 
  // Likewise for the ion mobility time.
 
  double ion_mobility_value = -1;
 
  MsDataFormat file_format = mcsp_msRunId->getMsDataFormat();
  bool ok                  = false;
 
  if(file_format == MsDataFormat::brukerTims)
    {
      // Start by looking if there is a OneOverK0 valid value, which
      // means we have effectively handled a genuine mobility scan spectrum.
      QVariant ion_mobility_variant_value =
        mass_spectrum_csp->getParameterValue(QualifiedMassSpectrumParameter::IonMobOneOverK0);
 
      if(ion_mobility_variant_value.isValid())
        {
          // Yes, genuine ion mobility scan handled here.
 
          ion_mobility_value = ion_mobility_variant_value.toDouble(&ok);
 
          if(!ok)
            {
              qFatal(
                "The data are Bruker timsTOF data but failed to convert valid "
                "QVariant 1/K0 value to double.");
            }
        }
      else
        {
          // We are not handling a genuine single ion mobility scan here.
          // We must be handling a mass spectrum that correspond to
          // the combination of all the ion mobility scans for a single
          // frame. This is when the user asks that ion mobility scans
          // be flattended. In this case, the OneOverK0 value is not valid
          // but there are two values that are set:
          // TimsFrameInvKoBegin and TimsFrameInvKoEnd.
          // See TimsFramesMsRunReader::readSpectrumCollection2.
 
          // Test one of these values as a sanity check. But
          // give the value of -1 for the ion mobility because we do not
          // have ion mobility data here.
 
          ion_mobility_variant_value = mass_spectrum_csp->getParameterValue(
            QualifiedMassSpectrumParameter::IonMobOneOverK0Begin);
 
          if(!ion_mobility_variant_value.isValid())
            {
              qFatal(
                "The data are Bruker timsTOF data but failed to get correct "
                "ion mobility data. Inconsistency found.");
            }
        }
    }
  else
    {
      ion_mobility_value = mass_spectrum_csp->getDtInMilliSeconds();
    }
 
  if(ion_mobility_value != -1)
    documentNodeInDtRtMap(ion_mobility_value, new_node_p, DataKind::dt);
 
  ++m_spectrumCount;
 
  // qDebug() << "New index/node map:"
  //<< mass_spectrum_csp->getMassSpectrumId().getSpectrumIndex() << "/"
  //<< new_node_p;
 
  return new_node_p;
}
 
 
const std::map<std::size_t, MsRunDataSetTreeNode *> &
MsRunDataSetTree::getIndexNodeMap() const
{
  return m_indexNodeMap;
}
 
 
std::size_t
MsRunDataSetTree::massSpectrumIndex(const MsRunDataSetTreeNode *node) const
{
  // We have a node and we want to get the matching mass spectrum index.
 
  if(node == nullptr)
    throw("Cannot be that the node pointer is nullptr");
 
  std::map<std::size_t, MsRunDataSetTreeNode *>::const_iterator iterator =
    std::find_if(m_indexNodeMap.begin(),
                 m_indexNodeMap.end(),
                 [node](const std::pair<std::size_t, MsRunDataSetTreeNode *> pair) {
                   return pair.second == node;
                 });
 
  if(iterator != m_indexNodeMap.end())
    return iterator->first;
 
  return std::numeric_limits<std::size_t>::max();
}
 
 
std::size_t
MsRunDataSetTree::massSpectrumIndex(QualifiedMassSpectrumCstSPtr qualified_mass_spectrum_csp) const
{
  MsRunDataSetTreeNode *node_p = findNode(qualified_mass_spectrum_csp);
 
  return massSpectrumIndex(node_p);
}
 
 
const std::vector<MsRunDataSetTreeNode *> &
MsRunDataSetTree::getRootNodes() const
{
  return m_rootNodes;
}
 
 
void
MsRunDataSetTree::accept(MsRunDataSetTreeNodeVisitorInterface &visitor)
{
  // qDebug() << "Going to call node->accept(visitor) for each root node.";
 
  for(auto &&node : m_rootNodes)
    {
      // qDebug() << "Calling accept for root node:" << node;
 
      if(visitor.shouldStop())
        break;
 
      node->accept(visitor);
    }
}
 
 
void
MsRunDataSetTree::accept(MsRunDataSetTreeNodeVisitorInterface &visitor,
                         std::vector<MsRunDataSetTreeNode *>::const_iterator nodes_begin_iterator,
                         std::vector<MsRunDataSetTreeNode *>::const_iterator nodes_end_iterator)
{
  // qDebug() << "Visitor:" << &visitor << "The distance is between iterators
  // is:"
  //<< std::distance(nodes_begin_iterator, nodes_end_iterator);
 
  using Iterator = std::vector<MsRunDataSetTreeNode *>::const_iterator;
 
  Iterator iter = nodes_begin_iterator;
 
  // Inform the visitor of the number of nodes to work on.
 
  std::size_t node_count = std::distance(nodes_begin_iterator, nodes_end_iterator);
 
  visitor.setNodesToProcessCount(node_count);
 
  while(iter != nodes_end_iterator)
    {
      // qDebug() << "Visitor:" << &visitor
      //<< "The distance is between iterators is:"
      //<< std::distance(nodes_begin_iterator, nodes_end_iterator);
 
      // qDebug() << "Node visited:" << (*iter)->toString();
 
      if(visitor.shouldStop())
        break;
 
      (*iter)->accept(visitor);
      ++iter;
    }
}
 
 
MsRunDataSetTreeNode *
MsRunDataSetTree::findNode(QualifiedMassSpectrumCstSPtr mass_spectrum_csp) const
{
  // qDebug();
 
  for(auto &node : m_rootNodes)
    {
      // qDebug() << "In one node of the root nodes.";
 
      MsRunDataSetTreeNode *iterNode = node->findNode(mass_spectrum_csp);
      if(iterNode != nullptr)
        return iterNode;
    }
 
  return nullptr;
}
 
 
MsRunDataSetTreeNode *
MsRunDataSetTree::findNode(std::size_t spectrum_index) const
{
  // qDebug();
 
  for(auto &node : m_rootNodes)
    {
      // qDebug() << "In one node of the root nodes.";
 
      MsRunDataSetTreeNode *iterNode = node->findNode(spectrum_index);
      if(iterNode != nullptr)
        return iterNode;
    }
 
  return nullptr;
}
 
 
std::vector<MsRunDataSetTreeNode *>
MsRunDataSetTree::flattenedView()
{
  // We want to push back all the nodes of the tree in a flat vector of nodes.
 
  std::vector<MsRunDataSetTreeNode *> nodes;
 
  for(auto &&node : m_rootNodes)
    {
      // The node will store itself and all of its children.
      node->flattenedView(nodes, true /* with_descendants */);
    }
 
  return nodes;
}
 
 
std::vector<MsRunDataSetTreeNode *>
MsRunDataSetTree::flattenedViewMsLevel(std::size_t ms_level, bool with_descendants)
{
  std::vector<MsRunDataSetTreeNode *> nodes;
 
  // Logically, ms_level cannot be 0.
 
  if(!ms_level)
    {
      throw ExceptionNotPossible("msrundatasettree.cpp -- ERROR the MS level cannot be 0.");
 
      return nodes;
    }
 
  // The depth of the tree at which we are right at this point is 0, we have not
  // gone into the children yet.
 
  std::size_t depth = 0;
 
  // If ms_level is 1, then that means that we want the nodes starting right at
  // the root nodes with or without the descendants.
 
  // std::cout << __FILE__ << " @ " << __LINE__ << " " << __FUNCTION__ << " () "
  //<< "ms_level: " << ms_level << " depth: " << depth << std::endl;
 
  if(ms_level == 1)
    {
      for(auto &&node : m_rootNodes)
        {
          // std::cout << __FILE__ << " @ " << __LINE__ << " " << __FUNCTION__
          //<< " () "
          //<< "Handling one of the root nodes at ms_level = 1."
          //<< std::endl;
 
          node->flattenedView(nodes, with_descendants);
        }
 
      return nodes;
    }
 
  // At this point, we know that we want the descendants of the root nodes since
  // we want ms_level > 1, so we need go to to the children of the root nodes.
 
  // Let depth to 0, because if we go to the children of the root nodes we will
  // still be at depth 0, that is MS level 1.
 
  for(auto &node : m_rootNodes)
    {
      // std::cout
      //<< __FILE__ << " @ " << __LINE__ << " " << __FUNCTION__ << " () "
      //<< std::setprecision(15)
      //<< "Requesting a flattened view of the root's child nodes with depth: "
      //<< depth << std::endl;
 
      node->flattenedViewMsLevelNodes(ms_level, depth, nodes, with_descendants);
    }
 
  return nodes;
}
 
 
MsRunDataSetTreeNode *
MsRunDataSetTree::precursorNodeByProductSpectrumIndex(std::size_t product_spectrum_index)
{
 
  // qDebug();
 
  // Find the node that holds the mass spectrum that was acquired as the
  // precursor that when fragmented gave a spectrum at spectrum_index;
 
  // Get the node that contains the product_spectrum_index first.
  MsRunDataSetTreeNode *node = nullptr;
  node                       = findNode(product_spectrum_index);
 
  // Now get the node that contains the precursor_spectrum_index.
 
  return findNode(node->mcsp_massSpectrum->getPrecursorSpectrumIndex());
}
 
 
std::vector<MsRunDataSetTreeNode *>
MsRunDataSetTree::productNodesByPrecursorSpectrumIndex(std::size_t precursor_spectrum_index)
{
  std::vector<MsRunDataSetTreeNode *> nodes;
 
  // First get the node of the precursor spectrum index.
 
  MsRunDataSetTreeNode *precursor_node = findNode(precursor_spectrum_index);
 
  if(precursor_node == nullptr)
    return nodes;
 
  nodes.assign(precursor_node->m_children.begin(), precursor_node->m_children.end());
 
  return nodes;
}
 
 
std::vector<MsRunDataSetTreeNode *>
MsRunDataSetTree::precursorNodesByPrecursorMz(pappso_double mz, PrecisionPtr precision_ptr)
{
 
  // Find all the precursor nodes holding a mass spectrum that contained a
  // precursor mz-value.
 
  if(precision_ptr == nullptr)
    throw ExceptionNotPossible("msrundatasettree.cpp -- ERROR precision_ptr cannot be nullptr.");
 
  std::vector<MsRunDataSetTreeNode *> product_nodes;
 
  // As a first step, find all the nodes that hold a mass spectrum that was
  // acquired as a fragmentation spectrum of an ion of mz,  that is, search all
  // the product ion nodes for which precursor was mz.
 
  for(auto &&node : m_rootNodes)
    {
      node->productNodesByPrecursorMz(mz, precision_ptr, product_nodes);
    }
 
  // Now, for each node found get the precursor node
 
  std::vector<MsRunDataSetTreeNode *> precursor_nodes;
 
  for(auto &&node : product_nodes)
    {
      precursor_nodes.push_back(findNode(node->mcsp_massSpectrum->getPrecursorSpectrumIndex()));
    }
 
  return precursor_nodes;
}
 
 
bool
MsRunDataSetTree::documentNodeInDtRtMap(double time,
                                        MsRunDataSetTreeNode *node_p,
                                        DataKind data_kind)
{
  // qDebug();
 
  using NodeVector          = std::vector<MsRunDataSetTreeNode *>;
  using DoubleNodeVectorMap = std::map<double, NodeVector>;
  using MapPair             = std::pair<double, NodeVector>;
  using MapIterator         = DoubleNodeVectorMap::iterator;
 
  DoubleNodeVectorMap *map_p;
 
  if(data_kind == DataKind::rt)
    {
      map_p = &m_rtDoubleNodeVectorMap;
    }
  else if(data_kind == DataKind::dt)
    {
      map_p = &m_dtDoubleNodeVectorMap;
    }
  else
    qFatal("Programming error.");
 
  // There are two possibilities:
  //
  // 1. The time was never encountered yet. We won't find it. We need to
  // allocate a vector of Node's and set it associated to time in the map.
  //
  // 2. The time was encountered already, we will find it in the maps, we'll
  // just push_back the Node in the vector of nodes.
 
  MapIterator found_iterator = map_p->find(time);
 
  if(found_iterator != map_p->end())
    {
      // The time value was encountered already.
 
      found_iterator->second.push_back(node_p);
 
      // qDebug() << "Found iterator for time:" << time;
    }
  else
    {
      // We need to create a new vector with the node.
 
      NodeVector node_vector = {node_p};
 
      map_p->insert(MapPair(time, node_vector));
 
      // qDebug() << "Inserted new time:node_vector pair.";
    }
 
  return true;
}
 
 
MsRunDataSetTreeNode *
MsRunDataSetTree::addMassSpectrum(QualifiedMassSpectrumCstSPtr mass_spectrum_csp,
                                  MsRunDataSetTreeNode *parent_p)
{
  // qDebug();
 
  // We want to add a mass spectrum. Either the parent_p argument is nullptr or
  // not. If it is nullptr, then we just append the mass spectrum to the vector
  // of root nodes. If it is not nullptr, we need to append the mass spectrum to
  // that node.
 
  MsRunDataSetTreeNode *new_node_p = new MsRunDataSetTreeNode(mass_spectrum_csp, parent_p);
 
  if(parent_p == nullptr)
    {
      m_rootNodes.push_back(new_node_p);
 
      // qDebug() << "Pushed back" << new_node << "to root nodes:" <<
      // &m_rootNodes;
    }
  else
    {
      parent_p->m_children.push_back(new_node_p);
 
      // qDebug() << "Pushed back" << new_node << "with parent:" << parent_p;
    }
 
  ++m_spectrumCount;
 
  // And now document that addition in the node index map.
  m_indexNodeMap.insert(std::pair<std::size_t, MsRunDataSetTreeNode *>(
    mass_spectrum_csp->getMassSpectrumId().getSpectrumIndex(), new_node_p));
 
  // We also want to document the new node relating to the
  // retention time.
 
  documentNodeInDtRtMap(mass_spectrum_csp->getRtInMinutes(), new_node_p, DataKind::rt);
 
  // Likewise for the ion mobility time.
 
  double ion_mobility_value = -1;
 
  MsDataFormat file_format = mcsp_msRunId->getMsDataFormat();
  bool ok                  = false;
 
  if(file_format == MsDataFormat::brukerTims)
    {
      // Start by looking if there is a OneOverK0 valid value, which
      // means we have effectively handled a genuine mobility scan spectrum.
      QVariant ion_mobility_variant_value =
        mass_spectrum_csp->getParameterValue(QualifiedMassSpectrumParameter::IonMobOneOverK0);
 
      if(ion_mobility_variant_value.isValid())
        {
          // Yes, genuine ion mobility scan handled here.
 
          ion_mobility_value = ion_mobility_variant_value.toDouble(&ok);
 
          if(!ok)
            {
              qFatal(
                "The data are Bruker timsTOF data but failed to convert valid "
                "QVariant 1/K0 value to double.");
            }
        }
      else
        {
          // We are not handling a genuine single ion mobility scan here.
          // We must be handling a mass spectrum that correspond to
          // the combination of all the ion mobility scans for a single
          // frame. This is when the user asks that ion mobility scans
          // be flattended. In this case, the OneOverK0 value is not valid
          // but there are two values that are set:
          // TimsFrameInvKoBegin and TimsFrameInvKoEnd.
          // See TimsFramesMsRunReader::readSpectrumCollection2.
 
          // Test one of these values as a sanity check. But
          // give the value of -1 for the ion mobility because we do not
          // have ion mobility data here.
 
          ion_mobility_variant_value = mass_spectrum_csp->getParameterValue(
            QualifiedMassSpectrumParameter::IonMobOneOverK0Begin);
 
          if(!ion_mobility_variant_value.isValid())
            {
              qFatal(
                "The data are Bruker timsTOF data but failed to get correct "
                "ion mobility data. Inconsistency found.");
            }
        }
    }
  else
    {
      ion_mobility_value = mass_spectrum_csp->getDtInMilliSeconds();
    }
 
  if(ion_mobility_value != -1)
    documentNodeInDtRtMap(ion_mobility_value, new_node_p, DataKind::dt);
 
  // qDebug() << "New index/node map:"
  //<< mass_spectrum_csp->getMassSpectrumId().getSpectrumIndex() << "/"
  //<< new_node;
 
  return new_node_p;
}
 
 
MsRunDataSetTreeNode *
MsRunDataSetTree::addMassSpectrum(QualifiedMassSpectrumCstSPtr mass_spectrum_csp,
                                  std::size_t precursor_spectrum_index)
{
  // qDebug();
 
  // First get the node containing the mass spectrum that was acquired at index
  // precursor_spectrum_index.
 
  // qDebug() << "Need to find the precursor's mass spectrum node for precursor
  // "
  //"spectrum index:"
  //<< precursor_spectrum_index;
 
  MsRunDataSetTreeNode *mass_spec_data_node_p = findNode(precursor_spectrum_index);
 
  // qDebug() << "Found node" << mass_spec_data_node_p
  //<< "for precursor index:" << precursor_spectrum_index;
 
  if(mass_spec_data_node_p == nullptr)
    {
      throw ExceptionNotPossible(
        "msrundatasettree.cpp -- ERROR could not find a a "
        "tree node matching the index.");
    }
 
  // qDebug() << "Calling addMassSpectrum with parent node:"
  //<< mass_spec_data_node_p;
 
  return addMassSpectrum(mass_spectrum_csp, mass_spec_data_node_p);
}
 
 
std::size_t
MsRunDataSetTree::addDataSetTreeNodesInsideDtRtRange(double start,
                                                     double end,
                                                     NodeVector &nodes,
                                                     DataKind data_kind) const
{
  using NodeVector          = std::vector<MsRunDataSetTreeNode *>;
  using DoubleNodeVectorMap = std::map<double, NodeVector>;
  using MapIterator         = DoubleNodeVectorMap::const_iterator;
 
  const DoubleNodeVectorMap *map_p;
 
  if(data_kind == DataKind::rt)
    {
      map_p = &m_rtDoubleNodeVectorMap;
    }
  else if(data_kind == DataKind::dt)
    {
      map_p = &m_dtDoubleNodeVectorMap;
    }
  else
    qFatal("Programming error.");
 
  std::size_t added_nodes = 0;
 
  // Get the iterator to the map item that has the key greater or equal to
  // start.
 
  MapIterator start_iterator = map_p->lower_bound(start);
 
  if(start_iterator == map_p->end())
    return 0;
 
  // Now get the end of the map useful range of items.
 
  MapIterator end_iterator = map_p->upper_bound(end);
 
  // Now that we have the iterator range, iterate in it and get the mass spectra
  // from each item's pair.second node vector.
 
  for(MapIterator iterator = start_iterator; iterator != end_iterator; ++iterator)
    {
      // We are iterating in MapPair items.
 
      NodeVector node_vector = iterator->second;
 
      // All the nodes in the node vector need to be copied to the mass_spectra
      // vector passed as parameter.
 
      for(auto &&node_p : node_vector)
        {
          nodes.push_back(node_p);
 
          ++added_nodes;
        }
    }
 
  return added_nodes;
}
 
 
std::size_t
MsRunDataSetTree::removeDataSetTreeNodesOutsideDtRtRange(double start,
                                                         double end,
                                                         NodeVector &nodes,
                                                         DataKind data_kind) const
{
  using NodeVector         = std::vector<MsRunDataSetTreeNode *>;
  using NodeVectorIterator = NodeVector::iterator;
 
  using DoubleNodeVectorMap = std::map<double, NodeVector>;
  using MapIterator         = DoubleNodeVectorMap::const_iterator;
 
  const DoubleNodeVectorMap *map_p;
 
  if(data_kind == DataKind::rt)
    {
      map_p = &m_rtDoubleNodeVectorMap;
    }
  else if(data_kind == DataKind::dt)
    {
      map_p = &m_dtDoubleNodeVectorMap;
    }
  else
    qFatal("Programming error.");
 
  std::size_t removed_vector_items = 0;
 
  // We want to remove from the nodes vector all the nodes that contain a mass
  // spectrum acquired at a time range outside of [ start-end ], that is, the
  // time values [begin() - start [ and ]end -- end()[.
 
  // Get the iterator to the map item that has the key less to
  // start (we want to keep the map item having key == start).
 
  MapIterator first_end_iterator = (*map_p).upper_bound(start);
 
  // Now that we have the first_end_iterator, we can iterate between [begin --
  // first_end_iterator[
 
  for(MapIterator iterator = map_p->begin(); iterator != first_end_iterator; ++iterator)
    {
      // Remove from the nodes vector the nodes.
 
      // We are iterating in MapPair items.
 
      NodeVector node_vector = iterator->second;
 
      // All the nodes in the node vector need to be removed from the
      // mass_spectra vector passed as parameter if found.
 
      for(auto &&node_p : node_vector)
        {
          NodeVectorIterator iterator = std::find(nodes.begin(), nodes.end(), node_p);
 
          if(iterator != nodes.end())
            {
              // We found the node: remove it.
 
              nodes.erase(iterator);
 
              ++removed_vector_items;
            }
        }
    }
 
  // Now the second begin iterator, so that we can remove all the items
  // contained in the second range, that is, ]end--end()[.
 
  // The second_first_iterator will point to the item having its time value less
  // or equal to end. But we do not want to get items having their time equal to
  // end, only < end. So, if the iterator is not begin(), we just need to
  // decrement it once.
  MapIterator second_first_iterator = map_p->upper_bound(end);
  if(second_first_iterator != map_p->begin())
    --second_first_iterator;
 
  for(MapIterator iterator = second_first_iterator; iterator != map_p->end(); ++iterator)
    {
      // We are iterating in MapPair items.
 
      NodeVector node_vector = iterator->second;
 
      // All the nodes in the node vector need to be removed from the
      // mass_spectra vector passed as parameter if found.
 
      for(auto &&node_p : node_vector)
        {
          NodeVectorIterator iterator = std::find(nodes.begin(), nodes.end(), node_p);
 
          if(iterator != nodes.end())
            {
              // We found the node: remove it.
 
              nodes.erase(iterator);
 
              ++removed_vector_items;
            }
        }
    }
 
  return removed_vector_items;
}
 
 
std::size_t
MsRunDataSetTree::addDataSetQualMassSpectraInsideDtRtRange(double start,
                                                           double end,
                                                           QualMassSpectraVector &mass_spectra,
                                                           DataKind data_kind) const
{
  // qDebug() << "With start:" << start << "and end:" << end;
 
  if(start == end)
    qDebug() << "Special case, start and end are equal:" << start;
 
  // We will use the maps that relate rt | dt to a vector of data tree nodes.
  // Indeed, we may have more than one mass spectrum acquired for a given rt, in
  // case of ion mobility mass spectrometry. Same for dt: we will have as many
  // spectra for each dt as there are retention time values...
 
  using DoubleNodeVectorMap = std::map<double, NodeVector>;
  using MapIterator         = DoubleNodeVectorMap::const_iterator;
 
  const DoubleNodeVectorMap *map_p;
 
  if(data_kind == DataKind::rt)
    {
      map_p = &m_rtDoubleNodeVectorMap;
 
      // qDebug() << "The RT map has size:" << map_p->size() << "start:" <<
      // start
      //<< "end:" << end;
    }
  else if(data_kind == DataKind::dt)
    {
      map_p = &m_dtDoubleNodeVectorMap;
 
      // qDebug() << "The DT map has size:" << map_p->size() << "start:" <<
      // start
      //<< "end:" << end;
    }
  else
    qFatal("Programming error.");
 
  // qDebug() << "The rt |dt / mass spectra map has size:" << map_p->size()
  //<< "The start:" << start << "the end:" << end;
 
  std::size_t added_mass_spectra = 0;
 
  // Get the iterator to the map item that has the key greater or equal to
  // start.
 
  MapIterator start_iterator = map_p->lower_bound(start);
 
  if(start_iterator == map_p->end())
    {
      qDebug() << "The start iterator is end()!";
      return 0;
    }
 
  // qDebug() << "The start_iterator points to:" << start_iterator->first
  //<< "as a rt|dt time.";
 
  // Now get the end of the map's useful range of items.
 
  // Returns an iterator pointing to the first element in the container whose
  // key is considered to go after 'end'.
 
  MapIterator end_iterator = map_p->upper_bound(end);
 
  // Immediately verify if there is no distance between start and end.
  if(!std::distance(start_iterator, end_iterator))
    {
      qDebug() << "No range of mass spectra could be selected.";
      return 0;
    }
 
  if(end_iterator == map_p->end())
    {
      // qDebug() << "The end_iterator points to the end of the map."
      //<< "The last map item is prev() at key value: "
      //<< std::prev(end_iterator)->first;
    }
  else
    {
      // qDebug() << "The end_iterator points to:" << end_iterator->first
      //<< "as a rt|dt time and the accounted key value is actually"
      //<< std::prev(end_iterator)->first;
    }
 
  // qDebug() << "The number of time values to iterate through:"
  //<< std::distance(start_iterator, end_iterator)
  //<< "with values: start: " << start_iterator->first
  //<< "and end: " << std::prev(end_iterator)->first;
 
  // Now that we have the iterator range, iterate in it and get the mass
  // spectra from each item's pair.second node vector.
 
  for(MapIterator iterator = start_iterator; iterator != end_iterator; ++iterator)
    {
      // We are iterating in DoubleNodeVectorMap MapPair items,
      // with double = rt or dt and NodeVector being just that.
 
      NodeVector node_vector = iterator->second;
 
      // All the nodes' mass spectra in the node vector need to be copied to
      // the mass_spectra vector passed as parameter.
 
      for(auto &&node_p : node_vector)
        {
          QualifiedMassSpectrumCstSPtr qualified_mass_spectrum_csp =
            node_p->getQualifiedMassSpectrum();
 
#if 0
          // Sanity check only for deep debugging.
 
          if(qualified_mass_spectrum_csp == nullptr ||
             qualified_mass_spectrum_csp.get() == nullptr)
            {
              throw ExceptionNotPossible(
                "The QualifiedMassSpectrumCstSPtr cannot be nullptr.");
            }
          else
            {
              //qDebug() << "Current mass spectrum is valid with rt:"
                       //<< qualified_mass_spectrum_csp->getRtInMinutes();
            }
#endif
 
          mass_spectra.push_back(qualified_mass_spectrum_csp);
 
          ++added_mass_spectra;
        }
    }
 
  // qDebug() << "Returning added_mass_spectra:" << added_mass_spectra;
 
  return added_mass_spectra;
}
 
 
std::size_t
MsRunDataSetTree::removeDataSetQualMassSpectraOutsideDtRtRange(double start,
                                                               double end,
                                                               QualMassSpectraVector &mass_spectra,
                                                               DataKind data_kind) const
{
  using QualMassSpectraVectorIterator = QualMassSpectraVector::iterator;
 
  using DoubleNodeVectorMap = std::map<double, NodeVector>;
  using MapIterator         = DoubleNodeVectorMap::const_iterator;
 
  const DoubleNodeVectorMap *map_p;
 
  if(data_kind == DataKind::rt)
    {
      map_p = &m_rtDoubleNodeVectorMap;
 
      // qDebug() << "The RT map has size:" << map_p->size() << "start:" <<
      // start
      //<< "end:" << end;
    }
  else if(data_kind == DataKind::dt)
    {
      map_p = &m_dtDoubleNodeVectorMap;
 
      // qDebug() << "The DT map has size:" << map_p->size() << "start:" <<
      // start
      //<< "end:" << end;
    }
  else
    qFatal("Programming error.");
 
  std::size_t removed_vector_items = 0;
 
  // We want to remove from the nodes vector all the nodes that contain a mass
  // spectrum acquired at a time range outside of [ start-end ], that is, the
  // time values [begin() - start [ and ]end -- end()[.
 
  // Looking for an iterator that points to an item having a time < start.
 
  // lower_bound returns an iterator pointing to the first element in the
  // range [first, last) that is not less than (i.e. greater or equal to)
  // value, or last if no such element is found.
 
  MapIterator first_end_iterator = (*map_p).lower_bound(start);
 
  // first_end_iterator points to the item that has the next time value with
  // respect to start. This is fine because we'll not remove that point
  // because the for loop below will stop one item short of
  // first_end_iterator. That means that we effectively remove all the items
  // [begin() -> start[ (start not include). Exactly what we want.
 
  // qDebug() << "lower_bound for start:" << first_end_iterator->first;
 
  // Now that we have the first_end_iterator, we can iterate between [begin --
  // first_end_iterator[
 
  for(MapIterator iterator = map_p->begin(); iterator != first_end_iterator; ++iterator)
    {
      // Remove from the nodes vector the nodes.
 
      // We are iterating in MapPair items.
 
      NodeVector node_vector = iterator->second;
 
      // All the nodes in the node vector need to be removed from the
      // mass_spectra vector passed as parameter if found.
 
      for(auto &&node_p : node_vector)
        {
          QualMassSpectraVectorIterator iterator =
            std::find(mass_spectra.begin(), mass_spectra.end(), node_p->getQualifiedMassSpectrum());
 
          if(iterator != mass_spectra.end())
            {
              // We found the mass spectrum: remove it.
 
              mass_spectra.erase(iterator);
 
              ++removed_vector_items;
            }
        }
    }
 
  // Now the second begin iterator, so that we can remove all the items
  // contained in the second range, that is, ]end--end()[.
 
  // The second_first_iterator will point to the item having its time value
  // less or equal to end. But we do not want to get items having their time
  // equal to end, only < end. So, if the iterator is not begin(), we just
  // need to decrement it once.
 
  MapIterator second_first_iterator = map_p->upper_bound(end);
 
  // second_first_iterator now points to the item after the one having time
  // end. Which is exactly what we want: we want to remove ]end--end()[ and
  // this is exactly what the loop starting a the point after end below.
 
  // qDebug() << "second_first_iterator for end:" <<
  // second_first_iterator->first;
 
  for(MapIterator iterator = second_first_iterator; iterator != map_p->end(); ++iterator)
    {
      // We are iterating in MapPair items.
 
      NodeVector node_vector = iterator->second;
 
      // All the nodes in the node vector need to be removed from the
      // mass_spectra vector passed as parameter if found.
 
      for(auto &&node_p : node_vector)
        {
          QualMassSpectraVectorIterator iterator =
            std::find(mass_spectra.begin(), mass_spectra.end(), node_p->getQualifiedMassSpectrum());
 
          if(iterator != mass_spectra.end())
            {
              // We found the node: remove it.
 
              mass_spectra.erase(iterator);
 
              ++removed_vector_items;
            }
        }
    }
 
  return removed_vector_items;
}
 
 
std::size_t
MsRunDataSetTree::depth() const
{
  // We want to know what is the depth of the tree, that is the highest level
  // of MSn, that is, n.
 
  if(!m_rootNodes.size())
    return 0;
 
  // qDebug() << "There are" << m_rootNodes.size() << "root nodes";
 
  // By essence, we are at MS0: only if we have at least one root node do we
  // know we have MS1 data. So we already know that we have at least one
  // child, so start with depth 1.
 
  std::size_t depth          = 1;
  std::size_t tmp_depth      = 0;
  std::size_t greatest_depth = 0;
 
  for(auto &node : m_rootNodes)
    {
      tmp_depth = node->depth(depth);
 
      // qDebug() << "Returned depth:" << tmp_depth;
 
      if(tmp_depth > greatest_depth)
        greatest_depth = tmp_depth;
    }
 
  return greatest_depth;
}
 
 
std::size_t
MsRunDataSetTree::size() const
{
 
  std::size_t cumulative_node_count = 0;
 
  for(auto &node : m_rootNodes)
    {
      node->size(cumulative_node_count);
 
      // qDebug() << "Returned node_count:" << node_count;
    }
 
  return cumulative_node_count;
}
 
 
std::size_t
MsRunDataSetTree::indexNodeMapSize() const
{
  return m_indexNodeMap.size();
}
 
 
std::size_t
MsRunDataSetTree::getSpectrumCount() const
{
  return m_spectrumCount;
}
 
 
} // namespace pappso