sliding_window_generator.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_WINDOW_SLIDING_WINDOW_GENERATOR_H_
#define GENESIS_POPULATION_WINDOW_SLIDING_WINDOW_GENERATOR_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 <cassert>
#include <deque>
#include <functional>
#include <stdexcept>
#include <string>
#include <vector>
 
#include "genesis/population/window/functions.hpp"
#include "genesis/population/window/window.hpp"
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Sliding Window Type
// =================================================================================================
 
enum class SlidingWindowType
{
    kInterval,
 
    kVariants,
 
    kChromosome
};
 
// =================================================================================================
//     Genomic Sliding Window Generator
// =================================================================================================
 
template<class D, class A = EmptyAccumulator>
class SlidingWindowGenerator
{
public:
 
    // -------------------------------------------------------------------------
    //     Typedefs and Enums
    // -------------------------------------------------------------------------
 
    using Data = D;
    using Accumulator = A;
    using Window = ::genesis::population::Window<D, A>;
    using Entry = typename Window::Entry;
 
    using self_type = SlidingWindowGenerator<Data, Accumulator>;
 
    using on_chromosome_start = std::function<void(
        std::string const& chromosome,
        typename Window::Accumulator& accumulator
    )>;
 
    using on_chromosome_finish = std::function<void(
        std::string const& chromosome,
        typename Window::Accumulator& accumulator
    )>;
 
    using on_enqueue = std::function<void(
        typename Window::Entry const& entry,
        typename Window::Accumulator& accumulator
    )>;
 
    using on_dequeue = std::function<void(
        typename Window::Entry const& entry,
        typename Window::Accumulator& accumulator
    )>;
 
    using on_emission = std::function<void(
        Window const& window
    )>;
 
    // -------------------------------------------------------------------------
    //     Constructors and Rule of Five
    // -------------------------------------------------------------------------
 
    // /**
    //  * @brief Create a default (empty) instance.
    //  *
    //  * Not really important for usage, as this cannot do much. Typically, one nees a window width
    //  * and stride for doing any work. But having a default constructor comes in handy at times.
    //  */
    // SlidingWindowGenerator() = default;
 
    SlidingWindowGenerator( SlidingWindowType type, size_t width, size_t stride = 0 )
        : window_type_(type)
        , width_(width)
        , stride_(stride)
    {
        if( width == 0 ) {
            throw std::runtime_error( "Cannot use SlidingWindowGenerator of width 0." );
        }
        if( stride == 0 ) {
            stride_ = width;
        }
        if( stride_ > width_ ) {
            throw std::runtime_error( "Cannot use SlidingWindowGenerator with stride > width." );
        }
    }
 
    ~SlidingWindowGenerator()
    {
        finish_chromosome();
    }
 
    SlidingWindowGenerator( SlidingWindowGenerator const& ) = default;
    SlidingWindowGenerator( SlidingWindowGenerator&& )      = default;
 
    SlidingWindowGenerator& operator= ( SlidingWindowGenerator const& ) = default;
    SlidingWindowGenerator& operator= ( SlidingWindowGenerator&& )      = default;
 
    // -------------------------------------------------------------------------
    //     Settings
    // -------------------------------------------------------------------------
 
    SlidingWindowType window_type() const
    {
        return window_type_;
    }
 
    size_t width() const
    {
        return width_;
    }
 
    size_t stride() const
    {
        return stride_;
    }
 
    bool emit_incomplete_windows() const
    {
        return emit_incomplete_windows_;
    }
 
    void emit_incomplete_windows( bool value )
    {
        emit_incomplete_windows_ = value;
    }
 
    // /* *
    //  * @brief Get the WindowAnchorType that we use for the emitted Window%s.
    //  */
    // WindowAnchorType anchor_type() const
    // {
    //     return window_.anchor_type();
    // }
    //
    // /* *
    //  * @brief Set the WindowAnchorType that we use for the emitted Window%s.
    //  */
    // void anchor_type( WindowAnchorType value )
    // {
    //     window_.anchor_type( value );
    // }
 
