From f6ad0319cc3c1593a44ba0620a2e5030d6044b0e Mon Sep 17 00:00:00 2001 From: molpopogen Date: Tue, 6 Jul 2021 14:48:23 -0700 Subject: [PATCH 01/19] Add types/tree_sequence.hpp --- fwdpp/ts/types/Makefile.am | 3 ++- fwdpp/ts/types/tree_sequence.hpp | 30 ++++++++++++++++++++++++++++++ 2 files changed, 32 insertions(+), 1 deletion(-) create mode 100644 fwdpp/ts/types/tree_sequence.hpp diff --git a/fwdpp/ts/types/Makefile.am b/fwdpp/ts/types/Makefile.am index d3b1752c2..df7fed654 100644 --- a/fwdpp/ts/types/Makefile.am +++ b/fwdpp/ts/types/Makefile.am @@ -5,4 +5,5 @@ pkginclude_HEADERS= node.hpp \ site.hpp \ mutation_record.hpp \ generate_null_id.hpp \ - table_collection.hpp + table_collection.hpp \ + tree_sequence.hpp diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp new file mode 100644 index 000000000..cf57fa1fd --- /dev/null +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -0,0 +1,30 @@ +#pragma once + +#include +#include +#include +#include "table_collection.hpp" + +namespace fwdpp +{ + namespace ts + { + namespace types + { + template class tree_sequence + { + private: + std::shared_ptr> tables_; + std::vector samples_; + std::size_t num_trees_; + + public: + tree_sequence( + std::shared_ptr> tables) + : tables_{std::move(tables)} + { + } + }; + } + } +} From 1eac491ff98eb4acce980ad8d383ae77785e8a2f Mon Sep 17 00:00:00 2001 From: molpopogen Date: Wed, 7 Jul 2021 09:29:47 -0700 Subject: [PATCH 02/19] Add initialization helper functions --- fwdpp/ts/types/tree_sequence.hpp | 69 +++++++++++++++++++++++++++++++- 1 file changed, 68 insertions(+), 1 deletion(-) diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index cf57fa1fd..6ff1d6c81 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -4,6 +4,7 @@ #include #include #include "table_collection.hpp" +#include "../node_flags.hpp" namespace fwdpp { @@ -18,13 +19,79 @@ namespace fwdpp std::vector samples_; std::size_t num_trees_; + static constexpr SignedInteger null + = types::table_collection::null; + + std::shared_ptr> + init_tables( + std::shared_ptr> tables) + { + if (tables == nullptr) + { + throw std::invalid_argument( + "input pointer to table_collection is nullptr"); + } + return tables; + } + + std::vector + init_samples(const types::table_collection& tables) + { + std::vector rv; + for (std::size_t i = 0; i < tables.nodes.size(); ++i) + { + if (tables.nodes[i].flags & node_flags::IS_SAMPLE) + { + rv.push_back(static_cast(i)); + } + } + return rv; + } + + std::vector + init_samples(const types::table_collection& tables, + std::vector samples) + { + std::vector is_sample(tables.nodes.size(), 0); + for (auto s : samples) + { + if (s == null) + { + throw std::invalid_argument( + "input sample ID is null"); + } + if (is_sample[static_cast(s)] == 1) + { + throw std::invalid_argument( + "node ID marked as sample multiple times"); + } + is_sample[static_cast(s)] = 1; + } + return samples; + } + public: tree_sequence( std::shared_ptr> tables) - : tables_{std::move(tables)} + : tables_{init_tables(std::move(tables))}, samples_{ + init_samples(*tables)} + { + } + + tree_sequence( + std::shared_ptr> tables, + std::vector samples) + : tables_{init_tables(std::move(tables))}, samples_{init_samples( + *tables, + std::move(samples))} { } }; + +#if __cplusplus < 201703L + template + constexpr SignedInteger tree_sequence::null; +#endif } } } From f0ac0da8a7cfcf1404a3605ae281e93b1e9c2918 Mon Sep 17 00:00:00 2001 From: molpopogen Date: Wed, 7 Jul 2021 09:36:06 -0700 Subject: [PATCH 03/19] Add test file. --- testsuite/Makefile.am | 3 ++- testsuite/tree_sequences/test_tree_sequence.cc | 12 ++++++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) create mode 100644 testsuite/tree_sequences/test_tree_sequence.cc diff --git a/testsuite/Makefile.am b/testsuite/Makefile.am index 8c042897d..3420a072e 100644 --- a/testsuite/Makefile.am +++ b/testsuite/Makefile.am @@ -96,7 +96,8 @@ tree_sequences_tree_sequence_tests_SOURCES=tree_sequences/tree_sequence_tests.cc tree_sequences/test_diploid_recording.cc \ tree_sequences/wfevolve_table_collection_fxns.cc \ tree_sequences/test_edge_buffering_std_table_collection.cc \ - tree_sequences/tskit_utils.cc + tree_sequences/tskit_utils.cc \ + tree_sequences/test_tree_sequence.cc tree_sequences_tree_sequence_tests_CFLAGS=-std=c99 diff --git a/testsuite/tree_sequences/test_tree_sequence.cc b/testsuite/tree_sequences/test_tree_sequence.cc new file mode 100644 index 000000000..0af150d81 --- /dev/null +++ b/testsuite/tree_sequences/test_tree_sequence.cc @@ -0,0 +1,12 @@ +#include +#include +#include +#include + +using tree_sequence = fwdpp::ts::types::tree_sequence; +using tree_sequence64 = fwdpp::ts::types::tree_sequence; + +//explicit instantiations +template class fwdpp::ts::types::tree_sequence; +template class fwdpp::ts::types::tree_sequence; + From 582480fabd28651d60d548e0bdab479c72c9bfa0 Mon Sep 17 00:00:00 2001 From: molpopogen Date: Wed, 7 Jul 2021 11:25:10 -0700 Subject: [PATCH 04/19] add init_num_trees --- fwdpp/ts/types/tree_sequence.hpp | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index 6ff1d6c81..1fac469e8 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -70,11 +70,17 @@ namespace fwdpp return samples; } + std::size_t + init_num_trees() + { + return 0; + } + public: tree_sequence( std::shared_ptr> tables) - : tables_{init_tables(std::move(tables))}, samples_{ - init_samples(*tables)} + : tables_{init_tables(std::move(tables))}, + samples_{init_samples(*tables)}, num_trees_{init_num_trees()} { } @@ -83,7 +89,8 @@ namespace fwdpp std::vector samples) : tables_{init_tables(std::move(tables))}, samples_{init_samples( *tables, - std::move(samples))} + std::move(samples))}, + num_trees_{init_num_trees()} { } }; From a7baf2af71ed73a533f62977d2a5bf24659047e2 Mon Sep 17 00:00:00 2001 From: molpopogen Date: Wed, 7 Jul 2021 12:58:37 -0700 Subject: [PATCH 05/19] Add ts/tree_sequence_flags.hpp --- fwdpp/ts/Makefile.am | 3 ++- fwdpp/ts/tree_sequence_flags.hpp | 15 +++++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) create mode 100644 fwdpp/ts/tree_sequence_flags.hpp diff --git a/fwdpp/ts/Makefile.am b/fwdpp/ts/Makefile.am index 736d42e7e..255f7c29b 100644 --- a/fwdpp/ts/Makefile.am +++ b/fwdpp/ts/Makefile.am @@ -34,6 +34,7 @@ pkginclude_HEADERS= definitions.hpp \ mutation_tools.hpp \ visit_sites.hpp \ site_visitor.hpp \ - recording.hpp + recording.hpp \ + tree_sequence_flags.hpp diff --git a/fwdpp/ts/tree_sequence_flags.hpp b/fwdpp/ts/tree_sequence_flags.hpp new file mode 100644 index 000000000..3bb1dc87f --- /dev/null +++ b/fwdpp/ts/tree_sequence_flags.hpp @@ -0,0 +1,15 @@ +#pragma once + +#include + +namespace fwdpp +{ + namespace ts + { + namespace tree_sequence_flags + { + static constexpr std::uint32_t NONE = 0; + static constexpr std::uint32_t TRACK_SAMPLES = 1 << 0; + } + } +} From d16a7a492af589ce2a79d71ae9f9c277514cca13 Mon Sep 17 00:00:00 2001 From: molpopogen Date: Wed, 7 Jul 2021 13:08:18 -0700 Subject: [PATCH 06/19] add support for treeseq flags --- fwdpp/ts/types/tree_sequence.hpp | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index 1fac469e8..21de1a8c8 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -4,7 +4,10 @@ #include #include #include "table_collection.hpp" +#include "../marginal_tree.hpp" +#include "../exceptions.hpp" #include "../node_flags.hpp" +#include "../tree_sequence_flags.hpp" namespace fwdpp { @@ -18,6 +21,8 @@ namespace fwdpp std::shared_ptr> tables_; std::vector samples_; std::size_t num_trees_; + std::uint32_t flags_; + marginal_tree tree_; static constexpr SignedInteger null = types::table_collection::null; @@ -31,6 +36,10 @@ namespace fwdpp throw std::invalid_argument( "input pointer to table_collection is nullptr"); } + if (!tables->indexed()) + { + throw tables_error("tables are not indexed"); + } return tables; } @@ -78,19 +87,27 @@ namespace fwdpp public: tree_sequence( - std::shared_ptr> tables) - : tables_{init_tables(std::move(tables))}, - samples_{init_samples(*tables)}, num_trees_{init_num_trees()} + std::shared_ptr> tables, + std::uint32_t flags) + : tables_{init_tables(std::move(tables))}, samples_{init_samples( + *tables)}, + num_trees_{init_num_trees()}, flags_{flags}, + tree_{ + tables->nodes.size(), samples_, + static_cast(flags_ & tree_sequence_flags::TRACK_SAMPLES)} { } tree_sequence( std::shared_ptr> tables, - std::vector samples) + std::vector samples, std::uint32_t flags) : tables_{init_tables(std::move(tables))}, samples_{init_samples( *tables, std::move(samples))}, - num_trees_{init_num_trees()} + num_trees_{init_num_trees()}, flags_{flags}, + tree_{ + tables->nodes.size(), samples_, + static_cast(flags_ & tree_sequence_flags::TRACK_SAMPLES)} { } }; From 4ed6aef768cfcce53f37097ae271f03c06d37a37 Mon Sep 17 00:00:00 2001 From: molpopogen Date: Wed, 7 Jul 2021 13:14:55 -0700 Subject: [PATCH 07/19] uncouple tree sequence from marginal_tree --- fwdpp/ts/Makefile.am | 2 +- .../ts/{tree_sequence_flags.hpp => tree_flags.hpp} | 2 +- fwdpp/ts/types/tree_sequence.hpp | 13 +++---------- 3 files changed, 5 insertions(+), 12 deletions(-) rename fwdpp/ts/{tree_sequence_flags.hpp => tree_flags.hpp} (85%) diff --git a/fwdpp/ts/Makefile.am b/fwdpp/ts/Makefile.am index 255f7c29b..6059418dc 100644 --- a/fwdpp/ts/Makefile.am +++ b/fwdpp/ts/Makefile.am @@ -35,6 +35,6 @@ pkginclude_HEADERS= definitions.hpp \ visit_sites.hpp \ site_visitor.hpp \ recording.hpp \ - tree_sequence_flags.hpp + tree_flags.hpp diff --git a/fwdpp/ts/tree_sequence_flags.hpp b/fwdpp/ts/tree_flags.hpp similarity index 85% rename from fwdpp/ts/tree_sequence_flags.hpp rename to fwdpp/ts/tree_flags.hpp index 3bb1dc87f..04f1ff266 100644 --- a/fwdpp/ts/tree_sequence_flags.hpp +++ b/fwdpp/ts/tree_flags.hpp @@ -6,7 +6,7 @@ namespace fwdpp { namespace ts { - namespace tree_sequence_flags + namespace tree_flags { static constexpr std::uint32_t NONE = 0; static constexpr std::uint32_t TRACK_SAMPLES = 1 << 0; diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index 21de1a8c8..3ba511c0a 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -7,7 +7,7 @@ #include "../marginal_tree.hpp" #include "../exceptions.hpp" #include "../node_flags.hpp" -#include "../tree_sequence_flags.hpp" +#include "../tree_flags.hpp" namespace fwdpp { @@ -22,7 +22,6 @@ namespace fwdpp std::vector samples_; std::size_t num_trees_; std::uint32_t flags_; - marginal_tree tree_; static constexpr SignedInteger null = types::table_collection::null; @@ -91,10 +90,7 @@ namespace fwdpp std::uint32_t flags) : tables_{init_tables(std::move(tables))}, samples_{init_samples( *tables)}, - num_trees_{init_num_trees()}, flags_{flags}, - tree_{ - tables->nodes.size(), samples_, - static_cast(flags_ & tree_sequence_flags::TRACK_SAMPLES)} + num_trees_{init_num_trees()}, flags_{flags} { } @@ -104,10 +100,7 @@ namespace fwdpp : tables_{init_tables(std::move(tables))}, samples_{init_samples( *tables, std::move(samples))}, - num_trees_{init_num_trees()}, flags_{flags}, - tree_{ - tables->nodes.size(), samples_, - static_cast(flags_ & tree_sequence_flags::TRACK_SAMPLES)} + num_trees_{init_num_trees()}, flags_{flags} { } }; From 3e57311e95410083ba7d51d27fd8aa9bdeabc577 Mon Sep 17 00:00:00 2001 From: molpopogen Date: Wed, 7 Jul 2021 13:18:56 -0700 Subject: [PATCH 08/19] remove flags from treeseq constructor --- fwdpp/ts/types/tree_sequence.hpp | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index 3ba511c0a..f332bb99a 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -21,7 +21,6 @@ namespace fwdpp std::shared_ptr> tables_; std::vector samples_; std::size_t num_trees_; - std::uint32_t flags_; static constexpr SignedInteger null = types::table_collection::null; @@ -86,21 +85,19 @@ namespace fwdpp public: tree_sequence( - std::shared_ptr> tables, - std::uint32_t flags) - : tables_{init_tables(std::move(tables))}, samples_{init_samples( - *tables)}, - num_trees_{init_num_trees()}, flags_{flags} + std::shared_ptr> tables) + : tables_{init_tables(std::move(tables))}, + samples_{init_samples(*tables)}, num_trees_{init_num_trees()} { } tree_sequence( std::shared_ptr> tables, - std::vector samples, std::uint32_t flags) + std::vector samples) : tables_{init_tables(std::move(tables))}, samples_{init_samples( *tables, std::move(samples))}, - num_trees_{init_num_trees()}, flags_{flags} + num_trees_{init_num_trees()} { } }; From 29ab1f274aa9bf5e9805518b88f044cc492e9e97 Mon Sep 17 00:00:00 2001 From: molpopogen Date: Wed, 7 Jul 2021 13:34:13 -0700 Subject: [PATCH 09/19] start thinking about tree iterators --- fwdpp/ts/types/tree_sequence.hpp | 39 ++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index f332bb99a..0c41c14b3 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -3,6 +3,9 @@ #include #include #include +#if __cplusplus >= 201703L +#include +#endif #include "table_collection.hpp" #include "../marginal_tree.hpp" #include "../exceptions.hpp" @@ -15,6 +18,42 @@ namespace fwdpp { namespace types { + template class tree_sequence; + + template class tree_iterator + { + private: + const tree_sequence* treeseq_; + marginal_tree tree_; + + public: + tree_iterator(const tree_sequence* ts, + const std::vector& samples, + std::uint32_t flags) + : treeseq_{ts}, tree_{ts->nodes.size(), samples, + static_cast(flags + & tree_flags::TRACK_SAMPLES)} + { + } + + const marginal_tree* + tree_ptr() const + { + // FIXME: should be nullptr if no more trees + return &tree_; + } + +#if __cplusplus >= 201703L + std::optional&> + tree() const + { + // FIXME: should be nullopt if no more trees + return std::optional&>{ + std::cref(tree_)}; + } +#endif + }; + template class tree_sequence { private: From cf96361cbfd779c6d01a11473d52f111296f2959 Mon Sep 17 00:00:00 2001 From: molpopogen Date: Wed, 7 Jul 2021 13:50:43 -0700 Subject: [PATCH 10/19] clean up ownership semantics --- fwdpp/ts/types/tree_sequence.hpp | 25 ++++++++++++++++--- .../tree_sequences/test_tree_sequence.cc | 2 ++ 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index 0c41c14b3..475fc9548 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -23,14 +23,13 @@ namespace fwdpp template class tree_iterator { private: - const tree_sequence* treeseq_; + std::reference_wrapper> treeseq_; marginal_tree tree_; public: - tree_iterator(const tree_sequence* ts, - const std::vector& samples, + tree_iterator(const tree_sequence& ts, std::uint32_t flags) - : treeseq_{ts}, tree_{ts->nodes.size(), samples, + : treeseq_{ts}, tree_{ts.num_nodes(), ts.samples(), static_cast(flags & tree_flags::TRACK_SAMPLES)} { @@ -139,6 +138,24 @@ namespace fwdpp num_trees_{init_num_trees()} { } + + std::size_t + num_nodes() const + { + return tables_->nodes.size(); + } + + const std::vector& + samples() const + { + return samples_; + } + + tree_iterator + trees(std::uint32_t flags) + { + return tree_iterator{*this, flags}; + } }; #if __cplusplus < 201703L diff --git a/testsuite/tree_sequences/test_tree_sequence.cc b/testsuite/tree_sequences/test_tree_sequence.cc index 0af150d81..bdb25e00b 100644 --- a/testsuite/tree_sequences/test_tree_sequence.cc +++ b/testsuite/tree_sequences/test_tree_sequence.cc @@ -9,4 +9,6 @@ using tree_sequence64 = fwdpp::ts::types::tree_sequence; //explicit instantiations template class fwdpp::ts::types::tree_sequence; template class fwdpp::ts::types::tree_sequence; +template class fwdpp::ts::types::tree_iterator; +template class fwdpp::ts::types::tree_iterator; From 38876049c9f1164d2dab5b33d3932eee304dd939 Mon Sep 17 00:00:00 2001 From: molpopogen Date: Wed, 7 Jul 2021 13:57:19 -0700 Subject: [PATCH 11/19] optional should wrap a reference_wrapper --- fwdpp/ts/types/tree_sequence.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index 475fc9548..6eeebca71 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -43,11 +43,12 @@ namespace fwdpp } #if __cplusplus >= 201703L - std::optional&> + std::optional>> tree() const { // FIXME: should be nullopt if no more trees - return std::optional&>{ + return std::optional< + std::reference_wrapper>>{ std::cref(tree_)}; } #endif From 8c4eb1aabde6f98efd49e19c3893dc8c788a74ec Mon Sep 17 00:00:00 2001 From: molpopogen Date: Wed, 7 Jul 2021 14:37:05 -0700 Subject: [PATCH 12/19] more api tweaks --- fwdpp/ts/types/tree_sequence.hpp | 40 +++++++++++++++++++++++++------- 1 file changed, 32 insertions(+), 8 deletions(-) diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index 6eeebca71..913be9992 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -25,33 +25,51 @@ namespace fwdpp private: std::reference_wrapper> treeseq_; marginal_tree tree_; + std::size_t input_edge_index, output_edge_index; + double position; + + bool + still_advancing() const + // assumes left-to-right iteration + { + return (input_edge_index < tables().input_left.size() + || position < tables().genome_length()); + } public: tree_iterator(const tree_sequence& ts, std::uint32_t flags) : treeseq_{ts}, tree_{ts.num_nodes(), ts.samples(), - static_cast(flags - & tree_flags::TRACK_SAMPLES)} + static_cast( + flags & tree_flags::TRACK_SAMPLES)}, + input_edge_index{0}, output_edge_index{0}, position{0.} { } const marginal_tree* tree_ptr() const { - // FIXME: should be nullptr if no more trees - return &tree_; + return still_advancing() ? &tree_ : nullptr; } #if __cplusplus >= 201703L std::optional>> tree() const { - // FIXME: should be nullopt if no more trees - return std::optional< - std::reference_wrapper>>{ - std::cref(tree_)}; + if (still_advancing()) + { + return std::optional>>{std::cref(tree_)}; + } + + return std::nullopt; } #endif + const types::table_collection& + tables() const + { + return treeseq_.get().tables(); + } }; template class tree_sequence @@ -157,6 +175,12 @@ namespace fwdpp { return tree_iterator{*this, flags}; } + + const types::table_collection& + tables() const + { + return *tables_; + } }; #if __cplusplus < 201703L From 447557a62170daaa1871c663f851f631176f2de7 Mon Sep 17 00:00:00 2001 From: molpopogen Date: Wed, 7 Jul 2021 14:42:48 -0700 Subject: [PATCH 13/19] add empty advance fn --- fwdpp/ts/types/tree_sequence.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index 913be9992..12c436e94 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -181,6 +181,12 @@ namespace fwdpp { return *tables_; } + + void + advance_right() + // Move 1 tree left-to-right + { + } }; #if __cplusplus < 201703L From 67001f43205c4a393dffc385ff37928704141f03 Mon Sep 17 00:00:00 2001 From: molpopogen Date: Thu, 8 Jul 2021 06:55:30 -0700 Subject: [PATCH 14/19] interface streamline/fix --- fwdpp/ts/types/tree_sequence.hpp | 29 ++++++++++++----------------- 1 file changed, 12 insertions(+), 17 deletions(-) diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index 12c436e94..02b257946 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -27,14 +27,7 @@ namespace fwdpp marginal_tree tree_; std::size_t input_edge_index, output_edge_index; double position; - - bool - still_advancing() const - // assumes left-to-right iteration - { - return (input_edge_index < tables().input_left.size() - || position < tables().genome_length()); - } + bool advanced_; public: tree_iterator(const tree_sequence& ts, @@ -42,21 +35,23 @@ namespace fwdpp : treeseq_{ts}, tree_{ts.num_nodes(), ts.samples(), static_cast( flags & tree_flags::TRACK_SAMPLES)}, - input_edge_index{0}, output_edge_index{0}, position{0.} + input_edge_index{0}, output_edge_index{0}, position{0.}, advanced_{ + false} { } const marginal_tree* tree_ptr() const + // tree() should be preferred for c++ >= 17 { - return still_advancing() ? &tree_ : nullptr; + return advanced_ ? &tree_ : nullptr; } #if __cplusplus >= 201703L std::optional>> tree() const { - if (still_advancing()) + if (advanced_) { return std::optional>>{std::cref(tree_)}; @@ -70,6 +65,12 @@ namespace fwdpp { return treeseq_.get().tables(); } + + void + advance_right() + // Move 1 tree left-to-right + { + } }; template class tree_sequence @@ -181,12 +182,6 @@ namespace fwdpp { return *tables_; } - - void - advance_right() - // Move 1 tree left-to-right - { - } }; #if __cplusplus < 201703L From fd5352add89c1ef7791a46148eeffe788a367fa9 Mon Sep 17 00:00:00 2001 From: molpopogen Date: Thu, 8 Jul 2021 08:11:14 -0700 Subject: [PATCH 15/19] NOTE to self --- fwdpp/ts/types/tree_sequence.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index 02b257946..8b2acb4cf 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -23,6 +23,9 @@ namespace fwdpp template class tree_iterator { private: + // NOTE: since tree stores the samples, + // we can probably just keep a shared_ptr to the tables here, + // adding a level of memory safety. std::reference_wrapper> treeseq_; marginal_tree tree_; std::size_t input_edge_index, output_edge_index; From 613d563c8f8426b4238a6b4ce42ca89d7026889f Mon Sep 17 00:00:00 2001 From: molpopogen Date: Thu, 8 Jul 2021 08:24:02 -0700 Subject: [PATCH 16/19] move marginal_tree definition into ts::types --- fwdpp/ts/marginal_tree.hpp | 367 ---------------------- fwdpp/ts/types/Makefile.am | 3 +- fwdpp/ts/types/marginal_tree.hpp | 366 +++++++++++++++++++++ testsuite/tree_sequences/preorder_adl.hpp | 7 +- 4 files changed, 372 insertions(+), 371 deletions(-) delete mode 100644 fwdpp/ts/marginal_tree.hpp create mode 100644 fwdpp/ts/types/marginal_tree.hpp diff --git a/fwdpp/ts/marginal_tree.hpp b/fwdpp/ts/marginal_tree.hpp deleted file mode 100644 index bbb89bf1c..000000000 --- a/fwdpp/ts/marginal_tree.hpp +++ /dev/null @@ -1,367 +0,0 @@ -#ifndef FWDPP_TS_MARGINAL_TREE_HPP -#define FWDPP_TS_MARGINAL_TREE_HPP - -#include -#include -#include -#include -#include -#include -#include -#include "exceptions.hpp" -#include "types/generate_null_id.hpp" - -namespace fwdpp -{ - namespace ts - { - template struct sample_group_map - /// \brief Maps a node id to a sample group - /// - /// When constructing a fwdpp::ts::tree_visitor, - /// vectors of this type may be used to mark - /// sample nodes as beloning to different "groups". - { - static_assert(std::is_integral::value, - "SignedInteger must be an integral type"); - static_assert(std::is_signed::value, - "SignedInteger must be a signed type"); - /// node id - SignedInteger node_id; - /// group label - std::int32_t group; - sample_group_map(SignedInteger n, std::int32_t g) : node_id(n), group(g) - /// \param n A node id - /// \param g A gropup label - { - if (n == types::generate_null_id()) - { - throw samples_error("null ID passed to sample_index_map"); - } - } - }; - - template class marginal_tree - /// \brief A non-recombining tree - /// - /// The tree is represented as a sparse tree - /// data structure where several linked lists - /// allow efficient retrieval of parent/child/siblings - /// of nodes. - /// - /// \version 0.7.0 Added to fwdpp - /// \version 0.7.1 Constructors throw exceptions when sample lists contain the - /// same node ID more than once. Changed from struct to class in order - /// to reuse some code in a private function. Initialization also - /// tracks the total sample size, which is number of nonzero elements - /// in sample_index_map. - /// \version 0.7.4 Update to include data structures for root tracking - /// \version 0.8.0 Now holds a list of samples. Samples may be assigned to groups. - { - private: - static_assert(std::is_integral::value, - "SignedInteger must be an integral type"); - static_assert(std::is_signed::value, - "SignedInteger must be a signed type"); - std::size_t num_nodes; - std::vector sample_groups; - std::vector samples_list; - bool advancing_sample_list_; - - static constexpr SignedInteger null - = types::generate_null_id(); - - std::vector - fill_sample_groups(const std::vector& samples) - { - std::vector rv(num_nodes, - std::numeric_limits::min()); - for (auto i : samples) - { - rv[i] = 0; - } - return rv; - } - - std::vector - fill_sample_groups( - const std::vector>& samples) - { - std::vector rv(num_nodes, - std::numeric_limits::min()); - for (auto i : samples) - { - rv[i.node_id] = i.group; - } - return rv; - } - - std::vector - fill_sample_groups(const std::vector& samples_a, - const std::vector& samples_b) - { - std::vector rv(num_nodes, - std::numeric_limits::min()); - for (auto i : samples_a) - { - rv[i] = 0; - } - for (auto i : samples_b) - { - rv[i] = 1; - } - return rv; - } - - std::vector - fill_sample_groups( - const std::vector>& samples_a, - const std::vector& samples_b) - { - std::vector rv(num_nodes, - std::numeric_limits::min()); - for (auto i : samples_a) - { - rv[i.node_id] = 0; - } - for (auto i : samples_b) - { - rv[i] = 1; - } - return rv; - } - - std::vector - init_samples_list(const std::vector& s) - { - return s; - } - - std::vector - init_samples_list(const std::vector>& s) - { - std::vector rv; - for (auto& i : s) - { - rv.push_back(i.node_id); - } - return rv; - } - - std::vector - init_samples_list(const std::vector& a, - const std::vector& b) - { - auto rv = a; - rv.insert(end(rv), begin(b), end(b)); - return rv; - } - - std::vector - init_samples_list(const std::vector>& a, - const std::vector& b) - { - auto rv = init_samples_list(a); - rv.insert(end(rv), begin(b), end(b)); - return rv; - } - - void - init_samples() - { - for (std::size_t i = 0; i < samples_list.size(); ++i) - { - auto s = samples_list[i]; - // See GitHub issue #158 for background - if (sample_index_map[s] != null) - { - throw samples_error("invalid sample list"); - } - sample_index_map[s] = static_cast(i); - left_sample[s] = right_sample[s] = sample_index_map[s]; - above_sample[s] = 1; - // Initialize roots - if (i < samples_list.size() - 1) - { - right_sib[s] = samples_list[i + 1]; - } - if (i > 0) - { - left_sib[s] = samples_list[i - 1]; - } - } - } - - public: - using id_type = SignedInteger; - std::vector parents, leaf_counts, preserved_leaf_counts, left_sib, - right_sib, left_child, right_child, left_sample, right_sample, - next_sample, sample_index_map; - std::vector above_sample; - double left, right; - id_type left_root; - - template - marginal_tree(std::size_t nnodes, const SAMPLES& samples, - bool advancing_sample_list) - : num_nodes(nnodes), sample_groups(fill_sample_groups(samples)), - samples_list(init_samples_list(samples)), - advancing_sample_list_(advancing_sample_list), parents(nnodes, null), - leaf_counts(nnodes, 0), preserved_leaf_counts(nnodes, 0), - left_sib(nnodes, null), right_sib(nnodes, null), - left_child(nnodes, null), right_child(nnodes, null), - left_sample(nnodes, null), right_sample(nnodes, null), - next_sample(nnodes, null), sample_index_map(nnodes, null), - above_sample(nnodes, 0), - left{std::numeric_limits::quiet_NaN()}, - right{std::numeric_limits::quiet_NaN()}, left_root(null) - /// \param nnodes Number of nodes in table_collection - /// \param samples The sample list - /// - /// \a samples may be either std::vector or - /// std::vector. For the former, all - /// sample nodes will be assigned group 0. - { - if (samples_list.empty()) - { - throw samples_error("marginal_tree: empty sample list"); - } - init_samples(); - for (auto s : samples_list) - { - leaf_counts[s] = 1; - } - left_root = samples_list[0]; - } - - template - marginal_tree(std::size_t nnodes, const SAMPLES& samples, - const std::vector preserved_nodes, - bool advancing_sample_list) - : num_nodes(nnodes), - sample_groups(fill_sample_groups(samples, preserved_nodes)), - samples_list(init_samples_list(samples, preserved_nodes)), - advancing_sample_list_(advancing_sample_list), parents(nnodes, null), - leaf_counts(nnodes, 0), preserved_leaf_counts(nnodes, 0), - left_sib(nnodes, null), right_sib(nnodes, null), - left_child(nnodes, null), right_child(nnodes, null), - left_sample(nnodes, null), right_sample(nnodes, null), - next_sample(nnodes, null), sample_index_map(nnodes, null), - above_sample(nnodes, 0), - left{std::numeric_limits::quiet_NaN()}, - right{std::numeric_limits::quiet_NaN()}, left_root(null) - /// Constructor - { - if (samples_list.empty()) - { - throw samples_error("marginal_tree: empty sample list"); - } - init_samples(); - left_root = samples_list[0]; - for (std::size_t i = 0; i < samples.size(); ++i) - { - leaf_counts[samples_list[i]] = 1; - } - for (auto s : preserved_nodes) - { - preserved_leaf_counts[s] = 1; - } - } - - marginal_tree(id_type nnodes) - : num_nodes(nnodes), sample_groups{}, samples_list{}, - advancing_sample_list_(false), - parents(nnodes, null), leaf_counts{}, preserved_leaf_counts{}, - left_sib(nnodes, null), right_sib(nnodes, null), - left_child(nnodes, null), right_child(nnodes, null), - left_sample(nnodes, null), right_sample(nnodes, null), - next_sample(nnodes, null), sample_index_map(nnodes, null), - above_sample(nnodes, 0), - left{std::numeric_limits::quiet_NaN()}, - right{std::numeric_limits::quiet_NaN()}, left_root(null) - /// Constructor - /// \todo Document - { - } - - int - num_roots() const - /// Return number of roots - { - if (left_root == null) - { - throw std::runtime_error("left_root is NULL"); - } - int nroots = 0; - auto lr = left_root; - while (lr != null) - { - ++nroots; - lr = right_sib[lr]; - } - return nroots; - } - - inline std::size_t - sample_size() const - /// Number of samples - { - return samples_list.size(); - } - - inline typename std::vector::const_iterator - samples_list_begin() const - /// Beginning of samples list - { - return begin(samples_list); - } - - inline typename std::vector::const_iterator - samples_list_end() const - /// End itertor for samples list - { - return end(samples_list); - } - - inline id_type - sample_group(id_type u) const - /// Return the sample group for node \u - { - if (static_cast(u) >= num_nodes) - { - throw std::invalid_argument("invalid node id"); - } - return sample_groups[u]; - } - - inline bool - advancing_sample_list() const - { - return advancing_sample_list_; - } - - inline std::size_t - size() const - /// Return the length of the internal vectors. - { - return num_nodes; - } - - inline id_type - sample_table_index_to_node(id_type u) const - /// If u is a sample index, return the associated - /// node id. - { - return samples_list[u]; - } - }; - -#if __cplusplus < 201703L - template - constexpr SignedInteger marginal_tree::null; -#endif - - } // namespace ts -} // namespace fwdpp - -#endif diff --git a/fwdpp/ts/types/Makefile.am b/fwdpp/ts/types/Makefile.am index df7fed654..a3f32c6b0 100644 --- a/fwdpp/ts/types/Makefile.am +++ b/fwdpp/ts/types/Makefile.am @@ -6,4 +6,5 @@ pkginclude_HEADERS= node.hpp \ mutation_record.hpp \ generate_null_id.hpp \ table_collection.hpp \ - tree_sequence.hpp + tree_sequence.hpp \ + marginal_tree.hpp diff --git a/fwdpp/ts/types/marginal_tree.hpp b/fwdpp/ts/types/marginal_tree.hpp new file mode 100644 index 000000000..39c30dcda --- /dev/null +++ b/fwdpp/ts/types/marginal_tree.hpp @@ -0,0 +1,366 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include "../exceptions.hpp" +#include "generate_null_id.hpp" + +namespace fwdpp +{ + namespace ts + { + namespace types + { + template struct sample_group_map + /// \brief Maps a node id to a sample group + /// + /// When constructing a fwdpp::ts::tree_visitor, + /// vectors of this type may be used to mark + /// sample nodes as beloning to different "groups". + { + static_assert(std::is_integral::value, + "SignedInteger must be an integral type"); + static_assert(std::is_signed::value, + "SignedInteger must be a signed type"); + /// node id + SignedInteger node_id; + /// group label + std::int32_t group; + sample_group_map(SignedInteger n, std::int32_t g) : node_id(n), group(g) + /// \param n A node id + /// \param g A gropup label + { + if (n == types::generate_null_id()) + { + throw samples_error("null ID passed to sample_index_map"); + } + } + }; + + template class marginal_tree + /// \brief A non-recombining tree + /// + /// The tree is represented as a sparse tree + /// data structure where several linked lists + /// allow efficient retrieval of parent/child/siblings + /// of nodes. + /// + /// \version 0.7.0 Added to fwdpp + /// \version 0.7.1 Constructors throw exceptions when sample lists contain the + /// same node ID more than once. Changed from struct to class in order + /// to reuse some code in a private function. Initialization also + /// tracks the total sample size, which is number of nonzero elements + /// in sample_index_map. + /// \version 0.7.4 Update to include data structures for root tracking + /// \version 0.8.0 Now holds a list of samples. Samples may be assigned to groups. + { + private: + static_assert(std::is_integral::value, + "SignedInteger must be an integral type"); + static_assert(std::is_signed::value, + "SignedInteger must be a signed type"); + std::size_t num_nodes; + std::vector sample_groups; + std::vector samples_list; + bool advancing_sample_list_; + + static constexpr SignedInteger null + = types::generate_null_id(); + + std::vector + fill_sample_groups(const std::vector& samples) + { + std::vector rv( + num_nodes, std::numeric_limits::min()); + for (auto i : samples) + { + rv[i] = 0; + } + return rv; + } + + std::vector + fill_sample_groups( + const std::vector>& samples) + { + std::vector rv( + num_nodes, std::numeric_limits::min()); + for (auto i : samples) + { + rv[i.node_id] = i.group; + } + return rv; + } + + std::vector + fill_sample_groups(const std::vector& samples_a, + const std::vector& samples_b) + { + std::vector rv( + num_nodes, std::numeric_limits::min()); + for (auto i : samples_a) + { + rv[i] = 0; + } + for (auto i : samples_b) + { + rv[i] = 1; + } + return rv; + } + + std::vector + fill_sample_groups( + const std::vector>& samples_a, + const std::vector& samples_b) + { + std::vector rv( + num_nodes, std::numeric_limits::min()); + for (auto i : samples_a) + { + rv[i.node_id] = 0; + } + for (auto i : samples_b) + { + rv[i] = 1; + } + return rv; + } + + std::vector + init_samples_list(const std::vector& s) + { + return s; + } + + std::vector + init_samples_list(const std::vector>& s) + { + std::vector rv; + for (auto& i : s) + { + rv.push_back(i.node_id); + } + return rv; + } + + std::vector + init_samples_list(const std::vector& a, + const std::vector& b) + { + auto rv = a; + rv.insert(end(rv), begin(b), end(b)); + return rv; + } + + std::vector + init_samples_list(const std::vector>& a, + const std::vector& b) + { + auto rv = init_samples_list(a); + rv.insert(end(rv), begin(b), end(b)); + return rv; + } + + void + init_samples() + { + for (std::size_t i = 0; i < samples_list.size(); ++i) + { + auto s = samples_list[i]; + // See GitHub issue #158 for background + if (sample_index_map[s] != null) + { + throw samples_error("invalid sample list"); + } + sample_index_map[s] = static_cast(i); + left_sample[s] = right_sample[s] = sample_index_map[s]; + above_sample[s] = 1; + // Initialize roots + if (i < samples_list.size() - 1) + { + right_sib[s] = samples_list[i + 1]; + } + if (i > 0) + { + left_sib[s] = samples_list[i - 1]; + } + } + } + + public: + using id_type = SignedInteger; + std::vector parents, leaf_counts, preserved_leaf_counts, + left_sib, right_sib, left_child, right_child, left_sample, + right_sample, next_sample, sample_index_map; + std::vector above_sample; + double left, right; + id_type left_root; + + template + marginal_tree(std::size_t nnodes, const SAMPLES& samples, + bool advancing_sample_list) + : num_nodes(nnodes), sample_groups(fill_sample_groups(samples)), + samples_list(init_samples_list(samples)), + advancing_sample_list_(advancing_sample_list), + parents(nnodes, null), leaf_counts(nnodes, 0), + preserved_leaf_counts(nnodes, 0), left_sib(nnodes, null), + right_sib(nnodes, null), left_child(nnodes, null), + right_child(nnodes, null), left_sample(nnodes, null), + right_sample(nnodes, null), next_sample(nnodes, null), + sample_index_map(nnodes, null), above_sample(nnodes, 0), + left{std::numeric_limits::quiet_NaN()}, + right{std::numeric_limits::quiet_NaN()}, left_root(null) + /// \param nnodes Number of nodes in table_collection + /// \param samples The sample list + /// + /// \a samples may be either std::vector or + /// std::vector. For the former, all + /// sample nodes will be assigned group 0. + { + if (samples_list.empty()) + { + throw samples_error("marginal_tree: empty sample list"); + } + init_samples(); + for (auto s : samples_list) + { + leaf_counts[s] = 1; + } + left_root = samples_list[0]; + } + + template + marginal_tree(std::size_t nnodes, const SAMPLES& samples, + const std::vector preserved_nodes, + bool advancing_sample_list) + : num_nodes(nnodes), + sample_groups(fill_sample_groups(samples, preserved_nodes)), + samples_list(init_samples_list(samples, preserved_nodes)), + advancing_sample_list_(advancing_sample_list), + parents(nnodes, null), leaf_counts(nnodes, 0), + preserved_leaf_counts(nnodes, 0), left_sib(nnodes, null), + right_sib(nnodes, null), left_child(nnodes, null), + right_child(nnodes, null), left_sample(nnodes, null), + right_sample(nnodes, null), next_sample(nnodes, null), + sample_index_map(nnodes, null), above_sample(nnodes, 0), + left{std::numeric_limits::quiet_NaN()}, + right{std::numeric_limits::quiet_NaN()}, left_root(null) + /// Constructor + { + if (samples_list.empty()) + { + throw samples_error("marginal_tree: empty sample list"); + } + init_samples(); + left_root = samples_list[0]; + for (std::size_t i = 0; i < samples.size(); ++i) + { + leaf_counts[samples_list[i]] = 1; + } + for (auto s : preserved_nodes) + { + preserved_leaf_counts[s] = 1; + } + } + + marginal_tree(id_type nnodes) + : num_nodes(nnodes), sample_groups{}, samples_list{}, + advancing_sample_list_(false), + parents(nnodes, null), leaf_counts{}, preserved_leaf_counts{}, + left_sib(nnodes, null), right_sib(nnodes, null), + left_child(nnodes, null), right_child(nnodes, null), + left_sample(nnodes, null), right_sample(nnodes, null), + next_sample(nnodes, null), sample_index_map(nnodes, null), + above_sample(nnodes, 0), + left{std::numeric_limits::quiet_NaN()}, + right{std::numeric_limits::quiet_NaN()}, left_root(null) + /// Constructor + /// \todo Document + { + } + + int + num_roots() const + /// Return number of roots + { + if (left_root == null) + { + throw std::runtime_error("left_root is NULL"); + } + int nroots = 0; + auto lr = left_root; + while (lr != null) + { + ++nroots; + lr = right_sib[lr]; + } + return nroots; + } + + inline std::size_t + sample_size() const + /// Number of samples + { + return samples_list.size(); + } + + inline typename std::vector::const_iterator + samples_list_begin() const + /// Beginning of samples list + { + return begin(samples_list); + } + + inline typename std::vector::const_iterator + samples_list_end() const + /// End itertor for samples list + { + return end(samples_list); + } + + inline id_type + sample_group(id_type u) const + /// Return the sample group for node \u + { + if (static_cast(u) >= num_nodes) + { + throw std::invalid_argument("invalid node id"); + } + return sample_groups[u]; + } + + inline bool + advancing_sample_list() const + { + return advancing_sample_list_; + } + + inline std::size_t + size() const + /// Return the length of the internal vectors. + { + return num_nodes; + } + + inline id_type + sample_table_index_to_node(id_type u) const + /// If u is a sample index, return the associated + /// node id. + { + return samples_list[u]; + } + }; + +#if __cplusplus < 201703L + template + constexpr SignedInteger marginal_tree::null; +#endif + } // namespace types + } // namespace ts +} // namespace fwdpp diff --git a/testsuite/tree_sequences/preorder_adl.hpp b/testsuite/tree_sequences/preorder_adl.hpp index 3c55184a4..79abba05c 100644 --- a/testsuite/tree_sequences/preorder_adl.hpp +++ b/testsuite/tree_sequences/preorder_adl.hpp @@ -27,10 +27,11 @@ class node_traversal_preorder_test : public fwdpp::ts::node_traversal_order Date: Thu, 8 Jul 2021 11:38:20 -0700 Subject: [PATCH 17/19] start moving tree advance functions over... --- fwdpp/ts/types/tree_sequence.hpp | 211 ++++++++++++++++++++++++++++++- 1 file changed, 207 insertions(+), 4 deletions(-) diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index 8b2acb4cf..cf932acdc 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -3,6 +3,7 @@ #include #include #include +#include #if __cplusplus >= 201703L #include #endif @@ -11,6 +12,7 @@ #include "../exceptions.hpp" #include "../node_flags.hpp" #include "../tree_flags.hpp" +#include "../detail/advance_marginal_tree_policies.hpp" namespace fwdpp { @@ -28,9 +30,13 @@ namespace fwdpp // adding a level of memory safety. std::reference_wrapper> treeseq_; marginal_tree tree_; - std::size_t input_edge_index, output_edge_index; - double position; + double position, maxpos; bool advanced_; + typename std::vector::const_iterator current_input_edge, + input_edge_end, current_output_edge, output_edge_end; + typename types::table_collection< + SignedInteger>::edge_table::const_iterator beg_edges, + end_edges; public: tree_iterator(const tree_sequence& ts, @@ -38,8 +44,12 @@ namespace fwdpp : treeseq_{ts}, tree_{ts.num_nodes(), ts.samples(), static_cast( flags & tree_flags::TRACK_SAMPLES)}, - input_edge_index{0}, output_edge_index{0}, position{0.}, advanced_{ - false} + position{0.}, maxpos{tables().genome_length()}, advanced_{false}, + current_input_edge{end(tables().input_left)}, + input_edge_end{end(tables().input_left)}, + current_output_edge{end(tables().output_right)}, + output_edge_end{end(tables().output_right)}, + beg_edges{begin(tables().edges)}, end_edges{end(tables().edges)} { } @@ -73,6 +83,199 @@ namespace fwdpp advance_right() // Move 1 tree left-to-right { + if (current_input_edge < input_edge_end || position < maxpos) + { + while (current_output_edge < output_edge_end + + && (beg_edges + *current_output_edge)->right + == position) // T4 + { + const auto p + = (beg_edges + *current_output_edge)->parent; + const auto c + = (beg_edges + *current_output_edge)->child; + const auto lsib = tree_.left_sib[c]; + const auto rsib = tree_.right_sib[c]; + if (lsib + == types::table_collection::null) + { + tree_.left_child[p] = rsib; + } + else + { + tree_.right_sib[lsib] = rsib; + } + if (rsib + == types::table_collection::null) + { + tree_.right_child[p] = lsib; + } + else + { + tree_.left_sib[rsib] = lsib; + } + tree_.parents[c] + = types::table_collection::null; + tree_.left_sib[c] + = types::table_collection::null; + tree_.right_sib[c] + = types::table_collection::null; + detail::outgoing_leaf_counts( + tree_, + (beg_edges + *current_output_edge)->parent, + (beg_edges + *current_output_edge)->child); + if (tree_.advancing_sample_list()) + { + detail::update_samples_list( + tree_, (beg_edges + *current_output_edge) + ->parent); + } + update_roots_outgoing(p, c, tree_); + ++current_input_edge; + } + while (current_input_edge < input_edge_end + && (beg_edges + *current_input_edge)->left + == position) // Step T2 + { + const auto p + = (beg_edges + *current_input_edge)->parent; + const auto c + = (beg_edges + *current_input_edge)->child; + const auto rchild = tree_.right_child[p]; + const auto lsib = tree_.left_sib[c]; + const auto rsib = tree_.right_sib[c]; + if (rchild + == types::table_collection::null) + { + tree_.left_child[p] = c; + tree_.left_sib[c] = types::table_collection< + SignedInteger>::null; + tree_.right_sib[c] = types::table_collection< + SignedInteger>::null; + } + else + { + tree_.right_sib[rchild] = c; + tree_.left_sib[c] = rchild; + tree_.right_sib[c] = types::table_collection< + SignedInteger>::null; + } + // The entry for the child refers to + // the parent's location in the node table. + tree_.parents[c] + = (beg_edges + *current_input_edge)->parent; + tree_.right_child[p] = c; + detail::incoming_leaf_counts( + tree_, (beg_edges + *current_input_edge)->parent, + (beg_edges + *current_input_edge)->child); + if (tree_.advancing_sample_list()) + { + detail::update_samples_list( + tree_, (beg_edges + *current_input_edge) + ->parent); + } + update_roots_incoming(p, c, lsib, rsib, tree_); + + ++current_input_edge; + } + + // This is a big "gotcha". + // The root tracking functions will sometimes + // result in left_root not actually being the left_root. + // We loop through the left_sibs to fix that. + if (tree_.left_root + != types::table_collection::null) + { + while ( + tree_.left_sib[tree_.left_root] + != types::table_collection::null) + { + tree_.left_root + = tree_.left_sib[tree_.left_root]; + } + } +#ifndef NDEBUG + // Validate the roots via brute-force. + auto lr = tree_.left_root; + if (lr == types::table_collection::null) + { + throw std::runtime_error( + "FWDPP DEBUG: left_root is null"); + } + std::vector is_root(tree_.sample_index_map.size(), 0); + std::vector processed(is_root.size(), 0); + for (std::size_t s = 0; s < tree_.sample_index_map.size(); + ++s) + { + if (tree_.sample_index_map[s] + != types::table_collection::null) + { + auto u = static_cast(s); + auto root = u; + bool early_exit = false; + while (u + != types::table_collection< + SignedInteger>::null) + { + if (processed[u]) + { + early_exit = true; + break; + } + processed[u] = 1; + root = u; + u = tree_.parents[u]; + } + if (early_exit == false) + { + is_root[root] = 1; + } + } + } + int nroots_brute = 0; + for (auto r : is_root) + { + nroots_brute += r; + } + if (nroots_brute != tree_.num_roots()) + { + throw std::runtime_error("FWDPP DEBUG: num_roots " + "disagreement"); + } + while (lr != types::table_collection::null) + { + if (is_root[lr] != 1) + { + throw std::runtime_error("FWDPP DEBUG: root " + "contents " + "disagreement"); + } + lr = tree_.right_sib[lr]; + } +#endif + double right = maxpos; + if (current_input_edge < input_edge_end) + { + right = std::min( + right, (beg_edges + *current_input_edge)->left); + } + if (current_output_edge < output_edge_end) + { + right = std::min( + right, + (beg_edges + *current_output_edge)->right); + } + tree_.left = position; + tree_.right = right; + // Must set return value before + // updating right, else the + // last tree will be skipped. + bool rv = current_input_edge < input_edge_end + || position < maxpos; + position = right; + advanced_ = rv; + } + advanced_ = false; } }; From 4fb333f7e606928cf1b1836fe544edaf7515af32 Mon Sep 17 00:00:00 2001 From: molpopogen Date: Thu, 8 Jul 2021 11:49:39 -0700 Subject: [PATCH 18/19] move root updating to ts::detail --- .../detail/advance_marginal_tree_policies.hpp | 170 ++++++++++++++++++ fwdpp/ts/types/tree_sequence.hpp | 5 +- 2 files changed, 173 insertions(+), 2 deletions(-) diff --git a/fwdpp/ts/detail/advance_marginal_tree_policies.hpp b/fwdpp/ts/detail/advance_marginal_tree_policies.hpp index 338678128..8c3ee63a4 100644 --- a/fwdpp/ts/detail/advance_marginal_tree_policies.hpp +++ b/fwdpp/ts/detail/advance_marginal_tree_policies.hpp @@ -6,6 +6,7 @@ #include #include "../types/generate_null_id.hpp" #include "../marginal_tree.hpp" +#include "../types/generate_null_id.hpp" namespace fwdpp { @@ -13,6 +14,175 @@ namespace fwdpp { namespace detail { + template + void + update_roots_outgoing(SignedInteger p, SignedInteger c, + marginal_tree& marginal) + // This is the algorithm used by tskit. + // The method is applied to nodes p and c + // AFTER the output index has been applied + // to a marginal tree. + // + // The criteria to mark something as a root + // are that its parent is NULL and it is + // above a sample. Here, c's parent is NULL + // because that is the result of processing + // "outgoing" nodes. Thus, c will be + // added to the list of nodes. + // + // The trick is to determine if the most + // ancient ancestor of p is NOT a root, + // which occurs if that node is not above any samples. + // + // While we are at it, we update above_sample. + // A node is above a sample if it is a sample + // or if any of its children are above a sample. + { + if (marginal.above_sample[c] == 1) + { + auto x = p; + auto root = x; + std::int8_t above_sample = 0; + while (x != types::generate_null_id() + && above_sample == 0) + { + // If this node is a sample, then + // it must descend from a root, + // which is why we OR this in the loop below. + above_sample + = (marginal.sample_index_map[x] + != types::generate_null_id()); + auto lc = marginal.left_child[x]; + while (lc != types::generate_null_id() + && above_sample == 0) + { + above_sample + = above_sample || marginal.above_sample[lc]; + lc = marginal.right_sib[lc]; + } + marginal.above_sample[x] = above_sample; + root = x; + x = marginal.parents[x]; + } + // Now, root refers to the most ancient + // ancestor of p found in the above loop + if (above_sample == 0) + { + // Remove root from list of roots + auto lroot = marginal.left_sib[root]; + auto rroot = marginal.right_sib[root]; + marginal.left_root + = types::generate_null_id(); + if (lroot != types::generate_null_id()) + { + marginal.right_sib[lroot] = rroot; + marginal.left_root = lroot; + } + if (rroot != types::generate_null_id()) + { + marginal.left_sib[rroot] = lroot; + marginal.left_root = rroot; + } + marginal.left_sib[root] + = types::generate_null_id(); + marginal.right_sib[root] + = types::generate_null_id(); + } + if (marginal.left_root + != types::generate_null_id()) + { + //Put c into a pre-existing root list + auto lroot = marginal.left_sib[marginal.left_root]; + if (lroot != types::generate_null_id()) + { + marginal.right_sib[lroot] = c; + } + marginal.left_sib[c] = lroot; + marginal.left_sib[marginal.left_root] = c; + } + marginal.right_sib[c] = marginal.left_root; + marginal.left_root = c; + } + } + + template + void + update_roots_incoming(SignedInteger p, SignedInteger c, SignedInteger lsib, + SignedInteger rsib, + marginal_tree& marginal) + // This is the algorithm used by tskit. It is applied + // AFTER the incoming node list has updated marginal. + // Thus, p is parent to c. + // + // lsib and rsib are with respect to c BEFORE the incoming + // node list has updated marginal. This is a bit confusing: + // these values are used to remove c from the root list if + // necessary. + // + // NOTE: all values of c for which above_sample[c] == 1 + // are in the root list. + { + if (marginal.above_sample[c]) + { + auto x = p; + auto root = x; + + std::int8_t above_sample = 0; + while (x != types::generate_null_id() + && above_sample == 0) + { + above_sample = marginal.above_sample[x]; + // c is above_sample and p is c's parent. + // Thus, all parents to p are above_sample, too. + marginal.above_sample[x] = 1; + root = x; + x = marginal.parents[x]; + } + if (above_sample == 0) + // If we are here, then the above loop terminated + // by encountering a null node, because above_sample[x] + // must have been 0 for all x. However, because c is + // above sample, all nodes encountered have been update + // to be above_sample as well. Thus, the new value of root + // replaces c in the root list. + { + // Replace c with root in root list + if (lsib != types::generate_null_id()) + { + marginal.right_sib[lsib] = root; + } + if (rsib != types::generate_null_id()) + { + marginal.left_sib[rsib] = root; + } + marginal.left_sib[root] = lsib; + marginal.right_sib[root] = rsib; + marginal.left_root = root; + } + else + // If we are here, then we encountered a node + // ancestral to c where above_sample == 1. + // Thus, c can no longer be a root. If the current + // p is also a c in a later call to this function, then + // it may also be removed, etc.. + { + // Remove c from root list + marginal.left_root + = types::generate_null_id(); + if (lsib != types::generate_null_id()) + { + marginal.right_sib[lsib] = rsib; + marginal.left_root = lsib; + } + if (rsib != types::generate_null_id()) + { + marginal.left_sib[rsib] = lsib; + marginal.left_root = rsib; + } + } + } + } + template inline void outgoing_leaf_counts(marginal_tree& marginal, diff --git a/fwdpp/ts/types/tree_sequence.hpp b/fwdpp/ts/types/tree_sequence.hpp index cf932acdc..2ceb354e5 100644 --- a/fwdpp/ts/types/tree_sequence.hpp +++ b/fwdpp/ts/types/tree_sequence.hpp @@ -130,7 +130,7 @@ namespace fwdpp tree_, (beg_edges + *current_output_edge) ->parent); } - update_roots_outgoing(p, c, tree_); + detail::update_roots_outgoing(p, c, tree_); ++current_input_edge; } while (current_input_edge < input_edge_end @@ -174,7 +174,8 @@ namespace fwdpp tree_, (beg_edges + *current_input_edge) ->parent); } - update_roots_incoming(p, c, lsib, rsib, tree_); + detail::update_roots_incoming(p, c, lsib, rsib, + tree_); ++current_input_edge; } From a8d7c47b0609b719f1736a08da9249cb2fdb99ea Mon Sep 17 00:00:00 2001 From: molpopogen Date: Thu, 8 Jul 2021 11:51:55 -0700 Subject: [PATCH 19/19] add ts/marginal_tree.hpp --- fwdpp/ts/marginal_tree.hpp | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 fwdpp/ts/marginal_tree.hpp diff --git a/fwdpp/ts/marginal_tree.hpp b/fwdpp/ts/marginal_tree.hpp new file mode 100644 index 000000000..2cc62e4fe --- /dev/null +++ b/fwdpp/ts/marginal_tree.hpp @@ -0,0 +1,15 @@ +#pragma once + +#include "types/marginal_tree.hpp" + +namespace fwdpp +{ + namespace ts + { + template + using marginal_tree = types::marginal_tree; + + template + using sample_group_map = types::sample_group_map; + } +}