variant_parallel_input_stream.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_STREAM_VARIANT_PARALLEL_INPUT_STREAM_H_
#define GENESIS_POPULATION_STREAM_VARIANT_PARALLEL_INPUT_STREAM_H_
 
/*
    Genesis - A toolkit for working with phylogenetic data.
    Copyright (C) 2014-2024 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 <lucas.czech@sund.ku.dk>
    University of Copenhagen, Globe Institute, Section for GeoGenetics
    Oster Voldgade 5-7, 1350 Copenhagen K, Denmark
*/
 
#include "genesis/population/sample_counts.hpp"
#include "genesis/population/function/genome_locus.hpp"
#include "genesis/population/genome_locus.hpp"
#include "genesis/population/stream/variant_input_stream.hpp"
#include "genesis/population/variant.hpp"
#include "genesis/sequence/sequence_dict.hpp"
#include "genesis/utils/containers/optional.hpp"
 
#include <cassert>
#include <functional>
#include <memory>
#include <set>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
namespace genesis {
namespace population {
 
// =================================================================================================
//      Variant Parallel Input Stream
// =================================================================================================
 
class VariantParallelInputStream
{
public:
 
    // -------------------------------------------------------------------------
    //     Typedefs and Enums
    // -------------------------------------------------------------------------
 
    enum class ContributionType
    {
        kCarrying,
 
        kFollowing
    };
 
    struct JoinedVariantParams
    {
        // Weird clang error here if the following is not defined the way it is now.
        // See https://stackoverflow.com/a/44693603/4184258
        JoinedVariantParams()
        {}
 
        bool allow_ref_base_mismatches = false;
        bool allow_alt_base_mismatches = true;
        bool treat_non_passing_variants_as_missing = false;
        bool move_samples = false;
    };
 
    using self_type  = VariantParallelInputStream;
    using value_type = Variant;
 
    // ======================================================================================
    //      Internal Iterator
    // ======================================================================================
 
    class Iterator
    {
    public:
 
        // -------------------------------------------------------------------------
        //     Constructors and Rule of Five
        // -------------------------------------------------------------------------
 
        using self_type         = VariantParallelInputStream::Iterator;
        using value_type        = std::vector<utils::Optional<Variant>>;
        using pointer           = value_type const*;
        using reference         = value_type const&;
        using iterator_category = std::input_iterator_tag;
 
    public:
 
        Iterator() = default;
 
    private:
 
        Iterator( VariantParallelInputStream* parent );
 
    public:
 
        ~Iterator() = default;
 
        Iterator( self_type const& ) = default;
        Iterator( self_type&& )      = default;
 
        Iterator& operator= ( self_type const& ) = default;
        Iterator& operator= ( self_type&& )      = default;
 
        friend VariantParallelInputStream;
 
        // -------------------------------------------------------------------------
        //     Accessors
        // -------------------------------------------------------------------------
 
        std::vector<utils::Optional<Variant>> const * operator->() const
        {
            return &variants_;
        }
 
        std::vector<utils::Optional<Variant>> * operator->()
        {
            return &variants_;
        }
 
        std::vector<utils::Optional<Variant>> const & operator*() const
        {
            return variants_;
        }
 
        std::vector<utils::Optional<Variant>> & operator*()
        {
            return variants_;
        }
 
        std::vector<utils::Optional<Variant>> const& variants() const
        {
            return variants_;
        }
 
        std::vector<utils::Optional<Variant>>& variants()
        {
            return variants_;
        }
 
        std::vector<VariantInputStream> const& inputs() const
        {
            // We assume that the user only does this when the iterator is not an end() iterator.
            assert( parent_ );
            return parent_->inputs_;
        }
 
        VariantInputStream const& input_at( size_t index ) const
        {
            return parent_->inputs_[index];
        }
 
        utils::Optional<Variant> const& variant_at( size_t index ) const
        {
            // Return with boundary check.
            return variants_.at( index );
        }
 
        utils::Optional<Variant>& variant_at( size_t index )
        {
            // Return with boundary check.
            return variants_.at( index );
        }
 
        Variant joined_variant( JoinedVariantParams const& params = JoinedVariantParams{} );
 
        GenomeLocus const& locus() const
        {
            return current_locus_;
        }
 
        // -------------------------------------------------------------------------
        //     Iteration
        // -------------------------------------------------------------------------
 
        self_type& operator ++ ()
        {
            advance_();
            return *this;
        }
 