    // -------------------------------------------------------------------------
    //     Accessors & Modifiers
    // -------------------------------------------------------------------------
 
    std::string const& chromosome() const
    {
        // We could keep our own chromosome here, but Window already as a member for this,
        // so we just re-use.
        return window_.chromosome();
    }
 
    bool empty() const
    {
        return next_index_ == 0;
    }
 
    void clear()
    {
        current_start_ = 0;
        next_index_ = 0;
        window_.clear();
    }
 
    // -------------------------------------------------------------------------
    //     Event Plugin Functionals
    // -------------------------------------------------------------------------
 
    self_type& add_chromosome_start_plugin( on_chromosome_start const& plugin )
    {
        chromosome_start_plugins_.push_back(plugin);
        return *this;
    }
 
    self_type& add_chromosome_finish_plugin( on_chromosome_finish const& plugin )
    {
        chromosome_finish_plugins_.push_back(plugin);
        return *this;
    }
 
    self_type& add_enqueue_plugin( on_enqueue const& plugin )
    {
        enqueue_plugins_.push_back(plugin);
        return *this;
    }
 
    self_type& add_dequeue_plugin( on_dequeue const& plugin )
    {
        dequeue_plugins_.push_back(plugin);
        return *this;
    }
 
    self_type& add_emission_plugin( on_emission const& plugin )
    {
        emission_plugins_.push_back(plugin);
        return *this;
    }
 
    void clear_plugins()
    {
        chromosome_start_plugins_.clear();
        chromosome_finish_plugins_.clear();
        enqueue_plugins_.clear();
        dequeue_plugins_.clear();
        emission_plugins_.clear();
    }
 
    // -------------------------------------------------------------------------
    //     Enqueue and Generate Windows
    // -------------------------------------------------------------------------
 
    void start_chromosome( std::string const& chromosome )
    {
        if( chromosome != window_.chromosome() ) {
            finish_chromosome();
            window_.chromosome( chromosome );
        }
    }
 
    void enqueue( std::string const& chromosome, size_t position, Data const& data )
    {
        start_chromosome( chromosome );
        enqueue_( position, Data{ data });
    }
 
    void enqueue( std::string const& chromosome, size_t position, Data&& data )
    {
        start_chromosome( chromosome );
        enqueue_( position, std::move( data ));
    }
 
    void enqueue( size_t position, Data const& data )
    {
        enqueue_( position, Data{ data });
    }
 
    void enqueue( size_t position, Data&& data )
    {
        enqueue_( position, std::move( data ));
    }
 
    void finish_chromosome( size_t last_position = 0 )
    {
        // If nothing was enqueued yet, there is nothing to finish.
        // This also makes sure that calling this function multiple times in a row does not
        // have any side effects.
        if( next_index_ == 0 ) {
            return;
        }
 
        // If we did not get a specific last position, we just finish the current interval.
        if( last_position == 0 ) {
            last_position = current_start_ + width_;
        }
 
        // Boundary check. We make sure that the given position is neither in front of the current
        // window, or, if there are entries in the list, also not in front of those.
        // See enqueue_() for details. Same here.
        size_t cur_end = window_.entries().empty() ? 0 : window_.entries().back().position;
        cur_end = current_start_ > cur_end ? current_start_ - 1 : cur_end;
        if( last_position <= cur_end ) {
            throw std::runtime_error(
                "Cannot call finish_chromosome() with position " + std::to_string(last_position) +
                ", as the current window/chromosome is already advanced up to position " +
                std::to_string( cur_end ) + "."
            );
        }
 
        // Emit the remaining data entries.
        if( window_type_ == SlidingWindowType::kInterval ) {
            synchronize_interval_( last_position );
            assert( current_start_ <= last_position );
            assert( last_position < current_start_ + width_ );
 
            // Special case for the emit_incomplete_windows_ setting. We have synchronized so that
            // the given last_position is within the current interval. Now we need to emit that
            // particular (incomplete) window and clean it up. The last p
            if( emit_incomplete_windows_ ) {
                emit_window_( current_start_, last_position + 1, last_position + 1 );
            }
 
        } else if( window_type_ == SlidingWindowType::kVariants ) {
            throw std::runtime_error( "Not yet implemented." );
        } else {
            assert( false );
        }
 
        // Wrap up the chromosome, and clear, so that we can start a new chromosome cleanly.
        for( auto const& chromosome_finish_plugin : chromosome_finish_plugins_ ) {
            if( chromosome_finish_plugin ) {
                chromosome_finish_plugin( window_.chromosome(), window_.accumulator() );
            }
        }
        clear();
    }
 
