variant_gapless_input_stream.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_STREAM_VARIANT_GAPLESS_INPUT_STREAM_H_
#define GENESIS_POPULATION_STREAM_VARIANT_GAPLESS_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/function/functions.hpp"
#include "genesis/population/function/genome_locus.hpp"
#include "genesis/population/genome_locus.hpp"
#include "genesis/population/genome_locus_set.hpp"
#include "genesis/population/stream/variant_input_stream.hpp"
#include "genesis/population/variant.hpp"
#include "genesis/sequence/reference_genome.hpp"
#include "genesis/sequence/sequence_dict.hpp"
 
#include <cassert>
#include <functional>
#include <memory>
#include <stdexcept>
#include <string>
#include <unordered_set>
#include <utility>
#include <vector>
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Variant Gapless Input Stream
// =================================================================================================
 
class VariantGaplessInputStream
{
public:
 
    // -------------------------------------------------------------------------
    //     Typedefs and Enums
    // -------------------------------------------------------------------------
 
    using self_type         = VariantGaplessInputStream;
    using value_type        = Variant;
    using pointer           = value_type const*;
    using reference         = value_type const&;
    using difference_type   = std::ptrdiff_t;
    using iterator_category = std::input_iterator_tag;
 
    // ======================================================================================
    //      Internal Iterator
    // ======================================================================================
 
    class Iterator
    {
    public:
 
        // -------------------------------------------------------------------------
        //     Constructors and Rule of Five
        // -------------------------------------------------------------------------
 
        using self_type         = VariantGaplessInputStream::Iterator;
        using value_type        = Variant;
        using pointer           = value_type const*;
        using reference         = value_type const&;
        using iterator_category = std::input_iterator_tag;
 
    public:
 
        Iterator() = default;
 
    private:
 
        Iterator( VariantGaplessInputStream* 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 VariantGaplessInputStream;
 
        // -------------------------------------------------------------------------
        //     Accessors
        // -------------------------------------------------------------------------
 
        Variant const * operator->() const
        {
            return &current_variant_();
        }
 
        // Variant * operator->()
        // {
        //     return &current_variant_();
        // }
 
        Variant const & operator*() const
        {
            return current_variant_();
        }
 
        // Variant & operator*()
        // {
        //     return current_variant_();
        // }
 
        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_();
 
        void init_chromosome_();
 
        void advance_current_locus_();
 
        void advance_current_locus_beyond_input_();
 
        bool has_more_ref_loci_on_current_chromosome_();
 
        std::string find_next_extra_chromosome_();
 
        void prepare_current_variant_();
 
        void prepare_current_variant_ref_base_();
 
        void check_input_iterator_();
 
        bool current_locus_is_covered_by_genome_locus_set_();
 
        Variant& current_variant_()
        {
            if( current_variant_is_missing_ ) {
                return missing_variant_;
            } else {
                assert( iterator_ );
                return *iterator_;
            }
        }
 
        Variant const& current_variant_() const
        {
            if( current_variant_is_missing_ ) {
                return missing_variant_;
            } else {
                assert( iterator_ );
                return *iterator_;
            }
        }
 
    private:
 
        // Parent
        VariantGaplessInputStream* parent_ = nullptr;
 
        // Keep track of the locus that the iterator currently is at.
        GenomeLocus current_locus_;
 
        // Is the current variant missing? If so, we are using the dummy missing_variant_,
        // otherwise the one of the input iterator.
        bool current_variant_is_missing_ = false;
 
        // Storage for the missing variants of the iterators. This serves as a dummy variant
        // for all positions of the input without data, so that we do not need to re-allocate
        // every time for this. It is initialized to zero counts for all base counts.
        // In each iteration step, the position and refrence base are updated.
        // To allow manupulation or move-from, we reset the whole Variant each time though.
        // To this end, we store the number of samples here as well.
        Variant missing_variant_;
        size_t num_samples_;
 
        // 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).
        VariantInputStream::Iterator iterator_;
 
        // Cache for the ref genome and seq dict sequences, as well as the genome locus set,
        // so that we do not have to find them on every iteration here.
        ::genesis::sequence::ReferenceGenome::const_iterator ref_genome_it_;
        ::genesis::sequence::SequenceDict::const_iterator seq_dict_it_;
        ::genesis::population::GenomeLocusSet::const_iterator genome_locus_set_it_;
 
        // We keep track of which chromosomes we have seen yet, in order to avoid errors due to
        // unordered input, and also to know which chromosomes to process later for the case that
        // iterate_extra_chromosomes() is set and there are unprocessed chromosomes left.
        std::unordered_set<std::string> processed_chromosomes_;
    };
 
