Divergence matrix tree-by-tree algorithms

Implement the basic version of the divergence matrix operation using tree-by-tree algorithms, and provide interface for parallelising along the genome.
jeromekelleher · Jul 11, 2023 · f84f5e3 · f84f5e3
1 parent 5eb173a
commit f84f5e3
Show file tree

Hide file tree

Showing 14 changed files with 2,230 additions and 29 deletions.
diff --git a/c/tests/test_stats.c b/c/tests/test_stats.c
diff --git a/c/tests/test_trees.c b/c/tests/test_trees.c
@@ -5395,7 +5395,6 @@ test_simplify_keep_input_roots_multi_tree(void)
 
     tsk_treeseq_from_text(&ts, 10, paper_ex_nodes, paper_ex_edges, NULL, paper_ex_sites,
         paper_ex_mutations, paper_ex_individuals, NULL, 0);
-    tsk_treeseq_dump(&ts, "tmp.trees", 0);
     ret = tsk_treeseq_simplify(
         &ts, samples, 2, TSK_SIMPLIFY_KEEP_INPUT_ROOTS, &simplified, NULL);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
@@ -7801,7 +7800,7 @@ test_time_uncalibrated(void)
     tsk_size_t sample_set_sizes[] = { 2, 2 };
     tsk_id_t samples[] = { 0, 1, 2, 3 };
     tsk_size_t num_samples;
-    double result[10];
+    double result[100];
     double *W;
     double *sigma;
 
@@ -7857,6 +7856,12 @@ test_time_uncalibrated(void)
         TSK_STAT_BRANCH | TSK_STAT_ALLOW_TIME_UNCALIBRATED, sigma);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
 
+    ret = tsk_treeseq_divergence_matrix(&ts2, 0, NULL, 0, NULL, TSK_STAT_BRANCH, result);
+    CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_TIME_UNCALIBRATED);
+    ret = tsk_treeseq_divergence_matrix(&ts2, 0, NULL, 0, NULL,
+        TSK_STAT_BRANCH | TSK_STAT_ALLOW_TIME_UNCALIBRATED, result);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+
     tsk_safe_free(W);
     tsk_safe_free(sigma);
     tsk_treeseq_free(&ts);

diff --git a/c/tests/testlib.c b/c/tests/testlib.c
@@ -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
@@ -966,6 +966,16 @@ tskit_suite_init(void)
     return CUE_SUCCESS;
 }
 
+void
+assert_arrays_almost_equal(tsk_size_t len, double *a, double *b)
+{
+    tsk_size_t j;
+
+    for (j = 0; j < len; j++) {
+        CU_ASSERT_DOUBLE_EQUAL(a[j], b[j], 1e-9);
+    }
+}
+
 static int
 tskit_suite_cleanup(void)
 {

diff --git a/c/tests/testlib.h b/c/tests/testlib.h
@@ -1,7 +1,7 @@
 /*
  * MIT License
  *
- * Copyright (c) 2019-2021 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
@@ -54,6 +54,8 @@ void parse_individuals(const char *text, tsk_individual_table_t *individual_tabl
 
 void unsort_edges(tsk_edge_table_t *edges, size_t start);
 
+void assert_arrays_almost_equal(tsk_size_t len, double *a, double *b);
+
 extern const char *single_tree_ex_nodes;
 extern const char *single_tree_ex_edges;
 extern const char *single_tree_ex_sites;

diff --git a/c/tskit/core.c b/c/tskit/core.c
@@ -466,6 +466,15 @@ tsk_strerror_internal(int err)
             ret = "Statistics using branch lengths cannot be calculated when time_units "
                   "is 'uncalibrated'. (TSK_ERR_TIME_UNCALIBRATED)";
             break;
+        case TSK_ERR_STAT_POLARISED_UNSUPPORTED:
+            ret = "The TSK_STAT_POLARISED option is not supported by this statistic. "
+                  "(TSK_ERR_STAT_POLARISED_UNSUPPORTED)";
+            break;
+        case TSK_ERR_STAT_SPAN_NORMALISE_UNSUPPORTED:
+            ret = "The TSK_STAT_SPAN_NORMALISE option is not supported by this "
+                  "statistic. "
+                  "(TSK_ERR_STAT_SPAN_NORMALISE_UNSUPPORTED)";
+            break;
 
         /* Mutation mapping errors */
         case TSK_ERR_GENOTYPES_ALL_MISSING:

diff --git a/c/tskit/core.h b/c/tskit/core.h
@@ -675,6 +675,16 @@ Statistics based on branch lengths were attempted when the ``time_units``
 were ``uncalibrated``.
 */
 #define TSK_ERR_TIME_UNCALIBRATED                                   -910
+/**
+The TSK_STAT_POLARISED option was passed to a statistic that does
+not support it.
+*/
+#define TSK_ERR_STAT_POLARISED_UNSUPPORTED                          -911
+/**
+The TSK_STAT_SPAN_NORMALISE option was passed to a statistic that does
+not support it.
+*/
+#define TSK_ERR_STAT_SPAN_NORMALISE_UNSUPPORTED                     -912
 /** @} */
 
 /**