Skip to content

Decapitate #2240

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

Closed
wants to merge 13 commits into from
Closed

Decapitate #2240

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
6a2c704
Implement decapitate function
jeromekelleher May 4, 2022
a0f960c
Add start of new operation
jeromekelleher May 16, 2022
a3c1569
Implement split_edges
jeromekelleher May 16, 2022
31367be
tests
jeromekelleher May 19, 2022
7b6c9c7
Implement semantics for extra values
jeromekelleher May 19, 2022
8f4fbf7
Update
jeromekelleher May 20, 2022
68deb71
C implementation of split_edges
jeromekelleher May 21, 2022
1168ec5
Apply suggestions from code review
jeromekelleher May 23, 2022
ab8a260
Review tweak
jeromekelleher May 23, 2022
b9c50dc
Add impute population option
jeromekelleher May 23, 2022
3c9ae2e
Implement py-c split edges
jeromekelleher May 23, 2022
e2de3c7
Split edges using C impl
jeromekelleher May 24, 2022
5100269
Start for delete older
jeromekelleher May 24, 2022
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
77 changes: 77 additions & 0 deletions c/tests/test_tables.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10155,6 +10155,80 @@ test_table_collection_clear(void)
| TSK_CLEAR_TS_METADATA_AND_SCHEMA);
}

static void
test_table_collection_decapitate(void)
{
int ret;
tsk_treeseq_t ts;
tsk_table_collection_t t;

tsk_treeseq_from_text(&ts, 10, paper_ex_nodes, paper_ex_edges, NULL, paper_ex_sites,
paper_ex_mutations, paper_ex_individuals, NULL, 0);
ret = tsk_treeseq_copy_tables(&ts, &t, 0);
CU_ASSERT_EQUAL_FATAL(ret, 0);
tsk_treeseq_free(&ts);

/* Add some migrations */
tsk_population_table_add_row(&t.populations, NULL, 0);
tsk_population_table_add_row(&t.populations, NULL, 0);
tsk_migration_table_add_row(&t.migrations, 0, 10, 0, 0, 1, 0.05, NULL, 0);
tsk_migration_table_add_row(&t.migrations, 0, 10, 0, 1, 0, 0.09, NULL, 0);
tsk_migration_table_add_row(&t.migrations, 0, 10, 0, 0, 1, 0.10, NULL, 0);
CU_ASSERT_EQUAL(t.migrations.num_rows, 3);

/* NOTE: haven't worked out the exact IDs on the branches here, just
* for illustration.
0.09┊ 9 5 10 ┊ 9 5 ┊11 5 ┊
┊ ┃ ┏┻┓ ┃ ┊ ┃ ┏━┻┓ ┊ ┃ ┏━┻┓ ┊
0.07┊ ┃ ┃ ┃ ┃ ┊ ┃ ┃ 4 ┊ ┃ ┃ 4 ┊
┊ ┃ ┃ ┃ ┃ ┊ ┃ ┃ ┏┻┓ ┊ ┃ ┃ ┏┻┓ ┊
0.00┊ 0 1 3 2 ┊ 0 1 2 3 ┊ 0 1 2 3 ┊
0.00 2.00 7.00 10.00
*/

ret = tsk_table_collection_decapitate(&t, 0.09, 0);
CU_ASSERT_EQUAL_FATAL(ret, 0);

ret = tsk_treeseq_init(&ts, &t, 0);
CU_ASSERT_EQUAL_FATAL(ret, 0);

CU_ASSERT_EQUAL(tsk_treeseq_get_num_trees(&ts), 3);
CU_ASSERT_EQUAL(tsk_treeseq_get_num_nodes(&ts), 12);
/* Lost the mutation over 5 */
CU_ASSERT_EQUAL(tsk_treeseq_get_num_mutations(&ts), 2);
/* We keep the migration at exactly 0.09. */
CU_ASSERT_EQUAL(tsk_treeseq_get_num_migrations(&ts), 2);

tsk_table_collection_free(&t);
tsk_treeseq_free(&ts);
}

