function/genome_locus.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_FUNCTION_GENOME_LOCUS_H_
#define GENESIS_POPULATION_FUNCTION_GENOME_LOCUS_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/genome_locus.hpp"
#include "genesis/sequence/sequence_dict.hpp"
#include "genesis/utils/math/common.hpp"
 
#include <cassert>
#include <cstring>
#include <iostream>
#include <memory>
#include <string>
#include <stdexcept>
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Output
// =================================================================================================
 
inline std::string to_string( GenomeLocus const& locus )
{
    // Error case.
    if( locus.chromosome.empty() ) {
        throw std::invalid_argument( "Invalid GenomeLocus with empty chromosome." );
    }
 
    // Special case.
    if( locus.position == 0 ) {
        return locus.chromosome;
    }
 
    // General case.
    return locus.chromosome + ":" + std::to_string( locus.position );
}
 
inline std::ostream& operator << ( std::ostream& os, GenomeLocus const& locus )
{
    os << to_string( locus );
    return os;
}
 
// =================================================================================================
//     Comparion Operators
// =================================================================================================
 
// Really good example of the claim that C++ needs a lot of boilerplate code...
// We here provide all comparison operators for GenomeLocus, but also comparing loci given
// as chromosome (`std::string`) and position (`size_t`) separately,
// and then also overload all of those for using references and shared pointers to a SequenceDict.
//
// Previously (before also offering overloads that take SequenceDict), we just spelled out all
// overloads manually, but that became too cumbersome with the extra overloads. So now we use a
// macro that defines all of them for us. It is not super nice for the doxygen code documentation,
// as the doxygen preprocessor does not fully seem to be able to expand the macro. Hence, we cannot
// use copydoc to refer to specific functions, and have re-written the doc blocks of the main
// functions (the ones that the macro-expanded functions call) so that they work for all overloads.
 
#define ADD_LOCUS_COMPARISON_OVERLOADS( CMP_FUNC ) \
 \
    inline int CMP_FUNC( \
        GenomeLocus const& l, \
        std::string const& r_chromosome, size_t r_position \
    ) { \
        return CMP_FUNC( l.chromosome, l.position, r_chromosome, r_position ); \
    } \
    \
 \
    inline int CMP_FUNC( \
        std::string const& l_chromosome, size_t l_position, \
        GenomeLocus const& r \
    ) { \
        return CMP_FUNC( l_chromosome, l_position, r.chromosome, r.position ); \
    } \
    \
 \
    inline int CMP_FUNC( \
        GenomeLocus const& l, \
        GenomeLocus const& r \
    ) { \
        return CMP_FUNC( l.chromosome, l.position, r.chromosome, r.position ); \
    } \
    \
 \
    inline int CMP_FUNC( \
        GenomeLocus const& l, \
        std::string const& r_chromosome, size_t r_position, \
        genesis::sequence::SequenceDict const& sequence_dict \
    ) { \
        return CMP_FUNC( l.chromosome, l.position, r_chromosome, r_position, sequence_dict ); \
    } \
    \
 \
    inline int CMP_FUNC( \
        std::string const& l_chromosome, size_t l_position, \
        GenomeLocus const& r, \
        genesis::sequence::SequenceDict const& sequence_dict \
    ) { \
        return CMP_FUNC( l_chromosome, l_position, r.chromosome, r.position, sequence_dict ); \
    } \
    \
 \
    inline int CMP_FUNC( \
        GenomeLocus const& l, \
        GenomeLocus const& r, \
        genesis::sequence::SequenceDict const& sequence_dict \
    ) { \
        return CMP_FUNC( l.chromosome, l.position, r.chromosome, r.position, sequence_dict ); \
    } \
    \
 \
    inline int CMP_FUNC( \
        std::string const& l_chromosome, size_t l_position, \
        std::string const& r_chromosome, size_t r_position, \
        std::shared_ptr<genesis::sequence::SequenceDict> const& sequence_dict \
    ) { \
        if( sequence_dict ) { \
            return CMP_FUNC( l_chromosome, l_position, r_chromosome, r_position, *sequence_dict ); \
        } else { \
            return CMP_FUNC( l_chromosome, l_position, r_chromosome, r_position ); \
        } \
    } \
    \
 \
    inline int CMP_FUNC( \
        GenomeLocus const& l, \
        std::string const& r_chromosome, size_t r_position, \
        std::shared_ptr<genesis::sequence::SequenceDict> const& sequence_dict \
    ) { \
        if( sequence_dict ) { \
            return CMP_FUNC( l.chromosome, l.position, r_chromosome, r_position, *sequence_dict ); \
        } else { \
            return CMP_FUNC( l.chromosome, l.position, r_chromosome, r_position ); \
        } \
    } \
    \
 \
    inline int CMP_FUNC( \
        std::string const& l_chromosome, size_t l_position, \
        GenomeLocus const& r, \
        std::shared_ptr<genesis::sequence::SequenceDict> const& sequence_dict \
    ) { \
        if( sequence_dict ) { \
            return CMP_FUNC( l_chromosome, l_position, r.chromosome, r.position, *sequence_dict ); \
        } else { \
            return CMP_FUNC( l_chromosome, l_position, r.chromosome, r.position ); \
        } \
    } \
    \
 \
    inline int CMP_FUNC( \
        GenomeLocus const& l, \
        GenomeLocus const& r, \
        std::shared_ptr<genesis::sequence::SequenceDict> const& sequence_dict \
    ) { \
        if( sequence_dict ) { \
            return CMP_FUNC( l.chromosome, l.position, r.chromosome, r.position, *sequence_dict ); \
        } else { \
            return CMP_FUNC( l.chromosome, l.position, r.chromosome, r.position ); \
        } \
    }
 
