Dks integration

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}")

src/dks/CMakeLists.txt

@@ -0,0 +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}
+  ${CMAKE_CURRENT_SOURCE_DIR}/../msa
+  )

src/dks/benchmark.cpp

@@ -0,0 +1,162 @@
+#include "benchmark.h"
+#include "partition.h"
+#include <algorithm>
+#include <cmath>
+#include <unordered_map>
+#include <utility>
+
+namespace dks {
+
+msa_t convert_pll_msa_t(const pll_msa_t *pll_msa) {
+  msa_t msa;
+  msa.reserve(pll_msa->count);
+  for (int i = 0; i < pll_msa->count; i++) {
+    msa.emplace_back(pll_msa->length);
+    for (int 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] +
+         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(double sites, double states, double taxa) {
+  kernel_weight_t kw;
+
+  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;
+}
+
+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,
+                                int taxa) {
+  return suggest_weights(sites, taxa, msa->states);
+}
+
+attributes_time_t
+select_kernel_verbose(const model_t &model,
+                      const std::vector<std::vector<double>> &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,
+                                        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, charmap);
+    times[attribs] = weight_kernel_times(kw, tc.benchmark(msa, weights, model));
+  }
+  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 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();
+}
+
+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_auto(msa, weights, charmap, states, rate_cats, gen);
+}
+
+unsigned int select_kernel_auto(const std::vector<std::vector<double>> &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<std::vector<double>> &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

src/dks/benchmark.h

@@ -0,0 +1,69 @@
+#ifndef DKS_BENCHMARK_H_
+#define DKS_BENCHMARK_H_
+#include "pll_msa.h"
+#include "test_case.h"
+#include <unordered_map>
+
+namespace dks {
+typedef std::unordered_map<test_kernel_t, double> kernel_weight_t;
+typedef std::unordered_map<attributes_t, benchmark_time_t> attributes_time_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 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<std::vector<double>> &clvs,
+                      const msa_weight_t &weights, const kernel_weight_t &kw,
+                      attributes_generator_t att_gen);
+
+attributes_t select_kernel_auto(const pll_partition_t *pll_partition,
+                                const pll_msa_t *pll_msa);
+
+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);
+
+unsigned int select_kernel_auto(const std::vector<std::vector<double>> &clvs,
+                                const msa_weight_t &weights,
+                                unsigned int states, unsigned int rate_cats,
+                                attributes_generator_t gen);
+
+unsigned int select_kernel_auto(const std::vector<std::vector<double>> &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
+#endif// DKS_BENCHMARK_H_

src/dks/dks.h

@@ -0,0 +1,6 @@
+#ifndef DKS_H_
+#define DKS_H_
+
+#include "benchmark.h"
+
+#endif

src/dks/model.cpp

@@ -0,0 +1,37 @@
+#include "model.h"
+#include <pll.h>
+
+namespace dks {
+
+unsigned int model_t::submodels() const { return 1; }
+
+unsigned int model_t::rate_categories() const { return _rate_categories; }
+
+uint64_t model_t::states() const { return _states; }
+
+const double *model_t::subst_params_raw() const { return _subst_params.data(); }
+
+const std::vector<double> &model_t::subst_params() const {
+  return _subst_params;
+}
+
+const double *model_t::frequencies_raw() const { return _frequencies.data(); }
+
+const std::vector<double> &model_t::frequencies() const { return _frequencies; }
+
+const tree_t &model_t::tree() const { return _tree; }
+
+std::vector<pll_operation_t> model_t::make_operations() const {
+  auto traversal_nodes = _tree.full_traverse();
+  std::vector<pll_operation_t> 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

src/dks/model.h

@@ -0,0 +1,72 @@
+#ifndef DKS_MODEL_H_
+#define DKS_MODEL_H_
+#include "tree.h"
+#include <memory>
+#include <pll.h>
+#include <vector>
+
+namespace dks {
+typedef std::vector<std::vector<char>> msa_t;
+typedef std::vector<unsigned int> msa_weight_t;
+class model_t {
+public:
+  model_t(const msa_t &msa, unsigned int states) : model_t{msa, states, 0} {};
+
+  model_t(size_t tip_count, unsigned int states)
+      : model_t{tip_count, states, 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, 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, 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],
+                      subst_params[model_index] +
+                          (_states - 1) * (_states - 2)},
+        _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;
+  const double *subst_params_raw() const;
+  const std::vector<double> &subst_params() const;
+  const double *frequencies_raw() const;
+  const std::vector<double> &frequencies() const;
+  const tree_t &tree() const;
+
+  void reset_tree();
+
+  std::vector<pll_operation_t> make_operations() const;
+
+private:
+  tree_t _tree;
+  unsigned int _states;
+  unsigned int _rate_categories;
+  double _prop_invar;
+  std::vector<double> _subst_params;
+  std::vector<double> _frequencies;
+};
+} // namespace dks
+#endif