        self_type operator ++(int)
        {
            auto cpy = *this;
            advance_();
            return cpy;
        }
 
        operator bool() const
        {
            return parent_ != nullptr;
        }
 
        bool operator==( self_type const& it ) const
        {
            return parent_ == it.parent_;
        }
 
        bool operator!=( self_type const& it ) const
        {
            return !(*this == it);
        }
 
        // -------------------------------------------------------------------------
        //     Internal Members
        // -------------------------------------------------------------------------
 
    private:
 
        void advance_()
        {
            // Some basic checks.
            assert( parent_ );
            assert( parent_->inputs_.size() == parent_->selections_.size() );
            assert( parent_->inputs_.size() == iterators_.size() );
 
            // Depending on what type of inputs we have, we need two different algorithms
            // to find the next position to iterate to.
            if( parent_->has_carrying_input_ ) {
                advance_using_carrying_();
            } else {
                advance_using_only_following_();
            }
        }
 
        void advance_using_carrying_();
 
        void advance_using_only_following_();
 
        void increment_iterator_( VariantInputStream::Iterator& iterator );
 
        void assert_correct_chr_and_pos_( VariantInputStream::Iterator const& iterator );
 
        void update_variants_();
 
    private:
 
        // Parent
        VariantParallelInputStream* parent_ = nullptr;
 
        // Keep track of the locus that the iterator currently is at.
        // Not all sources have to be there (if they don't have data for that locus), in which case
        // we want them to be at the next position in their data beyond the current locus.
        GenomeLocus current_locus_;
 
        // Keep the iterators that we want to traverse. We only need the begin() iterators,
        // as they are themselves able to tell us if they are still good (via their operator bool).
        std::vector<VariantInputStream::Iterator> iterators_;
 
        // We need to store how many samples (SampleCounts objects) the Variant of each iterator has,
        // in order to fill in the empty ones at the iterator positions where they don't have data.
        // We cannot always look that up from the iterators themselves, as they might already have
        // reached their end of the data while others are still having data, so we store it here.
        std::vector<size_t> variant_sizes_;
        size_t variant_size_sum_;
 
        // Storage for the variants of the iterators. We need these copies, as not all iterators
        // are expected to have all loci in the genome, so if we'd instead gave access to the
        // iterators directly to the user of this class, they'd have to check if the iterator is at
        // the correct locus, and so on. So instead, we offer a user-friendly interface that they
        // can simply iterator over and check if the optional is empty or not. Bit of copying,
        // but then again, each layer of abstraction comes at some cost...
        // At least, we move (not copy) data into here, for efficiency.
        std::vector<utils::Optional<Variant>> variants_;
 
        // Store the current additional carrying locus that we are at (if those have been
        // added; if not, we just store the end iterator here).
        std::set<GenomeLocus>::const_iterator carrying_locus_it_;
 
    };
 
    // ======================================================================================
    //      Main Class
    // ======================================================================================
 
    // -------------------------------------------------------------------------
    //     Constructors and Rule of Five
    // -------------------------------------------------------------------------
 
    VariantParallelInputStream() = default;
    ~VariantParallelInputStream() = default;
 
    VariantParallelInputStream( self_type const& ) = default;
    VariantParallelInputStream( self_type&& )      = default;
 
    self_type& operator= ( self_type const& ) = default;
    self_type& operator= ( self_type&& )      = default;
 
    friend Iterator;
 
    // -------------------------------------------------------------------------
    //     Iteration
    // -------------------------------------------------------------------------
 
    Iterator begin()
    {
        if( started_ ) {
            throw std::runtime_error(
                "VariantParallelInputStream implements an input iterator (single pass); "
                "begin() can hence not be called multiple times."
            );
        }
        started_ = true;
        return Iterator( this );
    }
 
    Iterator end()
    {
        return Iterator( nullptr );
    }
 
    // -------------------------------------------------------------------------
    //     Input Sources
    // -------------------------------------------------------------------------
 
    self_type& add_variant_input_stream(
        VariantInputStream const& input,
        ContributionType selection
    ) {
        if( started_ ) {
            throw std::runtime_error(
                "VariantParallelInputStream::add_variant_input_stream() cannot be called "
                "once the iteration has been started with begin()."
            );
        }
 
        inputs_.emplace_back( input );
        selections_.emplace_back( selection );
        assert( inputs_.size() == selections_.size() );
 
        if( selection == ContributionType::kCarrying ) {
            has_carrying_input_ = true;
        }
        return *this;
    }
 
