region_window_iterator.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_WINDOW_REGION_WINDOW_ITERATOR_H_
#define GENESIS_POPULATION_WINDOW_REGION_WINDOW_ITERATOR_H_
 
/*
    Genesis - A toolkit for working with phylogenetic data.
    Copyright (C) 2014-2022 Lucas Czech
 
    This program 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.
 
    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/>.
 
    Contact:
    Lucas Czech <lczech@carnegiescience.edu>
    Department of Plant Biology, Carnegie Institution For Science
    260 Panama Street, Stanford, CA 94305, USA
*/
 
#include "genesis/population/genome_locus.hpp"
#include "genesis/population/genome_region_list.hpp"
#include "genesis/population/genome_region.hpp"
#include "genesis/population/window/window.hpp"
#include "genesis/utils/containers/range.hpp"
 
#include <cassert>
#include <deque>
#include <functional>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Region Window Iterator
// =================================================================================================
 
template<class ForwardIterator, class InputType, class DataType = InputType>
class RegionWindowIterator
{
public:
 
    // -------------------------------------------------------------------------
    //     Typedefs and Enums
    // -------------------------------------------------------------------------
 
    // using Data              = typename ForwardIterator::value_type;
    // using Data              = typename std::result_of<decltype(EntryInputFunctor)()>::type;
 
    using Window            = ::genesis::population::Window<DataType>;
    using Entry             = typename Window::Entry;
 
    using settings_type     = RegionWindowIteratorSettings<InputType, DataType>;
    using self_type         = RegionWindowIterator<ForwardIterator, InputType, DataType>;
    using value_type        = Window;
    using reference         = value_type&;
    using const_reference   = value_type const&;
 
    // using difference_type   = typename Window::container::difference_type;
    // using size_type         = typename Window::container::size_type;
    using iterator_category = std::input_iterator_tag;
 
    // ======================================================================================
    //      Internal Iterator
    // ======================================================================================
 
    class Iterator
    {
    public:
 
        // -------------------------------------------------------------------------
        //     Constructors and Rule of Five
        // -------------------------------------------------------------------------
 
        using self_type         = RegionWindowIterator<T, D>::Iterator;
        using value_type        = T;
        using pointer           = value_type const*;
        using reference         = value_type const&;
        using difference_type   = std::ptrdiff_t;
        using iterator_category = std::input_iterator_tag;
 
        using Data              = D;
 
    private:
 
        Iterator() = default;
 
        Iterator(
            RegionWindowIterator const* parent,
            ForwardIterator begin, ForwardIterator end
        )
            : parent_(parent)
            , cur_locus_(begin)
            , end_locus_(end)
        {
            // Error checking.
            if( ! parent_->entry_input_function ) {
                throw std::runtime_error(
                    "Need to set RegionWindowIteratorSettings::entry_input_function before using it "
                    "to construct a RegionWindowIterator"
                );
            }
            if( ! parent_->chromosome_function ) {
                throw std::runtime_error(
                    "Need to set RegionWindowIteratorSettings::chromosome_function before using it "
                    "to construct a RegionWindowIterator"
                );
            }
            if( ! parent_->position_function ) {
                throw std::runtime_error(
                    "Need to set RegionWindowIteratorSettings::position_function before using it "
                    "to construct a RegionWindowIterator"
                );
            }
 
            // Let's get going.
            init_();
            increment_();
        }
 
    public:
 
        ~Iterator() = default;
 
        Iterator( self_type const& ) = default;
        Iterator( self_type&& )      = default;
 
        Iterator& operator= ( self_type const& ) = default;
        Iterator& operator= ( self_type&& )      = default;
 
        friend RegionWindowIterator;
 
        // -------------------------------------------------------------------------
        //     Accessors
        // -------------------------------------------------------------------------
 
        value_type const * operator->() const
        {
            return &(window_.front());
        }
 
        value_type * operator->()
        {
            return &(window_.front());
        }
 
        value_type const & operator*() const
        {
            return window_.front();
        }
 
        value_type & operator*()
        {
            return window_.front();
        }
 
        // -------------------------------------------------------------------------
        //     Iteration
        // -------------------------------------------------------------------------
 
        self_type& operator ++()
        {
            increment_();
            return *this;
        }
 
        // self_type operator ++(int)
        // {
        //     auto cpy = *this;
        //     increment_();
        //     return cpy;
        // }
 
        operator bool() const
        {
            return ! windows_.empty();
        }
 
        bool operator==( self_type const& other ) const
        {
            todo
            return cur_locus_ == it.cur_locus_;
        }
 
        bool operator!=( self_type const& other ) const
        {
            return !(*this == other);
        }
 
        // -------------------------------------------------------------------------
        //     Internal Members
        // -------------------------------------------------------------------------
 
    private:
 
        void init_()
        {
            // Saveguard. This might be called on an empty range, in which case we just do nothing.
            if( cur_locus_ == end_locus_ ) {
                return;
            }
 
            // Clear the window and prepare for new chromosome.
            window_.clear();
            window_.chromosome( settings_.chromosome_function( *cur_locus_ ));
            is_first_window_ = true;
            is_last_window_ = false;
            next_index_ = 0;
 
            if( settings_.emit_leading_empty_windows ) {
                cur_locus_start_ = 1;
            } else {
                // Set the start to the window position that we would get after going through all
                // the previous windows if they were emitted.
                auto const pos = settings_.position_function( *cur_locus_ );
                cur_locus_start_ = pos - (( pos - 1 ) % settings_.stride );
            }
        }
 
