Remove tree_pos iterator from KC distance and haplotype matching algo…

…rithms

c/tskit/haplotype_matching.c

@@ -209,7 +209,6 @@ int
 tsk_ls_hmm_free(tsk_ls_hmm_t *self)
 {
     tsk_tree_free(&self->tree);
-    tsk_tree_position_free(&self->tree_pos);
     tsk_safe_free(self->recombination_rate);
     tsk_safe_free(self->mutation_rate);
     tsk_safe_free(self->recombination_rate);
@@ -247,11 +246,6 @@ tsk_ls_hmm_reset(tsk_ls_hmm_t *self, double value)
     tsk_memset(self->transition_parent, 0xff,
         self->max_transitions * sizeof(*self->transition_parent));
 
-    tsk_tree_position_free(&self->tree_pos);
-    ret = tsk_tree_position_init(&self->tree_pos, self->tree_sequence, 0);
-    if (ret != 0) {
-        goto out;
-    }
     samples = tsk_treeseq_get_samples(self->tree_sequence);
     for (j = 0; j < self->num_samples; j++) {
         u = samples[j];
@@ -260,7 +254,6 @@ tsk_ls_hmm_reset(tsk_ls_hmm_t *self, double value)
         self->transition_index[u] = (tsk_id_t) j;
     }
     self->num_transitions = self->num_samples;
-out:
     return ret;
 }
 
@@ -310,12 +303,9 @@ tsk_ls_hmm_update_tree(tsk_ls_hmm_t *self, int direction)
     tsk_id_t u, c, p, j, e;
     tsk_value_transition_t *vt;
 
