microvariant_scanner.hpp

Go to the documentation of this file.

#ifndef GENESIS_SEQUENCE_KMER_MICROVARIANT_SCANNER_H_
#define GENESIS_SEQUENCE_KMER_MICROVARIANT_SCANNER_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/sequence/kmer/kmer.hpp"
#include "genesis/utils/containers/range.hpp"
 
#include <array>
#include <cassert>
#include <climits>
#include <cstdint>
#include <limits>
#include <stdexcept>
#include <string>
 
namespace genesis {
namespace sequence {
 
// =================================================================================================
//     Kmer Microvariant Scanner
// =================================================================================================
 
template<typename Tag>
class MicrovariantScanner
{
public:
 
    // -----------------------------------------------------
    //     Typedefs
    // -----------------------------------------------------
 
    using iterator_category = std::forward_iterator_tag;
    using self_type         = MicrovariantScanner;
    using value_type        = Kmer<Tag>;
 
    // -----------------------------------------------------
    //     Constructors and Rule of Five
    // -----------------------------------------------------
 
    MicrovariantScanner() // = default;
        : kmer_( 0 )
    {}
 
    explicit MicrovariantScanner( Kmer<Tag> const& kmer )
        : kmer_( kmer )
        , pos_( 0 )
        , cnt_( 0 )
    {
        // Iterate to the first microvariant.
        operator++();
    }
 
    ~MicrovariantScanner() = default;
 
    MicrovariantScanner( MicrovariantScanner const& ) = default;
    MicrovariantScanner( MicrovariantScanner&& )      = default;
 
    MicrovariantScanner& operator= ( MicrovariantScanner const& ) = default;
    MicrovariantScanner& operator= ( MicrovariantScanner&& )      = default;
 
    // -----------------------------------------------------
    //     Operators
    // -----------------------------------------------------
 
    value_type const& operator * ()
    {
        return kmer_;
    }
 
    value_type const* operator -> ()
    {
        return &kmer_;
    }
 
    self_type& operator ++ ()
    {
        // Check assumptions of this function.
        using WordType = typename Kmer<Tag>::WordType;
        static_assert( Kmer<Tag>::BITS_PER_CHAR == 2, "Kmer<Tag>::BITS_PER_CHAR != 2" );
        static_assert( std::is_same<WordType, std::uint64_t>::value, "Kmer::WordType != uint64_t" );
 
        // We use four xor's at the current position to cycle through the variants:
        // The first thee are the substitutions, the last one then restores the original value.
        // For this, we use the xor order 01 11 01 11.
        //
        // The table shows that this works for all four possible values.
        //
        //            | 00 01 10 11
        //     ---------------------
        //     0 | 01 | 01 00 11 10
        //     1 | 11 | 10 11 00 01
        //     2 | 01 | 11 10 01 00
        //     3 | 11 | 00 01 10 11
 
        // Helper function that cycles the value at a position, using the above table.
        auto cycle_ = [&]( size_t pos, size_t& cnt )
        {
            // Move the needed xor value to the position in the word and apply it.
            WordType const xor_val = ( cnt % 2 == 0 ? 0x1 : 0x3 );
            kmer_.value() ^= ( xor_val << ( 2 * pos ));
            ++cnt;
        };
 
        // Do at least one cycle at the current position.
        cycle_( pos_, cnt_ );
 
        // Check if we already did all three possible substitutions at the current position.
        if( cnt_ == 4 ) {
 
            // If this is not the last possible position, move to the next one.
            // Otherwise, we are done, indicated by setting everything to the default.
            assert( kmer_.k() > 0 );
            if( pos_ < kmer_.k() - 1 ) {
                ++pos_;
                cnt_ = 0;
                cycle_( pos_, cnt_ );
            } else {
                pos_ = std::numeric_limits<size_t>::max();
                cnt_ = std::numeric_limits<size_t>::max();
            }
        }
 
        return *this;
    }
 
    // self_type operator ++ (int)
    // {
    //     self_type tmp = *this;
    //     ++(*this);
    //     return tmp;
    // }
 
    bool operator == ( self_type const& other ) const
    {
        return ( pos_ == other.pos_ ) && ( cnt_ == other.cnt_ );
    }
 
    bool operator != ( self_type const& other ) const
    {
        return !( other == *this );
    }
 
    // -----------------------------------------------------
    //     Members
    // -----------------------------------------------------
 
    size_t position() const
    {
        return pos_;
    }
 
    Kmer<Tag> const& kmer() const
    {
        return kmer_;
    }
 
private:
 
    // The current k-mer, which always has an additional value (compared to the original k-mer).
    Kmer<Tag> kmer_;
 
    // The position where currently a value is inserted, and the counter for the possible
    // microvariant (substitution) possibilities per position.
    size_t pos_ = std::numeric_limits<size_t>::max();
    size_t cnt_ = std::numeric_limits<size_t>::max();
 
};
 
// =================================================================================================
//     Range Wrapper
// =================================================================================================
 
template<typename Tag>
inline utils::Range<MicrovariantScanner<Tag>> iterate_microvariants( Kmer<Tag> const& kmer )
{
    return {
        MicrovariantScanner<Tag>( kmer ),
        MicrovariantScanner<Tag>()
    };
}
 
} // namespace sequence
} // namespace genesis
 
#endif // include guard