Skip to content

Commit

Permalink
compare to python version and other tidying
Browse files Browse the repository at this point in the history
  • Loading branch information
@petrelharp
petrelharp committed Aug 3, 2023
1 parent be44dc1 commit 33ff451
Show file tree
Hide file tree
Showing 5 changed files with 109 additions and 136 deletions.
26 changes: 11 additions & 15 deletions c/tests/test_tables.c
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
/*
* MIT License
*
* Copyright (c) 2019-2022 Tskit Developers
* Copyright (c) 2019-2023 Tskit Developers
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
Expand Down Expand Up @@ -841,7 +841,7 @@ test_table_collection_metadata(void)
takeset_metadata = tsk_malloc(example_metadata_length * sizeof(char));
CU_ASSERT_FATAL(takeset_metadata != NULL);
memcpy(takeset_metadata, &example_metadata,
(size_t)(example_metadata_length * sizeof(char)));
(size_t) (example_metadata_length * sizeof(char)));

ret = tsk_table_collection_init(&tc1, 0);
CU_ASSERT_EQUAL_FATAL(ret, 0);
Expand Down Expand Up @@ -959,7 +959,7 @@ test_node_table(void)
CU_ASSERT_EQUAL(table.individual[j], j);
CU_ASSERT_EQUAL(table.num_rows, (tsk_size_t) j + 1);
CU_ASSERT_EQUAL(
table.metadata_length, (tsk_size_t)(j + 1) * test_metadata_length);
table.metadata_length, (tsk_size_t) (j + 1) * test_metadata_length);
CU_ASSERT_EQUAL(table.metadata_offset[j + 1], table.metadata_length);
/* check the metadata */
tsk_memcpy(metadata_copy, table.metadata + table.metadata_offset[j],
Expand Down Expand Up @@ -1518,7 +1518,7 @@ test_edge_table_with_options(tsk_flags_t options)
CU_ASSERT_EQUAL(table.metadata_offset, NULL);
} else {
CU_ASSERT_EQUAL(
table.metadata_length, (tsk_size_t)(j + 1) * test_metadata_length);
table.metadata_length, (tsk_size_t) (j + 1) * test_metadata_length);
CU_ASSERT_EQUAL(table.metadata_offset[j + 1], table.metadata_length);
/* check the metadata */
tsk_memcpy(metadata_copy, table.metadata + table.metadata_offset[j],
Expand Down Expand Up @@ -2213,8 +2213,7 @@ test_extend_edges(void)
ret = tsk_table_collection_extend_edges(&tables, 10, 0);
CU_ASSERT_EQUAL_FATAL(ret, 0);
CU_ASSERT_TRUE(tsk_table_collection_equals(&tables, &tables_copy, 0));

// Free things.

tsk_table_collection_free(&tables);
tsk_table_collection_free(&tables_copy);
}
Expand All @@ -2224,9 +2223,6 @@ test_edge_table_squash(void)
{
int ret;
tsk_table_collection_t tables;
ret = tsk_table_collection_init(&tables, 0);
CU_ASSERT_EQUAL_FATAL(ret, 0);
tables.sequence_length = 10;

const char *nodes_ex = "1 0 -1 -1\n"
"1 0 -1 -1\n"
Expand Down Expand Up @@ -3695,7 +3691,7 @@ test_migration_table(void)
CU_ASSERT_EQUAL(table.time[j], j);
CU_ASSERT_EQUAL(table.num_rows, (tsk_size_t) j + 1);
CU_ASSERT_EQUAL(
table.metadata_length, (tsk_size_t)(j + 1) * test_metadata_length);
table.metadata_length, (tsk_size_t) (j + 1) * test_metadata_length);
CU_ASSERT_EQUAL(table.metadata_offset[j + 1], table.metadata_length);
/* check the metadata */
tsk_memcpy(metadata_copy, table.metadata + table.metadata_offset[j],
Expand Down Expand Up @@ -4314,7 +4310,7 @@ test_individual_table(void)
table.location[spatial_dimension * (size_t) j + k], test_location[k]);
}
CU_ASSERT_EQUAL(
table.metadata_length, (tsk_size_t)(j + 1) * test_metadata_length);
table.metadata_length, (tsk_size_t) (j + 1) * test_metadata_length);
CU_ASSERT_EQUAL(table.metadata_offset[j + 1], table.metadata_length);
/* check the metadata */
tsk_memcpy(metadata_copy, table.metadata + table.metadata_offset[j],
Expand Down Expand Up @@ -4368,7 +4364,7 @@ test_individual_table(void)
flags = tsk_malloc(num_rows * sizeof(tsk_flags_t));
CU_ASSERT_FATAL(flags != NULL);
for (k = 0; k < num_rows; k++) {
flags[k] = (tsk_flags_t)(k + num_rows);
flags[k] = (tsk_flags_t) (k + num_rows);
}
location = tsk_malloc(spatial_dimension * num_rows * sizeof(double));
CU_ASSERT_FATAL(location != NULL);
Expand All @@ -4383,7 +4379,7 @@ test_individual_table(void)
parents = tsk_malloc(num_parents * num_rows * sizeof(tsk_id_t));
CU_ASSERT_FATAL(parents != NULL);
for (k = 0; k < num_parents * num_rows; k++) {
parents[k] = (tsk_id_t)(k + (num_rows * 4));
parents[k] = (tsk_id_t) (k + (num_rows * 4));
}
parents_offset = tsk_malloc((num_rows + 1) * sizeof(tsk_size_t));
CU_ASSERT_FATAL(parents_offset != NULL);
Expand Down Expand Up @@ -10352,9 +10348,9 @@ test_table_collection_takeset_indexes(void)
rem = tsk_malloc(t1.edges.num_rows * sizeof(*rem));
CU_ASSERT_FATAL(rem != NULL);
memcpy(ins, t1.indexes.edge_insertion_order,
(size_t)(t1.edges.num_rows * sizeof(*ins)));
(size_t) (t1.edges.num_rows * sizeof(*ins)));
memcpy(
rem, t1.indexes.edge_removal_order, (size_t)(t1.edges.num_rows * sizeof(*rem)));
rem, t1.indexes.edge_removal_order, (size_t) (t1.edges.num_rows * sizeof(*rem)));

