Skip to content

add tree sequence #349

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
molpopgen wants to merge 19 commits into
base: dev
Choose a base branch
from
Draft

add tree sequence #349

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
f6ad031
Add types/tree_sequence.hpp
molpopgen Jul 6, 2021
1eac491
Add initialization helper functions
molpopgen Jul 7, 2021
f0ac0da
Add test file.
molpopgen Jul 7, 2021
582480f
add init_num_trees
molpopgen Jul 7, 2021
a7baf2a
Add ts/tree_sequence_flags.hpp
molpopgen Jul 7, 2021
d16a7a4
add support for treeseq flags
molpopgen Jul 7, 2021
4ed6aef
uncouple tree sequence from marginal_tree
molpopgen Jul 7, 2021
3e57311
remove flags from treeseq constructor
molpopgen Jul 7, 2021
29ab1f2
start thinking about tree iterators
molpopgen Jul 7, 2021
cf96361
clean up ownership semantics
molpopgen Jul 7, 2021
3887604
optional should wrap a reference_wrapper
molpopgen Jul 7, 2021
8c4eb1a
more api tweaks
molpopgen Jul 7, 2021
447557a
add empty advance fn
molpopgen Jul 7, 2021
67001f4
interface streamline/fix
molpopgen Jul 8, 2021
fd5352a
NOTE to self
molpopgen Jul 8, 2021
613d563
move marginal_tree definition into ts::types
molpopgen Jul 8, 2021
fc093ef
start moving tree advance functions over...
molpopgen Jul 8, 2021
4fb333f
move root updating to ts::detail
molpopgen Jul 8, 2021
a8d7c47
add ts/marginal_tree.hpp
molpopgen Jul 8, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion fwdpp/ts/Makefile.am
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ pkginclude_HEADERS= definitions.hpp \
mutation_tools.hpp \
visit_sites.hpp \
site_visitor.hpp \
recording.hpp
recording.hpp \
tree_flags.hpp


170 changes: 170 additions & 0 deletions fwdpp/ts/detail/advance_marginal_tree_policies.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,183 @@
#include <cassert>
#include "../types/generate_null_id.hpp"
#include "../marginal_tree.hpp"
#include "../types/generate_null_id.hpp"

namespace fwdpp
{
namespace ts
{
namespace detail
{
template <typename SignedInteger>
void
update_roots_outgoing(SignedInteger p, SignedInteger c,
marginal_tree<SignedInteger>& 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<SignedInteger>()
&& 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<SignedInteger>());
auto lc = marginal.left_child[x];
while (lc != types::generate_null_id<SignedInteger>()
&& 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<SignedInteger>();
if (lroot != types::generate_null_id<SignedInteger>())
{
marginal.right_sib[lroot] = rroot;
marginal.left_root = lroot;
}
if (rroot != types::generate_null_id<SignedInteger>())
{
marginal.left_sib[rroot] = lroot;
marginal.left_root = rroot;
}
marginal.left_sib[root]
= types::generate_null_id<SignedInteger>();
marginal.right_sib[root]
= types::generate_null_id<SignedInteger>();
}
if (marginal.left_root
!= types::generate_null_id<SignedInteger>())
{
//Put c into a pre-existing root list
auto lroot = marginal.left_sib[marginal.left_root];
if (lroot != types::generate_null_id<SignedInteger>())
{
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 <typename SignedInteger>
void
update_roots_incoming(SignedInteger p, SignedInteger c, SignedInteger lsib,
SignedInteger rsib,
marginal_tree<SignedInteger>& 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<SignedInteger>()
&& 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<SignedInteger>())
{
marginal.right_sib[lsib] = root;
}
if (rsib != types::generate_null_id<SignedInteger>())
{
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<SignedInteger>();
if (lsib != types::generate_null_id<SignedInteger>())
{
marginal.right_sib[lsib] = rsib;
marginal.left_root = lsib;
}
if (rsib != types::generate_null_id<SignedInteger>())
{
marginal.left_sib[rsib] = lsib;
marginal.left_root = rsib;
}
}
}
}

template <typename SignedInteger>
inline void
outgoing_leaf_counts(marginal_tree<SignedInteger>& marginal,
Expand Down
Loading