-    tsk_tree_position_step(&self->tree_pos, direction);
-    tsk_bug_assert(self->tree_pos.index != -1);
-    tsk_bug_assert(self->tree_pos.index == self->tree.index);
-
-    for (j = self->tree_pos.out.start; j != self->tree_pos.out.stop; j += direction) {
-        e = self->tree_pos.out.order[j];
+    for (j = self->tree.tree_pos.out.start; j != self->tree.tree_pos.out.stop;
+         j += direction) {
+        e = self->tree.tree_pos.out.order[j];
         c = edges_child[e];
         u = c;
         if (T_index[u] == TSK_NULL) {
@@ -333,8 +323,9 @@ tsk_ls_hmm_update_tree(tsk_ls_hmm_t *self, int direction)
         parent[c] = TSK_NULL;
     }
 
-    for (j = self->tree_pos.in.start; j != self->tree_pos.in.stop; j += direction) {
-        e = self->tree_pos.in.order[j];
+    for (j = self->tree.tree_pos.in.start; j != self->tree.tree_pos.in.stop;
+         j += direction) {
+        e = self->tree.tree_pos.in.order[j];
         c = edges_child[e];
         p = edges_parent[e];
         parent[c] = p;

c/tskit/haplotype_matching.h

@@ -1,7 +1,7 @@
 /*
  * MIT License
  *
- * Copyright (c) 2019-2023 Tskit Developers
+ * Copyright (c) 2019-2024 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
@@ -98,10 +98,6 @@ typedef struct _tsk_ls_hmm_t {
     tsk_size_t num_nodes;
     /* state */
     tsk_tree_t tree;
-    /* NOTE: this tree_position will be redundant once we integrate the top-level
-     * tree class with this.
-     */
-    tsk_tree_position_t tree_pos;
     tsk_id_t *parent;
     /* The probability value transitions on the tree */
     tsk_value_transition_t *transitions;

c/tskit/trees.c

@@ -6974,8 +6974,7 @@ update_kc_subtree_state(
 }
 
 static int
-update_kc_incremental(
-    tsk_tree_t *tree, kc_vectors *kc, tsk_tree_position_t *tree_pos, tsk_size_t *depths)
+update_kc_incremental(tsk_tree_t *tree, kc_vectors *kc, tsk_size_t *depths)
 {
     int ret = 0;
     tsk_id_t u, v, e, j;
@@ -6984,13 +6983,9 @@ update_kc_incremental(
     const tsk_id_t *restrict edges_child = tree->tree_sequence->tables->edges.child;
     const tsk_id_t *restrict edges_parent = tree->tree_sequence->tables->edges.parent;
 
-    tsk_bug_assert(tree_pos->index == tree->index);
-    tsk_bug_assert(tree_pos->interval.left == tree->interval.left);
-    tsk_bug_assert(tree_pos->interval.right == tree->interval.right);
-
     /* Update state of detached subtrees */
-    for (j = tree_pos->out.stop - 1; j >= tree_pos->out.start; j--) {
-        e = tree_pos->out.order[j];
+    for (j = tree->tree_pos.out.stop - 1; j >= tree->tree_pos.out.start; j--) {
+        e = tree->tree_pos.out.order[j];
         u = edges_child[e];
         depths[u] = 0;
 
@@ -7004,8 +6999,8 @@ update_kc_incremental(
     }
 
     /* Propagate state change down into reattached subtrees. */
-    for (j = tree_pos->in.stop - 1; j >= tree_pos->in.start; j--) {
-        e = tree_pos->in.order[j];
+    for (j = tree->tree_pos.in.stop - 1; j >= tree->tree_pos.in.start; j--) {
+        e = tree->tree_pos.in.order[j];
         u = edges_child[e];
         v = edges_parent[e];
 
@@ -7038,17 +7033,11 @@ tsk_treeseq_kc_distance(const tsk_treeseq_t *self, const tsk_treeseq_t *other,
     const tsk_treeseq_t *treeseqs[2] = { self, other };
     tsk_tree_t trees[2];
     kc_vectors kcs[2];
-    /* TODO the tree_pos here is redundant because we should be using this interally
-     * in the trees to do the advancing. Once we have converted the tree over to using
-     * tree_pos internally, we can get rid of these tree_pos variables and use
-     * the values stored in the trees themselves */
-    tsk_tree_position_t tree_pos[2];
     tsk_size_t *depths[2];
     int ret = 0;
 
     for (i = 0; i < 2; i++) {
         tsk_memset(&trees[i], 0, sizeof(trees[i]));
-        tsk_memset(&tree_pos[i], 0, sizeof(tree_pos[i]));
         tsk_memset(&kcs[i], 0, sizeof(kcs[i]));
         depths[i] = NULL;
     }
@@ -7064,10 +7053,6 @@ tsk_treeseq_kc_distance(const tsk_treeseq_t *self, const tsk_treeseq_t *other,
         if (ret != 0) {
             goto out;
         }
-        ret = tsk_tree_position_init(&tree_pos[i], treeseqs[i], 0);
-        if (ret != 0) {
-            goto out;
-        }
         ret = kc_vectors_alloc(&kcs[i], n);
         if (ret != 0) {
             goto out;
@@ -7091,10 +7076,8 @@ tsk_treeseq_kc_distance(const tsk_treeseq_t *self, const tsk_treeseq_t *other,
     if (ret != 0) {
         goto out;
     }
-    tsk_tree_position_next(&tree_pos[0]);
-    tsk_bug_assert(tree_pos[0].index == 0);
 
-    ret = update_kc_incremental(&trees[0], &kcs[0], &tree_pos[0], depths[0]);
+    ret = update_kc_incremental(&trees[0], &kcs[0], depths[0]);
     if (ret != 0) {
         goto out;
     }
@@ -7103,17 +7086,11 @@ tsk_treeseq_kc_distance(const tsk_treeseq_t *self, const tsk_treeseq_t *other,
         if (ret != 0) {
             goto out;
         }
-        tsk_tree_position_next(&tree_pos[1]);
-        tsk_bug_assert(tree_pos[1].index != -1);
 
-        ret = update_kc_incremental(&trees[1], &kcs[1], &tree_pos[1], depths[1]);
+        ret = update_kc_incremental(&trees[1], &kcs[1], depths[1]);
         if (ret != 0) {
             goto out;
         }
-        tsk_bug_assert(trees[0].interval.left == tree_pos[0].interval.left);
-        tsk_bug_assert(trees[0].interval.right == tree_pos[0].interval.right);
-        tsk_bug_assert(trees[1].interval.left == tree_pos[1].interval.left);
-        tsk_bug_assert(trees[1].interval.right == tree_pos[1].interval.right);
         while (trees[0].interval.right < trees[1].interval.right) {
             span = trees[0].interval.right - left;
             total += norm_kc_vectors(&kcs[0], &kcs[1], lambda_) * span;
@@ -7125,9 +7102,7 @@ tsk_treeseq_kc_distance(const tsk_treeseq_t *self, const tsk_treeseq_t *other,
             if (ret != 0) {
                 goto out;
             }
-            tsk_tree_position_next(&tree_pos[0]);
-            tsk_bug_assert(tree_pos[0].index != -1);
-            ret = update_kc_incremental(&trees[0], &kcs[0], &tree_pos[0], depths[0]);
+            ret = update_kc_incremental(&trees[0], &kcs[0], depths[0]);
             if (ret != 0) {
                 goto out;
             }
@@ -7144,7 +7119,6 @@ tsk_treeseq_kc_distance(const tsk_treeseq_t *self, const tsk_treeseq_t *other,
 out:
     for (i = 0; i < 2; i++) {
         tsk_tree_free(&trees[i]);
-        tsk_tree_position_free(&tree_pos[i]);
         kc_vectors_free(&kcs[i]);
         tsk_safe_free(depths[i]);
     }