fst_pool_calculator.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_FUNCTION_FST_POOL_CALCULATOR_H_
#define GENESIS_POPULATION_FUNCTION_FST_POOL_CALCULATOR_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/filter/sample_counts_filter.hpp"
#include "genesis/population/sample_counts.hpp"
#include "genesis/population/variant.hpp"
 
#include <algorithm>
#include <cassert>
#include <vector>
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Fst Pool Calculator
// =================================================================================================
 
class BaseFstPoolCalculator
{
public:
 
    // -------------------------------------------------------------------------
    //     Constructors and Rule of Five
    // -------------------------------------------------------------------------
 
    BaseFstPoolCalculator() = default;
    virtual ~BaseFstPoolCalculator() = default;
 
    BaseFstPoolCalculator( BaseFstPoolCalculator const& ) = default;
    BaseFstPoolCalculator( BaseFstPoolCalculator&& )      = default;
 
    BaseFstPoolCalculator& operator= ( BaseFstPoolCalculator const& ) = default;
    BaseFstPoolCalculator& operator= ( BaseFstPoolCalculator&& )      = default;
 
    // -------------------------------------------------------------------------
    //     Calculator Functions
    // -------------------------------------------------------------------------
 
    void reset()
    {
        reset_();
        filter_stats_.clear();
    }
 
    void process( SampleCounts const& p1, SampleCounts const& p2 )
    {
        // For now, we accept two populations as input here.
        // If we ever implement an FST equation that also works for multiple populations,
        // see https://stackoverflow.com/a/9377363/4184258 for a way to implement this.
        // Or simply make this function accept the Variant instead, and store the needed sample
        // indices here, instead of the list of sample pairs stored in the FstPoolProcessor.
 
        // We only want to process samples that are both passing.
        // Otherwise, if both fail, we want to record that, but it does not really matter which
        // of the two reasons for failing we record. However, as we generally have our filter status
        // tags numerically sorted in the order in which they are applied, we here use the first
        // tag (smaller number) of the two samples that failed, as that is usually a better and
        // more general idicator of what failed. Probably does not really matter much anyway,
        // but why not. Lastly, if only one of the samples failed, we just record that one.
        // This can be done as the max of the two values, as one of them will be zero.
        if( p1.status.passing() && p2.status.passing() ) {
            process_( p1, p2 );
            ++filter_stats_[SampleCountsFilterTag::kPassed];
        } else if( ! p1.status.passing() && ! p2.status.passing() ) {
            ++filter_stats_[ std::min( p1.status.get(), p2.status.get() )];
        } else {
            ++filter_stats_[ std::max( p1.status.get(), p2.status.get() )];
        }
    }
 
    double get_result() const
    {
        return get_result_();
    }
 
    SampleCountsFilterStats get_filter_stats() const
    {
        return filter_stats_;
    }
 
    // -------------------------------------------------------------------------
    //     Abstract Members
    // -------------------------------------------------------------------------
 
protected:
 
    // Non-virtual interface.
    virtual void reset_() = 0;
    virtual void process_( SampleCounts const& p1, SampleCounts const& p2 ) = 0;
    virtual double get_result_() const = 0;
 
    // -------------------------------------------------------------------------
    //     Private Members
    // -------------------------------------------------------------------------
 
private:
 
    // Count how many sample pairs were processed here,
    // and how many of them passed or failed the filters.
    SampleCountsFilterStats filter_stats_;
 
};
 
} // namespace population
} // namespace genesis
 
#endif // include guard