static void
test_table_collection_decapitate_errors(void)
{
int ret;
tsk_treeseq_t ts;
tsk_table_collection_t t;

tsk_treeseq_from_text(&ts, 10, paper_ex_nodes, paper_ex_edges, NULL, paper_ex_sites,
paper_ex_mutations, paper_ex_individuals, NULL, 0);
ret = tsk_treeseq_copy_tables(&ts, &t, 0);
CU_ASSERT_EQUAL_FATAL(ret, 0);
tsk_treeseq_free(&ts);

/* This should be caught later when we try to index */
reverse_edges(&t);
ret = tsk_table_collection_decapitate(&t, 0.09, 0);
CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_EDGES_NOT_SORTED_CHILD);

/* This should be caught immediately on entry to the function */
t.sequence_length = -1;
ret = tsk_table_collection_decapitate(&t, 0.09, 0);
CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_BAD_SEQUENCE_LENGTH);

tsk_table_collection_free(&t);
}

static void
test_table_collection_takeset_indexes(void)
{
Expand Down Expand Up @@ -10317,6 +10391,9 @@ main(int argc, char **argv)
test_table_collection_union_middle_merge },
{ "test_table_collection_union_errors", test_table_collection_union_errors },
{ "test_table_collection_clear", test_table_collection_clear },
{ "test_table_collection_decapitate", test_table_collection_decapitate },
{ "test_table_collection_decapitate_errors",
test_table_collection_decapitate_errors },
{ "test_table_collection_takeset_indexes",
test_table_collection_takeset_indexes },
{ NULL, NULL },
Expand Down
199 changes: 199 additions & 0 deletions c/tests/test_trees.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7249,6 +7249,202 @@ test_reference_sequence(void)
tsk_table_collection_free(&tables);
}

static void
test_split_edges_no_populations(void)
{
int ret;
tsk_treeseq_t ts, split_ts;
tsk_table_collection_t tables;
tsk_id_t new_nodes[] = { 9, 10, 11 };
tsk_size_t num_new_nodes = 3;
const char *metadata = "some metadata";
tsk_size_t j;
tsk_node_t node;
double time = 0.09;
tsk_id_t ret_id;

tsk_treeseq_from_text(&ts, 10, paper_ex_nodes, paper_ex_edges, NULL, paper_ex_sites,
paper_ex_mutations, paper_ex_individuals, NULL, 0);

ret_id = tsk_table_collection_copy(ts.tables, &tables, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 0);
tsk_treeseq_free(&ts);
ret_id = tsk_population_table_add_row(&tables.populations, NULL, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 0);
ret = tsk_table_collection_compute_mutation_times(&tables, NULL, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 0);
ret_id = tsk_treeseq_init(&ts, &tables, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 0);

/* NOTE: haven't worked out the exact IDs on the branches here, just
* for illustration.

0.25┊ 8 ┊ ┊ ┊
┊ ┏━┻━┓ ┊ ┊ ┊
0.20┊ ┃ ┃ ┊ ┊ 7 ┊
┊ ┃ ┃ ┊ ┊ ┏━┻━┓ ┊
0.17┊ 6 ┃ ┊ 6 ┊ ┃ ┃ ┊
┊ ┏━┻┓ ┃ ┊ ┏━┻━┓ ┊ ┃ ┃ ┊
0.09┊ 9 5 10┊ 9 5 ┊ 11 5 ┊
┊ ┃ ┏┻┓ ┃ ┊ ┃ ┏━┻┓ ┊ ┃ ┏━┻┓ ┊
0.07┊ ┃ ┃ ┃ ┃ ┊ ┃ ┃ 4 ┊ ┃ ┃ 4 ┊
┊ ┃ ┃ ┃ ┃ ┊ ┃ ┃ ┏┻┓ ┊ ┃ ┃ ┏┻┓ ┊
0.00┊ 0 1 3 2 ┊ 0 1 2 3 ┊ 0 1 2 3 ┊
0.00 2.00 7.00 10.00
*/
ret = tsk_treeseq_split_edges(
&ts, time, 1234, 0, metadata, strlen(metadata), 0, &split_ts);
CU_ASSERT_EQUAL_FATAL(ret, 0);
CU_ASSERT_EQUAL(tsk_treeseq_get_num_trees(&split_ts), 3);
CU_ASSERT_EQUAL(tsk_treeseq_get_num_nodes(&split_ts), 12);

for (j = 0; j < num_new_nodes; j++) {
ret = tsk_treeseq_get_node(&split_ts, new_nodes[j], &node);
CU_ASSERT_EQUAL_FATAL(ret, 0);
CU_ASSERT_EQUAL(node.time, time);
CU_ASSERT_EQUAL(node.flags, 1234);
CU_ASSERT_EQUAL(node.individual, TSK_NULL);
CU_ASSERT_EQUAL(node.population, 0);
CU_ASSERT_EQUAL(node.metadata_length, strlen(metadata));
CU_ASSERT_EQUAL(strncmp(node.metadata, metadata, strlen(metadata)), 0);
}
tsk_treeseq_free(&split_ts);

/* And again with imputed population value */
ret = tsk_treeseq_split_edges(&ts, time, 1234, 0, metadata, strlen(metadata),
TSK_SPLIT_EDGES_IMPUTE_POPULATION, &split_ts);
CU_ASSERT_EQUAL_FATAL(ret, 0);
CU_ASSERT_EQUAL(tsk_treeseq_get_num_trees(&split_ts), 3);
CU_ASSERT_EQUAL(tsk_treeseq_get_num_nodes(&split_ts), 12);

for (j = 0; j < num_new_nodes; j++) {
ret = tsk_treeseq_get_node(&split_ts, new_nodes[j], &node);
CU_ASSERT_EQUAL_FATAL(ret, 0);
CU_ASSERT_EQUAL(node.time, time);
CU_ASSERT_EQUAL(node.flags, 1234);
CU_ASSERT_EQUAL(node.individual, TSK_NULL);
CU_ASSERT_EQUAL(node.population, TSK_NULL);
CU_ASSERT_EQUAL(node.metadata_length, strlen(metadata));
CU_ASSERT_EQUAL(strncmp(node.metadata, metadata, strlen(metadata)), 0);
}
tsk_treeseq_free(&split_ts);

tsk_table_collection_free(&tables);
tsk_treeseq_free(&ts);
}

