cathedral_plot.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_PLOTTING_CATHEDRAL_PLOT_H_
#define GENESIS_POPULATION_PLOTTING_CATHEDRAL_PLOT_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/utils/color/color.hpp"
#include "genesis/utils/color/heat_map.hpp"
#include "genesis/utils/containers/matrix.hpp"
#include "genesis/utils/io/base_output_target.hpp"
 
#include <cassert>
#include <cmath>
#include <cstdint>
#include <deque>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
// =================================================================================================
//     Forward declacations
// =================================================================================================
 
namespace genesis {
namespace utils {
 
    class JsonDocument;
    class SvgDocument;
 
} // namespace utils
} // namespace genesis
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Cathedral Plot Parameters
// =================================================================================================
 
enum class CathedralWindowWidthMethod
{
    kExponential,
    kGeometric,
    kLinear
};
 
struct CathedralPlotParameters
{
    // Plot parameters
    size_t width = 0;
    size_t height = 0;
    CathedralWindowWidthMethod window_width_method = CathedralWindowWidthMethod::kExponential;
};
 
struct CathedralPlotRecord
{
    // Virtual destructor, just in case. We are problably not going to delete the derived class
    // through a pointer to base in our use case, but one never knows...
    virtual ~CathedralPlotRecord() = default;
 
    // Data-derived parameters from the initial input.
    std::string title;
    std::string plot_name;
    std::string chromosome_name;
    size_t chromosome_length = 0;
 
    // User-provided parameters, added here to keep track of the record.
    CathedralPlotParameters parameters;
 
    // Data matrix containing per-pixel values.
    genesis::utils::Matrix<double> value_matrix;
    std::vector<double> window_widths;
};
 
void validate_cathedral_plot_record( CathedralPlotRecord const& record );
 
double cathedral_window_width(
    CathedralPlotRecord const& record, size_t row
);
 
std::string cathedral_window_width_method_to_string( CathedralWindowWidthMethod method );
 
CathedralWindowWidthMethod cathedral_window_width_method_from_string( std::string const& method );
 
// =================================================================================================
//     Compute Matrix Functions
// =================================================================================================
 
