genome_locus_set.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_GENOME_LOCUS_SET_H_
#define GENESIS_POPULATION_GENOME_LOCUS_SET_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 <cassert>
#include <cstdint>
#include <limits>
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <vector>
 
#include "genesis/population/genome_locus.hpp"
#include "genesis/population/genome_region.hpp"
#include "genesis/population/genome_region_list.hpp"
#include "genesis/sequence/sequence_dict.hpp"
#include "genesis/utils/math/bitvector.hpp"
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Genome Locus Set
// =================================================================================================
 
class GenomeLocusSet
{
public:
 
    // -------------------------------------------------------------------------
    //     Typedefs and Enums
    // -------------------------------------------------------------------------
 
    using key_type    = std::string;
    using mapped_type = utils::Bitvector;
    using value_type  = std::pair<std::string const, utils::Bitvector>;
 
    using reference       = value_type&;
    using const_reference = value_type const&;
 
    using iterator       = typename std::unordered_map<std::string, utils::Bitvector>::iterator;
    using const_iterator = typename std::unordered_map<std::string, utils::Bitvector>::const_iterator;
 
    static const size_t npos = std::numeric_limits<size_t>::max();
    static_assert(
        npos == utils::Bitvector::npos,
        "Differing definitions of GenomeLocusSet::npos and Bitvector::npos"
    );
 
    // -------------------------------------------------------------------------
    //     Constructors and Rule of Five
    // -------------------------------------------------------------------------
 
    GenomeLocusSet() = default;
    ~GenomeLocusSet() = default;
 
    GenomeLocusSet( GenomeLocusSet const& ) = default;
    GenomeLocusSet( GenomeLocusSet&& )      = default;
 
    GenomeLocusSet& operator= ( GenomeLocusSet const& ) = default;
    GenomeLocusSet& operator= ( GenomeLocusSet&& )      = default;
 
    // -------------------------------------------------------------------------
    //     Modifiers
    // -------------------------------------------------------------------------
 
    // -------------------------------------------
    //         Add Manually
    // -------------------------------------------
 
    void add( std::string const& chromosome )
    {
        // We use the special value 0 to denote that we want the whole chromosome.
        add( chromosome, 0, 0 );
    }
 
    void add( std::string const& chromosome, size_t position )
    {
        add( chromosome, position, position );
    }
 
    void add(
        std::string const& chromosome,
        size_t start,
        size_t end
    );
 
    // -------------------------------------------
    //         Add Locus
    // -------------------------------------------
 
    void add( GenomeLocus const& locus )
    {
        add( locus.chromosome, locus.position, locus.position );
    }
 
    void add( GenomeLocus const& start, GenomeLocus const& end )
    {
        if( start.chromosome != end.chromosome ) {
            throw std::invalid_argument(
                "Cannot use two GenomeLocus instances with different chromosomes ( start == \"" +
                start.chromosome + "\", end == \"" + end.chromosome + "\") as an entry in a "
                "GenomeLocusSet."
            );
        }
        add( start.chromosome, start.position, end.position );
    }
 
    // -------------------------------------------
    //         Add Region
    // -------------------------------------------
 
    void add( GenomeRegion const& region )
    {
        add( region.chromosome, region.start, region.end );
    }
 
    // -------------------------------------------
    //         Add Region List
    // -------------------------------------------
 
    void add( GenomeRegionList const& list )
    {
        for( auto const& chr : list.chromosome_map() ) {
            for( auto const& interval : chr.second ) {
                add( chr.first, interval.low(), interval.high() );
            }
        }
    }
 
    // -------------------------------------------
    //         Add Bitvector
    // -------------------------------------------
 
    void add( std::string const& chromosome, utils::Bitvector const& values );
 
    void add( std::string const& chromosome, utils::Bitvector&& values );
 
    // -------------------------------------------
    //         Clear
    // -------------------------------------------
 
    void clear()
    {
        locus_map_.clear();
    }
 
    void clear( std::string const& chromosome )
    {
        if( locus_map_.count( chromosome ) == 0 ) {
            throw std::invalid_argument(
                "Chromosome name \"" + chromosome + "\" not found in GenomeLocusSet"
            );
        }
        locus_map_.erase( chromosome );
    }
 
    // -------------------------------------------------------------------------
    //     Operations
    // -------------------------------------------------------------------------
 
    void set_intersect( GenomeLocusSet const& rhs );
 
    void set_union( GenomeLocusSet const& rhs );
 
    void invert( sequence::SequenceDict const& sequence_dict );
 
    // -------------------------------------------------------------------------
    //     Locus Covered
    // -------------------------------------------------------------------------
 
