fst_cathedral.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_FUNCTION_FST_CATHEDRAL_H_
#define GENESIS_POPULATION_FUNCTION_FST_CATHEDRAL_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/stream/variant_input_stream.hpp"
#include "genesis/population/function/fst_pool_processor.hpp"
#include "genesis/population/function/fst_pool_unbiased.hpp"
#include "genesis/population/plotting/cathedral_plot.hpp"
#include "genesis/sequence/sequence_dict.hpp"
#include "genesis/utils/containers/matrix.hpp"
#include "genesis/utils/math/compensated_sum.hpp"
 
#include <cassert>
#include <cmath>
#include <stdexcept>
#include <string>
#include <vector>
 
// =================================================================================================
//     Forward declacations
// =================================================================================================
 
namespace genesis {
namespace utils {
 
    class JsonDocument;
 
} // namespace utils
} // namespace genesis
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Fst Cathedral Plot Record
// =================================================================================================
 
struct FstCathedralPlotRecord final : public CathedralPlotRecord
{
    struct Entry
    {
        Entry( size_t pos_, double pi_w_, double pi_b_, double pi_t_ )
            : position( pos_ )
            , pi_within( pi_w_ )
            , pi_between( pi_b_ )
            , pi_total( pi_t_ )
        {}
 
        size_t position   = 0;
        double pi_within  = 0.0;
        double pi_between = 0.0;
        double pi_total   = 0.0;
    };
 
    // The actual components of FST values per position.
    std::vector<Entry> entries;
 
    // Data-derived properties.
    std::string sample_name_1;
    std::string sample_name_2;
 
    // User-provided properties.
    // Type of accumulator. We store all three pi values here independently though,
    // to keep it simple, but use this to know what estimator was used for the data.
    FstPoolCalculatorUnbiased::Estimator fst_estimator;
};
 
// =================================================================================================
//     Fst Cathedral Accumulator
// =================================================================================================
 
class FstCathedralAccumulator
{
public:
 
    // -------------------------------------------------------------------------
    //     Constructors and Rule of Five
    // -------------------------------------------------------------------------
 
    FstCathedralAccumulator(
        FstPoolCalculatorUnbiased::Estimator fst_estimator
    )
        : fst_estimator_( fst_estimator )
    {}
 
    ~FstCathedralAccumulator() = default;
 
    FstCathedralAccumulator( FstCathedralAccumulator const& ) = default;
    FstCathedralAccumulator( FstCathedralAccumulator&& )      = default;
 
    FstCathedralAccumulator& operator= ( FstCathedralAccumulator const& ) = default;
    FstCathedralAccumulator& operator= ( FstCathedralAccumulator&& )      = default;
 
    // -------------------------------------------------------------------------
    //     Accumulator Functions
    // -------------------------------------------------------------------------
 
    void accumulate( FstCathedralPlotRecord::Entry const& entry );
    void dissipate( FstCathedralPlotRecord::Entry const& entry );
    double aggregate() const;
    void reset();
 
    // -------------------------------------------------------------------------
    //     Private Member Variables
    // -------------------------------------------------------------------------
 
private:
 
    // Type of accumulator.
    FstPoolCalculatorUnbiased::Estimator fst_estimator_;
 
    // Store our accumualted values. We are using a Neumaier summation here,
    // as we might be adding and subtracting values of different orders of magnitude,
    // which would lead to large errors with the standard Kahan sum.
    utils::NeumaierSum pi_within_sum_  = 0.0;
    utils::NeumaierSum pi_between_sum_ = 0.0;
    utils::NeumaierSum pi_total_sum_   = 0.0;
    size_t value_count_ = 0;
};
 
// =================================================================================================
//     Fst Cathedral Functions
// =================================================================================================
 
std::vector<FstCathedralPlotRecord> compute_fst_cathedral_records_for_chromosome(
    VariantInputStream::Iterator& iterator,
    FstPoolProcessor& processor,
    FstPoolCalculatorUnbiased::Estimator fst_estimator,
    std::vector<std::string> const& sample_names = std::vector<std::string>{},
    std::shared_ptr<genesis::sequence::SequenceDict> const& sequence_dict = nullptr
);
 
std::vector<FstCathedralPlotRecord> compute_fst_cathedral_records(
    VariantInputStream& iterator,
    FstPoolProcessor& processor,
    FstPoolCalculatorUnbiased::Estimator fst_estimator,
    std::vector<std::string> const& sample_names = std::vector<std::string>{},
    std::shared_ptr<genesis::sequence::SequenceDict> const& sequence_dict = nullptr
);
 
inline void compute_fst_cathedral_matrix(
    CathedralPlotParameters const& parameters,
    FstCathedralPlotRecord& record
) {
    auto accu = FstCathedralAccumulator( record.fst_estimator );
    return compute_cathedral_matrix<FstCathedralPlotRecord, FstCathedralAccumulator>(
        parameters, record, accu
    );
}
 
// =================================================================================================
//     Storage Functions
// =================================================================================================
 
genesis::utils::JsonDocument fst_cathedral_plot_record_to_json_document(
    FstCathedralPlotRecord const& record
);
 
} // namespace population
} // namespace genesis
 
#endif // include guard