static void
test_split_edges_populations(void)
{
int ret;
tsk_treeseq_t ts, split_ts;
tsk_table_collection_t tables;
double time = 0.5;
tsk_node_t node;
tsk_id_t valid_pops[] = { -1, 0, 1 };
tsk_id_t num_valid_pops = 3;
tsk_id_t j, population, ret_id;

ret = tsk_table_collection_init(&tables, 0);
CU_ASSERT_EQUAL_FATAL(ret, 0);
tables.sequence_length = 1;

ret_id = tsk_population_table_add_row(&tables.populations, NULL, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 0);
ret_id = tsk_population_table_add_row(&tables.populations, NULL, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 1);
ret_id = tsk_node_table_add_row(&tables.nodes, 0, 0, 0, TSK_NULL, NULL, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 0);
ret_id = tsk_node_table_add_row(&tables.nodes, 0, 1, 1, TSK_NULL, NULL, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 1);
ret_id = tsk_edge_table_add_row(&tables.edges, 0, 1, 1, 0, NULL, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 0);

ret = tsk_treeseq_init(&ts, &tables, TSK_TS_INIT_BUILD_INDEXES);
CU_ASSERT_EQUAL_FATAL(ret, 0);

for (j = 0; j < num_valid_pops; j++) {
population = valid_pops[j];
ret = tsk_treeseq_split_edges(&ts, time, 0, population, NULL, 0, 0, &split_ts);
CU_ASSERT_EQUAL_FATAL(ret, 0);
CU_ASSERT_EQUAL(tsk_treeseq_get_num_trees(&split_ts), 1);
CU_ASSERT_EQUAL(tsk_treeseq_get_num_nodes(&split_ts), 3);
ret = tsk_treeseq_get_node(&split_ts, 2, &node);
CU_ASSERT_EQUAL_FATAL(ret, 0);
CU_ASSERT_EQUAL(node.population, population);
tsk_treeseq_free(&split_ts);

ret = tsk_treeseq_split_edges(&ts, time, 0, population, NULL, 0,
TSK_SPLIT_EDGES_IMPUTE_POPULATION, &split_ts);
CU_ASSERT_EQUAL_FATAL(ret, 0);
CU_ASSERT_EQUAL(tsk_treeseq_get_num_trees(&split_ts), 1);
CU_ASSERT_EQUAL(tsk_treeseq_get_num_nodes(&split_ts), 3);
ret = tsk_treeseq_get_node(&split_ts, 2, &node);
CU_ASSERT_EQUAL_FATAL(ret, 0);
CU_ASSERT_EQUAL(node.population, 0);
tsk_treeseq_free(&split_ts);
}

tsk_table_collection_free(&tables);
tsk_treeseq_free(&ts);
}