    self_type& add_variant_input(
        std::function<bool(Variant&)> input_element_generator,
        ContributionType selection
    ) {
        add_variant_input_stream( VariantInputStream( input_element_generator ), selection );
        return *this;
    }
 
    std::vector<VariantInputStream> const& inputs() const
    {
        return inputs_;
    }
 
    std::vector<VariantInputStream>& inputs()
    {
        return inputs_;
    }
 
    VariantInputStream const& input_at( size_t index ) const
    {
        return inputs_[index];
    }
 
    VariantInputStream& input_at( size_t index )
    {
        return inputs_[index];
    }
 
    size_t input_size() const
    {
        assert( inputs_.size() == selections_.size() );
        return inputs_.size();
    }
 
    // -------------------------------------------------------------------------
    //     Input Loci
    // -------------------------------------------------------------------------
 
    self_type& add_carrying_locus( GenomeLocus const& locus )
    {
        if( started_ ) {
            throw std::runtime_error(
                "VariantParallelInputStream::add_carrying_locus() cannot be called "
                "once the iteration has been started with begin()."
            );
        }
 
        // Error check.
        if( locus.chromosome.empty() || locus.position == 0 ) {
            throw std::invalid_argument(
                "Cannot add a carrying locus with empty chromosome or position 0 "
                "to VariantParallelInputStream"
            );
        }
 
        // Add to the list. Also, if loci are added with this function, these serve as carrying loci,
        // and so we can always use advance_using_carrying_() to find the next locus;
        // mark this by setting has_carrying_input_.
        carrying_loci_.insert( locus );
        has_carrying_input_ = true;
        return *this;
    }
 
    self_type& add_carrying_loci( std::vector<GenomeLocus> const& loci )
    {
        add_carrying_loci( loci.begin(), loci.end() );
        return *this;
    }
 
    template<class ForwardIterator>
    self_type& add_carrying_loci( ForwardIterator first, ForwardIterator last )
    {
        while( first != last ) {
            add_carrying_locus( *first );
            ++first;
        }
 
        // Version for if we wanted to switch the set for a vector.
        // Sort the list of loci. All this is so inefficient, as we store the chromosome names
        // again and again for each locus. The sorting is okay though, we need to have that
        // complexity somewhere - using a std::set for example would just shift the place where
        // we do the sorting, but would make iteration a bit more tricky, and would need even more
        // memory.
        // std::sort( carrying_loci_.begin(), carrying_loci_.end() );
        // carrying_loci_.erase(
        //     std::unique( carrying_loci_.begin(), carrying_loci_.end() ),
        //     carrying_loci_.end()
        // );
 
        return *this;
    }
 
    // -------------------------------------------------------------------------
    //     Sequence Dict
    // -------------------------------------------------------------------------
 
    std::shared_ptr<genesis::sequence::SequenceDict> sequence_dict() const
    {
        return sequence_dict_;
    }
 
    self_type& sequence_dict( std::shared_ptr<genesis::sequence::SequenceDict> value )
    {
        if( started_ ) {
            throw std::runtime_error(
                "VariantParallelInputStream::sequence_dict() cannot be called "
                "once the iteration has been started with begin()."
            );
        }
        sequence_dict_ = value;
        return *this;
    }
 
    // -------------------------------------------------------------------------
    //     Data Members
    // -------------------------------------------------------------------------
 
private:
 
    // Store all input sources, as well as the type (carrying or following) of how we want
    // to traverse them. We keep track whether at least one of them is of type carrying.
    // If not (all following), the advance function of the iterator needs to be special.
    std::vector<VariantInputStream> inputs_;
    std::vector<ContributionType> selections_;
    bool has_carrying_input_ = false;
    bool started_ = false;
 
    // Store all additional loci that we want to include as stops in the iterator.
    // Memory-wise, this is highly inefficient, as we store the chromosome name for each of them.
    // But for now, this is easiest and fastest. We use a set, so that adding loci one after another
    // always results in a sorted container, without having to re-sort every time.
    // This again has a bit of a higher memory impact, but that should be okay for now.
    std::set<GenomeLocus> carrying_loci_;
 
    // Keep a sequence dictionary for the order of chromosomes. If not provided, we assume
    // chromosomes are sorted lexicographically. If provided however, the order as specified
    // by the dictionary is expected of the input sources instead.
    std::shared_ptr<genesis::sequence::SequenceDict> sequence_dict_;
 
};
 
} // namespace population
} // namespace genesis
 
#endif // include guard