// -------------------------------------------------------------------------
//     Spaceship <=>
// -------------------------------------------------------------------------
 
inline int locus_compare(
    std::string const& l_chromosome, size_t l_position,
    std::string const& r_chromosome, size_t r_position
) {
    // We compare the chromosomes with strcmp, so that we only run the relatively expensive string
    // comparison once, and then check the result - if equal, do the threeway position comparison.
    auto const chr_cmp = std::strcmp( l_chromosome.c_str(), r_chromosome.c_str() );
    if( chr_cmp != 0 ) {
        return chr_cmp;
    }
    assert( l_chromosome == r_chromosome );
    return utils::compare_threeway( l_position, r_position );
}
 
inline int locus_compare(
    std::string const& l_chromosome, size_t l_position,
    std::string const& r_chromosome, size_t r_position,
    ::genesis::sequence::SequenceDict const& sequence_dict
) {
    // Here, we want to compare chromosome names order based on the given dict.
    // However, string lookup is a bit expensive, so we first do a quick check for equality,
    // and only if they are not equal, we get their indices, and compare those.
    if( l_chromosome == r_chromosome ) {
        // For identical chromosomes, we compare the positions.
        return utils::compare_threeway( l_position, r_position );
    }
 
    // Here, we know the chromosomes are different, so we compare their indices.
    // No need to compare the positions here again. We assert that they are indeed different.
    auto const l_chr_idx = sequence_dict.index_of( l_chromosome );
    auto const r_chr_idx = sequence_dict.index_of( r_chromosome );
    auto const chr_cmp = utils::compare_threeway( l_chr_idx, r_chr_idx );
    assert( chr_cmp != 0 );
    return chr_cmp;
}
 
// Add all other overloads for GenomeLocus and SequenceDict combinations.
ADD_LOCUS_COMPARISON_OVERLOADS( locus_compare )
 
// Let's be cool and add the actual spaceship, even if genesis is currently using C++11.
// Users might compile with later versions, so this might be useful to have.
 
#if __cplusplus >= 202002L
 
inline int operator <=> ( GenomeLocus const& l, GenomeLocus const& r )
{
    return locus_compare( l.chromosome, l.position, r.chromosome, r.position );
}
 