static void
test_split_edges_errors(void)
{
int ret;
tsk_treeseq_t ts, split_ts;
tsk_table_collection_t tables;
double time = 0.5;
tsk_id_t invalid_pops[] = { -2, 2, 3 };
tsk_id_t num_invalid_pops = 3;
tsk_id_t j, population, ret_id;

ret = tsk_table_collection_init(&tables, 0);
CU_ASSERT_EQUAL_FATAL(ret, 0);
tables.sequence_length = 1;

ret_id = tsk_population_table_add_row(&tables.populations, NULL, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 0);
ret_id = tsk_population_table_add_row(&tables.populations, NULL, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 1);
ret_id = tsk_node_table_add_row(&tables.nodes, 0, 0, 0, TSK_NULL, NULL, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 0);
ret_id = tsk_node_table_add_row(&tables.nodes, 0, 1, 1, TSK_NULL, NULL, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 1);
ret_id = tsk_edge_table_add_row(&tables.edges, 0, 1, 1, 0, NULL, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 0);

ret = tsk_treeseq_init(&ts, &tables, TSK_TS_INIT_BUILD_INDEXES);
CU_ASSERT_EQUAL_FATAL(ret, 0);

for (j = 0; j < num_invalid_pops; j++) {
population = invalid_pops[j];
ret = tsk_treeseq_split_edges(&ts, time, 0, population, NULL, 0, 0, &split_ts);
CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_POPULATION_OUT_OF_BOUNDS);
tsk_treeseq_free(&split_ts);

/* We always check population values, even if they aren't used */
ret = tsk_treeseq_split_edges(&ts, time, 0, population, NULL, 0,
TSK_SPLIT_EDGES_IMPUTE_POPULATION, &split_ts);
CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_POPULATION_OUT_OF_BOUNDS);
tsk_treeseq_free(&split_ts);
}
tsk_treeseq_free(&ts);

ret_id
= tsk_migration_table_add_row(&tables.migrations, 0, 1, 0, 0, 1, 1.0, NULL, 0);
CU_ASSERT_EQUAL_FATAL(ret_id, 0);
ret = tsk_treeseq_init(&ts, &tables, TSK_TS_INIT_BUILD_INDEXES);
CU_ASSERT_EQUAL_FATAL(ret, 0);
ret = tsk_treeseq_split_edges(&ts, time, 0, population, NULL, 0, 0, &split_ts);
CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_MIGRATIONS_NOT_SUPPORTED);
tsk_treeseq_free(&split_ts);

tsk_table_collection_free(&tables);
tsk_treeseq_free(&ts);
}

static void
test_init_take_ownership_no_edge_metadata(void)
{
Expand Down Expand Up @@ -7449,6 +7645,9 @@ main(int argc, char **argv)
{ "test_tree_sequence_metadata", test_tree_sequence_metadata },
{ "test_time_uncalibrated", test_time_uncalibrated },
{ "test_reference_sequence", test_reference_sequence },
{ "test_split_edges_no_populations", test_split_edges_no_populations },
{ "test_split_edges_populations", test_split_edges_populations },
{ "test_split_edges_errors", test_split_edges_errors },
{ "test_init_take_ownership_no_edge_metadata",
test_init_take_ownership_no_edge_metadata },
{ NULL, NULL },
Expand Down
Loading