    static bool is_covered( utils::Bitvector const& bitvector, size_t position )
    {
        // Boundary check.
        if( bitvector.empty() ) {
            throw std::invalid_argument(
                "GenomeLocusSet::is_covered( Bitvector const&, size_t ) called with empty Bitvector"
            );
        }
        // if( position == 0 ) {
        //     throw std::invalid_argument(
        //         "GenomeLocusSet::is_covered( Bitvector const&, size_t ) called with position==0"
        //     );
        // }
 
        // If the chromosome has the 0th bit set, the whole chromosome is covered.
        if( bitvector.get( 0 ) ) {
            return true;
        }
 
        // If the above is not the case, check the actual position.
        // If the position is outside of the bitvector, it is not covered, obviously.
        if( position >= bitvector.size() ) {
            return false;
        }
        return bitvector.get( position );
    }
 
    static bool is_covered( const_iterator const& it, size_t position )
    {
        return is_covered( it->second, position );
    }
 
    bool is_covered( std::string const& chromosome, size_t position ) const
    {
        // Using find(), so we only have to search in the map once, for speed.
        auto const it = locus_map_.find( chromosome );
        if( it == locus_map_.end() ) {
            return false;
        }
        auto const& bv = it->second;
        assert( ! bv.empty() );
        return is_covered( bv, position );
    }
 
    bool is_covered( std::string const& chromosome ) const
    {
        auto const it = locus_map_.find( chromosome );
        if( it == locus_map_.end() ) {
            return false;
        }
        auto const& bv = it->second;
        assert( ! bv.empty() );
        return bv.get( 0 );
    }
 
    // -------------------------------------------------------------------------
    //     Any Covered Locus
    // -------------------------------------------------------------------------
 
    bool any_covered( std::string const& chromosome ) const
    {
        auto const it = locus_map_.find( chromosome );
        if( it == locus_map_.end() ) {
            return false;
        }
        auto const& bv = it->second;
        assert( ! bv.empty() );
 
        // We do not need to to an extra check for position 0 here.
        // If it is true, then so will be the result.
        return bv.any_set();
    }
 
    // -------------------------------------------------------------------------
    //     Next Covered Locus
    // -------------------------------------------------------------------------
 
    static size_t next_covered( utils::Bitvector const& bitvector, size_t start_position )
    {
        // Boundary check.
        if( bitvector.empty() ) {
            throw std::invalid_argument(
                "GenomeLocusSet::next_covered() called with empty Bitvector"
            );
        }
        if( start_position == 0 ) {
            throw std::invalid_argument(
                "GenomeLocusSet::next_covered() called with start_position==0"
            );
        }
 
        // If the chromosome has the 0th bit set, the whole chromosome is covered,
        // so that the start_position we are at is also covered.
        if( bitvector.get( 0 ) ) {
            return start_position;
        }
 
        // If the above is not the case, check the actual start_position.
        // If the start_position is outside of the bitvector, it is not covered, obviously.
        return bitvector.find_next_set( start_position );
    }
 
    static size_t next_covered( const_iterator const& it, size_t start_position )
    {
        return next_covered( it->second, start_position );
    }
 
    size_t next_covered( std::string const& chromosome, size_t start_position ) const
    {
        auto const it = locus_map_.find( chromosome );
        if( it == locus_map_.end() ) {
            return GenomeLocusSet::npos;
        }
        auto const& bv = it->second;
        assert( ! bv.empty() );
        return next_covered( bv, start_position );
    }
 
    // -------------------------------------------------------------------------
    //     Chromosome Iterators
    // -------------------------------------------------------------------------
 
    const_iterator begin() const
    {
        return locus_map_.begin();
    }
 
    const_iterator end() const
    {
        return locus_map_.end();
    }
 
    // -------------------------------------------------------------------------
    //     Chromosome Accessors
    // -------------------------------------------------------------------------
 
    bool empty() const
    {
        return locus_map_.empty();
    }
 
    size_t chromosome_count() const
    {
        return locus_map_.size();
    }
 
    std::vector<std::string> chromosome_names() const
    {
        std::vector<std::string> result;
        for( auto const& p : locus_map_ ) {
            result.push_back( p.first );
        }
        return result;
    }
 
    bool has_chromosome( std::string const& chromosome ) const
    {
        return locus_map_.count( chromosome ) > 0;
    }
 
    const_iterator find( std::string const& chromosome )
    {
        return locus_map_.find( chromosome );
    }
 
    utils::Bitvector const& chromosome_positions( std::string const& chromosome ) const
    {
        return locus_map_.at( chromosome );
    }
 
    utils::Bitvector& chromosome_positions( std::string const& chromosome )
    {
        return locus_map_.at( chromosome );
    }
 
    // -------------------------------------------------------------------------
    //     Data Members
    // -------------------------------------------------------------------------
 
private:
 
    // Map from chromosome names to bitvectors representing which positions are in (true)
    // and out (false). Note that position 0 is special; if set, it means that we consider
    // the whole chromosome as covered.
    std::unordered_map<std::string, utils::Bitvector> locus_map_;
 
};
 
} // namespace population
} // namespace genesis
 
#endif // include guard