    // -------------------------------------------------------------------------
    //     General Internal Members
    // -------------------------------------------------------------------------
 
private:
 
    void enqueue_( size_t position, Data&& data )
    {
        // If this is the first enqueuing of the window or the chromosome,
        // we need to call the start plugin.
        if( next_index_ == 0 ) {
            for( auto const& chromosome_start_plugin : chromosome_start_plugins_ ) {
                if( chromosome_start_plugin ) {
                    chromosome_start_plugin( window_.chromosome(), window_.accumulator() );
                }
            }
        }
 
        // Boundary check. We make sure that the given position is neither in front of the current
        // window, or, if there are entries in the list, also not in front of those.
        // (There might be cases were we are already in the middle of the chromosome, but the
        // entries list is empty. Not entirely sure when this can occurr, but it feels like it can,
        // and just checking this doesn't cost us much. If anyone wants to think this through,
        // feel free.)
        size_t cur_end = window_.entries().empty() ? 0 : window_.entries().back().position;
        cur_end = current_start_ > cur_end ? current_start_ - 1 : cur_end;
        if( position <= cur_end ) {
            throw std::invalid_argument(
                "Cannot enqueue at position " + std::to_string( position ) +
                ", as the current window/chromosome is already advanced up to position " +
                std::to_string( cur_end ) +
                ". Either start a new window or a new chromosome within the window. "
                "Note that this error might be caused by a VCF file that is not sorted by "
                "chromosome and position."
            );
        }
        assert( position >= current_start_ );
        assert( window_.entries().empty() || position > window_.entries().back().position );
 
        // Do the correct type of enqueuing.
        if( window_type_ == SlidingWindowType::kInterval ) {
            enqueue_interval_( position, std::move( data ));
        } else if( window_type_ == SlidingWindowType::kVariants ) {
            throw std::runtime_error( "Not yet implemented." );
        } else {
            assert( false );
        }
    }
 
    // -------------------------------------------------------------------------
    //     Interval Internal Members
    // -------------------------------------------------------------------------
 
private:
 
    void enqueue_interval_( size_t position, Data&& data )
    {
        assert( window_type_ == SlidingWindowType::kInterval );
 
        // Make sure that we move to the interval where our position needs to be added to.
        synchronize_interval_( position );
        assert( current_start_ <= position );
        assert( position < current_start_ + width_ );
 
        // Add the new data to our entry queue.
        window_.entries().emplace_back( next_index_, position, std::move( data ));
        ++next_index_;
 
        // Run the enqueue event plugins. We do not emit anything here. That will be done once
        // the interval is finished, that is, above, when a new position outside of the interval
        // is added (or we finish the whole iteration).
        for( auto const& enqueue_plugin : enqueue_plugins_ ) {
            if( enqueue_plugin ) {
                assert( window_.entries().size() > 0 );
                enqueue_plugin( window_.entries().back(), window_.accumulator() );
            }
        }
 
        // Make sure that all entries in the queue are within our current bounds,
        // and are in the correct order. We use a lambda that we immediately call.
        assert( [&](){
            size_t cur_pos = 0;
            for( auto it : window_.entries() ) {
                if( it.position < current_start_ || it.position >= current_start_ + width_ ) {
                    return false;
                }
                if( it.position < cur_pos ) {
                    return false;
                }
                cur_pos = it.position;
            }
            return true;
        }() );
    }
 
