From 9b42c03786eaee987f93d10da79c99c2e4b75faf Mon Sep 17 00:00:00 2001 From: computations Date: Thu, 25 Apr 2019 09:52:20 +0200 Subject: [PATCH 01/22] Initial commit with DKS Currently builds, but is not tested to _actually_ work. --- src/CMakeLists.txt | 2 +- src/dks/CMakeLists.txt | 3 + src/dks/benchmark.cpp | 160 ++++++++++++++++++++++++++++++++ src/dks/benchmark.h | 27 ++++++ src/dks/dks.h | 2 + src/dks/model.cpp | 37 ++++++++ src/dks/model.h | 54 +++++++++++ src/dks/msa.cpp | 199 ++++++++++++++++++++++++++++++++++++++++ src/dks/msa.h | 45 +++++++++ src/dks/partition.cpp | 169 ++++++++++++++++++++++++++++++++++ src/dks/partition.h | 58 ++++++++++++ src/dks/test_case.cpp | 150 ++++++++++++++++++++++++++++++ src/dks/test_case.h | 112 ++++++++++++++++++++++ src/dks/tree.cpp | 164 +++++++++++++++++++++++++++++++++ src/dks/tree.h | 41 +++++++++ src/util/CMakeLists.txt | 1 + 16 files changed, 1223 insertions(+), 1 deletion(-) create mode 100644 src/dks/CMakeLists.txt create mode 100644 src/dks/benchmark.cpp create mode 100644 src/dks/benchmark.h create mode 100644 src/dks/dks.h create mode 100644 src/dks/model.cpp create mode 100644 src/dks/model.h create mode 100644 src/dks/msa.cpp create mode 100644 src/dks/msa.h create mode 100644 src/dks/partition.cpp create mode 100644 src/dks/partition.h create mode 100644 src/dks/test_case.cpp create mode 100644 src/dks/test_case.h create mode 100644 src/dks/tree.cpp create mode 100644 src/dks/tree.h diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2ebf181..23f2408 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -8,7 +8,7 @@ endif() set(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} -Wl,-undefined -Wl,dynamic_lookup") set(PLLMOD_INCLUDE_PATH ${PLLMODULES_LIBPLL_PATH}/src/) -set(PLLMODULES_COMPONENTS "optimize;algorithm;binary;msa;tree;util" CACHE STRING "pll-modules components to build") +set(PLLMODULES_COMPONENTS "optimize;algorithm;binary;msa;tree;util;dks" CACHE STRING "pll-modules components to build") foreach(module ${PLLMODULES_COMPONENTS}) MESSAGE(STATUS "Will compile pll-module ${module}") diff --git a/src/dks/CMakeLists.txt b/src/dks/CMakeLists.txt new file mode 100644 index 0000000..61fc030 --- /dev/null +++ b/src/dks/CMakeLists.txt @@ -0,0 +1,3 @@ +file(GLOB DKS_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp) +add_pllmodules_lib(pllmoddks "${DKS_SOURCES}") +target_include_directories(pllmoddks_obj PRIVATE ${PLLMOD_INCLUD_PATH}) diff --git a/src/dks/benchmark.cpp b/src/dks/benchmark.cpp new file mode 100644 index 0000000..b6caff7 --- /dev/null +++ b/src/dks/benchmark.cpp @@ -0,0 +1,160 @@ +#include "benchmark.h" +#include "msa.h" +#include "partition.h" +#include +#include +#include +#include + +namespace dks { + +inline benchmark_time_t weight_kernel_times(kernel_weight_t kw, + benchmark_result_t bmr) { + return kw[test_kernel_t::partial] * bmr[test_kernel_t::partial] + + kw[test_kernel_t::likelihood] * bmr[test_kernel_t::likelihood] + + kw[test_kernel_t::derivative] * bmr[test_kernel_t::derivative] + + kw[test_kernel_t::pmatrix] * bmr[test_kernel_t::pmatrix]; +} + +inline attributes_t best_attrib_time(const attributes_time_t &at) { + return std::max_element(at.begin(), at.end(), + [](const attributes_time_t::value_type &a, + const attributes_time_t::value_type &b) { + return a.second > b.second; + }) + ->first; +} + +kernel_weight_t suggest_weights(const msa_t &msa) { + kernel_weight_t kw; + double taxa = msa.count(); + double sites = msa.length(); + double states = msa.states(); + + kw[test_kernel_t::partial] = + 0.007156765 * taxa + -2.719444e-05 * sites + 1.328822 + states + 37.70555; + kw[test_kernel_t::likelihood] = 0.0001289232 * sites + -0.004789004 * taxa + + -0.9559354 * states + 35.20843; + kw[test_kernel_t::derivative] = 2.240574e-06 * sites + -0.003947451 * taxa + + -0.6615589 * states + 25.72586; + kw[test_kernel_t::pmatrix] = -0.0001090363 * sites + 0.001506259 * sites + + 0.2866474 * states + 1.448459; + + for (auto &kv : kw) { + kv.second = kv.second < 0.0 ? 0.0 : kv.second; + } + return kw; +} + +kernel_weight_t suggest_weights_2(const msa_t &msa) { + kernel_weight_t kw; + double taxa = msa.count(); + double sites = msa.length(); + double states = msa.states(); + + kw[test_kernel_t::partial] = + 0.4866 * sites + 437.1470 * states + 6.5094 * taxa - 3557.8645; + kw[test_kernel_t::likelihood] = + 0.327 * sites + 28.952 * states + 1.147 * taxa + -43.042; + kw[test_kernel_t::derivative] = + 0.2174 * sites + 26.1509 * states + 0.7898 * taxa + -108.6298; + kw[test_kernel_t::pmatrix] = + 3.221e-03 * sites + 2.672e+01 * states + 7.741e-01 * taxa + -2.195e+02; + + for (auto &kv : kw) { + kv.second = kv.second < 0.0 ? 0.1 : kv.second; + } + return kw; +} + +attributes_t select_kernel(const pll_partition_t *pll_partition, + const pll_msa_t *pll_msa, const kernel_weight_t &kw, + bool fast) { + msa_t msa{pll_msa}; + model_t model{pll_partition}; + return select_kernel(model, msa, kw, fast); +} + +attributes_time_t select_kernel_verbose(const model_t &model, const msa_t &msa, + const kernel_weight_t &kw, bool fast) { + return fast ? select_kernel_fast_verbose(model, msa, kw) + : select_kernel_slow_verbose(model, msa, kw); +} + +/* + * Select the attributes in pairs, so that we can eliminate some of the + * searching. The order is + * 1. SIMD + * 2. (Site repeats, tip pattern, none) + * 3. Rate scalers + * for the others, we assume a default set of parameters that are likely too be + * correct. + */ +attributes_time_t select_kernel_fast_verbose(const model_t &model, + const msa_t &msa, + const kernel_weight_t &kw) { + attributes_time_t times; + msa_compressed_t cmsa(msa); + + attributes_t attribs(false, true, true, test_cpu_t::avx); + + for (int i = test_cpu_t::avx2; i >= test_cpu_t::none; --i) { + attribs.simd = static_cast(i); + test_case_t tc(attribs); + times[attribs] = weight_kernel_times(kw, tc.benchmark(msa, model)); + std::cout << "timing " << attribs << ":" << times[attribs].count() + << std::endl; + } + + attribs = best_attrib_time(times); + + for (size_t i = 0; i < 2; ++i) { + attribs.site_repeats = false; + attribs.pattern_tip = static_cast(i); + test_case_t tc(attribs); + times[attribs] = weight_kernel_times(kw, tc.benchmark(msa, model)); + std::cout << "timing " << attribs << ":" << times[attribs].count() + << std::endl; + } + + attribs = best_attrib_time(times); + + for (size_t i = 0; i < 1; ++i) { + attribs.rate_scalers = static_cast(i); + test_case_t tc(attribs); + times[attribs] = weight_kernel_times(kw, tc.benchmark(msa, model)); + std::cout << "timing " << attribs << ":" << times[attribs].count() + << std::endl; + } + + return times; +} + +attributes_time_t select_kernel_slow_verbose(const model_t &model, + const msa_t &msa, + const kernel_weight_t &kw) { + attributes_time_t times; + msa_compressed_t cmsa(msa); + for (uint8_t bit_attribs = 0; bit_attribs < 0x4; ++bit_attribs) { + for (uint8_t simd = test_cpu_t::none; simd <= test_cpu_t::avx2; ++simd) { + if (static_cast(bit_attribs & (1 << 0)) && + static_cast(bit_attribs & (1 << 1))) { + continue; + } + + attributes_t attribs(static_cast(bit_attribs & (1 << 0)), + static_cast(bit_attribs & (1 << 1)), + static_cast(0), static_cast(simd)); + test_case_t tc(attribs); + times[attribs] = weight_kernel_times(kw, tc.benchmark(msa, model)); + } + } + return times; +} + +attributes_t select_kernel(const model_t &model, const msa_t &msa, + const kernel_weight_t &kw, bool fast) { + auto times = select_kernel_verbose(model, msa, kw, fast); + return best_attrib_time(times); +} +} // namespace dks diff --git a/src/dks/benchmark.h b/src/dks/benchmark.h new file mode 100644 index 0000000..16a365b --- /dev/null +++ b/src/dks/benchmark.h @@ -0,0 +1,27 @@ +#pragma once +#include "test_case.h" +#include + +namespace dks { +typedef std::unordered_map kernel_weight_t; +typedef std::unordered_map attributes_time_t; + +kernel_weight_t suggest_weights(const msa_t& msa); +kernel_weight_t suggest_weights_2(const msa_t& msa); + +attributes_time_t select_kernel_fast_verbose(const model_t &model, const msa_t &msa, + const kernel_weight_t &); +attributes_time_t select_kernel_slow_verbose(const model_t &model, const msa_t &msa, + const kernel_weight_t &); +attributes_t select_kernel(const pll_partition_t *, const pll_msa_t *, + const kernel_weight_t &, bool fast); + +attributes_t select_kernel(const model_t &, const msa_t &, + const kernel_weight_t &, bool fast); + +attributes_time_t select_kernel_verbose(const model_t &, const msa_t &, + const kernel_weight_t &, bool fast); + +int physical_cpu_count(); + +} // namespace dks diff --git a/src/dks/dks.h b/src/dks/dks.h new file mode 100644 index 0000000..191ab50 --- /dev/null +++ b/src/dks/dks.h @@ -0,0 +1,2 @@ +#pragma once +#include diff --git a/src/dks/model.cpp b/src/dks/model.cpp new file mode 100644 index 0000000..3532481 --- /dev/null +++ b/src/dks/model.cpp @@ -0,0 +1,37 @@ +#include "model.h" +#include + +namespace dks { + +unsigned int model_t::submodels() const { return 1; } + +unsigned int model_t::rate_categories() const { return 1; } + +uint64_t model_t::states() const { return _states; } + +const double *model_t::subst_params_raw() const { return _subst_params.data(); } + +const std::vector &model_t::subst_params() const { + return _subst_params; +} + +const double *model_t::frequencies_raw() const { return _frequencies.data(); } + +const std::vector &model_t::frequencies() const { return _frequencies; } + +const tree_t &model_t::tree() const { return _tree; } + +std::vector model_t::make_operations() const { + auto traversal_nodes = _tree.full_traverse(); + std::vector operations(traversal_nodes.size()); + unsigned int operations_count = 0; + pll_utree_create_operations(traversal_nodes.data(), traversal_nodes.size(), + nullptr, nullptr, operations.data(), nullptr, + &operations_count); + + operations.resize(operations_count); + return operations; +} + +void model_t::reset_tree() { _tree = tree_t(_tree.tip_count()); } +} // namespace dks diff --git a/src/dks/model.h b/src/dks/model.h new file mode 100644 index 0000000..7d236ca --- /dev/null +++ b/src/dks/model.h @@ -0,0 +1,54 @@ +#pragma once +#include "tree.h" +#include +#include +#include + +namespace dks { +class model_t { +public: + model_t(const msa_t &msa) : model_t{msa, 0} {}; + + model_t(size_t tip_count) : model_t{tip_count, 0} {}; + + model_t(const msa_t &msa, uint64_t seed) + : _tree{msa.count(), seed}, _states{msa.states()}, + _subst_params((_states - 1) * (_states - 2), 1.0), + _frequencies(_states, 1.0 / _states){}; + + model_t(size_t tip_count, uint64_t seed) + : _tree{tip_count, seed}, _subst_params{6, 1.0}, _frequencies{4, .25} {}; + + model_t(const pll_partition_t *pll_partition) + : model_t(pll_partition->tips, pll_partition->states, 0, + pll_partition->subst_params, pll_partition->frequencies){}; + + model_t(size_t tip_count, size_t states, size_t model_index, + double **subst_params, double **frequencies) + : _tree{tree_t(tip_count)}, _states{states}, + _subst_params{subst_params[model_index], + subst_params[model_index] + + (_states - 1) * (_states - 2)}, + _frequencies{frequencies[model_index], + frequencies[model_index] + _states} {}; + + unsigned int submodels() const; + unsigned int rate_categories() const; + uint64_t states() const; + const double *subst_params_raw() const; + const std::vector &subst_params() const; + const double *frequencies_raw() const; + const std::vector &frequencies() const; + const tree_t &tree() const; + + void reset_tree(); + + std::vector make_operations() const; + +private: + tree_t _tree; + size_t _states; + std::vector _subst_params; + std::vector _frequencies; +}; +} // namespace dks diff --git a/src/dks/msa.cpp b/src/dks/msa.cpp new file mode 100644 index 0000000..2568e59 --- /dev/null +++ b/src/dks/msa.cpp @@ -0,0 +1,199 @@ +#include "msa.h" +#include +#include +#include +#include +#include + +namespace dks { +msa_t::msa_t(const pll_msa_t *msa) { + init(msa); + _states = 4; +} + + +msa_t::msa_t(const pll_msa_t *msa, size_t states) : _states(states) { + init(msa); +} + +msa_t::msa_t(const std::string &filename, size_t states) { + if (pll_phylip_t *fd = pll_phylip_open(filename.c_str(), pll_map_generic)) { + pll_msa_t *pll_msa = nullptr; + if ((pll_msa = pll_phylip_parse_interleaved(fd)) || + (pll_msa = pll_phylip_parse_sequential(fd))) { + init(pll_msa); + pll_msa_destroy(pll_msa); + pll_phylip_close(fd); + _states = states; + return; + } else { + pll_phylip_close(fd); + } + } + if (pll_fasta_t *fd = pll_fasta_open(filename.c_str(), pll_map_fasta)) { + char *header = nullptr; + char *sequence = nullptr; + long sequence_len = 0; + long header_len = 0; + long sequence_number = 0; + long expected_sequence_length = 0; + + while (pll_fasta_getnext(fd, &header, &header_len, &sequence, &sequence_len, + &sequence_number)) { + if (expected_sequence_length == 0) { + expected_sequence_length = sequence_len; + } else if (expected_sequence_length != sequence_len) { + throw std::runtime_error("Sequences don't match in size"); + } + _sequences.emplace_back(sequence, sequence + sequence_len); + free(sequence); + _labels.emplace_back(header); + free(header); + } + _states = states; + pll_fasta_close(fd); + } + if (!valid()) { + throw std::runtime_error("failed to parse the msa file"); + } +} + +void msa_t::init(const pll_msa_t *msa) { + for (int i = 0; i < msa->count; i++) { + _labels.emplace_back(msa->label[i]); + _sequences.emplace_back(msa->sequence[i], msa->sequence[i] + msa->length); + } +} + +size_t msa_t::count() const { return _labels.size(); } + +size_t msa_t::length() const { + if (_sequences.empty()) { + return 0; + } + return _sequences.front().size(); +} + +const char *msa_t::label(size_t i) const { return _labels[i].data(); } + +const char *msa_t::sequence(size_t i) const { return _sequences[i].data(); } + +const pll_state_t *msa_t::char_map() const { + if (states() == 4) + return pll_map_nt; + if (states() == 20) + return pll_map_aa; + throw std::runtime_error("no state map found for msa"); +} + +bool msa_t::valid() const { return _sequences.size() > 0; } + +double msa_t::column_entropy() const { + uint64_t taxa_count = _sequences.size(); + if (taxa_count == 0) { + return 0.0; + } + + double entropy = 0.0; + uint64_t max_states = 0; + + size_t sequence_len = length(); + + for (size_t i = 0; i < sequence_len; i++) { + uint64_t counts[256] = {0}; + for (const auto &s : _sequences) { + counts[static_cast(s[i])] += 1; + } + + uint64_t states = 0; + double site_entropy = 0.0; + for (size_t j = 0; j < 256; j++) { + if (counts[j] == 0) + continue; + double px = static_cast(counts[j]) / taxa_count; + states += 1; + site_entropy += px * std::log2(px); + } + + entropy -= site_entropy; + max_states = std::max(max_states, states); + } + + return 1 - (entropy / sequence_len) / -std::log2(1.0 / max_states); +} + +double msa_t::row_entropy() const { + uint64_t taxa_count = _sequences.size(); + if (taxa_count == 0) { + return 0.0; + } + + double entropy = 0.0; + uint64_t max_states = 0; + + size_t sequence_len = length(); + + for (const auto &s : _sequences) { + uint64_t counts[256] = {0}; + for (size_t i = 0; i < sequence_len; i++) { + counts[static_cast(s[i])] += 1; + } + + uint64_t states = 0; + double site_entropy = 0.0; + for (size_t j = 0; j < 256; j++) { + if (counts[j] == 0) + continue; + double px = static_cast(counts[j]) / sequence_len; + states += 1; + site_entropy += px * std::log2(px); + } + + entropy -= site_entropy; + max_states = std::max(max_states, states); + } + + return 1 - (entropy / taxa_count) / -std::log2(1.0 / max_states); +} + +inline size_t msa_t::states() const { return _states; } +void msa_t::set_states(size_t s) { _states = s; } + +msa_compressed_t::msa_compressed_t(const msa_t &msa) { + char **tmp_sequences = (char **)malloc(msa.count() * sizeof(char *)); + for (size_t i = 0; i < msa.count(); i++) { + tmp_sequences[i] = (char *)malloc(msa.length() * sizeof(char)); + for (size_t j = 0; j < msa.length(); j++) { + tmp_sequences[i][j] = msa.sequence(i)[j]; + } + } + + int tmp_length = msa.length(); + unsigned int *tmp_weights = pll_compress_site_patterns( + tmp_sequences, msa.char_map(), msa.count(), &tmp_length); + + if (!tmp_weights) { + throw std::runtime_error(std::string("failed to compress the msa: ") + + pll_errmsg); + } + _weights = std::vector(tmp_weights, tmp_weights + tmp_length); + + for (size_t i = 0; i < msa.count(); i++) { + _sequences.emplace_back( + sequence_t(tmp_sequences[i], tmp_sequences[i] + tmp_length)); + + _labels.emplace_back(msa.label(i)); + } + + free(tmp_weights); + + for (size_t i = 0; i < msa.count(); i++) { + free(tmp_sequences[i]); + } + free(tmp_sequences); +} + +const unsigned int *msa_compressed_t::weights() const { + return _weights.data(); +} +} // namespace dks diff --git a/src/dks/msa.h b/src/dks/msa.h new file mode 100644 index 0000000..1585dae --- /dev/null +++ b/src/dks/msa.h @@ -0,0 +1,45 @@ +#pragma once +#include +#include +#include +#include + +namespace dks { +typedef std::vector sequence_t; +typedef std::string label_t; + +class msa_t { +public: + msa_t() = default; + msa_t(const pll_msa_t *); + msa_t(const pll_msa_t *, size_t); + msa_t(const std::string & f) : msa_t(f, 4) {}; + msa_t(const std::string &, size_t); + void init(const pll_msa_t *msa); + size_t count() const; + size_t length() const; + const char *label(size_t i) const; + const char *sequence(size_t i) const; + const pll_state_t *char_map() const; + bool valid() const; + size_t states() const; + void set_states(size_t); + + double column_entropy() const; + double row_entropy() const; + +protected: + std::vector _sequences; + std::vector _labels; + size_t _states; +}; + +class msa_compressed_t : public msa_t { +public: + msa_compressed_t(const msa_t &); + const unsigned int *weights() const; + +private: + std::vector _weights; +}; +} // namespace dks diff --git a/src/dks/partition.cpp b/src/dks/partition.cpp new file mode 100644 index 0000000..c0547ca --- /dev/null +++ b/src/dks/partition.cpp @@ -0,0 +1,169 @@ +#include "partition.h" +#include + +namespace dks { + +constexpr unsigned int partition_t::_params_indices[]; +constexpr double partition_t::_rate_cats[]; + +partition_t::partition_t(const msa_t &msa, const model_t &model, + unsigned int attributes) { + + unsigned int tip_count = msa.count(); + unsigned int inner_count = tip_count - 2; + + _partition = pll_partition_create(tip_count, // tips + inner_count, // clv_buffers + msa.states(), // states + msa.length(), // sites + model.submodels(), // rate_matrices + 2 * tip_count - 3, // prob_matrices + model.rate_categories(), // rate_cats + inner_count, // scale_buffes + attributes // attributes + ); + initialize_tips(msa); + initialize_rates(model); + update_probability_matrices(model.tree()); + alloc_sumtable(attributes); +} + +void partition_t::alloc_sumtable(unsigned int attribs) { + unsigned int align = PLL_ALIGNMENT_CPU; + switch (attribs & PLL_ATTRIB_ARCH_MASK) { + case PLL_ATTRIB_ARCH_SSE: + align = PLL_ALIGNMENT_SSE; + break; + case PLL_ATTRIB_ARCH_AVX: + case PLL_ATTRIB_ARCH_AVX2: + case PLL_ATTRIB_ARCH_AVX512: + align = PLL_ALIGNMENT_AVX; + break; + default: + break; + } + _sumtable = (double *)pll_aligned_alloc( + _partition->sites * _partition->rate_cats * _partition->states_padded * + sizeof(double), + align); +} + +void partition_t::initialize_tips(const msa_t &msa) { + for (size_t tip_id = 0; tip_id < msa.count(); tip_id++) { + pll_set_tip_states(_partition, tip_id, msa.char_map(), + msa.sequence(tip_id)); + } +} + +void partition_t::initialize_rates(const model_t &model) { + for (size_t i = 0; i < model.submodels(); i++) { + pll_set_subst_params(_partition, i, model.subst_params_raw()); + pll_set_frequencies(_partition, i, model.frequencies_raw()); + } + pll_set_category_rates(_partition, _rate_cats); +} + +void partition_t::update_probability_matrices(const tree_t &tree) { + pll_update_prob_matrices( + _partition, _params_indices, tree.matrix_indices().data(), + tree.branch_lengths().data(), tree.branch_lengths().size()); +} + +void partition_t::update_partials(const std::vector &ops) { + pll_update_partials(_partition, ops.data(), ops.size()); +} + +void partition_t::update_site_repeats(const tree_t &tree) { + update_site_repeats(tree.make_operations()); +} + +void partition_t::update_site_repeats(const model_t &model) { + update_site_repeats(model.make_operations()); +} + +void partition_t::update_site_repeats(const std::vector &ops) { + for (size_t i = 0; i < ops.size(); i++) { + pll_update_repeats(_partition, &ops[i]); + } +} + +void partition_t::set_pattern_weights(const msa_compressed_t &msa) { + pll_set_pattern_weights(_partition, msa.weights()); +} + +void partition_t::set_pattern_weights(const msa_t &) {} + +void partition_t::update_partials(const tree_t &tree) { + update_partials(tree.make_operations()); +} + +void partition_t::update_partials(const model_t &model) { + update_partials(model.tree()); +} + +void partition_t::invalidate_prob_matrix() { + *_partition->eigen_decomp_valid = 0; +} + +double partition_t::loglh(const tree_t &tree) { + pll_unode_t *parent = tree.vroot(); + pll_unode_t *child = parent->back; + return pll_compute_edge_loglikelihood( + _partition, parent->clv_index, parent->scaler_index, child->clv_index, + child->scaler_index, parent->pmatrix_index, _params_indices, nullptr); +} + +double partition_t::loglh(const model_t &model) { return loglh(model.tree()); } + +std::vector partition_t::loglh_persite(const model_t &model, + size_t sites) { + pll_unode_t *parent = model.tree().vroot(); + pll_unode_t *child = parent->back; + + std::vector persites(sites, 0.0); + + pll_compute_edge_loglikelihood(_partition, parent->clv_index, + parent->scaler_index, child->clv_index, + child->scaler_index, parent->pmatrix_index, + _params_indices, persites.data()); + return persites; +} + +unsigned int partition_t::attributes() const { return _partition->attributes; } + +void partition_t::update_sumtable(const tree_t &tree) { + pll_unode_t *parent = tree.vroot(); + pll_unode_t *child = parent->back; + if (site_repeats()) { + update_partials(tree); + } + + pll_update_sumtable(_partition, parent->clv_index, child->clv_index, + parent->scaler_index, child->scaler_index, + _params_indices, _sumtable); +} + +derivative_t partition_t::compute_derivative(const tree_t &tree, double brlen) { + pll_unode_t *parent = tree.vroot(); + pll_unode_t *child = parent->back; + + derivative_t d; + pll_compute_likelihood_derivatives( + _partition, parent->scaler_index, child->scaler_index, brlen, + _params_indices, _sumtable, &d.first, &d.second); + + return d; +} + +bool partition_t::site_repeats() const { + return attributes() & PLL_ATTRIB_SITE_REPEATS; +} + +partition_t::~partition_t() { + pll_partition_destroy(_partition); + pll_aligned_free(_sumtable); +} + +const pll_partition_t *partition_t::partition_raw() const { return _partition; } + +} // namespace dks diff --git a/src/dks/partition.h b/src/dks/partition.h new file mode 100644 index 0000000..ba01c4e --- /dev/null +++ b/src/dks/partition.h @@ -0,0 +1,58 @@ +#pragma once +#include "model.h" +#include "msa.h" +#include +#include +#include + +namespace dks { +typedef std::pair derivative_t; +class partition_t { +public: + partition_t(unsigned int tips, unsigned int clv_buffers, unsigned int states, + unsigned int sites, unsigned int rate_matrices, + unsigned int prob_matrices, unsigned int rate_cats, + unsigned int scale_buffers, unsigned int attributes) + : _partition{pll_partition_create(tips, clv_buffers, states, sites, + rate_matrices, prob_matrices, rate_cats, + scale_buffers, attributes)} {}; + partition_t(const msa_t &, const model_t &, unsigned int); + ~partition_t(); + + void initialize_tips(const msa_t &); + void initialize_rates(const model_t &model); + void set_pattern_weights(const msa_compressed_t &); + void set_pattern_weights(const msa_t &); + void update_probability_matrices(const tree_t &tree); + + void update_partials(const std::vector &); + void update_partials(const tree_t &tree); + void update_partials(const model_t &model); + + void update_site_repeats(const std::vector &); + void update_site_repeats(const tree_t &tree); + void update_site_repeats(const model_t &model); + + void invalidate_prob_matrix(); + + void update_sumtable(const tree_t &tree); + derivative_t compute_derivative(const tree_t &tree, double brlen = 1.0); + + double loglh(const tree_t &tree); + double loglh(const model_t &model); + std::vector loglh_persite(const model_t &model, size_t sites); + + bool site_repeats() const; + unsigned int attributes() const; + + const pll_partition_t *partition_raw() const; + +private: + void alloc_sumtable(unsigned int attribs); + + pll_partition_t *_partition; + double *_sumtable; + constexpr static unsigned int _params_indices[] = {0}; + constexpr static double _rate_cats[] = {1.0}; +}; +} // namespace dks diff --git a/src/dks/test_case.cpp b/src/dks/test_case.cpp new file mode 100644 index 0000000..71c08aa --- /dev/null +++ b/src/dks/test_case.cpp @@ -0,0 +1,150 @@ +#include "partition.h" +#include "test_case.h" +#include +#include +#include +#include + +namespace dks { +benchmark_result_t test_case_t::benchmark(const msa_t &msa, + const model_t &model) { + partition_t partition(msa, model, attributes()); + partition.set_pattern_weights(msa); + benchmark_result_t br; + + br[test_kernel_t::partial] = benchmark_partials(partition, model); + br[test_kernel_t::likelihood] = benchmark_likelihood(partition, model); + br[test_kernel_t::pmatrix] = benchmark_pmatrix(partition, model); + br[test_kernel_t::derivative] = benchmark_derivative(partition, model); + + return br; +} + +benchmark_time_t test_case_t::benchmark_partials(partition_t &partition, + const model_t &model) { + for (size_t i = 0; i < _trials; i++) { + partition.update_partials(model); + } + auto t1 = std::chrono::high_resolution_clock::now(); + for (size_t i = 0; i < _trials; i++) { + partition.update_partials(model); + } + auto t2 = std::chrono::high_resolution_clock::now(); + return (t2 - t1) / _trials; +} + +benchmark_time_t test_case_t::benchmark_likelihood(partition_t &partition, + const model_t &model) { + partition.update_partials(model); + for (size_t i = 0; i < _trials; i++) { + partition.loglh(model); + } + auto t1 = std::chrono::high_resolution_clock::now(); + for (size_t i = 0; i < _trials; i++) { + partition.loglh(model); + } + auto t2 = std::chrono::high_resolution_clock::now(); + return (t2 - t1) / _trials; +} + +benchmark_time_t test_case_t::benchmark_pmatrix(partition_t &partition, + const model_t &model) { + for (size_t i = 0; i < _trials; i++) { + partition.update_probability_matrices(model.tree()); + } + auto t1 = std::chrono::high_resolution_clock::now(); + for (size_t i = 0; i < _trials; i++) { + partition.update_probability_matrices(model.tree()); + } + auto t2 = std::chrono::high_resolution_clock::now(); + return (t2 - t1) / _trials; +} + +benchmark_time_t test_case_t::benchmark_derivative(partition_t &partition, + const model_t &model) { + for (size_t i = 0; i < _trials; i++) { + partition.update_sumtable(model.tree()); + partition.compute_derivative(model.tree()); + } + auto t1 = std::chrono::high_resolution_clock::now(); + for (size_t i = 0; i < _trials; i++) { + partition.update_sumtable(model.tree()); + partition.compute_derivative(model.tree()); + } + auto t2 = std::chrono::high_resolution_clock::now(); + return (t2 - t1) / _trials; +} + +benchmark_time_t +test_case_t::benchmark_update_site_repeats(partition_t &partition, + const model_t &model) { + for (size_t i = 0; i < _trials; i++) { + partition.update_site_repeats(model); + } + auto t1 = std::chrono::high_resolution_clock::now(); + for (size_t i = 0; i < _trials; i++) { + partition.update_site_repeats(model); + } + auto t2 = std::chrono::high_resolution_clock::now(); + return (t2 - t1) / _trials; +} + +attributes_t test_case_t::attributes_struct() const { + return attributes_t(_pattern_tip, _site_repeats, _rate_scalers, _cpu); +} + +unsigned int test_case_t::attributes() const { + return cpu_attributes() | misc_attributes(); +} + +unsigned int test_case_t::misc_attributes() const { + return _pattern_tip ? PLL_ATTRIB_PATTERN_TIP + : 0 | _site_repeats + ? PLL_ATTRIB_SITE_REPEATS + : 0 | _rate_scalers ? PLL_ATTRIB_RATE_SCALERS : 0; +} + +unsigned int test_case_t::cpu_attributes() const { + if (test_cpu_t::none == _cpu) { + return PLL_ATTRIB_ARCH_CPU; + } + if (test_cpu_t::sse == _cpu) { + return PLL_ATTRIB_ARCH_SSE; + } + if (test_cpu_t::avx == _cpu) { + return PLL_ATTRIB_ARCH_AVX; + } + if (test_cpu_t::avx2 == _cpu) { + return PLL_ATTRIB_ARCH_AVX2; + } + if (test_cpu_t::avx512 == _cpu) { + return PLL_ATTRIB_ARCH_AVX512; + } + throw std::runtime_error("Unrecognized CPU type"); +} + +} // namespace dks + +std::ostream &operator<<(std::ostream &stream, const dks::test_cpu_t &cpu) { + if (cpu == dks::test_cpu_t::none) { + stream << "none"; + } else if (cpu == dks::test_cpu_t::sse) { + stream << "sse"; + } else if (cpu == dks::test_cpu_t::avx) { + stream << "avx"; + } else if (cpu == dks::test_cpu_t::avx2) { + stream << "avx2"; + } else if (cpu == dks::test_cpu_t::avx512) { + stream << "avx512"; + } + return stream; +} + +std::ostream &operator<<(std::ostream &stream, + const dks::attributes_t &attribs) { + stream << "{\"cpu\": \"" << attribs.simd << "\"" << std::boolalpha + << ", \"pattern tip\": \"" << attribs.pattern_tip << "\"" + << ", \"rate scalers\": \"" << attribs.rate_scalers << "\"" + << ", \"site repeats\": \"" << attribs.site_repeats << "\"}"; + return stream; +} diff --git a/src/dks/test_case.h b/src/dks/test_case.h new file mode 100644 index 0000000..a430c97 --- /dev/null +++ b/src/dks/test_case.h @@ -0,0 +1,112 @@ +#pragma once +#include "model.h" +#include "msa.h" +#include "partition.h" +#include "tree.h" +#include +#include +#include +#include +#include + +namespace dks { +enum test_cpu_t { + none, + sse, + avx, + avx2, + avx512, +}; + +enum test_kernel_t { + partial, + likelihood, + derivative, + pmatrix, +}; + +typedef std::chrono::duration benchmark_time_t; + +typedef std::unordered_map benchmark_result_t; + +struct attributes_t { + bool pattern_tip; + bool site_repeats; + bool rate_scalers; + test_cpu_t simd; + + attributes_t() = default; + + attributes_t(bool pt, bool sr, bool rs, test_cpu_t simd) + : pattern_tip{pt}, site_repeats{sr}, rate_scalers{rs}, simd{simd} {}; + + bool operator==(const attributes_t &other) const { + return pattern_tip == other.pattern_tip && + site_repeats == other.site_repeats && + rate_scalers == other.rate_scalers && simd == other.simd; + } +}; + +class test_case_t { +public: + test_case_t() + : _cpu{test_cpu_t::none}, _trials{30}, _random_seed{0}, + _pattern_tip{false}, _site_repeats{false}, _rate_scalers{false} {} + + test_case_t(test_cpu_t cpu, bool pt, bool sr, bool rs, uint64_t seed) + : _cpu{cpu}, _trials{30}, _random_seed{seed}, _pattern_tip{pt}, + _site_repeats{sr}, _rate_scalers{rs} {} + + test_case_t(test_cpu_t cpu) : test_case_t{cpu, 0, 0, 0, 0} {} + + test_case_t(const attributes_t &attribs) + : test_case_t(attribs.simd, attribs.pattern_tip, attribs.site_repeats, + attribs.rate_scalers, 0){}; + + benchmark_result_t benchmark(const msa_t &, const model_t &); + benchmark_time_t benchmark_partials(partition_t &partition, + const model_t &model); + benchmark_time_t benchmark_likelihood(partition_t &partition, + const model_t &model); + benchmark_time_t benchmark_pmatrix(partition_t &partition, + const model_t &model); + benchmark_time_t benchmark_derivative(partition_t &partition, + const model_t &model); + benchmark_time_t benchmark_update_site_repeats(partition_t &partition, + const model_t &model); + attributes_t attributes_struct() const; + unsigned int attributes() const; + unsigned int cpu_attributes() const; + unsigned int misc_attributes() const; + +private: + test_cpu_t _cpu; + size_t _trials; + uint64_t _random_seed; + bool _pattern_tip; + bool _site_repeats; + bool _rate_scalers; +}; +} // namespace dks + +namespace std { +template <> struct hash { + typedef dks::attributes_t argument_type; + typedef size_t result_type; + result_type operator()(const argument_type &s) const noexcept { + return (s.pattern_tip << 0) ^ (s.site_repeats << 1) ^ + (s.rate_scalers << 2) ^ (s.simd << 3); + } +}; +template <> struct hash { + typedef dks::test_kernel_t argument_type; + typedef size_t result_type; + result_type operator()(const argument_type &s) const noexcept { + return static_cast(s); + } +}; +} // namespace std + +std::ostream &operator<<(std::ostream &stream, const dks::test_cpu_t &cpu); +std::ostream &operator<<(std::ostream &stream, + const dks::attributes_t &attribs); diff --git a/src/dks/tree.cpp b/src/dks/tree.cpp new file mode 100644 index 0000000..f6262d3 --- /dev/null +++ b/src/dks/tree.cpp @@ -0,0 +1,164 @@ +#include "pll.h" +#include "tree.h" +#include + +namespace dks { +int full_traverse_cb(pll_unode_t *n) { + // silence a warning + do { + (void)(n); + } while (0); + return PLL_SUCCESS; +} + +tree_t::tree_t(size_t tip_count, uint64_t random_seed) { + size_t inner_count = tip_count - 2; + + pll_unode_t *vroot = make_triplet(make_tip(), make_tip(), make_tip()); + + // generate a list of nodes + size_t node_buffer_size = tip_count + 3 * inner_count; + pll_unode_t **nodes = + (pll_unode_t **)malloc(sizeof(pll_unode_t *) * node_buffer_size); + + std::minstd_rand random_engine(random_seed); + + for (size_t i = 3; i < tip_count; i++) { + unsigned int node_count; + pll_utree_traverse(vroot, PLL_TREE_TRAVERSE_POSTORDER, full_traverse_cb, + nodes, &node_count); + + std::uniform_int_distribution<> roller(0, node_count - 1); + pll_unode_t *insert_node = nodes[roller(random_engine)]; + insert_tip(insert_node); + } + _tree = pll_utree_wraptree(vroot, tip_count); + pll_utree_reset_template_indices(vroot, tip_count); + update_internal_lists(); + + free(nodes); + if (!pll_utree_check_integrity(_tree)) { + pll_utree_destroy(_tree, nullptr); + throw std::runtime_error("Tree validation check failed"); + } +} + +tree_t::~tree_t() { pll_utree_destroy(_tree, nullptr); } + +pll_unode_t *tree_t::vroot() const { return _tree->vroot; } + +size_t tree_t::node_count() const { + return _tree->tip_count + _tree->inner_count; +} + +size_t tree_t::edge_count() const { return _tree->tip_count * 2 - 3; } + +size_t tree_t::tip_count() const { return _tree->tip_count; } + +void tree_t::insert_tip(pll_unode_t *insert_node) { + pll_unode_t *new_a = make_node(); + pll_unode_t *new_b = make_node(); + pll_unode_t *new_c = make_node(); + pll_unode_t *new_tip = make_tip(); + + pll_unode_t *saved_back = insert_node->back; + + pair_nodes(insert_node, new_a); + pair_nodes(saved_back, new_b); + pair_nodes(new_tip, new_c); + + make_circle(new_a, new_b, new_c); +} + +void tree_t::pair_nodes(pll_unode_t *a, pll_unode_t *b) { + a->back = b; + b->back = a; +} + +void tree_t::make_circle(pll_unode_t *a, pll_unode_t *b, pll_unode_t *c) { + a->next = b; + b->next = c; + c->next = a; +} + +pll_unode_t *tree_t::make_triplet(pll_unode_t *a, pll_unode_t *b, + pll_unode_t *c) { + pll_unode_t *inner_a = make_node(); + pll_unode_t *inner_b = make_node(); + pll_unode_t *inner_c = make_node(); + + pair_nodes(inner_a, a); + pair_nodes(inner_b, b); + pair_nodes(inner_c, c); + + make_circle(inner_a, inner_b, inner_c); + return inner_a; +} + +pll_unode_t *tree_t::make_node() { + pll_unode_t *a = (pll_unode_t *)malloc(sizeof(pll_unode_t)); + a->label = nullptr; + a->next = nullptr; + a->back = nullptr; + a->length = 0.1; + return a; +} + +pll_unode_t *tree_t::make_tip() { return make_node(); } + +const std::vector &tree_t::branch_lengths() const { + return _branch_lengths; +} + +const std::vector &tree_t::matrix_indices() const { + return _matrix_indices; +} + +std::vector tree_t::full_traverse() const { + unsigned int node_number = 0; + std::vector trav_nodes(node_count()); + if (!(pll_utree_traverse(_tree->vroot, PLL_TREE_TRAVERSE_POSTORDER, + full_traverse_cb, trav_nodes.data(), + &node_number))) { + throw std::runtime_error("Generic failure to traverse the tree"); + } + + return trav_nodes; +} + +void tree_t::fill_branch_lengths(const std::vector &nodes) { + if (_branch_lengths.size() != nodes.size()) { + _branch_lengths.resize(nodes.size(), 0.1); + } + for (size_t i = 0; i < nodes.size(); i++) { + _branch_lengths[i] = nodes[i]->length; + } +} + +void tree_t::fill_matrix_indices(const std::vector &nodes) { + if (_matrix_indices.size() != nodes.size()) { + _matrix_indices.resize(nodes.size(), 0.1); + } + for (size_t i = 0; i < nodes.size(); i++) { + _matrix_indices[i] = nodes[i]->pmatrix_index; + } +} + +void tree_t::update_internal_lists() { + auto nodes = full_traverse(); + fill_branch_lengths(nodes); + fill_matrix_indices(nodes); +} + +std::vector tree_t::make_operations() const { + auto traversal_nodes = full_traverse(); + std::vector operations(traversal_nodes.size()); + unsigned int operations_count = 0; + pll_utree_create_operations(traversal_nodes.data(), traversal_nodes.size(), + nullptr, nullptr, operations.data(), nullptr, + &operations_count); + + operations.resize(operations_count); + return operations; +} +} // namespace dks diff --git a/src/dks/tree.h b/src/dks/tree.h new file mode 100644 index 0000000..d1d4dd4 --- /dev/null +++ b/src/dks/tree.h @@ -0,0 +1,41 @@ +#pragma once +#include "msa.h" +#include +#include + +namespace dks { +class tree_t { +public: + tree_t(const msa_t &msa) : tree_t{msa.count(), 0} {} + tree_t(const msa_t &msa, uint64_t random_seed) + : tree_t{msa.count(), random_seed} {} + tree_t(size_t tip_count) : tree_t{tip_count, 0} {}; + tree_t(size_t tip_count, uint64_t random_seed); + pll_unode_t *vroot() const; + size_t node_count() const; + size_t edge_count() const; + size_t tip_count() const; + const std::vector &branch_lengths() const; + const std::vector &matrix_indices() const; + std::vector full_traverse() const; + std::vector make_operations() const; + ~tree_t(); + +private: + static void insert_tip(pll_unode_t *); + static void pair_nodes(pll_unode_t *, pll_unode_t *); + static void make_circle(pll_unode_t *, pll_unode_t *, pll_unode_t *); + static pll_unode_t *make_triplet(pll_unode_t *, pll_unode_t *, pll_unode_t *); + static pll_unode_t *make_node(); + static pll_unode_t *make_tip(); + + void fill_branch_lengths(const std::vector &); + void fill_matrix_indices(const std::vector &); + + void update_internal_lists(); + + pll_utree_t *_tree; + std::vector _branch_lengths; + std::vector _matrix_indices; +}; +} // namespace dks diff --git a/src/util/CMakeLists.txt b/src/util/CMakeLists.txt index 61bc93a..d56626f 100644 --- a/src/util/CMakeLists.txt +++ b/src/util/CMakeLists.txt @@ -5,3 +5,4 @@ set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${PLLMOD_CFLAGS}") add_pllmodules_lib(pllmodutil "${UTIL_SOURCES}") + From 02c13a5e47102fd2c1934243227b85a1335f661a Mon Sep 17 00:00:00 2001 From: computations Date: Mon, 29 Apr 2019 16:10:07 +0200 Subject: [PATCH 02/22] removed the old model --- src/dks/benchmark.cpp | 21 --------------------- src/dks/benchmark.h | 1 - 2 files changed, 22 deletions(-) diff --git a/src/dks/benchmark.cpp b/src/dks/benchmark.cpp index b6caff7..dba25e6 100644 --- a/src/dks/benchmark.cpp +++ b/src/dks/benchmark.cpp @@ -31,27 +31,6 @@ kernel_weight_t suggest_weights(const msa_t &msa) { double sites = msa.length(); double states = msa.states(); - kw[test_kernel_t::partial] = - 0.007156765 * taxa + -2.719444e-05 * sites + 1.328822 + states + 37.70555; - kw[test_kernel_t::likelihood] = 0.0001289232 * sites + -0.004789004 * taxa + - -0.9559354 * states + 35.20843; - kw[test_kernel_t::derivative] = 2.240574e-06 * sites + -0.003947451 * taxa + - -0.6615589 * states + 25.72586; - kw[test_kernel_t::pmatrix] = -0.0001090363 * sites + 0.001506259 * sites + - 0.2866474 * states + 1.448459; - - for (auto &kv : kw) { - kv.second = kv.second < 0.0 ? 0.0 : kv.second; - } - return kw; -} - -kernel_weight_t suggest_weights_2(const msa_t &msa) { - kernel_weight_t kw; - double taxa = msa.count(); - double sites = msa.length(); - double states = msa.states(); - kw[test_kernel_t::partial] = 0.4866 * sites + 437.1470 * states + 6.5094 * taxa - 3557.8645; kw[test_kernel_t::likelihood] = diff --git a/src/dks/benchmark.h b/src/dks/benchmark.h index 16a365b..5688dab 100644 --- a/src/dks/benchmark.h +++ b/src/dks/benchmark.h @@ -7,7 +7,6 @@ typedef std::unordered_map kernel_weight_t; typedef std::unordered_map attributes_time_t; kernel_weight_t suggest_weights(const msa_t& msa); -kernel_weight_t suggest_weights_2(const msa_t& msa); attributes_time_t select_kernel_fast_verbose(const model_t &model, const msa_t &msa, const kernel_weight_t &); From 56ebfeaa98b0c4a2a1930b3ee0aace808f1ea0c9 Mon Sep 17 00:00:00 2001 From: computations Date: Tue, 30 Apr 2019 08:54:05 +0200 Subject: [PATCH 03/22] updated dks --- src/dks/CMakeLists.txt | 20 +++++++++++-- src/dks/benchmark.cpp | 66 +++++++++++++++++++++++++++++++++++++++++- src/dks/benchmark.h | 3 +- src/dks/main.cpp | 60 ++++++++++++++++++++++++++++++++++++++ src/dks/model.cpp | 2 +- src/dks/model.h | 14 ++++++--- src/dks/test_case.h | 6 ++-- src/dks/tree.cpp | 2 +- 8 files changed, 159 insertions(+), 14 deletions(-) create mode 100644 src/dks/main.cpp diff --git a/src/dks/CMakeLists.txt b/src/dks/CMakeLists.txt index 61fc030..cb2cf9a 100644 --- a/src/dks/CMakeLists.txt +++ b/src/dks/CMakeLists.txt @@ -1,3 +1,17 @@ -file(GLOB DKS_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp) -add_pllmodules_lib(pllmoddks "${DKS_SOURCES}") -target_include_directories(pllmoddks_obj PRIVATE ${PLLMOD_INCLUD_PATH}) +add_executable(raxml-dks + main.cpp + msa.cpp + test_case.cpp + tree.cpp + partition.cpp + model.cpp + benchmark.cpp +) +set_target_properties(raxml-dks PROPERTIES + CXX_STANDARD 11 + CXX_STANDARD_REQUIRED YES +) +set(${CMAKE_BUILD_TYPE} DEBUG) +target_include_directories(raxml-dks PRIVATE ${CMAKE_SOURCE_DIR}/lib/libpll/src/) +target_link_libraries(raxml-dks pll_static) +target_compile_options(raxml-dks PRIVATE -Wall -Wextra -pedantic) diff --git a/src/dks/benchmark.cpp b/src/dks/benchmark.cpp index dba25e6..0361a77 100644 --- a/src/dks/benchmark.cpp +++ b/src/dks/benchmark.cpp @@ -6,6 +6,13 @@ #include #include +#ifdef __linux__ +#include +#include +#elif _WIN32 +// nothing +#endif + namespace dks { inline benchmark_time_t weight_kernel_times(kernel_weight_t kw, @@ -31,6 +38,27 @@ kernel_weight_t suggest_weights(const msa_t &msa) { double sites = msa.length(); double states = msa.states(); + kw[test_kernel_t::partial] = + 0.007156765 * taxa + -2.719444e-05 * sites + 1.328822 + states + 37.70555; + kw[test_kernel_t::likelihood] = 0.0001289232 * sites + -0.004789004 * taxa + + -0.9559354 * states + 35.20843; + kw[test_kernel_t::derivative] = 2.240574e-06 * sites + -0.003947451 * taxa + + -0.6615589 * states + 25.72586; + kw[test_kernel_t::pmatrix] = -0.0001090363 * sites + 0.001506259 * sites + + 0.2866474 * states + 1.448459; + + for (auto &kv : kw) { + kv.second = kv.second < 0.0 ? 0.0 : kv.second; + } + return kw; +} + +kernel_weight_t suggest_weights_2(const msa_t &msa) { + kernel_weight_t kw; + double taxa = msa.count(); + double sites = msa.length(); + double states = msa.states(); + kw[test_kernel_t::partial] = 0.4866 * sites + 437.1470 * states + 6.5094 * taxa - 3557.8645; kw[test_kernel_t::likelihood] = @@ -123,7 +151,8 @@ attributes_time_t select_kernel_slow_verbose(const model_t &model, attributes_t attribs(static_cast(bit_attribs & (1 << 0)), static_cast(bit_attribs & (1 << 1)), - static_cast(0), static_cast(simd)); + static_cast(0), + static_cast(simd)); test_case_t tc(attribs); times[attribs] = weight_kernel_times(kw, tc.benchmark(msa, model)); } @@ -136,4 +165,39 @@ attributes_t select_kernel(const model_t &model, const msa_t &msa, auto times = select_kernel_verbose(model, msa, kw, fast); return best_attrib_time(times); } + +#ifdef __linux__ +std::string build_path(size_t cpu_number) { + std::string s; + s += "/sys/devices/system/cpu/cpu"; + int leading_zeros = static_cast(std::floor(std::log10(cpu_number))); + for (int i = 0; i < leading_zeros; ++i) { + s += "0"; + } + s += std::to_string(cpu_number); + s += "/topology/core_id"; + return s; +} + +size_t get_core_id(size_t cpu_number) { + std::string filename = build_path(cpu_number); + std::ifstream f(filename.c_str()); + char buf[32]; + f.getline(buf, 32); + return std::stoul(buf); +} + +size_t physical_cpu_count() { + size_t cpu_max = 0; + size_t n_cpu = static_cast(sysconf(_SC_NPROCESSORS_ONLN)); + for (size_t i = 0; i < n_cpu; ++i) { + size_t core_id = get_core_id(i); + cpu_max = core_id > cpu_max ? core_id : cpu_max; + } + return cpu_max + 1; +} + +#elif _WIN32 +size_t physical_cpu_count() { return 0; } +#endif } // namespace dks diff --git a/src/dks/benchmark.h b/src/dks/benchmark.h index 5688dab..63e0f63 100644 --- a/src/dks/benchmark.h +++ b/src/dks/benchmark.h @@ -7,6 +7,7 @@ typedef std::unordered_map kernel_weight_t; typedef std::unordered_map attributes_time_t; kernel_weight_t suggest_weights(const msa_t& msa); +kernel_weight_t suggest_weights_2(const msa_t& msa); attributes_time_t select_kernel_fast_verbose(const model_t &model, const msa_t &msa, const kernel_weight_t &); @@ -21,6 +22,6 @@ attributes_t select_kernel(const model_t &, const msa_t &, attributes_time_t select_kernel_verbose(const model_t &, const msa_t &, const kernel_weight_t &, bool fast); -int physical_cpu_count(); +size_t physical_cpu_count(); } // namespace dks diff --git a/src/dks/main.cpp b/src/dks/main.cpp new file mode 100644 index 0000000..f311dd7 --- /dev/null +++ b/src/dks/main.cpp @@ -0,0 +1,60 @@ +#include "benchmark.h" +#include +#include +#include +#include +#include +#include + +#include + +using namespace std; + +option cli_options[] = { + {"msa", required_argument, 0, 0}, + {"states", required_argument, 0, 0}, + {0, 0, 0, 0}, +}; + +int main(int argc, char **argv) { + int opt_index; + string filename; + char c; + size_t states = 4; + while ((c = getopt_long(argc, argv, "", cli_options, &opt_index)) == 0) { + switch (opt_index) { + case 0: + filename = string(optarg); + break; + case 1: + states = std::stoi(optarg); + break; + } + } + if (filename.empty()) { + cout << "Please give an msa file with the --msa switch" << endl; + return 1; + } + + dks::msa_t msa(filename, states); + dks::model_t model(msa); + + dks::kernel_weight_t kw = dks::suggest_weights_2(msa); + + auto results = dks::select_kernel_verbose(model, msa, kw, false); + + vector> sorted_times{ + results.begin(), results.end()}; + std::sort(sorted_times.begin(), sorted_times.end(), + [](const decltype(sorted_times)::value_type &a, + const decltype(sorted_times)::value_type &b) { + return a.second < b.second; + }); + + for (const auto &kv : sorted_times) { + cout << "{\"attributes\":" << kv.first + << ", \"score\": " << kv.second.count() << "}" << endl; + } + + return 0; +} diff --git a/src/dks/model.cpp b/src/dks/model.cpp index 3532481..cea5f46 100644 --- a/src/dks/model.cpp +++ b/src/dks/model.cpp @@ -5,7 +5,7 @@ namespace dks { unsigned int model_t::submodels() const { return 1; } -unsigned int model_t::rate_categories() const { return 1; } +unsigned int model_t::rate_categories() const { return _rate_categories; } uint64_t model_t::states() const { return _states; } diff --git a/src/dks/model.h b/src/dks/model.h index 7d236ca..dd8f12c 100644 --- a/src/dks/model.h +++ b/src/dks/model.h @@ -12,20 +12,24 @@ class model_t { model_t(size_t tip_count) : model_t{tip_count, 0} {}; model_t(const msa_t &msa, uint64_t seed) - : _tree{msa.count(), seed}, _states{msa.states()}, + : _tree{msa.count(), seed}, _states{msa.states()}, _rate_categories{1}, _subst_params((_states - 1) * (_states - 2), 1.0), _frequencies(_states, 1.0 / _states){}; model_t(size_t tip_count, uint64_t seed) - : _tree{tip_count, seed}, _subst_params{6, 1.0}, _frequencies{4, .25} {}; + : _tree{tip_count, seed}, _rate_categories{1}, _subst_params{6, 1.0}, + _frequencies{4, .25} {}; model_t(const pll_partition_t *pll_partition) : model_t(pll_partition->tips, pll_partition->states, 0, - pll_partition->subst_params, pll_partition->frequencies){}; + pll_partition->rate_cats, pll_partition->subst_params, + pll_partition->frequencies, *(pll_partition->prop_invar)){}; model_t(size_t tip_count, size_t states, size_t model_index, - double **subst_params, double **frequencies) + size_t rate_categories, double **subst_params, double **frequencies, + double pinv) : _tree{tree_t(tip_count)}, _states{states}, + _rate_categories{rate_categories}, _prop_invar{pinv}, _subst_params{subst_params[model_index], subst_params[model_index] + (_states - 1) * (_states - 2)}, @@ -48,6 +52,8 @@ class model_t { private: tree_t _tree; size_t _states; + size_t _rate_categories; + double _prop_invar; std::vector _subst_params; std::vector _frequencies; }; diff --git a/src/dks/test_case.h b/src/dks/test_case.h index a430c97..9858869 100644 --- a/src/dks/test_case.h +++ b/src/dks/test_case.h @@ -37,8 +37,8 @@ struct attributes_t { attributes_t() = default; - attributes_t(bool pt, bool sr, bool rs, test_cpu_t simd) - : pattern_tip{pt}, site_repeats{sr}, rate_scalers{rs}, simd{simd} {}; + attributes_t(bool pt, bool sr, bool rs, test_cpu_t s) + : pattern_tip{pt}, site_repeats{sr}, rate_scalers{rs}, simd{s} {}; bool operator==(const attributes_t &other) const { return pattern_tip == other.pattern_tip && @@ -102,7 +102,7 @@ template <> struct hash { typedef dks::test_kernel_t argument_type; typedef size_t result_type; result_type operator()(const argument_type &s) const noexcept { - return static_cast(s); + return static_cast(s); } }; } // namespace std diff --git a/src/dks/tree.cpp b/src/dks/tree.cpp index f6262d3..397c14e 100644 --- a/src/dks/tree.cpp +++ b/src/dks/tree.cpp @@ -137,7 +137,7 @@ void tree_t::fill_branch_lengths(const std::vector &nodes) { void tree_t::fill_matrix_indices(const std::vector &nodes) { if (_matrix_indices.size() != nodes.size()) { - _matrix_indices.resize(nodes.size(), 0.1); + _matrix_indices.resize(nodes.size(), 0); } for (size_t i = 0; i < nodes.size(); i++) { _matrix_indices[i] = nodes[i]->pmatrix_index; From 6422a0e4927d0eb37e590e2dcd531bb6ac01154b Mon Sep 17 00:00:00 2001 From: computations Date: Tue, 30 Apr 2019 08:56:05 +0200 Subject: [PATCH 04/22] removed reference to physical cpu count --- src/dks/benchmark.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/dks/benchmark.h b/src/dks/benchmark.h index 63e0f63..51a9138 100644 --- a/src/dks/benchmark.h +++ b/src/dks/benchmark.h @@ -21,7 +21,4 @@ attributes_t select_kernel(const model_t &, const msa_t &, attributes_time_t select_kernel_verbose(const model_t &, const msa_t &, const kernel_weight_t &, bool fast); - -size_t physical_cpu_count(); - } // namespace dks From 6d0c3842e97befb62cdf4023346d1155394f94b9 Mon Sep 17 00:00:00 2001 From: computations Date: Tue, 30 Apr 2019 08:58:01 +0200 Subject: [PATCH 05/22] updated formatting, removed main --- src/dks/CMakeLists.txt | 20 +++----------- src/dks/benchmark.cpp | 3 +-- src/dks/benchmark.h | 14 +++++----- src/dks/main.cpp | 60 ------------------------------------------ src/dks/msa.cpp | 1 - src/dks/msa.h | 2 +- 6 files changed, 13 insertions(+), 87 deletions(-) delete mode 100644 src/dks/main.cpp diff --git a/src/dks/CMakeLists.txt b/src/dks/CMakeLists.txt index cb2cf9a..61fc030 100644 --- a/src/dks/CMakeLists.txt +++ b/src/dks/CMakeLists.txt @@ -1,17 +1,3 @@ -add_executable(raxml-dks - main.cpp - msa.cpp - test_case.cpp - tree.cpp - partition.cpp - model.cpp - benchmark.cpp -) -set_target_properties(raxml-dks PROPERTIES - CXX_STANDARD 11 - CXX_STANDARD_REQUIRED YES -) -set(${CMAKE_BUILD_TYPE} DEBUG) -target_include_directories(raxml-dks PRIVATE ${CMAKE_SOURCE_DIR}/lib/libpll/src/) -target_link_libraries(raxml-dks pll_static) -target_compile_options(raxml-dks PRIVATE -Wall -Wextra -pedantic) +file(GLOB DKS_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp) +add_pllmodules_lib(pllmoddks "${DKS_SOURCES}") +target_include_directories(pllmoddks_obj PRIVATE ${PLLMOD_INCLUD_PATH}) diff --git a/src/dks/benchmark.cpp b/src/dks/benchmark.cpp index 0361a77..fdf1106 100644 --- a/src/dks/benchmark.cpp +++ b/src/dks/benchmark.cpp @@ -151,8 +151,7 @@ attributes_time_t select_kernel_slow_verbose(const model_t &model, attributes_t attribs(static_cast(bit_attribs & (1 << 0)), static_cast(bit_attribs & (1 << 1)), - static_cast(0), - static_cast(simd)); + static_cast(0), static_cast(simd)); test_case_t tc(attribs); times[attribs] = weight_kernel_times(kw, tc.benchmark(msa, model)); } diff --git a/src/dks/benchmark.h b/src/dks/benchmark.h index 51a9138..e1220b0 100644 --- a/src/dks/benchmark.h +++ b/src/dks/benchmark.h @@ -6,13 +6,15 @@ namespace dks { typedef std::unordered_map kernel_weight_t; typedef std::unordered_map attributes_time_t; -kernel_weight_t suggest_weights(const msa_t& msa); -kernel_weight_t suggest_weights_2(const msa_t& msa); +kernel_weight_t suggest_weights(const msa_t &msa); +kernel_weight_t suggest_weights_2(const msa_t &msa); -attributes_time_t select_kernel_fast_verbose(const model_t &model, const msa_t &msa, - const kernel_weight_t &); -attributes_time_t select_kernel_slow_verbose(const model_t &model, const msa_t &msa, - const kernel_weight_t &); +attributes_time_t select_kernel_fast_verbose(const model_t &model, + const msa_t &msa, + const kernel_weight_t &); +attributes_time_t select_kernel_slow_verbose(const model_t &model, + const msa_t &msa, + const kernel_weight_t &); attributes_t select_kernel(const pll_partition_t *, const pll_msa_t *, const kernel_weight_t &, bool fast); diff --git a/src/dks/main.cpp b/src/dks/main.cpp deleted file mode 100644 index f311dd7..0000000 --- a/src/dks/main.cpp +++ /dev/null @@ -1,60 +0,0 @@ -#include "benchmark.h" -#include -#include -#include -#include -#include -#include - -#include - -using namespace std; - -option cli_options[] = { - {"msa", required_argument, 0, 0}, - {"states", required_argument, 0, 0}, - {0, 0, 0, 0}, -}; - -int main(int argc, char **argv) { - int opt_index; - string filename; - char c; - size_t states = 4; - while ((c = getopt_long(argc, argv, "", cli_options, &opt_index)) == 0) { - switch (opt_index) { - case 0: - filename = string(optarg); - break; - case 1: - states = std::stoi(optarg); - break; - } - } - if (filename.empty()) { - cout << "Please give an msa file with the --msa switch" << endl; - return 1; - } - - dks::msa_t msa(filename, states); - dks::model_t model(msa); - - dks::kernel_weight_t kw = dks::suggest_weights_2(msa); - - auto results = dks::select_kernel_verbose(model, msa, kw, false); - - vector> sorted_times{ - results.begin(), results.end()}; - std::sort(sorted_times.begin(), sorted_times.end(), - [](const decltype(sorted_times)::value_type &a, - const decltype(sorted_times)::value_type &b) { - return a.second < b.second; - }); - - for (const auto &kv : sorted_times) { - cout << "{\"attributes\":" << kv.first - << ", \"score\": " << kv.second.count() << "}" << endl; - } - - return 0; -} diff --git a/src/dks/msa.cpp b/src/dks/msa.cpp index 2568e59..aadc709 100644 --- a/src/dks/msa.cpp +++ b/src/dks/msa.cpp @@ -11,7 +11,6 @@ msa_t::msa_t(const pll_msa_t *msa) { _states = 4; } - msa_t::msa_t(const pll_msa_t *msa, size_t states) : _states(states) { init(msa); } diff --git a/src/dks/msa.h b/src/dks/msa.h index 1585dae..296ca9c 100644 --- a/src/dks/msa.h +++ b/src/dks/msa.h @@ -13,7 +13,7 @@ class msa_t { msa_t() = default; msa_t(const pll_msa_t *); msa_t(const pll_msa_t *, size_t); - msa_t(const std::string & f) : msa_t(f, 4) {}; + msa_t(const std::string &f) : msa_t(f, 4){}; msa_t(const std::string &, size_t); void init(const pll_msa_t *msa); size_t count() const; From 1f1f7c95ee7617e02a66e37e1078597916b70f50 Mon Sep 17 00:00:00 2001 From: computations Date: Tue, 30 Apr 2019 09:28:08 +0200 Subject: [PATCH 06/22] changed the dks include to be a local reference --- src/dks/dks.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dks/dks.h b/src/dks/dks.h index 191ab50..cf0a557 100644 --- a/src/dks/dks.h +++ b/src/dks/dks.h @@ -1,2 +1,2 @@ #pragma once -#include +#include "benchmark.h" From 840e6df20bd05475559bc1679408a1a43e7e2ff9 Mon Sep 17 00:00:00 2001 From: computations Date: Tue, 30 Apr 2019 09:34:30 +0200 Subject: [PATCH 07/22] fixed a typo --- src/dks/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dks/CMakeLists.txt b/src/dks/CMakeLists.txt index 61fc030..fb5b27b 100644 --- a/src/dks/CMakeLists.txt +++ b/src/dks/CMakeLists.txt @@ -1,3 +1,3 @@ file(GLOB DKS_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp) add_pllmodules_lib(pllmoddks "${DKS_SOURCES}") -target_include_directories(pllmoddks_obj PRIVATE ${PLLMOD_INCLUD_PATH}) +target_include_directories(pllmoddks_obj PRIVATE ${PLLMOD_INCLUDE_PATH}) From 798f6c812f1a5667b569bc5085caf94625e1e74a Mon Sep 17 00:00:00 2001 From: computations Date: Tue, 30 Apr 2019 09:38:25 +0200 Subject: [PATCH 08/22] removed an inline that was causing linking errors --- src/dks/msa.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dks/msa.cpp b/src/dks/msa.cpp index aadc709..7c10ffa 100644 --- a/src/dks/msa.cpp +++ b/src/dks/msa.cpp @@ -155,7 +155,7 @@ double msa_t::row_entropy() const { return 1 - (entropy / taxa_count) / -std::log2(1.0 / max_states); } -inline size_t msa_t::states() const { return _states; } +size_t msa_t::states() const { return _states; } void msa_t::set_states(size_t s) { _states = s; } msa_compressed_t::msa_compressed_t(const msa_t &msa) { From 9bb64cf48708a2ad30ac22d5986e1021dd86479a Mon Sep 17 00:00:00 2001 From: computations Date: Tue, 30 Apr 2019 17:41:53 +0200 Subject: [PATCH 09/22] started integration, added support for constraints --- src/dks/CMakeLists.txt | 4 +- src/dks/benchmark.cpp | 155 +++++++---------------------------------- src/dks/benchmark.h | 26 +++---- src/dks/test_case.cpp | 9 +++ src/dks/test_case.h | 92 ++++++++++++++++++++++++ 5 files changed, 142 insertions(+), 144 deletions(-) diff --git a/src/dks/CMakeLists.txt b/src/dks/CMakeLists.txt index fb5b27b..ef6556e 100644 --- a/src/dks/CMakeLists.txt +++ b/src/dks/CMakeLists.txt @@ -1,3 +1,5 @@ file(GLOB DKS_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp) add_pllmodules_lib(pllmoddks "${DKS_SOURCES}") -target_include_directories(pllmoddks_obj PRIVATE ${PLLMOD_INCLUDE_PATH}) +target_include_directories(pllmoddks_obj PRIVATE ${PLLMOD_INCLUDE_PATH} + ${CMAKE_CURRENT_SOURCE_DIR}/../msa + ) diff --git a/src/dks/benchmark.cpp b/src/dks/benchmark.cpp index fdf1106..bfe588e 100644 --- a/src/dks/benchmark.cpp +++ b/src/dks/benchmark.cpp @@ -32,32 +32,8 @@ inline attributes_t best_attrib_time(const attributes_time_t &at) { ->first; } -kernel_weight_t suggest_weights(const msa_t &msa) { +kernel_weight_t suggest_weights(double sites, double states, double taxa) { kernel_weight_t kw; - double taxa = msa.count(); - double sites = msa.length(); - double states = msa.states(); - - kw[test_kernel_t::partial] = - 0.007156765 * taxa + -2.719444e-05 * sites + 1.328822 + states + 37.70555; - kw[test_kernel_t::likelihood] = 0.0001289232 * sites + -0.004789004 * taxa + - -0.9559354 * states + 35.20843; - kw[test_kernel_t::derivative] = 2.240574e-06 * sites + -0.003947451 * taxa + - -0.6615589 * states + 25.72586; - kw[test_kernel_t::pmatrix] = -0.0001090363 * sites + 0.001506259 * sites + - 0.2866474 * states + 1.448459; - - for (auto &kv : kw) { - kv.second = kv.second < 0.0 ? 0.0 : kv.second; - } - return kw; -} - -kernel_weight_t suggest_weights_2(const msa_t &msa) { - kernel_weight_t kw; - double taxa = msa.count(); - double sites = msa.length(); - double states = msa.states(); kw[test_kernel_t::partial] = 0.4866 * sites + 437.1470 * states + 6.5094 * taxa - 3557.8645; @@ -74,129 +50,46 @@ kernel_weight_t suggest_weights_2(const msa_t &msa) { return kw; } +kernel_weight_t suggest_weights(const msa_t &msa) { + return suggest_weights(msa.count(), msa.length(), msa.states()); +} + +kernel_weight_t suggest_weights(const pllmod_msa_stats_t *msa, int sites, + int taxa) { + return suggest_weights(sites, taxa, msa->states); +} + attributes_t select_kernel(const pll_partition_t *pll_partition, - const pll_msa_t *pll_msa, const kernel_weight_t &kw, - bool fast) { + const pll_msa_t *pll_msa, + const kernel_weight_t &kw) { msa_t msa{pll_msa}; model_t model{pll_partition}; - return select_kernel(model, msa, kw, fast); + return select_kernel(model, msa, kw); } -attributes_time_t select_kernel_verbose(const model_t &model, const msa_t &msa, - const kernel_weight_t &kw, bool fast) { - return fast ? select_kernel_fast_verbose(model, msa, kw) - : select_kernel_slow_verbose(model, msa, kw); +attributes_t select_kernel_auto(const pll_partition_t *pll_partition, + const pll_msa_t *pll_msa) { + auto kw = suggest_weights(pll_msa); + return select_kernel(pll_partition, pll_msa, kw); } -/* - * Select the attributes in pairs, so that we can eliminate some of the - * searching. The order is - * 1. SIMD - * 2. (Site repeats, tip pattern, none) - * 3. Rate scalers - * for the others, we assume a default set of parameters that are likely too be - * correct. - */ -attributes_time_t select_kernel_fast_verbose(const model_t &model, - const msa_t &msa, - const kernel_weight_t &kw) { +attributes_time_t select_kernel_verbose(const model_t &model, const msa_t &msa, + const kernel_weight_t &kw, + attributes_generator_t att_gen) { attributes_time_t times; msa_compressed_t cmsa(msa); + for (attributes_t attribs = att_gen.next(); attribs != att_gen.end(); + attribs = att_gen.next()) { - attributes_t attribs(false, true, true, test_cpu_t::avx); - - for (int i = test_cpu_t::avx2; i >= test_cpu_t::none; --i) { - attribs.simd = static_cast(i); - test_case_t tc(attribs); - times[attribs] = weight_kernel_times(kw, tc.benchmark(msa, model)); - std::cout << "timing " << attribs << ":" << times[attribs].count() - << std::endl; - } - - attribs = best_attrib_time(times); - - for (size_t i = 0; i < 2; ++i) { - attribs.site_repeats = false; - attribs.pattern_tip = static_cast(i); - test_case_t tc(attribs); - times[attribs] = weight_kernel_times(kw, tc.benchmark(msa, model)); - std::cout << "timing " << attribs << ":" << times[attribs].count() - << std::endl; - } - - attribs = best_attrib_time(times); - - for (size_t i = 0; i < 1; ++i) { - attribs.rate_scalers = static_cast(i); test_case_t tc(attribs); times[attribs] = weight_kernel_times(kw, tc.benchmark(msa, model)); - std::cout << "timing " << attribs << ":" << times[attribs].count() - << std::endl; } - return times; -} - -attributes_time_t select_kernel_slow_verbose(const model_t &model, - const msa_t &msa, - const kernel_weight_t &kw) { - attributes_time_t times; - msa_compressed_t cmsa(msa); - for (uint8_t bit_attribs = 0; bit_attribs < 0x4; ++bit_attribs) { - for (uint8_t simd = test_cpu_t::none; simd <= test_cpu_t::avx2; ++simd) { - if (static_cast(bit_attribs & (1 << 0)) && - static_cast(bit_attribs & (1 << 1))) { - continue; - } - - attributes_t attribs(static_cast(bit_attribs & (1 << 0)), - static_cast(bit_attribs & (1 << 1)), - static_cast(0), static_cast(simd)); - test_case_t tc(attribs); - times[attribs] = weight_kernel_times(kw, tc.benchmark(msa, model)); - } - } - return times; -} +} // namespace dks attributes_t select_kernel(const model_t &model, const msa_t &msa, const kernel_weight_t &kw, bool fast) { - auto times = select_kernel_verbose(model, msa, kw, fast); + auto times = select_kernel_verbose(model, msa, kw); return best_attrib_time(times); } - -#ifdef __linux__ -std::string build_path(size_t cpu_number) { - std::string s; - s += "/sys/devices/system/cpu/cpu"; - int leading_zeros = static_cast(std::floor(std::log10(cpu_number))); - for (int i = 0; i < leading_zeros; ++i) { - s += "0"; - } - s += std::to_string(cpu_number); - s += "/topology/core_id"; - return s; -} - -size_t get_core_id(size_t cpu_number) { - std::string filename = build_path(cpu_number); - std::ifstream f(filename.c_str()); - char buf[32]; - f.getline(buf, 32); - return std::stoul(buf); -} - -size_t physical_cpu_count() { - size_t cpu_max = 0; - size_t n_cpu = static_cast(sysconf(_SC_NPROCESSORS_ONLN)); - for (size_t i = 0; i < n_cpu; ++i) { - size_t core_id = get_core_id(i); - cpu_max = core_id > cpu_max ? core_id : cpu_max; - } - return cpu_max + 1; -} - -#elif _WIN32 -size_t physical_cpu_count() { return 0; } -#endif } // namespace dks diff --git a/src/dks/benchmark.h b/src/dks/benchmark.h index e1220b0..9ab7511 100644 --- a/src/dks/benchmark.h +++ b/src/dks/benchmark.h @@ -1,4 +1,5 @@ #pragma once +#include "pll_msa.h" #include "test_case.h" #include @@ -7,20 +8,21 @@ typedef std::unordered_map kernel_weight_t; typedef std::unordered_map attributes_time_t; kernel_weight_t suggest_weights(const msa_t &msa); -kernel_weight_t suggest_weights_2(const msa_t &msa); - -attributes_time_t select_kernel_fast_verbose(const model_t &model, - const msa_t &msa, - const kernel_weight_t &); -attributes_time_t select_kernel_slow_verbose(const model_t &model, - const msa_t &msa, - const kernel_weight_t &); +kernel_weight_t suggest_weights(const pllmod_msa_stats_t *msa, int sites, + int taxa); + +attributes_time_t select_kernel_verbose(const model_t &model, const msa_t &msa, + const kernel_weight_t &); + attributes_t select_kernel(const pll_partition_t *, const pll_msa_t *, - const kernel_weight_t &, bool fast); + const kernel_weight_t &); attributes_t select_kernel(const model_t &, const msa_t &, - const kernel_weight_t &, bool fast); + const kernel_weight_t &); + +attributes_t select_kernel_auto(const pll_partition_t *pll_partition, + const pll_msa_t *pll_msa); + +attributes_t select_kernel(attributes_generator_t); -attributes_time_t select_kernel_verbose(const model_t &, const msa_t &, - const kernel_weight_t &, bool fast); } // namespace dks diff --git a/src/dks/test_case.cpp b/src/dks/test_case.cpp index 71c08aa..146fb61 100644 --- a/src/dks/test_case.cpp +++ b/src/dks/test_case.cpp @@ -6,6 +6,15 @@ #include namespace dks { + +test_cpu_t test_cpu_from_attribs(uint32_t attribs) { + if ((attribs & PLL_ATTRIB_ARCH_MASK) == 0){ + return dks::test_cpu_t::none; + } + return static_cast(32 - + __builtin_clz(attribs & PLL_ATTRIB_ARCH_MASK)); +} + benchmark_result_t test_case_t::benchmark(const msa_t &msa, const model_t &model) { partition_t partition(msa, model, attributes()); diff --git a/src/dks/test_case.h b/src/dks/test_case.h index 9858869..de43a48 100644 --- a/src/dks/test_case.h +++ b/src/dks/test_case.h @@ -16,8 +16,11 @@ enum test_cpu_t { avx, avx2, avx512, + invalid, }; +test_cpu_t test_cpu_from_attribs(uint32_t attribs); + enum test_kernel_t { partial, likelihood, @@ -45,6 +48,95 @@ struct attributes_t { site_repeats == other.site_repeats && rate_scalers == other.rate_scalers && simd == other.simd; } + + bool operator!=(const attributes_t &other) const { return !(*this == other); } + + int cpu_attrib() const { + if (simd == dks::test_cpu_t::avx2) { + return PLL_ATTRIB_ARCH_AVX2; + } + if (simd == dks::test_cpu_t::avx) { + return PLL_ATTRIB_ARCH_AVX; + } + if (simd == dks::test_cpu_t::sse) { + return PLL_ATTRIB_ARCH_SSE; + } + if (simd == dks::test_cpu_t::none) { + return PLL_ATTRIB_ARCH_CPU; + } + if (simd == dks::test_cpu_t::avx512) { + return PLL_ATTRIB_ARCH_AVX512; + } + } + + int siterepeat_attrib() const { + return site_repeats ? PLL_ATTRIB_SITE_REPEATS : 0; + } + + int patterntip_attrib() const { + return pattern_tip ? PLL_ATTRIB_PATTERN_TIP : 0; + } + + int pll_attributes() const { + return cpu_attrib() | siterepeat_attrib() | patterntip_attrib(); + } +}; + +class attributes_generator_t { +public: + attributes_generator_t() + : _off_flags{PLL_ATTRIB_AB_MASK | PLL_ATTRIB_AB_FLAG | + PLL_ATTRIB_RATE_SCALERS | PLL_ATTRIB_ARCH_AVX512}, + _on_flags{0}, _max{(1 << 11) - 1}, _state{0} {}; + + attributes_t next() { + while (!valid()) { + _state++; + if (_state > _max) { + return end(); + } + } + attributes_t attrib(_state & PLL_ATTRIB_PATTERN_TIP, + _state & PLL_ATTRIB_SITE_REPEATS, 0, + test_cpu_from_attribs(_state)); + _state++; + return attrib; + } + + attributes_t end() { + return attributes_t(false, false, false, dks::test_cpu_t::invalid); + } + + void disable(int attribs) { _off_flags |= attribs; } + void enable(int attribs) { _on_flags |= attribs; } + +private: + inline bool check_cases() const { + // do some bit math to check that a bit is only set if the corrisponding + // flag is set + + return ((_on_flags & _off_flags) | (~_on_flags & _off_flags & _state) | + (_on_flags & ~_off_flags & ~_state)) & 0x3ff; + } + + bool check_xor(int attrib) const { + return __builtin_popcount(attrib & _state) <= 1; + } + + bool valid() const { + if (!check_xor(PLL_ATTRIB_SITE_REPEATS | PLL_ATTRIB_PATTERN_TIP)) { + return false; + } + if (!check_xor(PLL_ATTRIB_ARCH_MASK)) { + return false; + } + return !check_cases(); + } + + uint32_t _off_flags; + uint32_t _on_flags; + uint32_t _max; + uint32_t _state; }; class test_case_t { From 04c86c3ea2bc363ea7555bb289e3303c67b44afb Mon Sep 17 00:00:00 2001 From: computations Date: Tue, 30 Apr 2019 17:53:36 +0200 Subject: [PATCH 10/22] removed some extra libraries --- src/dks/msa.cpp | 1 - src/dks/test_case.cpp | 1 - src/dks/test_case.h | 1 - 3 files changed, 3 deletions(-) diff --git a/src/dks/msa.cpp b/src/dks/msa.cpp index 7c10ffa..71fca67 100644 --- a/src/dks/msa.cpp +++ b/src/dks/msa.cpp @@ -2,7 +2,6 @@ #include #include #include -#include #include namespace dks { diff --git a/src/dks/test_case.cpp b/src/dks/test_case.cpp index 146fb61..fd530f5 100644 --- a/src/dks/test_case.cpp +++ b/src/dks/test_case.cpp @@ -2,7 +2,6 @@ #include "test_case.h" #include #include -#include #include namespace dks { diff --git a/src/dks/test_case.h b/src/dks/test_case.h index de43a48..57601df 100644 --- a/src/dks/test_case.h +++ b/src/dks/test_case.h @@ -4,7 +4,6 @@ #include "partition.h" #include "tree.h" #include -#include #include #include #include From d468992d6c8e891d89c19f48efd96deb4e20010c Mon Sep 17 00:00:00 2001 From: computations Date: Tue, 30 Apr 2019 18:09:14 +0200 Subject: [PATCH 11/22] changed the interface of some functions to incoperate the generator --- src/dks/benchmark.cpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/dks/benchmark.cpp b/src/dks/benchmark.cpp index bfe588e..1cdf787 100644 --- a/src/dks/benchmark.cpp +++ b/src/dks/benchmark.cpp @@ -88,8 +88,15 @@ attributes_time_t select_kernel_verbose(const model_t &model, const msa_t &msa, } // namespace dks attributes_t select_kernel(const model_t &model, const msa_t &msa, - const kernel_weight_t &kw, bool fast) { - auto times = select_kernel_verbose(model, msa, kw); + const kernel_weight_t &kw, + attributes_generator_t gen) { + auto times = select_kernel_verbose(model, msa, kw, gen); return best_attrib_time(times); } + +attributes_t select_kernel(const model_t &model, const msa_t &msa, + const kernel_weight_t &kw) { + attributes_generator_t gen; + return select_kernel(model, msa, kw, gen); +} } // namespace dks From a7fb33ba262e1f6a45856bac2a29b48116a39ec0 Mon Sep 17 00:00:00 2001 From: computations Date: Thu, 2 May 2019 14:51:01 +0200 Subject: [PATCH 12/22] cleaned up some code --- src/dks/benchmark.cpp | 1 - src/dks/test_case.h | 5 +++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/dks/benchmark.cpp b/src/dks/benchmark.cpp index 1cdf787..ef4eeaa 100644 --- a/src/dks/benchmark.cpp +++ b/src/dks/benchmark.cpp @@ -77,7 +77,6 @@ attributes_time_t select_kernel_verbose(const model_t &model, const msa_t &msa, const kernel_weight_t &kw, attributes_generator_t att_gen) { attributes_time_t times; - msa_compressed_t cmsa(msa); for (attributes_t attribs = att_gen.next(); attribs != att_gen.end(); attribs = att_gen.next()) { diff --git a/src/dks/test_case.h b/src/dks/test_case.h index 57601df..8d978e7 100644 --- a/src/dks/test_case.h +++ b/src/dks/test_case.h @@ -8,6 +8,8 @@ #include #include +#define DKS_ATTRIB_MASK 0x3ff + namespace dks { enum test_cpu_t { none, @@ -114,8 +116,7 @@ class attributes_generator_t { // do some bit math to check that a bit is only set if the corrisponding // flag is set - return ((_on_flags & _off_flags) | (~_on_flags & _off_flags & _state) | - (_on_flags & ~_off_flags & ~_state)) & 0x3ff; + return ((_off_flags & _state) | (_on_flags & ~_state)) & DKS_ATTRIB_MASK; } bool check_xor(int attrib) const { From 7a726bcbd58cc868ac95e6b16023a7e1f47eebaf Mon Sep 17 00:00:00 2001 From: computations Date: Thu, 2 May 2019 16:45:13 +0200 Subject: [PATCH 13/22] Converted the interface to be more "friendly" to pll-moduels --- src/dks/benchmark.cpp | 55 +++++------- src/dks/benchmark.h | 14 ++- src/dks/model.h | 36 +++++--- src/dks/msa.cpp | 197 ------------------------------------------ src/dks/msa.h | 45 ---------- src/dks/partition.cpp | 26 +++--- src/dks/partition.h | 16 +--- src/dks/test_case.cpp | 56 +++++++++++- src/dks/test_case.h | 79 +++++------------ src/dks/tree.h | 6 +- 10 files changed, 147 insertions(+), 383 deletions(-) delete mode 100644 src/dks/msa.cpp delete mode 100644 src/dks/msa.h diff --git a/src/dks/benchmark.cpp b/src/dks/benchmark.cpp index ef4eeaa..d34fcee 100644 --- a/src/dks/benchmark.cpp +++ b/src/dks/benchmark.cpp @@ -1,18 +1,10 @@ #include "benchmark.h" -#include "msa.h" #include "partition.h" #include #include #include #include -#ifdef __linux__ -#include -#include -#elif _WIN32 -// nothing -#endif - namespace dks { inline benchmark_time_t weight_kernel_times(kernel_weight_t kw, @@ -50,8 +42,8 @@ kernel_weight_t suggest_weights(double sites, double states, double taxa) { return kw; } -kernel_weight_t suggest_weights(const msa_t &msa) { - return suggest_weights(msa.count(), msa.length(), msa.states()); +kernel_weight_t suggest_weights(const msa_t &msa, unsigned int states) { + return suggest_weights(msa.size(), msa[0].size(), states); } kernel_weight_t suggest_weights(const pllmod_msa_stats_t *msa, int sites, @@ -59,43 +51,36 @@ kernel_weight_t suggest_weights(const pllmod_msa_stats_t *msa, int sites, return suggest_weights(sites, taxa, msa->states); } -attributes_t select_kernel(const pll_partition_t *pll_partition, - const pll_msa_t *pll_msa, - const kernel_weight_t &kw) { - msa_t msa{pll_msa}; - model_t model{pll_partition}; - return select_kernel(model, msa, kw); -} - -attributes_t select_kernel_auto(const pll_partition_t *pll_partition, - const pll_msa_t *pll_msa) { - auto kw = suggest_weights(pll_msa); - return select_kernel(pll_partition, pll_msa, kw); -} - attributes_time_t select_kernel_verbose(const model_t &model, const msa_t &msa, + const msa_weight_t &weights, + const pll_state_t *charmap, const kernel_weight_t &kw, attributes_generator_t att_gen) { attributes_time_t times; for (attributes_t attribs = att_gen.next(); attribs != att_gen.end(); attribs = att_gen.next()) { - test_case_t tc(attribs); - times[attribs] = weight_kernel_times(kw, tc.benchmark(msa, model)); + test_case_t tc(attribs, charmap); + times[attribs] = weight_kernel_times(kw, tc.benchmark(msa, weights, model)); } return times; -} // namespace dks +} -attributes_t select_kernel(const model_t &model, const msa_t &msa, - const kernel_weight_t &kw, - attributes_generator_t gen) { - auto times = select_kernel_verbose(model, msa, kw, gen); - return best_attrib_time(times); +unsigned int select_kernel_auto(const msa_t &msa, const msa_weight_t &weights, + const pll_state_t *charmap, unsigned int states, + unsigned int rate_cats, + attributes_generator_t gen) { + auto kw = suggest_weights(weights.size(), msa.size(), states); + model_t model{msa, states, rate_cats}; + auto result = select_kernel_verbose(model, msa, weights, charmap, kw, gen); + return best_attrib_time(result).pll_attributes(); } -attributes_t select_kernel(const model_t &model, const msa_t &msa, - const kernel_weight_t &kw) { +unsigned int select_kernel_auto(const msa_t &msa, const msa_weight_t &weights, + const pll_state_t *charmap, unsigned int states, + unsigned int rate_cats) { attributes_generator_t gen; - return select_kernel(model, msa, kw, gen); + return select_kernel_auto(msa, weights, charmap, states, rate_cats, gen); } + } // namespace dks diff --git a/src/dks/benchmark.h b/src/dks/benchmark.h index 9ab7511..ad93723 100644 --- a/src/dks/benchmark.h +++ b/src/dks/benchmark.h @@ -7,12 +7,15 @@ namespace dks { typedef std::unordered_map kernel_weight_t; typedef std::unordered_map attributes_time_t; +typedef std::vector> msa_t; + kernel_weight_t suggest_weights(const msa_t &msa); kernel_weight_t suggest_weights(const pllmod_msa_stats_t *msa, int sites, int taxa); attributes_time_t select_kernel_verbose(const model_t &model, const msa_t &msa, - const kernel_weight_t &); + const kernel_weight_t &, + attributes_generator_t); attributes_t select_kernel(const pll_partition_t *, const pll_msa_t *, const kernel_weight_t &); @@ -23,6 +26,13 @@ attributes_t select_kernel(const model_t &, const msa_t &, attributes_t select_kernel_auto(const pll_partition_t *pll_partition, const pll_msa_t *pll_msa); -attributes_t select_kernel(attributes_generator_t); +unsigned int select_kernel_auto(const msa_t &msa, const msa_weight_t &weights, + const pll_state_t *charmap, unsigned int states, + unsigned int rate_cats, + attributes_generator_t gen); + +unsigned int select_kernel_auto(const msa_t &msa, const msa_weight_t &weights, + const pll_state_t *charmap, unsigned int states, + unsigned int rate_cats); } // namespace dks diff --git a/src/dks/model.h b/src/dks/model.h index dd8f12c..e15fdcb 100644 --- a/src/dks/model.h +++ b/src/dks/model.h @@ -5,29 +5,32 @@ #include namespace dks { +typedef std::vector> msa_t; +typedef std::vector msa_weight_t; class model_t { public: - model_t(const msa_t &msa) : model_t{msa, 0} {}; + model_t(const msa_t &msa, unsigned int states) : model_t{msa, states, 0} {}; - model_t(size_t tip_count) : model_t{tip_count, 0} {}; + model_t(size_t tip_count, unsigned int states) + : model_t{tip_count, states, 0} {}; - model_t(const msa_t &msa, uint64_t seed) - : _tree{msa.count(), seed}, _states{msa.states()}, _rate_categories{1}, - _subst_params((_states - 1) * (_states - 2), 1.0), + model_t(const msa_t &msa, unsigned int states, uint64_t seed) + : _tree{msa.size(), seed}, _states{states}, _rate_categories{1}, + _prop_invar{0.0}, _subst_params((_states - 1) * (_states - 2), 1.0), _frequencies(_states, 1.0 / _states){}; - model_t(size_t tip_count, uint64_t seed) - : _tree{tip_count, seed}, _rate_categories{1}, _subst_params{6, 1.0}, - _frequencies{4, .25} {}; + model_t(size_t tip_count, unsigned int states, uint64_t seed) + : _tree{tip_count, seed}, _states{states}, _rate_categories{1}, + _prop_invar{0.0}, _subst_params{6, 1.0}, _frequencies{4, .25} {}; model_t(const pll_partition_t *pll_partition) : model_t(pll_partition->tips, pll_partition->states, 0, pll_partition->rate_cats, pll_partition->subst_params, pll_partition->frequencies, *(pll_partition->prop_invar)){}; - model_t(size_t tip_count, size_t states, size_t model_index, - size_t rate_categories, double **subst_params, double **frequencies, - double pinv) + model_t(size_t tip_count, unsigned int states, size_t model_index, + unsigned int rate_categories, double **subst_params, + double **frequencies, double pinv) : _tree{tree_t(tip_count)}, _states{states}, _rate_categories{rate_categories}, _prop_invar{pinv}, _subst_params{subst_params[model_index], @@ -36,6 +39,13 @@ class model_t { _frequencies{frequencies[model_index], frequencies[model_index] + _states} {}; + model_t(size_t tip_count, unsigned int states, unsigned int rate_categories, + unsigned int seed) + : _tree{tip_count, seed}, _states{states}, + _rate_categories{rate_categories}, _prop_invar{0.0}, + _subst_params((_states - 1) * (_states - 2) / 2, 1.0), + _frequencies(_states, 1.0 / _states) {} + unsigned int submodels() const; unsigned int rate_categories() const; uint64_t states() const; @@ -51,8 +61,8 @@ class model_t { private: tree_t _tree; - size_t _states; - size_t _rate_categories; + unsigned int _states; + unsigned int _rate_categories; double _prop_invar; std::vector _subst_params; std::vector _frequencies; diff --git a/src/dks/msa.cpp b/src/dks/msa.cpp deleted file mode 100644 index 71fca67..0000000 --- a/src/dks/msa.cpp +++ /dev/null @@ -1,197 +0,0 @@ -#include "msa.h" -#include -#include -#include -#include - -namespace dks { -msa_t::msa_t(const pll_msa_t *msa) { - init(msa); - _states = 4; -} - -msa_t::msa_t(const pll_msa_t *msa, size_t states) : _states(states) { - init(msa); -} - -msa_t::msa_t(const std::string &filename, size_t states) { - if (pll_phylip_t *fd = pll_phylip_open(filename.c_str(), pll_map_generic)) { - pll_msa_t *pll_msa = nullptr; - if ((pll_msa = pll_phylip_parse_interleaved(fd)) || - (pll_msa = pll_phylip_parse_sequential(fd))) { - init(pll_msa); - pll_msa_destroy(pll_msa); - pll_phylip_close(fd); - _states = states; - return; - } else { - pll_phylip_close(fd); - } - } - if (pll_fasta_t *fd = pll_fasta_open(filename.c_str(), pll_map_fasta)) { - char *header = nullptr; - char *sequence = nullptr; - long sequence_len = 0; - long header_len = 0; - long sequence_number = 0; - long expected_sequence_length = 0; - - while (pll_fasta_getnext(fd, &header, &header_len, &sequence, &sequence_len, - &sequence_number)) { - if (expected_sequence_length == 0) { - expected_sequence_length = sequence_len; - } else if (expected_sequence_length != sequence_len) { - throw std::runtime_error("Sequences don't match in size"); - } - _sequences.emplace_back(sequence, sequence + sequence_len); - free(sequence); - _labels.emplace_back(header); - free(header); - } - _states = states; - pll_fasta_close(fd); - } - if (!valid()) { - throw std::runtime_error("failed to parse the msa file"); - } -} - -void msa_t::init(const pll_msa_t *msa) { - for (int i = 0; i < msa->count; i++) { - _labels.emplace_back(msa->label[i]); - _sequences.emplace_back(msa->sequence[i], msa->sequence[i] + msa->length); - } -} - -size_t msa_t::count() const { return _labels.size(); } - -size_t msa_t::length() const { - if (_sequences.empty()) { - return 0; - } - return _sequences.front().size(); -} - -const char *msa_t::label(size_t i) const { return _labels[i].data(); } - -const char *msa_t::sequence(size_t i) const { return _sequences[i].data(); } - -const pll_state_t *msa_t::char_map() const { - if (states() == 4) - return pll_map_nt; - if (states() == 20) - return pll_map_aa; - throw std::runtime_error("no state map found for msa"); -} - -bool msa_t::valid() const { return _sequences.size() > 0; } - -double msa_t::column_entropy() const { - uint64_t taxa_count = _sequences.size(); - if (taxa_count == 0) { - return 0.0; - } - - double entropy = 0.0; - uint64_t max_states = 0; - - size_t sequence_len = length(); - - for (size_t i = 0; i < sequence_len; i++) { - uint64_t counts[256] = {0}; - for (const auto &s : _sequences) { - counts[static_cast(s[i])] += 1; - } - - uint64_t states = 0; - double site_entropy = 0.0; - for (size_t j = 0; j < 256; j++) { - if (counts[j] == 0) - continue; - double px = static_cast(counts[j]) / taxa_count; - states += 1; - site_entropy += px * std::log2(px); - } - - entropy -= site_entropy; - max_states = std::max(max_states, states); - } - - return 1 - (entropy / sequence_len) / -std::log2(1.0 / max_states); -} - -double msa_t::row_entropy() const { - uint64_t taxa_count = _sequences.size(); - if (taxa_count == 0) { - return 0.0; - } - - double entropy = 0.0; - uint64_t max_states = 0; - - size_t sequence_len = length(); - - for (const auto &s : _sequences) { - uint64_t counts[256] = {0}; - for (size_t i = 0; i < sequence_len; i++) { - counts[static_cast(s[i])] += 1; - } - - uint64_t states = 0; - double site_entropy = 0.0; - for (size_t j = 0; j < 256; j++) { - if (counts[j] == 0) - continue; - double px = static_cast(counts[j]) / sequence_len; - states += 1; - site_entropy += px * std::log2(px); - } - - entropy -= site_entropy; - max_states = std::max(max_states, states); - } - - return 1 - (entropy / taxa_count) / -std::log2(1.0 / max_states); -} - -size_t msa_t::states() const { return _states; } -void msa_t::set_states(size_t s) { _states = s; } - -msa_compressed_t::msa_compressed_t(const msa_t &msa) { - char **tmp_sequences = (char **)malloc(msa.count() * sizeof(char *)); - for (size_t i = 0; i < msa.count(); i++) { - tmp_sequences[i] = (char *)malloc(msa.length() * sizeof(char)); - for (size_t j = 0; j < msa.length(); j++) { - tmp_sequences[i][j] = msa.sequence(i)[j]; - } - } - - int tmp_length = msa.length(); - unsigned int *tmp_weights = pll_compress_site_patterns( - tmp_sequences, msa.char_map(), msa.count(), &tmp_length); - - if (!tmp_weights) { - throw std::runtime_error(std::string("failed to compress the msa: ") + - pll_errmsg); - } - _weights = std::vector(tmp_weights, tmp_weights + tmp_length); - - for (size_t i = 0; i < msa.count(); i++) { - _sequences.emplace_back( - sequence_t(tmp_sequences[i], tmp_sequences[i] + tmp_length)); - - _labels.emplace_back(msa.label(i)); - } - - free(tmp_weights); - - for (size_t i = 0; i < msa.count(); i++) { - free(tmp_sequences[i]); - } - free(tmp_sequences); -} - -const unsigned int *msa_compressed_t::weights() const { - return _weights.data(); -} -} // namespace dks diff --git a/src/dks/msa.h b/src/dks/msa.h deleted file mode 100644 index 296ca9c..0000000 --- a/src/dks/msa.h +++ /dev/null @@ -1,45 +0,0 @@ -#pragma once -#include -#include -#include -#include - -namespace dks { -typedef std::vector sequence_t; -typedef std::string label_t; - -class msa_t { -public: - msa_t() = default; - msa_t(const pll_msa_t *); - msa_t(const pll_msa_t *, size_t); - msa_t(const std::string &f) : msa_t(f, 4){}; - msa_t(const std::string &, size_t); - void init(const pll_msa_t *msa); - size_t count() const; - size_t length() const; - const char *label(size_t i) const; - const char *sequence(size_t i) const; - const pll_state_t *char_map() const; - bool valid() const; - size_t states() const; - void set_states(size_t); - - double column_entropy() const; - double row_entropy() const; - -protected: - std::vector _sequences; - std::vector _labels; - size_t _states; -}; - -class msa_compressed_t : public msa_t { -public: - msa_compressed_t(const msa_t &); - const unsigned int *weights() const; - -private: - std::vector _weights; -}; -} // namespace dks diff --git a/src/dks/partition.cpp b/src/dks/partition.cpp index c0547ca..6268f4a 100644 --- a/src/dks/partition.cpp +++ b/src/dks/partition.cpp @@ -7,23 +7,25 @@ constexpr unsigned int partition_t::_params_indices[]; constexpr double partition_t::_rate_cats[]; partition_t::partition_t(const msa_t &msa, const model_t &model, - unsigned int attributes) { + const msa_weight_t &weights, + const pll_state_t *charmap, unsigned int attributes) { - unsigned int tip_count = msa.count(); + unsigned int tip_count = msa.size(); unsigned int inner_count = tip_count - 2; _partition = pll_partition_create(tip_count, // tips inner_count, // clv_buffers - msa.states(), // states - msa.length(), // sites + model.states(), // states + msa[0].size(), // sites model.submodels(), // rate_matrices 2 * tip_count - 3, // prob_matrices model.rate_categories(), // rate_cats inner_count, // scale_buffes attributes // attributes ); - initialize_tips(msa); + initialize_tips(msa, charmap); initialize_rates(model); + set_pattern_weights(weights); update_probability_matrices(model.tree()); alloc_sumtable(attributes); } @@ -48,10 +50,10 @@ void partition_t::alloc_sumtable(unsigned int attribs) { align); } -void partition_t::initialize_tips(const msa_t &msa) { - for (size_t tip_id = 0; tip_id < msa.count(); tip_id++) { - pll_set_tip_states(_partition, tip_id, msa.char_map(), - msa.sequence(tip_id)); +void partition_t::initialize_tips(const msa_t &msa, + const pll_state_t *charmap) { + for (size_t tip_id = 0; tip_id < msa.size(); tip_id++) { + pll_set_tip_states(_partition, tip_id, charmap, msa[tip_id].data()); } } @@ -87,12 +89,10 @@ void partition_t::update_site_repeats(const std::vector &ops) { } } -void partition_t::set_pattern_weights(const msa_compressed_t &msa) { - pll_set_pattern_weights(_partition, msa.weights()); +void partition_t::set_pattern_weights(const msa_weight_t &weights) { + pll_set_pattern_weights(_partition, weights.data()); } -void partition_t::set_pattern_weights(const msa_t &) {} - void partition_t::update_partials(const tree_t &tree) { update_partials(tree.make_operations()); } diff --git a/src/dks/partition.h b/src/dks/partition.h index ba01c4e..ca18a11 100644 --- a/src/dks/partition.h +++ b/src/dks/partition.h @@ -1,6 +1,5 @@ #pragma once #include "model.h" -#include "msa.h" #include #include #include @@ -9,20 +8,13 @@ namespace dks { typedef std::pair derivative_t; class partition_t { public: - partition_t(unsigned int tips, unsigned int clv_buffers, unsigned int states, - unsigned int sites, unsigned int rate_matrices, - unsigned int prob_matrices, unsigned int rate_cats, - unsigned int scale_buffers, unsigned int attributes) - : _partition{pll_partition_create(tips, clv_buffers, states, sites, - rate_matrices, prob_matrices, rate_cats, - scale_buffers, attributes)} {}; - partition_t(const msa_t &, const model_t &, unsigned int); + partition_t(const msa_t &, const model_t &, const msa_weight_t &, + const pll_state_t *, unsigned int); ~partition_t(); - void initialize_tips(const msa_t &); + void initialize_tips(const msa_t &, const pll_state_t *); void initialize_rates(const model_t &model); - void set_pattern_weights(const msa_compressed_t &); - void set_pattern_weights(const msa_t &); + void set_pattern_weights(const msa_weight_t &); void update_probability_matrices(const tree_t &tree); void update_partials(const std::vector &); diff --git a/src/dks/test_case.cpp b/src/dks/test_case.cpp index fd530f5..9fc5323 100644 --- a/src/dks/test_case.cpp +++ b/src/dks/test_case.cpp @@ -2,22 +2,72 @@ #include "test_case.h" #include #include +#include #include namespace dks { test_cpu_t test_cpu_from_attribs(uint32_t attribs) { - if ((attribs & PLL_ATTRIB_ARCH_MASK) == 0){ + if ((attribs & PLL_ATTRIB_ARCH_MASK) == 0) { return dks::test_cpu_t::none; } return static_cast(32 - __builtin_clz(attribs & PLL_ATTRIB_ARCH_MASK)); } +int attributes_t::cpu_attrib() const { + if (simd == dks::test_cpu_t::avx2) { + return PLL_ATTRIB_ARCH_AVX2; + } + if (simd == dks::test_cpu_t::avx) { + return PLL_ATTRIB_ARCH_AVX; + } + if (simd == dks::test_cpu_t::sse) { + return PLL_ATTRIB_ARCH_SSE; + } + if (simd == dks::test_cpu_t::none) { + return PLL_ATTRIB_ARCH_CPU; + } + if (simd == dks::test_cpu_t::avx512) { + return PLL_ATTRIB_ARCH_AVX512; + } + throw std::runtime_error("Unrecognized CPU type"); +} + +int attributes_t::siterepeat_attrib() const { + return site_repeats ? PLL_ATTRIB_SITE_REPEATS : 0; +} + +int attributes_t::patterntip_attrib() const { + return pattern_tip ? PLL_ATTRIB_PATTERN_TIP : 0; +} + +int attributes_t::pll_attributes() const { + return cpu_attrib() | siterepeat_attrib() | patterntip_attrib(); +} + +attributes_t attributes_generator_t::next() { + while (!valid()) { + _state++; + if (_state > _max) { + return end(); + } + } + attributes_t attrib(_state & PLL_ATTRIB_PATTERN_TIP, + _state & PLL_ATTRIB_SITE_REPEATS, 0, + test_cpu_from_attribs(_state)); + _state++; + return attrib; +} + +attributes_t attributes_generator_t::end() { + return attributes_t(false, false, false, dks::test_cpu_t::invalid); +} + benchmark_result_t test_case_t::benchmark(const msa_t &msa, + const msa_weight_t &weights, const model_t &model) { - partition_t partition(msa, model, attributes()); - partition.set_pattern_weights(msa); + partition_t partition(msa, model, weights, _charmap, attributes()); benchmark_result_t br; br[test_kernel_t::partial] = benchmark_partials(partition, model); diff --git a/src/dks/test_case.h b/src/dks/test_case.h index 8d978e7..4170919 100644 --- a/src/dks/test_case.h +++ b/src/dks/test_case.h @@ -1,9 +1,9 @@ #pragma once #include "model.h" -#include "msa.h" #include "partition.h" #include "tree.h" #include +#include #include #include #include @@ -52,35 +52,13 @@ struct attributes_t { bool operator!=(const attributes_t &other) const { return !(*this == other); } - int cpu_attrib() const { - if (simd == dks::test_cpu_t::avx2) { - return PLL_ATTRIB_ARCH_AVX2; - } - if (simd == dks::test_cpu_t::avx) { - return PLL_ATTRIB_ARCH_AVX; - } - if (simd == dks::test_cpu_t::sse) { - return PLL_ATTRIB_ARCH_SSE; - } - if (simd == dks::test_cpu_t::none) { - return PLL_ATTRIB_ARCH_CPU; - } - if (simd == dks::test_cpu_t::avx512) { - return PLL_ATTRIB_ARCH_AVX512; - } - } + int cpu_attrib() const; - int siterepeat_attrib() const { - return site_repeats ? PLL_ATTRIB_SITE_REPEATS : 0; - } + int siterepeat_attrib() const; - int patterntip_attrib() const { - return pattern_tip ? PLL_ATTRIB_PATTERN_TIP : 0; - } + int patterntip_attrib() const; - int pll_attributes() const { - return cpu_attrib() | siterepeat_attrib() | patterntip_attrib(); - } + int pll_attributes() const; }; class attributes_generator_t { @@ -88,25 +66,11 @@ class attributes_generator_t { attributes_generator_t() : _off_flags{PLL_ATTRIB_AB_MASK | PLL_ATTRIB_AB_FLAG | PLL_ATTRIB_RATE_SCALERS | PLL_ATTRIB_ARCH_AVX512}, - _on_flags{0}, _max{(1 << 11) - 1}, _state{0} {}; - - attributes_t next() { - while (!valid()) { - _state++; - if (_state > _max) { - return end(); - } - } - attributes_t attrib(_state & PLL_ATTRIB_PATTERN_TIP, - _state & PLL_ATTRIB_SITE_REPEATS, 0, - test_cpu_from_attribs(_state)); - _state++; - return attrib; - } + _on_flags{0}, _max{DKS_ATTRIB_MASK + 1}, _state{0} {}; - attributes_t end() { - return attributes_t(false, false, false, dks::test_cpu_t::invalid); - } + attributes_t next(); + + attributes_t end(); void disable(int attribs) { _off_flags |= attribs; } void enable(int attribs) { _on_flags |= attribs; } @@ -119,11 +83,11 @@ class attributes_generator_t { return ((_off_flags & _state) | (_on_flags & ~_state)) & DKS_ATTRIB_MASK; } - bool check_xor(int attrib) const { + inline bool check_xor(int attrib) const { return __builtin_popcount(attrib & _state) <= 1; } - bool valid() const { + inline bool valid() const { if (!check_xor(PLL_ATTRIB_SITE_REPEATS | PLL_ATTRIB_PATTERN_TIP)) { return false; } @@ -141,21 +105,17 @@ class attributes_generator_t { class test_case_t { public: - test_case_t() - : _cpu{test_cpu_t::none}, _trials{30}, _random_seed{0}, - _pattern_tip{false}, _site_repeats{false}, _rate_scalers{false} {} - - test_case_t(test_cpu_t cpu, bool pt, bool sr, bool rs, uint64_t seed) - : _cpu{cpu}, _trials{30}, _random_seed{seed}, _pattern_tip{pt}, - _site_repeats{sr}, _rate_scalers{rs} {} - - test_case_t(test_cpu_t cpu) : test_case_t{cpu, 0, 0, 0, 0} {} + test_case_t(test_cpu_t cpu, bool pt, bool sr, bool rs, uint64_t seed, + const pll_state_t *charmap) + : _charmap{charmap}, _cpu{cpu}, _trials{30}, _random_seed{seed}, + _pattern_tip{pt}, _site_repeats{sr}, _rate_scalers{rs} {} - test_case_t(const attributes_t &attribs) + test_case_t(const attributes_t &attribs, const pll_state_t *charmap) : test_case_t(attribs.simd, attribs.pattern_tip, attribs.site_repeats, - attribs.rate_scalers, 0){}; + attribs.rate_scalers, 0, charmap){}; - benchmark_result_t benchmark(const msa_t &, const model_t &); + benchmark_result_t benchmark(const msa_t &, const msa_weight_t &, + const model_t &); benchmark_time_t benchmark_partials(partition_t &partition, const model_t &model); benchmark_time_t benchmark_likelihood(partition_t &partition, @@ -172,6 +132,7 @@ class test_case_t { unsigned int misc_attributes() const; private: + const pll_state_t *_charmap; test_cpu_t _cpu; size_t _trials; uint64_t _random_seed; diff --git a/src/dks/tree.h b/src/dks/tree.h index d1d4dd4..39c534a 100644 --- a/src/dks/tree.h +++ b/src/dks/tree.h @@ -1,14 +1,12 @@ #pragma once -#include "msa.h" +#include "pll.h" +#include #include #include namespace dks { class tree_t { public: - tree_t(const msa_t &msa) : tree_t{msa.count(), 0} {} - tree_t(const msa_t &msa, uint64_t random_seed) - : tree_t{msa.count(), random_seed} {} tree_t(size_t tip_count) : tree_t{tip_count, 0} {}; tree_t(size_t tip_count, uint64_t random_seed); pll_unode_t *vroot() const; From 79b1f114f2845ed2c6cae248b3bf0eb6ed1b4bac Mon Sep 17 00:00:00 2001 From: computations Date: Fri, 3 May 2019 11:31:24 +0200 Subject: [PATCH 14/22] reverting change relating to valid states in attrib generator --- src/dks/test_case.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/dks/test_case.h b/src/dks/test_case.h index 4170919..0b056b1 100644 --- a/src/dks/test_case.h +++ b/src/dks/test_case.h @@ -80,7 +80,8 @@ class attributes_generator_t { // do some bit math to check that a bit is only set if the corrisponding // flag is set - return ((_off_flags & _state) | (_on_flags & ~_state)) & DKS_ATTRIB_MASK; + return ((_on_flags & _off_flags) | (~_on_flags & _off_flags & _state) | + (_on_flags & ~_off_flags & ~_state)) & DKS_ATTRIB_MASK; } inline bool check_xor(int attrib) const { From 1c441185e5bb0e1936525eca3c2fff9a2b838c35 Mon Sep 17 00:00:00 2001 From: computations Date: Fri, 3 May 2019 14:02:23 +0200 Subject: [PATCH 15/22] added a 'raw' clv option --- src/dks/benchmark.cpp | 31 +++++++++++++++++++++++++++++++ src/dks/benchmark.h | 20 +++++++++++++++++--- src/dks/partition.cpp | 30 ++++++++++++++++++++++++++++++ src/dks/partition.h | 4 ++++ src/dks/test_case.cpp | 22 ++++++++++++++++++---- src/dks/test_case.h | 13 ++++++++++++- 6 files changed, 112 insertions(+), 8 deletions(-) diff --git a/src/dks/benchmark.cpp b/src/dks/benchmark.cpp index d34fcee..588aba0 100644 --- a/src/dks/benchmark.cpp +++ b/src/dks/benchmark.cpp @@ -51,6 +51,21 @@ kernel_weight_t suggest_weights(const pllmod_msa_stats_t *msa, int sites, return suggest_weights(sites, taxa, msa->states); } +attributes_time_t +select_kernel_verbose(const model_t &model, + const std::vector> &clvs, + const msa_weight_t &weights, const kernel_weight_t &kw, + attributes_generator_t att_gen) { + attributes_time_t times; + for (attributes_t attribs = att_gen.next(); attribs != att_gen.end(); + attribs = att_gen.next()) { + + test_case_t tc(attribs); + times[attribs] = weight_kernel_times(kw, tc.benchmark(clvs, weights, model)); + } + return times; +} + attributes_time_t select_kernel_verbose(const model_t &model, const msa_t &msa, const msa_weight_t &weights, const pll_state_t *charmap, @@ -83,4 +98,20 @@ unsigned int select_kernel_auto(const msa_t &msa, const msa_weight_t &weights, return select_kernel_auto(msa, weights, charmap, states, rate_cats, gen); } +unsigned int select_kernel_auto(const std::vector> &clvs, + const msa_weight_t &weights, + unsigned int states, unsigned int rate_cats, + attributes_generator_t gen) { + auto kw = suggest_weights(weights.size(), clvs.size(), states); + model_t model{clvs.size(), states, rate_cats}; + auto result = select_kernel_verbose(model, clvs, weights, kw, gen); + return best_attrib_time(result).pll_attributes(); +} + +unsigned int select_kernel_auto(const std::vector> &clvs, + const msa_weight_t &weights, + unsigned int states, unsigned int rate_cats) { + attributes_generator_t gen; + return select_kernel_auto(clvs, weights, states, rate_cats, gen); +} } // namespace dks diff --git a/src/dks/benchmark.h b/src/dks/benchmark.h index ad93723..c03326b 100644 --- a/src/dks/benchmark.h +++ b/src/dks/benchmark.h @@ -13,9 +13,12 @@ kernel_weight_t suggest_weights(const msa_t &msa); kernel_weight_t suggest_weights(const pllmod_msa_stats_t *msa, int sites, int taxa); -attributes_time_t select_kernel_verbose(const model_t &model, const msa_t &msa, - const kernel_weight_t &, - attributes_generator_t); +template +attributes_time_t select_kernel_verbose(const model_t &model, const T &msa, + const msa_weight_t &weights, + const pll_state_t *charmap, + const kernel_weight_t &kw, + attributes_generator_t att_gen); attributes_t select_kernel(const pll_partition_t *, const pll_msa_t *, const kernel_weight_t &); @@ -35,4 +38,15 @@ unsigned int select_kernel_auto(const msa_t &msa, const msa_weight_t &weights, const pll_state_t *charmap, unsigned int states, unsigned int rate_cats); +unsigned int select_kernel_auto(const std::vector> &clvs, + const msa_weight_t &weights, + const pll_state_t *charmap, unsigned int states, + unsigned int rate_cats, + attributes_generator_t gen); + +unsigned int select_kernel_auto(const std::vector> &clvs, + const msa_weight_t &weights, + const pll_state_t *charmap, unsigned int states, + unsigned int rate_cats); + } // namespace dks diff --git a/src/dks/partition.cpp b/src/dks/partition.cpp index 6268f4a..dab3640 100644 --- a/src/dks/partition.cpp +++ b/src/dks/partition.cpp @@ -30,6 +30,30 @@ partition_t::partition_t(const msa_t &msa, const model_t &model, alloc_sumtable(attributes); } +partition_t::partition_t(const std::vector> &clvs, + const model_t &model, const msa_weight_t &weights, + const pll_state_t *charmap, unsigned int attributes) { + + unsigned int tip_count = clvs.size(); + unsigned int inner_count = tip_count - 2; + + _partition = pll_partition_create(tip_count, // tips + inner_count, // clv_buffers + model.states(), // states + clvs[0].size(), // sites + model.submodels(), // rate_matrices + 2 * tip_count - 3, // prob_matrices + model.rate_categories(), // rate_cats + inner_count, // scale_buffes + attributes // attributes + ); + initialize_tips(clvs); + initialize_rates(model); + set_pattern_weights(weights); + update_probability_matrices(model.tree()); + alloc_sumtable(attributes); +} + void partition_t::alloc_sumtable(unsigned int attribs) { unsigned int align = PLL_ALIGNMENT_CPU; switch (attribs & PLL_ATTRIB_ARCH_MASK) { @@ -57,6 +81,12 @@ void partition_t::initialize_tips(const msa_t &msa, } } +void partition_t::initialize_tips(const std::vector> &clvs){ + for (size_t tip_id = 0; tip_id < clvs.size(); tip_id++) { + pll_set_tip_clv(_partition, tip_id, clvs[tip_id].data(), PLL_FALSE); + } +} + void partition_t::initialize_rates(const model_t &model) { for (size_t i = 0; i < model.submodels(); i++) { pll_set_subst_params(_partition, i, model.subst_params_raw()); diff --git a/src/dks/partition.h b/src/dks/partition.h index ca18a11..51bd620 100644 --- a/src/dks/partition.h +++ b/src/dks/partition.h @@ -10,9 +10,13 @@ class partition_t { public: partition_t(const msa_t &, const model_t &, const msa_weight_t &, const pll_state_t *, unsigned int); + partition_t(const std::vector> &, const model_t &, + const msa_weight_t &, const pll_state_t *, unsigned int); ~partition_t(); void initialize_tips(const msa_t &, const pll_state_t *); + void initialize_tips(const std::vector> &); + void initialize_rates(const model_t &model); void set_pattern_weights(const msa_weight_t &); void update_probability_matrices(const tree_t &tree); diff --git a/src/dks/test_case.cpp b/src/dks/test_case.cpp index 9fc5323..179e99b 100644 --- a/src/dks/test_case.cpp +++ b/src/dks/test_case.cpp @@ -64,10 +64,8 @@ attributes_t attributes_generator_t::end() { return attributes_t(false, false, false, dks::test_cpu_t::invalid); } -benchmark_result_t test_case_t::benchmark(const msa_t &msa, - const msa_weight_t &weights, - const model_t &model) { - partition_t partition(msa, model, weights, _charmap, attributes()); +benchmark_result_t test_case_t::run_benchmarks(partition_t &partition, + const model_t &model) { benchmark_result_t br; br[test_kernel_t::partial] = benchmark_partials(partition, model); @@ -78,6 +76,22 @@ benchmark_result_t test_case_t::benchmark(const msa_t &msa, return br; } +benchmark_result_t +test_case_t::benchmark(const std::vector> &clvs , + const msa_weight_t &weights, const model_t &model) { + partition_t partition(clvs, model, weights, _charmap, attributes()); + + return run_benchmarks(partition, model); +} + +benchmark_result_t test_case_t::benchmark(const msa_t &msa, + const msa_weight_t &weights, + const model_t &model) { + partition_t partition(msa, model, weights, _charmap, attributes()); + + return run_benchmarks(partition, model); +} + benchmark_time_t test_case_t::benchmark_partials(partition_t &partition, const model_t &model) { for (size_t i = 0; i < _trials; i++) { diff --git a/src/dks/test_case.h b/src/dks/test_case.h index 0b056b1..ec50d95 100644 --- a/src/dks/test_case.h +++ b/src/dks/test_case.h @@ -81,7 +81,8 @@ class attributes_generator_t { // flag is set return ((_on_flags & _off_flags) | (~_on_flags & _off_flags & _state) | - (_on_flags & ~_off_flags & ~_state)) & DKS_ATTRIB_MASK; + (_on_flags & ~_off_flags & ~_state)) & + DKS_ATTRIB_MASK; } inline bool check_xor(int attrib) const { @@ -111,12 +112,22 @@ class test_case_t { : _charmap{charmap}, _cpu{cpu}, _trials{30}, _random_seed{seed}, _pattern_tip{pt}, _site_repeats{sr}, _rate_scalers{rs} {} + test_case_t(const attributes_t &attribs) + : test_case_t(attribs.simd, attribs.pattern_tip, attribs.site_repeats, + attribs.rate_scalers, 0, nullptr){}; + test_case_t(const attributes_t &attribs, const pll_state_t *charmap) : test_case_t(attribs.simd, attribs.pattern_tip, attribs.site_repeats, attribs.rate_scalers, 0, charmap){}; benchmark_result_t benchmark(const msa_t &, const msa_weight_t &, const model_t &); + benchmark_result_t benchmark(const std::vector> &, + const msa_weight_t &, const model_t &); + + benchmark_result_t run_benchmarks(partition_t &partition, + const model_t &model); + benchmark_time_t benchmark_partials(partition_t &partition, const model_t &model); benchmark_time_t benchmark_likelihood(partition_t &partition, From 9d84e9f828a1a332ae1a7bb0f6c698bad5a875b1 Mon Sep 17 00:00:00 2001 From: computations Date: Fri, 3 May 2019 14:04:45 +0200 Subject: [PATCH 16/22] fixed the benchmark header --- src/dks/benchmark.h | 6 ++---- src/dks/partition.cpp | 2 +- src/dks/partition.h | 2 +- src/dks/test_case.cpp | 2 +- 4 files changed, 5 insertions(+), 7 deletions(-) diff --git a/src/dks/benchmark.h b/src/dks/benchmark.h index c03326b..07bb03d 100644 --- a/src/dks/benchmark.h +++ b/src/dks/benchmark.h @@ -40,13 +40,11 @@ unsigned int select_kernel_auto(const msa_t &msa, const msa_weight_t &weights, unsigned int select_kernel_auto(const std::vector> &clvs, const msa_weight_t &weights, - const pll_state_t *charmap, unsigned int states, - unsigned int rate_cats, + unsigned int states, unsigned int rate_cats, attributes_generator_t gen); unsigned int select_kernel_auto(const std::vector> &clvs, const msa_weight_t &weights, - const pll_state_t *charmap, unsigned int states, - unsigned int rate_cats); + unsigned int states, unsigned int rate_cats); } // namespace dks diff --git a/src/dks/partition.cpp b/src/dks/partition.cpp index dab3640..ea11290 100644 --- a/src/dks/partition.cpp +++ b/src/dks/partition.cpp @@ -32,7 +32,7 @@ partition_t::partition_t(const msa_t &msa, const model_t &model, partition_t::partition_t(const std::vector> &clvs, const model_t &model, const msa_weight_t &weights, - const pll_state_t *charmap, unsigned int attributes) { + unsigned int attributes) { unsigned int tip_count = clvs.size(); unsigned int inner_count = tip_count - 2; diff --git a/src/dks/partition.h b/src/dks/partition.h index 51bd620..472f9bb 100644 --- a/src/dks/partition.h +++ b/src/dks/partition.h @@ -11,7 +11,7 @@ class partition_t { partition_t(const msa_t &, const model_t &, const msa_weight_t &, const pll_state_t *, unsigned int); partition_t(const std::vector> &, const model_t &, - const msa_weight_t &, const pll_state_t *, unsigned int); + const msa_weight_t &, unsigned int); ~partition_t(); void initialize_tips(const msa_t &, const pll_state_t *); diff --git a/src/dks/test_case.cpp b/src/dks/test_case.cpp index 179e99b..96ae3f8 100644 --- a/src/dks/test_case.cpp +++ b/src/dks/test_case.cpp @@ -79,7 +79,7 @@ benchmark_result_t test_case_t::run_benchmarks(partition_t &partition, benchmark_result_t test_case_t::benchmark(const std::vector> &clvs , const msa_weight_t &weights, const model_t &model) { - partition_t partition(clvs, model, weights, _charmap, attributes()); + partition_t partition(clvs, model, weights, attributes()); return run_benchmarks(partition, model); } From 8795a2c0b557e737eac21f389c477baca6b705fe Mon Sep 17 00:00:00 2001 From: computations Date: Fri, 3 May 2019 14:26:07 +0200 Subject: [PATCH 17/22] reverting a reversion, formatting --- src/dks/benchmark.cpp | 3 ++- src/dks/benchmark.h | 9 +++++++-- src/dks/partition.cpp | 3 ++- src/dks/test_case.cpp | 2 +- src/dks/test_case.h | 6 ++---- 5 files changed, 14 insertions(+), 9 deletions(-) diff --git a/src/dks/benchmark.cpp b/src/dks/benchmark.cpp index 588aba0..85c3468 100644 --- a/src/dks/benchmark.cpp +++ b/src/dks/benchmark.cpp @@ -61,7 +61,8 @@ select_kernel_verbose(const model_t &model, attribs = att_gen.next()) { test_case_t tc(attribs); - times[attribs] = weight_kernel_times(kw, tc.benchmark(clvs, weights, model)); + times[attribs] = + weight_kernel_times(kw, tc.benchmark(clvs, weights, model)); } return times; } diff --git a/src/dks/benchmark.h b/src/dks/benchmark.h index 07bb03d..378539b 100644 --- a/src/dks/benchmark.h +++ b/src/dks/benchmark.h @@ -13,13 +13,18 @@ kernel_weight_t suggest_weights(const msa_t &msa); kernel_weight_t suggest_weights(const pllmod_msa_stats_t *msa, int sites, int taxa); -template -attributes_time_t select_kernel_verbose(const model_t &model, const T &msa, +attributes_time_t select_kernel_verbose(const model_t &model, const msa_t &msa, const msa_weight_t &weights, const pll_state_t *charmap, const kernel_weight_t &kw, attributes_generator_t att_gen); +attributes_time_t +select_kernel_verbose(const model_t &model, + const std::vector> &clvs, + const msa_weight_t &weights, const kernel_weight_t &kw, + attributes_generator_t att_gen); + attributes_t select_kernel(const pll_partition_t *, const pll_msa_t *, const kernel_weight_t &); diff --git a/src/dks/partition.cpp b/src/dks/partition.cpp index ea11290..fd750a4 100644 --- a/src/dks/partition.cpp +++ b/src/dks/partition.cpp @@ -81,7 +81,8 @@ void partition_t::initialize_tips(const msa_t &msa, } } -void partition_t::initialize_tips(const std::vector> &clvs){ +void partition_t::initialize_tips( + const std::vector> &clvs) { for (size_t tip_id = 0; tip_id < clvs.size(); tip_id++) { pll_set_tip_clv(_partition, tip_id, clvs[tip_id].data(), PLL_FALSE); } diff --git a/src/dks/test_case.cpp b/src/dks/test_case.cpp index 96ae3f8..1d1585e 100644 --- a/src/dks/test_case.cpp +++ b/src/dks/test_case.cpp @@ -77,7 +77,7 @@ benchmark_result_t test_case_t::run_benchmarks(partition_t &partition, } benchmark_result_t -test_case_t::benchmark(const std::vector> &clvs , +test_case_t::benchmark(const std::vector> &clvs, const msa_weight_t &weights, const model_t &model) { partition_t partition(clvs, model, weights, attributes()); diff --git a/src/dks/test_case.h b/src/dks/test_case.h index ec50d95..62e4ee5 100644 --- a/src/dks/test_case.h +++ b/src/dks/test_case.h @@ -80,9 +80,7 @@ class attributes_generator_t { // do some bit math to check that a bit is only set if the corrisponding // flag is set - return ((_on_flags & _off_flags) | (~_on_flags & _off_flags & _state) | - (_on_flags & ~_off_flags & ~_state)) & - DKS_ATTRIB_MASK; + return ((_off_flags & _state) | (_on_flags & ~_state)) & DKS_ATTRIB_MASK; } inline bool check_xor(int attrib) const { @@ -126,7 +124,7 @@ class test_case_t { const msa_weight_t &, const model_t &); benchmark_result_t run_benchmarks(partition_t &partition, - const model_t &model); + const model_t &model); benchmark_time_t benchmark_partials(partition_t &partition, const model_t &model); From 503893143bc138d8e8a13f9c59ebc2681381ce39 Mon Sep 17 00:00:00 2001 From: computations Date: Tue, 7 May 2019 15:25:13 +0200 Subject: [PATCH 18/22] adds more interafaces into dks --- src/dks/benchmark.cpp | 52 +++++++++++++++++++++++++++++++++++++++++++ src/dks/benchmark.h | 28 ++++++++++++++++------- src/dks/test_case.h | 6 ++--- 3 files changed, 74 insertions(+), 12 deletions(-) diff --git a/src/dks/benchmark.cpp b/src/dks/benchmark.cpp index 85c3468..e1cfcd5 100644 --- a/src/dks/benchmark.cpp +++ b/src/dks/benchmark.cpp @@ -7,6 +7,18 @@ namespace dks { +msa_t convert_pll_msa_t(const pll_msa_t *pll_msa) { + msa_t msa; + msa.reserve(pll_msa->count); + for (size_t i = 0; i < pll_msa->count; i++) { + msa.emplace_back(pll_msa->length); + for (size_t j = 0; j < pll_msa->length; j++) { + msa[i][j] = pll_msa->sequence[i][j]; + } + } + return msa; +} + inline benchmark_time_t weight_kernel_times(kernel_weight_t kw, benchmark_result_t bmr) { return kw[test_kernel_t::partial] * bmr[test_kernel_t::partial] + @@ -82,6 +94,46 @@ attributes_time_t select_kernel_verbose(const model_t &model, const msa_t &msa, return times; } +unsigned int select_kernel_auto(const pll_partition_t *pll_partition, + const pll_msa_t *pll_msa, + const pll_state_t *charmap, + const kernel_weight_t &kw, + attributes_generator_t gen) { + auto msa = convert_pll_msa_t(pll_msa); + model_t model{pll_partition}; + msa_weight_t weights(pll_partition->pattern_weights, + pll_partition->pattern_weights + pll_msa->length); + return best_attrib_time( + select_kernel_verbose(model, msa, weights, charmap, kw, gen)) + .pll_attributes(); +} + +unsigned int select_kernel_auto(const pll_partition_t *pll_partition, + const pll_msa_t *pll_msa, + const pll_state_t *charmap) { + auto kw = + suggest_weights(pll_msa->length, pll_msa->count, pll_partition->states); + attributes_generator_t gen; + return select_kernel_auto(pll_partition, pll_msa, charmap, kw, gen); +} + +unsigned int select_kernel_auto(const pll_partition_t *pll_partition, + const pll_msa_t *pll_msa, + const pll_state_t *charmap, + attributes_generator_t gen) { + auto kw = + suggest_weights(pll_msa->length, pll_msa->count, pll_partition->states); + return select_kernel_auto(pll_partition, pll_msa, charmap, kw, gen); +} + +unsigned int select_kernel_auto(const pll_partition_t *pll_partition, + const pll_msa_t *pll_msa, + const pll_state_t *charmap, + const msa_weight_t &weights) { + attributes_generator_t gen; + return select_kernel_auto(pll_partition, pll_msa, charmap, weights, gen); +} + unsigned int select_kernel_auto(const msa_t &msa, const msa_weight_t &weights, const pll_state_t *charmap, unsigned int states, unsigned int rate_cats, diff --git a/src/dks/benchmark.h b/src/dks/benchmark.h index 378539b..8ae47f7 100644 --- a/src/dks/benchmark.h +++ b/src/dks/benchmark.h @@ -7,8 +7,6 @@ namespace dks { typedef std::unordered_map kernel_weight_t; typedef std::unordered_map attributes_time_t; -typedef std::vector> msa_t; - kernel_weight_t suggest_weights(const msa_t &msa); kernel_weight_t suggest_weights(const pllmod_msa_stats_t *msa, int sites, int taxa); @@ -25,12 +23,6 @@ select_kernel_verbose(const model_t &model, const msa_weight_t &weights, const kernel_weight_t &kw, attributes_generator_t att_gen); -attributes_t select_kernel(const pll_partition_t *, const pll_msa_t *, - const kernel_weight_t &); - -attributes_t select_kernel(const model_t &, const msa_t &, - const kernel_weight_t &); - attributes_t select_kernel_auto(const pll_partition_t *pll_partition, const pll_msa_t *pll_msa); @@ -52,4 +44,24 @@ unsigned int select_kernel_auto(const std::vector> &clvs, const msa_weight_t &weights, unsigned int states, unsigned int rate_cats); +unsigned int select_kernel_auto(const pll_partition_t *pll_partition, + const pll_msa_t *pll_msa, + const pll_state_t *charmap, + const msa_weight_t &weights, + attributes_generator_t gen); + +unsigned int select_kernel_auto(const pll_partition_t *pll_partition, + const pll_msa_t *pll_msa, + const pll_state_t *charmap); + +unsigned int select_kernel_auto(const pll_partition_t *pll_partition, + const pll_msa_t *pll_msa, + const pll_state_t *charmap, + attributes_generator_t gen); + +unsigned int select_kernel_auto(const pll_partition_t *pll_partition, + const pll_msa_t *pll_msa, + const pll_state_t *charmap, + const msa_weight_t &weights); + } // namespace dks diff --git a/src/dks/test_case.h b/src/dks/test_case.h index 62e4ee5..179ba8c 100644 --- a/src/dks/test_case.h +++ b/src/dks/test_case.h @@ -8,8 +8,6 @@ #include #include -#define DKS_ATTRIB_MASK 0x3ff - namespace dks { enum test_cpu_t { none, @@ -66,7 +64,7 @@ class attributes_generator_t { attributes_generator_t() : _off_flags{PLL_ATTRIB_AB_MASK | PLL_ATTRIB_AB_FLAG | PLL_ATTRIB_RATE_SCALERS | PLL_ATTRIB_ARCH_AVX512}, - _on_flags{0}, _max{DKS_ATTRIB_MASK + 1}, _state{0} {}; + _on_flags{0}, _max{PLL_ATTRIB_MASK + 1}, _state{0} {}; attributes_t next(); @@ -80,7 +78,7 @@ class attributes_generator_t { // do some bit math to check that a bit is only set if the corrisponding // flag is set - return ((_off_flags & _state) | (_on_flags & ~_state)) & DKS_ATTRIB_MASK; + return ((_off_flags & _state) | (_on_flags & ~_state)) & PLL_ATTRIB_MASK; } inline bool check_xor(int attrib) const { From 3794627e502b460013a819834a277a78fcb6b5ba Mon Sep 17 00:00:00 2001 From: computations Date: Wed, 8 May 2019 09:50:36 +0200 Subject: [PATCH 19/22] fixs a signed compare warning --- src/dks/benchmark.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dks/benchmark.cpp b/src/dks/benchmark.cpp index e1cfcd5..e74605c 100644 --- a/src/dks/benchmark.cpp +++ b/src/dks/benchmark.cpp @@ -10,9 +10,9 @@ namespace dks { msa_t convert_pll_msa_t(const pll_msa_t *pll_msa) { msa_t msa; msa.reserve(pll_msa->count); - for (size_t i = 0; i < pll_msa->count; i++) { + for (int i = 0; i < pll_msa->count; i++) { msa.emplace_back(pll_msa->length); - for (size_t j = 0; j < pll_msa->length; j++) { + for (int j = 0; j < pll_msa->length; j++) { msa[i][j] = pll_msa->sequence[i][j]; } } From 6ac86e1dfd003cc05e5ef3014d6e399beeccc628 Mon Sep 17 00:00:00 2001 From: computations Date: Wed, 8 May 2019 09:50:54 +0200 Subject: [PATCH 20/22] bumps libpll version up --- libs/libpll | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libs/libpll b/libs/libpll index fb40c17..1ac5ceb 160000 --- a/libs/libpll +++ b/libs/libpll @@ -1 +1 @@ -Subproject commit fb40c1766ed24039969ae2367ba007c58e0ee2c7 +Subproject commit 1ac5cebae551b95bbe903f5be0014c1ecada4a58 From 186f5dde03947eee1f2df7efa8f0ad457abc0b57 Mon Sep 17 00:00:00 2001 From: computations Date: Wed, 8 May 2019 09:58:25 +0200 Subject: [PATCH 21/22] Eliminates some dead code --- src/dks/benchmark.cpp | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/dks/benchmark.cpp b/src/dks/benchmark.cpp index e74605c..913acbe 100644 --- a/src/dks/benchmark.cpp +++ b/src/dks/benchmark.cpp @@ -126,14 +126,6 @@ unsigned int select_kernel_auto(const pll_partition_t *pll_partition, return select_kernel_auto(pll_partition, pll_msa, charmap, kw, gen); } -unsigned int select_kernel_auto(const pll_partition_t *pll_partition, - const pll_msa_t *pll_msa, - const pll_state_t *charmap, - const msa_weight_t &weights) { - attributes_generator_t gen; - return select_kernel_auto(pll_partition, pll_msa, charmap, weights, gen); -} - unsigned int select_kernel_auto(const msa_t &msa, const msa_weight_t &weights, const pll_state_t *charmap, unsigned int states, unsigned int rate_cats, From d7b03e12bceacb79f5b836c58ad2cd0625ae3860 Mon Sep 17 00:00:00 2001 From: computations Date: Wed, 8 May 2019 10:50:48 +0200 Subject: [PATCH 22/22] Replace `#pragma once` with real include guards --- src/dks/benchmark.h | 4 +++- src/dks/dks.h | 6 +++++- src/dks/model.h | 4 +++- src/dks/partition.h | 4 +++- src/dks/test_case.h | 4 +++- src/dks/tree.h | 4 +++- 6 files changed, 20 insertions(+), 6 deletions(-) diff --git a/src/dks/benchmark.h b/src/dks/benchmark.h index 8ae47f7..701f995 100644 --- a/src/dks/benchmark.h +++ b/src/dks/benchmark.h @@ -1,4 +1,5 @@ -#pragma once +#ifndef DKS_BENCHMARK_H_ +#define DKS_BENCHMARK_H_ #include "pll_msa.h" #include "test_case.h" #include @@ -65,3 +66,4 @@ unsigned int select_kernel_auto(const pll_partition_t *pll_partition, const msa_weight_t &weights); } // namespace dks +#endif// DKS_BENCHMARK_H_ diff --git a/src/dks/dks.h b/src/dks/dks.h index cf0a557..3a81feb 100644 --- a/src/dks/dks.h +++ b/src/dks/dks.h @@ -1,2 +1,6 @@ -#pragma once +#ifndef DKS_H_ +#define DKS_H_ + #include "benchmark.h" + +#endif diff --git a/src/dks/model.h b/src/dks/model.h index e15fdcb..71b1598 100644 --- a/src/dks/model.h +++ b/src/dks/model.h @@ -1,4 +1,5 @@ -#pragma once +#ifndef DKS_MODEL_H_ +#define DKS_MODEL_H_ #include "tree.h" #include #include @@ -68,3 +69,4 @@ class model_t { std::vector _frequencies; }; } // namespace dks +#endif diff --git a/src/dks/partition.h b/src/dks/partition.h index 472f9bb..a663a2f 100644 --- a/src/dks/partition.h +++ b/src/dks/partition.h @@ -1,4 +1,5 @@ -#pragma once +#ifndef DKS_PARTITION_H_ +#define DKS_PARTITION_H_ #include "model.h" #include #include @@ -52,3 +53,4 @@ class partition_t { constexpr static double _rate_cats[] = {1.0}; }; } // namespace dks +#endif diff --git a/src/dks/test_case.h b/src/dks/test_case.h index 179ba8c..9ed66ac 100644 --- a/src/dks/test_case.h +++ b/src/dks/test_case.h @@ -1,4 +1,5 @@ -#pragma once +#ifndef DKS_TESTCASE_H_ +#define DKS_TESTCASE_H_ #include "model.h" #include "partition.h" #include "tree.h" @@ -171,3 +172,4 @@ template <> struct hash { std::ostream &operator<<(std::ostream &stream, const dks::test_cpu_t &cpu); std::ostream &operator<<(std::ostream &stream, const dks::attributes_t &attribs); +#endif diff --git a/src/dks/tree.h b/src/dks/tree.h index 39c534a..d2e194f 100644 --- a/src/dks/tree.h +++ b/src/dks/tree.h @@ -1,4 +1,5 @@ -#pragma once +#ifndef DKS_TREE_H_ +#define DKS_TREE_H_ #include "pll.h" #include #include @@ -37,3 +38,4 @@ class tree_t { std::vector _matrix_indices; }; } // namespace dks +#endif