template<class Record, class Accumulator>
void compute_cathedral_matrix(
    CathedralPlotParameters const& parameters,
    Record& record,
    Accumulator accumulator = Accumulator{}
) {
    // Prepare a result matrix for the values, of the desired dimensions.
    record.value_matrix = genesis::utils::Matrix<double>( parameters.height, parameters.width );
    record.window_widths = std::vector<double>( parameters.height );
 
    // Also store the parameters in the record, for later reference to have them in one place.
    record.parameters = parameters;
 
    // How far (in genome coordinates) do we advance between windows?
    double const chr_len  = static_cast<double>( record.chromosome_length );
    auto const advance = chr_len / static_cast<double>( parameters.width );
 
    // Compute each cell of the result. We experimented with parallelizing this loop across threads,
    // but the computation seems to be memory bound, and even when trying to avoid false sharing
    // (of writing to individual cells of the matrix in each iteration), the result was never faster
    // (and often way slower) than the single threaded code here. So let's keep it simple.
    for( size_t row = 0; row < parameters.height; ++row ) {
 
        // How wide (in genome coordinates) is each window in the current row?
        auto const window_width = cathedral_window_width( record, row );
        assert( std::isfinite( window_width ) && window_width > 0.0 );
        record.window_widths[row] = window_width;
 
        // Per row, we have a lot of overlap between the windows, up until the very last few
        // rows where windows tend to overlap less. Using this gives massive speedup,
        // as we only need to add entries once, and then remove them again once,
        // instead of computing their accumulated sums over and over again.
        // We use a deque for the entries that are in the current window,
        // and keep track of the next index in the entry vector that needs to be enqueued.
        std::deque<typename Record::Entry> queue;
        size_t entry_idx = 0;
 
        // Start a new accumulation of values for the row.
        accumulator.reset();
 
        // Assert that for one row, we accumulate and dissipate each value exactly once.
        size_t accu_cnt = 0;
        size_t diss_cnt = 0;
        (void) accu_cnt;
        (void) diss_cnt;
 
        // We move along the windows using the cur_gen_pos (in genome coordinates) to indicate
        // where in the data positions we are. As we center the windows around their pixel positions,
        // we start at half the width.
        double cur_gen_pos = - window_width / 2.0;
 
        // Iterate the pixels of the columns, computing a window for each of them,
        for( size_t col = 0; col < parameters.width; ++col ) {
            assert( cur_gen_pos + window_width >= 0.0 );
            assert( cur_gen_pos <= static_cast<double>( record.chromosome_length ));
 
            // Find the genome positions that correspond to the boundaries of the current window,
            // limited to their possible ranges, and making sure to include the last one
            // (can be a bit off due to rounding?! might need to check).
            auto l_gen_pos = static_cast<size_t>( std::max( cur_gen_pos, 0.0 ));
            auto r_gen_pos = static_cast<size_t>( std::min( cur_gen_pos + window_width, chr_len ));
            if( col == parameters.width - 1 ) {
                r_gen_pos = record.chromosome_length;
            }
            cur_gen_pos += advance;
 
            // Some checks that hold true if this function is called with correct data.
            assert( l_gen_pos <= r_gen_pos );
            assert( r_gen_pos <= record.chromosome_length );
 
            // Remove entries in the beginning of the queue that are not part of the window any more.
            while( ! queue.empty() && queue.front().position < l_gen_pos ) {
                ++diss_cnt;
                accumulator.dissipate( queue.front() );
                queue.pop_front();
            }
 
            // Now accumulate entries that need to be added for the current window.
            // In cases where due to rounding the windows do not overlap and leave entries
            // between their boundaries, this here also incluse those entries in the later window.
            // This is good, as that way, every entry is used at least once.
            while(
                entry_idx < record.entries.size() &&
                record.entries[ entry_idx ].position <= r_gen_pos
            ) {
 
                // Assert that the entries are in order.
                assert( queue.empty() || queue.back().position < record.entries[entry_idx].position );
                assert( record.entries[entry_idx].position <= record.chromosome_length );
 
                // Accumulate the values and add the entry to the queue.
                ++accu_cnt;
                accumulator.accumulate( record.entries[entry_idx] );
                queue.push_back( record.entries[entry_idx] );
 
                // Move to the next entry to be enqueued.
                ++entry_idx;
            }
            assert(
                entry_idx == record.entries.size() ||
                record.entries[entry_idx].position > r_gen_pos
            );
 
            // The queue contains entries that are exactly within the window.
            assert( queue.empty() || queue.front().position >= l_gen_pos );
            assert( queue.empty() || queue.back().position  <= r_gen_pos );
 
            // Now we have processed everything for this pixel, and can store the result.
            record.value_matrix( row, col ) = accumulator.aggregate();
        }
 
        // Assert that for one row, we accumulate and dissipate each value exactly once.
        // The dissipated ones are not including the remainder of the queue in the last window,
        // so we need to account for those here.
        assert( record.entries.size() == accu_cnt );
        assert( record.entries.size() == diss_cnt + queue.size() );
    }
}
 
// =================================================================================================
//     Storage Functions
// =================================================================================================
 
genesis::utils::JsonDocument cathedral_plot_parameters_to_json_document(
    CathedralPlotParameters const& parameters
);
 
genesis::utils::JsonDocument cathedral_plot_record_to_json_document(
    CathedralPlotRecord const& record
);
 
void save_cathedral_plot_record_to_targets(
    genesis::utils::JsonDocument const& record_document,
    genesis::utils::Matrix<double> const& record_value_matrix,
    std::shared_ptr<genesis::utils::BaseOutputTarget> json_target,
    std::shared_ptr<genesis::utils::BaseOutputTarget> csv_target
);
 
void save_cathedral_plot_record_to_files(
    genesis::utils::JsonDocument const& record_document,
    genesis::utils::Matrix<double> const& record_value_matrix,
    std::string const& base_path
);
 
void save_cathedral_plot_record_to_files(
    CathedralPlotRecord const& record,
    std::string const& base_path
);
 
std::pair<genesis::utils::JsonDocument, genesis::utils::Matrix<double>>
load_cathedral_plot_record_components_from_files(
    std::string const& base_path
);
 
CathedralPlotRecord load_cathedral_plot_record_from_files(
    std::string const& base_path
);
 
// =================================================================================================
//     Plotting Functions
// =================================================================================================
 
genesis::utils::Matrix<genesis::utils::Color> make_cathedral_plot_heatmap(
    CathedralPlotRecord const& record,
    genesis::utils::HeatmapParameters const& heatmap_parameters
);
 
genesis::utils::SvgDocument make_cathedral_plot_svg(
    CathedralPlotRecord const& record,
    genesis::utils::HeatmapParameters const& heatmap_parameters,
    genesis::utils::Matrix<genesis::utils::Color> const& image
);
 
genesis::utils::SvgDocument make_cathedral_plot_svg(
    CathedralPlotRecord const& record,
    genesis::utils::HeatmapParameters const& heatmap_parameters
);
 
} // namespace population
} // namespace genesis
 
#endif // include guard