    // ======================================================================================
    //      Main Class
    // ======================================================================================
 
    // -------------------------------------------------------------------------
    //     Constructors and Rule of Five
    // -------------------------------------------------------------------------
 
    VariantGaplessInputStream() = default;
 
    VariantGaplessInputStream(
        VariantInputStream const& input
    )
        : input_( input )
    {}
 
    VariantGaplessInputStream(
        VariantInputStream&& input
    )
        : input_( std::move( input ))
    {}
 
    ~VariantGaplessInputStream() = default;
 
    VariantGaplessInputStream( VariantGaplessInputStream const& ) = default;
    VariantGaplessInputStream( VariantGaplessInputStream&& )      = default;
 
    VariantGaplessInputStream& operator= ( VariantGaplessInputStream const& ) = default;
    VariantGaplessInputStream& operator= ( VariantGaplessInputStream&& )      = default;
 
    friend Iterator;
 
    // -------------------------------------------------------------------------
    //     Input
    // -------------------------------------------------------------------------
 
    VariantInputStream const& input() const
    {
        return input_;
    }
 
    VariantInputStream& input()
    {
        return input_;
    }
 
    // -------------------------------------------------------------------------
    //     Iteration
    // -------------------------------------------------------------------------
 
    Iterator begin()
    {
        if( started_ ) {
            throw std::runtime_error(
                "VariantGaplessInputStream implements an input iterator (single pass); "
                "begin() can hence not be called multiple times."
            );
        }
        started_ = true;
        return Iterator( this );
    }
 
    Iterator end()
    {
        return Iterator( nullptr );
    }
 
    // ---------------------------------------------------------------------
    //     Settings
    // ---------------------------------------------------------------------
 
    bool iterate_extra_chromosomes() const
    {
        return iterate_extra_chromosomes_;
    }
 
    self_type& iterate_extra_chromosomes( bool value )
    {
        if( started_ ) {
            throw std::runtime_error(
                "VariantGaplessInputStream::iterate_extra_chromosomes() cannot be called "
                "once the iteration has been started with begin()."
            );
        }
        iterate_extra_chromosomes_ = value;
        return *this;
    }
 
    std::shared_ptr<::genesis::sequence::ReferenceGenome> reference_genome() const
    {
        return ref_genome_;
    }
 
    self_type& reference_genome( std::shared_ptr<::genesis::sequence::ReferenceGenome> value )
    {
        if( started_ ) {
            throw std::runtime_error(
                "VariantGaplessInputStream::reference_genome() cannot be called "
                "once the iteration has been started with begin()."
            );
        }
        if( value && seq_dict_ ) {
            throw std::invalid_argument(
                "Cannot set reference_genome() in VariantGaplessInputStream "
                "when sequence_dict() is already provided."
            );
        }
        ref_genome_ = value;
        return *this;
    }
 
    std::shared_ptr<genesis::sequence::SequenceDict> sequence_dict() const
    {
        return seq_dict_;
    }
 
    self_type& sequence_dict( std::shared_ptr<genesis::sequence::SequenceDict> value )
    {
        if( started_ ) {
            throw std::runtime_error(
                "VariantGaplessInputStream::sequence_dict() cannot be called "
                "once the iteration has been started with begin()."
            );
        }
        if( value && ref_genome_ ) {
            throw std::invalid_argument(
                "Cannot set sequence_dict() in VariantGaplessInputStream "
                "when reference_genome() is already provided."
            );
        }
        seq_dict_ = value;
        return *this;
    }
 
    std::shared_ptr<GenomeLocusSet> genome_locus_set() const
    {
        return genome_locus_set_;
    }
 
    self_type& genome_locus_set( std::shared_ptr<GenomeLocusSet> value )
    {
        if( started_ ) {
            throw std::runtime_error(
                "VariantGaplessInputStream::genome_locus_set() cannot be called "
                "once the iteration has been started with begin()."
            );
        }
        genome_locus_set_ = value;
        return *this;
    }
 
    // -------------------------------------------------------------------------
    //     Internal Members
    // -------------------------------------------------------------------------
 
private:
 
    VariantInputStream input_;
    bool iterate_extra_chromosomes_ = true;
    bool started_ = false;
 
    // We offer two ways of specifying chromosome lengths.
    // With ref genome, we additionally gain access to the bases.
    // Also, we here subset to regions if needed, to avoid unnecessary work later.
    std::shared_ptr<::genesis::sequence::ReferenceGenome> ref_genome_;
    std::shared_ptr<::genesis::sequence::SequenceDict> seq_dict_;
    std::shared_ptr<GenomeLocusSet> genome_locus_set_;
 
};
 
} // namespace population
} // namespace genesis
 
#endif // include guard