        void increment_()
        {
            assert( parent_ );
 
            auto const cur_chr = windows_.front().chromosome();
 
            // We go to the next window, meaning we are done with this one.
            windows_.pop_front();
 
            // If there are no windows left, that means we have processed all regions of
            // the current chromosome. Hence, move the data input forward to the next chromosome.
            if( windows_.empty() ) {
                do {
 
                    // Move the input data to the next chromosome.
                    nope not if we are already at the next chr...
                    auto const cur_chr = parent_->chromosome_function( *cur_locus_ );
                    while (
                        cur_locus_ != end_locus_ &&
                        cur_chr == parent_->chromosome_function( *cur_locus_ )
                    ) {
                        ++cur_locus_;
                    }
 
                    // We might have reached the end. If so, we are done with the whole thing.
                    if( cur_locus_ == end_locus_ ) {
                        parent_ = nullptr;
                        return;
                    }
 
                    // If we are here, our input has still data. Let's see if the chromosome that
                    // we just moved to also is present in our regions list. If so, we are done
                    // with this part - we are at a new chromosome. Let's store the iterators
                    // to its regions then.
                    // However, if the chromosome is not present in the region, we continue with
                    // another iteration of the while loop, to find the next chromosome of the
                    // input data.
                    auto const new_chr = parent_->chromosome_function( *cur_locus_ );
                    if( parent_->region_list_->has_chromosome( new_chr )) {
 
                        // Get the interval tree iterators for the new chromosome.
                        // If there are regions in there, we use them!
                        auto const& tree = parent_->region_list_->get_chromosome_regions( new_chr );
                        if( tree.begin() != tree.end() ) {
                            cur_region_ = tree.begin();
                            end_region_ = tree.end();
                            cur_index_ = 0;
                            break;
                        }
                    }
                } while ( true );
            }
 
            while( cur_locus_ != end_locus_ ) {
                // We need those often, so let's store them here.
                auto const cur_chr = parent_->chromosome_function( *cur_locus_ );
                auto const cur_pos = parent_->position_function( *cur_locus_ )
 
                // Add all regions as windows that are in the region list and span the current locus.
                while( cur_region_ != end_region_ && cur_region_->interval().within( cur_pos )) {
 
                    // Add the region to our current window list.
                    windows_.emplace_back();
                    auto& window = windows_.back();
                    window.chromosome( cur_chr );
                    window.first_position( cur_region_.->interval().low() );
                    window.last_position( cur_region_->interval().high() );
 
                    // Move to the next region, and assert that (if we have not reached the end of
                    // the region list already) its low is actually not before the previous region,
                    // i.e., that the regions are sorted by their low ends (instead of high ends).
                    // This should be a property of the IntervalTree already, but let's make sure.
                    ++cur_region_;
                    assert(
                        cur_region_ == end_region_ ||
                        cur_region_->interval().low() >= window.first_position()
                    );
                }
 
                // Add the current locus to all windows that span it.
                for( auto& window : windows_ ) {
 
                    // We only keep windows in the list that are of the current chromosome,
                    // and are already started at the current position.
                    assert( cur_chr == window.chromosome() );
                    assert( cur_pos >= window.first_position() );
 
                    // If the window spans the current locus, add it.
                    if( cur_pos <= window.last_position() ) {
                        window.entries().emplace_back(
                            cur_index_,
                            cur_pos,
                            parent_->entry_input_function( *cur_locus_ )
                        );
                    }
                }
 
                // Move to next input.
                ++cur_index_;
                ++cur_locus_;
 
 
                if( parent_->chromosome_function( *cur_locus_ ) != window_.chromosome() ) {
 
                }
            }
        }
 
    private:
 
        // Parent
        RegionWindowIterator<> const* parent_ = nullptr;
 
        // Underlying data iterator
        ForwardIterator cur_locus_;
        ForwardIterator end_locus_;
        size_t cur_index_ = 0;
 
        // Region list
        GenomeRegionList::iterator cur_region_;
        GenomeRegionList::iterator end_region_;
 
        // Store the windows that we are filling
        std::deque<Window> windows_;
 
    };
 
    // ======================================================================================
    //      Main Class
    // ======================================================================================
 
    // -------------------------------------------------------------------------
    //     Constructors and Rule of Five
    // -------------------------------------------------------------------------
 
    RegionWindowIterator(
        ForwardIterator begin, ForwardIterator end
    )
        : cur_locus_(begin)
        , end_locus_(end)
    {
 
    }
 
    ~RegionWindowIterator() = default;
 
    RegionWindowIterator( RegionWindowIterator const& ) = default;
    RegionWindowIterator( RegionWindowIterator&& )      = default;
 
    RegionWindowIterator& operator= ( RegionWindowIterator const& ) = default;
    RegionWindowIterator& operator= ( RegionWindowIterator&& )      = default;
 
    // -------------------------------------------------------------------------
    //     Functors
    // -------------------------------------------------------------------------
 
    std::function<DataType( InputType const& )> entry_input_function;
 
    std::function<std::string( InputType const& )> chromosome_function;
 
    std::function<size_t( InputType const& )> position_function;
 
    // -------------------------------------------------------------------------
    //     Data Members
    // -------------------------------------------------------------------------
 
private:
 
    // Settings
    bool skip_empty_regions = true;
 
    GenomeRegionList const* region_list_;
 
};
 
} // namespace population
} // namespace genesis
 
#endif // include guard