ret = tsk_table_collection_copy(&t1, &t2, 0);
CU_ASSERT_EQUAL_FATAL(ret, 0);
Expand Down
89 changes: 43 additions & 46 deletions c/tskit/tables.c
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
/*
* MIT License
*
* Copyright (c) 2019-2022 Tskit Developers
* Copyright (c) 2019-2023 Tskit Developers
* Copyright (c) 2017-2018 University of Oxford
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
Expand Down Expand Up @@ -12340,14 +12340,15 @@ tsk_table_collection_add_and_remap_node(tsk_table_collection_t *self,

typedef struct _edge_list_t {
tsk_id_t edge;
bool extended; // have we decided to extend this one on the current tree?
// the `extended` flags records whether we have decided to extend
// this entry to the current tree?
bool extended;
struct _edge_list_t *next;
} edge_list_t;

static edge_list_t *TSK_WARN_UNUSED
extend_edges_alloc_entry(tsk_blkalloc_t *heap, tsk_id_t edge)
{
// see ancestor_mapper_alloc_interval_list
edge_list_t *x = NULL;

x = tsk_blkalloc_get(heap, sizeof(*x));
Expand Down Expand Up @@ -12378,7 +12379,6 @@ remove_unextended(edge_list_t **head, edge_list_t **tail)
x = px->next;
while (x != NULL) {
if (x->extended) {
// keep it
x->extended = false;
px->next = x;
px = x;
Expand All @@ -12393,14 +12393,12 @@ static int
forward_extend(tsk_table_collection_t *self, int direction)
{
int ret = 0;
double *new_left, *new_right;
tsk_id_t *num_edges;
tsk_id_t tj, tk, ret_id;
tsk_id_t e1, e2, e_in;
tsk_id_t *I, *O;
const tsk_id_t M = (tsk_id_t) self->edges.num_rows;
double left, right;
double *near_edge, *far_edge;
double *near_side, *far_side;
tsk_blkalloc_t edge_list_heap;
edge_list_t *edges_in_head, *edges_in_tail;
edge_list_t *edges_out_head, *edges_out_tail;
Expand All @@ -12409,20 +12407,19 @@ forward_extend(tsk_table_collection_t *self, int direction)
tsk_edge_t edge;
double sign, here, there;
tsk_id_t sign_int;
bool forwards;

forwards = (direction == TSK_DIR_FORWARD);
bool forwards = (direction == TSK_DIR_FORWARD);

num_edges = tsk_malloc(self->nodes.num_rows * sizeof(*num_edges));
new_left = tsk_malloc(self->edges.num_rows * sizeof(*new_left));
new_right = tsk_malloc(self->edges.num_rows * sizeof(*new_right));
if (num_edges == NULL || new_left == NULL || new_right == NULL) {
if (num_edges == NULL) {
ret = TSK_ERR_NO_MEMORY;
goto out;
}
tsk_memset(num_edges, 0x00, self->nodes.num_rows * sizeof(*num_edges));
memcpy(new_left, self->edges.left, self->edges.num_rows * sizeof(*new_left));
memcpy(new_right, self->edges.right, self->edges.num_rows * sizeof(*new_right));

ret = tsk_edge_table_copy(&self->edges, &edges, 0);
if (ret != 0) {
goto out;
}

ret = tsk_blkalloc_init(&edge_list_heap, 8192);
if (ret != 0) {
Expand All @@ -12434,8 +12431,8 @@ forward_extend(tsk_table_collection_t *self, int direction)
here = 0;
sign = 1;
sign_int = 1;
near_edge = self->edges.left;
far_edge = self->edges.right;
near_side = self->edges.left;
far_side = self->edges.right;
tj = 0;
tk = 0;
} else {
Expand All @@ -12444,10 +12441,10 @@ forward_extend(tsk_table_collection_t *self, int direction)
here = self->sequence_length;
sign = -1;
sign_int = -1;
near_edge = self->edges.right;
far_edge = self->edges.left;
tj = M-1;
tk = M-1;
near_side = self->edges.right;
far_side = self->edges.left;
tj = M - 1;
tk = M - 1;
}
left = 0;
edges_in_head = NULL;
Expand All @@ -12459,9 +12456,13 @@ forward_extend(tsk_table_collection_t *self, int direction)
remove_unextended(&edges_in_head, &edges_in_tail);
remove_unextended(&edges_out_head, &edges_out_tail);

while (((tk < M) && (tk >= 0)) && (far_edge[O[tk]] == here)) {
while (((tk < M) && (tk >= 0)) && (far_side[O[tk]] == here)) {
// add edge tk to pending_out
x = extend_edges_alloc_entry(&edge_list_heap, O[tk]);
if (x == NULL) {
ret = TSK_ERR_NO_MEMORY;
goto out;
}
if (edges_out_tail == NULL) {
edges_out_head = x;
} else {
Expand All @@ -12473,9 +12474,13 @@ forward_extend(tsk_table_collection_t *self, int direction)
num_edges[self->edges.child[O[tk]]] -= 1;
tk += sign_int;
}
while (((tj < M) && (tj >= 0)) && (near_edge[I[tj]] == here)) {
while (((tj < M) && (tj >= 0)) && (near_side[I[tj]] == here)) {
// add edge tj to pending_in
x = extend_edges_alloc_entry(&edge_list_heap, I[tj]);
if (x == NULL) {
ret = TSK_ERR_NO_MEMORY;
goto out;
}
if (edges_in_tail == NULL) {
edges_in_head = x;
} else {
Expand All @@ -12490,17 +12495,17 @@ forward_extend(tsk_table_collection_t *self, int direction)
there = forwards ? self->sequence_length : 0;
if (forwards) {
if (tk < M) {
there = TSK_MIN(there, far_edge[O[tk]]);
there = TSK_MIN(there, far_side[O[tk]]);
}
if (tj < M) {
there = TSK_MIN(there, near_edge[I[tj]]);
there = TSK_MIN(there, near_side[I[tj]]);
}
} else {
if (tk >= 0) {
there = TSK_MAX(there, far_edge[O[tk]]);
there = TSK_MAX(there, far_side[O[tk]]);
}
if (tj >= 0) {
there = TSK_MAX(there, near_edge[I[tj]]);
there = TSK_MAX(there, near_side[I[tj]]);
}
}

Expand All @@ -12516,7 +12521,7 @@ forward_extend(tsk_table_collection_t *self, int direction)
for (ex_in = edges_in_head; ex_in != NULL;
ex_in = ex_in->next) {
e_in = ex_in->edge;
if (sign * far_edge[e_in] > sign * here) {
if (sign * far_side[e_in] > sign * here) {
if ((self->edges.child[e1]
== self->edges.child[e_in])
&& (self->edges.parent[e2]
Expand All @@ -12525,13 +12530,13 @@ forward_extend(tsk_table_collection_t *self, int direction)
ex2->extended = true;
ex_in->extended = true;
if (forwards) {
new_right[e1] = there;
new_right[e2] = there;
new_left[e_in] = there;
edges.right[e1] = there;
edges.right[e2] = there;
edges.left[e_in] = there;
} else {
new_left[e1] = there;
new_left[e2] = there;
new_right[e_in] = there;
edges.left[e1] = there;
edges.left[e2] = there;
edges.right[e_in] = there;
}
num_edges[self->edges.parent[e1]] += 2;
}
Expand All @@ -12547,21 +12552,15 @@ forward_extend(tsk_table_collection_t *self, int direction)
}

// done! write out new edge tables
ret = tsk_edge_table_copy(&self->edges, &edges, 0);
if (ret != 0) {
goto out;
}
ret = tsk_edge_table_clear(&self->edges);
if (ret != 0) {
goto out;
}
for (tj = 0; tj < (tsk_id_t) edges.num_rows; tj++) {
left = new_left[tj];
right = new_right[tj];
if (left < right) {
tsk_edge_table_get_row_unsafe(&edges, tj, &edge);
ret_id = tsk_edge_table_add_row(&self->edges, left, right, edge.parent,
edge.child, edge.metadata, edge.metadata_length);
tsk_edge_table_get_row_unsafe(&edges, tj, &edge);
if (edge.left < edge.right) {
ret_id = tsk_edge_table_add_row(&self->edges, edge.left, edge.right,
edge.parent, edge.child, edge.metadata, edge.metadata_length);
if (ret_id < 0) {
ret = (int) ret_id;
goto out;
Expand All @@ -12576,8 +12575,6 @@ forward_extend(tsk_table_collection_t *self, int direction)
out:
tsk_blkalloc_free(&edge_list_heap);
tsk_safe_free(num_edges);
tsk_safe_free(new_left);
tsk_safe_free(new_right);
tsk_edge_table_free(&edges);
return ret;
}
Expand Down
Loading

0 comments on commit 33ff451

Please sign in to comment.