BamPairIterator2.h

#ifndef theplu_yat_omic_bam_pair_iterator2
#define theplu_yat_omic_bam_pair_iterator2

// $Id: BamPairIterator2.h 4057 2021-03-31 05:46:00Z peter $

/*
  Copyright (C) 2020, 2021 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it 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.

  The yat library 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 yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include "BamPair.h"
#include "BamRead.h"

#include "yat/utility/config_public.h"
#include "yat/utility/yat_assert.h"

#include <boost/concept/assert.hpp>
#include <boost/iterator/iterator_concepts.hpp>
#include <boost/iterator/iterator_facade.hpp>

#include <iterator>
#include <map>
#include <memory>
#include <utility>

namespace theplu {
namespace yat {
namespace omic {

  template<typename Base>
  class BamPairIterator2
    : public boost::iterator_facade<
    BamPairIterator2<Base>, const BamPair, std::input_iterator_tag,
    const BamPairProxy
    >
  {
  public:
    BamPairIterator2(void);

    explicit BamPairIterator2(Base begin, Base end);

#ifndef YAT_HAVE_BOOST_ITERATOR_FACADE_PROXY_PTR

    typename BamPairIterator2::pointer operator->(void) const
    {
      // older versions of boost mandates pointer to be
      // value_type*. To accomplish that we implement this slow thing
      // which keeps a copy of a value_type as member.
      // See https://svn.boost.org/trac/boost/ticket/5697 for discussion.

      // Possibly stupid to assign each time, but why bother optimize
      // for old bugs in boost. Users who are eager for speed, should
      // either upgrade their boost or use (*iterator).function()
      // rather than iterator->function().
      YAT_ASSERT(first_);
      YAT_ASSERT(second_);
      dummie_.first() = first_->first;
      dummie_.second() = second_->second;
      return &dummie_;
    }
  private:
    mutable BamPair dummie_;
#endif

  private:
    Base iter_;
    Base end_;
    typedef std::pair<int32_t, int32_t> Position;
    // multimap with reads' position as key
    typedef std::multimap<Position, BamRead> MultiMap;
    std::shared_ptr<MultiMap> reads_;
    MultiMap::iterator first_;
    MultiMap::iterator second_;
    friend class boost::iterator_core_access;

    void assign_pair(void);
    const BamPairProxy dereference(void) const;
    bool equal(const BamPairIterator2& other) const;
    //
    std::multimap<Position, BamRead>::iterator
    find_mate(const yat::omic::BamRead&);
    void increment(void);
    void pop(void);
    bool push(void);
  };


  template<typename Base>
  BamPairIterator2<Base> bam_pair_iterator2(Base base, Base end)
  { return BamPairIterator2<Base>(base, end); }


  template<typename Base>
  BamPairIterator2<Base> bam_pair_iterator2(Base end)
  { return BamPairIterator2<Base>(end, end); }


  //             template implementations

  template<typename Base>
  BamPairIterator2<Base>::BamPairIterator2(Base base, Base end)
    : iter_(base), end_(end),
      reads_(new std::multimap<Position, BamRead>),
      first_(reads_->end()), second_(reads_->end())
  {
    BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<Base>));
    BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<Base>));
    using boost::Convertible;
    typedef typename std::iterator_traits<Base>::value_type value_type;
    BOOST_CONCEPT_ASSERT((Convertible<value_type, BamRead>));
    assign_pair();
  }


  template<typename Base>
  BamPairIterator2<Base>::BamPairIterator2(void)
  {
    BOOST_CONCEPT_ASSERT((boost::InputIterator<Base>));
    using boost::Convertible;
    typedef typename std::iterator_traits<Base>::value_type value_type;
    BOOST_CONCEPT_ASSERT((Convertible<value_type, BamRead>));
  }


  template<typename Base>
  void BamPairIterator2<Base>::assign_pair(void)
  {
    while (true) {
      if (reads_->empty())
        if (!push()) // if no reads, we are done
          return;
      YAT_ASSERT(reads_->size());
      first_ = reads_->begin();
      second_ = find_mate(first_->second);
      if (second_ == reads_->end())
        pop();
      else
        break;
    }
  }


  template<typename Base>
  const BamPairProxy
  BamPairIterator2<Base>::dereference(void) const
  {
    return BamPairProxy(&first_->second, &second_->second);
  }


  template<typename Base>
  bool BamPairIterator2<Base>::equal(const BamPairIterator2& other) const
  {
    return iter_ == other.iter_ && reads_->size() == other.reads_->size();
  }


  template<typename Base>
  void BamPairIterator2<Base>::increment(void)
  {
    pop();
    assign_pair();
  }


  template<typename Base>
  BamPairIterator2<Base>::MultiMap::iterator
  BamPairIterator2<Base>::find_mate(const yat::omic::BamRead& bam)
  {
    Position mate_pos(bam.mtid(), bam.mpos());
    YAT_ASSERT(reads_->size());
    // push in reads, ensure all reads with pos==mate_pos has been
    // pushed in by pushing until pos of last one is one-past mate_pos
    while (reads_->rbegin()->first <= mate_pos && iter_!=end_)
      push();

    typedef MultiMap::iterator iterator;
    std::pair<iterator, iterator> range = reads_->equal_range(mate_pos);
    for (; range.first != range.second; ++range.first) {
      if (!same_query_name(bam, range.first->second))
        continue;
      // for siamese pair avoid pairing up with self
      if (bam.pos()==bam.mpos() && bam.tid()==bam.mtid() &&
          bam.flag() == range.first->second.flag())
        continue;
      return range.first;
    }
    return reads_->end();
  }


  template<typename Base>
  void BamPairIterator2<Base>::pop(void)
  {
    YAT_ASSERT(reads_->size());
    reads_->erase(first_);
  }


  template<typename Base>
  bool BamPairIterator2<Base>::push(void)
  {
    if (iter_ == end_)
      return false;
    Position position(iter_->tid(), iter_->pos());
    reads_->insert(reads_->end(), std::make_pair(position, *iter_));
    ++iter_;
    return true;
  }

}}}
#endif