BamPairIterator.h

#ifndef theplu_yat_omic_bam_pair_iterator
#define theplu_yat_omic_bam_pair_iterator

// $Id: BamPairIterator.h 3550 2017-01-03 05:41:02Z peter $

/*
  Copyright (C) 2014, 2016 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 this program. 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 <boost/shared_ptr.hpp>

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

namespace theplu {
namespace yat {
namespace omic {

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

    explicit BamPairIterator(Base begin, Base end);

#ifndef YAT_HAVE_BOOST_ITERATOR_FACADE_PROXY_PTR

    typename BamPairIterator::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(mate_);
      YAT_ASSERT(iter_ != end_);
      dummie_.first() = *mate_;
      dummie_.second() = *iter_;
      return &dummie_;
    }
  private:
    mutable BamPair dummie_;
#endif

  private:
    Base iter_;
    Base end_;
    const BamRead* mate_;
    boost::shared_ptr<std::map<std::string, BamRead> > siam_reads_;
    typedef std::pair<int32_t, int32_t> Position;
    boost::shared_ptr<std::multimap<Position, BamRead> > reads_;
    friend class boost::iterator_core_access;

    const BamPairProxy dereference(void) const;
    bool equal(const BamPairIterator& other) const;
    void increment(void);
    void find_next(void);
  };


  template<typename Base>
  BamPairIterator<Base> bam_pair_iterator(Base base, Base end)
  { return BamPairIterator<Base>(base, end); }


  template<typename Base>
  BamPairIterator<Base> bam_pair_iterator(Base end)
  { return BamPairIterator<Base>(end, end); }


  //             template implementations

  template<typename Base>
  BamPairIterator<Base>::BamPairIterator(Base base, Base end)
    : iter_(base), end_(end),
      siam_reads_(new std::map<std::string, BamRead>),
      reads_(new std::multimap<Position, BamRead>)
  {
    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>));
    find_next();
  }


  template<typename Base>
  BamPairIterator<Base>::BamPairIterator(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>
  const BamPairProxy
  BamPairIterator<Base>::dereference(void) const
  {
    return BamPairProxy(mate_, &*iter_);
  }


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


  template<typename Base>
  void BamPairIterator<Base>::increment(void)
  {
    ++iter_;
    find_next();
  }


  template<typename Base>
  void BamPairIterator<Base>::find_next(void)
  {
    BamLessPos less;
    for (; iter_!=end_; ++iter_) {
      Position position(iter_->tid(), iter_->pos());
      Position mate_position(iter_->mtid(), iter_->mpos());
      // clear siam reads if iter is more advanced than siam reads
      if (siam_reads_->size() && less(siam_reads_->begin()->second, *iter_))
        siam_reads_->clear();

      // have not seen mate yet
      if (mate_position > position)
        reads_->insert(std::make_pair(mate_position, *iter_));
      else if (position > mate_position) {
        std::multimap<Position, BamRead>::iterator
          lower = reads_->lower_bound(position);
        // erase all reads with mate position less than current
        // position (assuming input range is sorted)
        reads_->erase(reads_->begin(), lower);

        // find mate and assign pair
        for (; lower!=reads_->end() && lower->first == position; ++lower)
          if (same_query_name(lower->second, *iter_)) {
            mate_ = &lower->second;
            return;
          }
      }
      else { // siam read, i.e., iter_ and mate have same position
        // check if we have already seen mate and stored it in map
        std::map<std::string, BamRead>::iterator
          mate = siam_reads_->lower_bound(iter_->name());
        if (mate!=siam_reads_->end() && same_query_name(mate->second, *iter_)) {
          mate_ = &mate->second;
          return;
        }
        else
          // insert with hint
          siam_reads_->insert(mate, std::make_pair(iter_->name(), *iter_));
      }
    }
  }

}}}
#endif