#endif // __cplusplus >= 202002L
 
// -------------------------------------------------------------------------
//     Equality ==
// -------------------------------------------------------------------------
 
inline bool locus_equal(
    std::string const& l_chromosome, size_t l_position,
    std::string const& r_chromosome, size_t r_position
) {
    return l_chromosome == r_chromosome && l_position == r_position;
}
 
inline bool locus_equal(
    GenomeLocus const& l,
    std::string const& r_chromosome, size_t r_position
) {
    return locus_equal( l.chromosome, l.position, r_chromosome, r_position );
}
 
inline bool locus_equal(
    std::string const& l_chromosome, size_t l_position,
    GenomeLocus const& r
) {
    return locus_equal( l_chromosome, l_position, r.chromosome, r.position );
}
 
inline bool locus_equal(
    GenomeLocus const& l,
    GenomeLocus const& r
) {
    return locus_equal( l.chromosome, l.position, r.chromosome, r.position );
}
 
inline bool operator == ( GenomeLocus const& l, GenomeLocus const& r )
{
    return locus_equal( l.chromosome, l.position, r.chromosome, r.position );
}
 
// -------------------------------------------------------------------------
//     Inequality !=
// -------------------------------------------------------------------------
 
inline bool locus_inequal(
    std::string const& l_chromosome, size_t l_position,
    std::string const& r_chromosome, size_t r_position
) {
    return ! locus_equal( l_chromosome, l_position, r_chromosome, r_position );
}
 
inline bool locus_inequal(
    GenomeLocus const& l,
    std::string const& r_chromosome, size_t r_position
) {
    return locus_inequal( l.chromosome, l.position, r_chromosome, r_position );
}
 
inline bool locus_inequal(
    std::string const& l_chromosome, size_t l_position,
    GenomeLocus const& r
) {
    return locus_inequal( l_chromosome, l_position, r.chromosome, r.position );
}
 
inline bool locus_inequal(
    GenomeLocus const& l,
    GenomeLocus const& r
) {
    return locus_inequal( l.chromosome, l.position, r.chromosome, r.position );
}
 
inline bool operator != ( GenomeLocus const& l, GenomeLocus const& r )
{
    return locus_inequal( l.chromosome, l.position, r.chromosome, r.position );
}
 
// -------------------------------------------------------------------------
//     Less than <
// -------------------------------------------------------------------------
 
inline bool locus_less(
    std::string const& l_chromosome, size_t l_position,
    std::string const& r_chromosome, size_t r_position
) {
    return l_chromosome < r_chromosome || ( l_chromosome == r_chromosome && l_position < r_position );
}
 
inline bool locus_less(
    std::string const& l_chromosome, size_t l_position,
    std::string const& r_chromosome, size_t r_position,
    ::genesis::sequence::SequenceDict const& sequence_dict
) {
    // Same logic as above, but using chromosome indices in the dict, instead of their names.
    // We also apply the speedup of the locus_compare(), by first checking the strings for
    // equality, before doing the expensive index lookup in the dict.
    if( l_chromosome == r_chromosome ) {
        // For identical chromosomes, we compare the positions.
        return l_position < r_position;
    }
 
    // Here, we know the chromosomes are different, so we compare their indices.
    // No need to compare the positions here again. We assert that they are indeed different.
    auto const l_chr_idx = sequence_dict.index_of( l_chromosome );
    auto const r_chr_idx = sequence_dict.index_of( r_chromosome );
    assert( l_chr_idx != r_chr_idx );
    return l_chr_idx < r_chr_idx;
}
 
// Add all other overloads for GenomeLocus and SequenceDict combinations.
ADD_LOCUS_COMPARISON_OVERLOADS( locus_less )
 
 
inline bool operator < ( GenomeLocus const& l, GenomeLocus const& r )
{
    return locus_less( l.chromosome, l.position, r.chromosome, r.position );
}
 