    void synchronize_interval_( size_t position )
    {
        assert( window_type_ == SlidingWindowType::kInterval );
 
        // This function is only called internally, and only if we are sure that the position is
        // valid. Let's assert this.
        assert( position >= current_start_ );
        assert( window_.entries().empty() || window_.entries().back().position < position );
 
        // Either there are no entries, or they are all within the current interval.
        // That has to be the case, because we emit if we finish an interval, and remove the data.
        // So, there should never be data that is from an old interval at this point here.
        assert(
            window_.entries().empty() || window_.entries().front().position >= current_start_
        );
        assert(
            window_.entries().empty() || window_.entries().front().position < current_start_ + width_
        );
        assert(
            window_.entries().empty() || window_.entries().back().position >= current_start_
        );
        assert(
            window_.entries().empty() || window_.entries().back().position < current_start_ + width_
        );
 
        // Emit the windows up to the position where we want to enqueue the new data entry.
        // As we slide over intervals of fixed size along the genome, this can mean that we
        // have to emit multiple (potentially empty) windows along the way, until we are at the
        // interval that contains our new position.
        while( current_start_ + width_ <= position ) {
 
            // Emit and move to next interval.
            emit_window_( current_start_, current_start_ + width_, current_start_ + stride_ );
            current_start_ += stride_;
        }
 
        // We are now within the exact interval where we need to be.
        assert( current_start_ <= position );
        assert( position < current_start_ + width_ );
    }
 
    void emit_window_( size_t first_position, size_t last_position, size_t dequeue_until )
    {
        assert( window_type_ == SlidingWindowType::kInterval );
 
        // Make sure that all entries in the queue are within our current bounds,
        // and are in the correct order. We use a lambda that we immediately call.
        assert( [&](){
            size_t cur_pos = 0;
            for( auto const& it : window_.entries() ) {
                if( it.position < first_position || it.position >= last_position ) {
                    return false;
                }
                if( it.position < cur_pos ) {
                    return false;
                }
                cur_pos = it.position;
            }
            return true;
        }() );
 
        // Prepare the window properties.
        assert( last_position > first_position );
        window_.first_position( first_position );
        window_.last_position( last_position );
 
        // Now emit all plugin functions.
        for( auto const& emission_plugin : emission_plugins_ ) {
            if( emission_plugin ) {
                emission_plugin( window_ );
            }
        }
 
        // Dequeue everything that just moved out of the current interval.
        while(
            window_.entries().size() > 0 && window_.entries().front().position < dequeue_until
        ) {
            for( auto const& dequeue_plugin : dequeue_plugins_ ) {
                if( dequeue_plugin ) {
                    dequeue_plugin( window_.entries().front(), window_.accumulator() );
                }
            }
            window_.entries().pop_front();
        }
    }
 
    // -------------------------------------------------------------------------
    //     Data Members
    // -------------------------------------------------------------------------
 
private:
 
    // Fixed settings
    SlidingWindowType window_type_;
    size_t width_;
    size_t stride_;
 
    // Other settings
    bool emit_incomplete_windows_ = false;
 
    // Current window and its position
    size_t current_start_ = 0;
    size_t next_index_ = 0;
    Window window_;
 
    // Plugin functions
    std::vector<on_chromosome_start>  chromosome_start_plugins_;
    std::vector<on_chromosome_finish> chromosome_finish_plugins_;
    std::vector<on_enqueue>           enqueue_plugins_;
    std::vector<on_dequeue>           dequeue_plugins_;
    std::vector<on_emission>          emission_plugins_;
 
};
 
} // namespace population
} // namespace genesis
 
#endif // include guard