// -------------------------------------------------------------------------
//     Greater than >
// -------------------------------------------------------------------------
 
inline bool locus_greater(
    std::string const& l_chromosome, size_t l_position,
    std::string const& r_chromosome, size_t r_position
) {
    // Just use the existing function, but with reverse l and r.
    return locus_less( r_chromosome, r_position, l_chromosome, l_position );
}
 
inline bool locus_greater(
    std::string const& l_chromosome, size_t l_position,
    std::string const& r_chromosome, size_t r_position,
    ::genesis::sequence::SequenceDict const& sequence_dict
) {
    // Just use the existing function, but with reverse l and r.
    return locus_less( r_chromosome, r_position, l_chromosome, l_position, sequence_dict );
}
 
// Add all other overloads for GenomeLocus and SequenceDict combinations.
ADD_LOCUS_COMPARISON_OVERLOADS( locus_greater )
 
 
inline bool operator > ( GenomeLocus const& l, GenomeLocus const& r )
{
    return locus_greater( l.chromosome, l.position, r.chromosome, r.position );
}
 
// -------------------------------------------------------------------------
//     Less than or equal <=
// -------------------------------------------------------------------------
 
inline bool locus_less_or_equal(
    std::string const& l_chromosome, size_t l_position,
    std::string const& r_chromosome, size_t r_position
) {
    // We could do the simple default way of implementing this as `a == b || a < b`,
    // but this seems wasteful; in this case, we can do with fewer comparisons!
    return l_chromosome < r_chromosome || ( l_chromosome == r_chromosome && l_position <= r_position );
}
 
inline bool locus_less_or_equal(
    std::string const& l_chromosome, size_t l_position,
    std::string const& r_chromosome, size_t r_position,
    ::genesis::sequence::SequenceDict const& sequence_dict
) {
    // Same logic as in locus_less. See there for details.
    if( l_chromosome == r_chromosome ) {
        return l_position <= r_position;
    }
    auto const l_chr_idx = sequence_dict.index_of( l_chromosome );
    auto const r_chr_idx = sequence_dict.index_of( r_chromosome );
    assert( l_chr_idx != r_chr_idx );
    return l_chr_idx < r_chr_idx;
}
 
// Add all other overloads for GenomeLocus and SequenceDict combinations.
ADD_LOCUS_COMPARISON_OVERLOADS( locus_less_or_equal )
 
 
inline bool operator <= ( GenomeLocus const& l, GenomeLocus const& r )
{
    return locus_less_or_equal( l.chromosome, l.position, r.chromosome, r.position );
}
 
// -------------------------------------------------------------------------
//     Greater than or equal >=
// -------------------------------------------------------------------------
 
inline bool locus_greater_or_equal(
    std::string const& l_chromosome, size_t l_position,
    std::string const& r_chromosome, size_t r_position
) {
    // Just use the existing function, but with reverse l and r.
    return locus_less_or_equal( r_chromosome, r_position, l_chromosome, l_position );
}
 
inline bool locus_greater_or_equal(
    std::string const& l_chromosome, size_t l_position,
    std::string const& r_chromosome, size_t r_position,
    ::genesis::sequence::SequenceDict const& sequence_dict
) {
    // Just use the existing function, but with reverse l and r.
    return locus_less_or_equal( r_chromosome, r_position, l_chromosome, l_position, sequence_dict );
}
 
// Add all other overloads for GenomeLocus and SequenceDict combinations.
ADD_LOCUS_COMPARISON_OVERLOADS( locus_greater_or_equal )
 
 
inline bool operator >= ( GenomeLocus const& l, GenomeLocus const& r )
{
    return locus_greater_or_equal( l.chromosome, l.position, r.chromosome, r.position );
}
 
// -------------------------------------------------------------------------
//     Clean Up
// -------------------------------------------------------------------------
 
// We do not need to propagate this macro throughout our code base.
#undef ADD_LOCUS_COMPARISON_OVERLOADS
 
} // namespace population
} // namespace genesis
 
#endif // include guard