Merge pull request #312 from jeromekelleher/remove-pairwise-diversity

Removed old pairwise_diversity code.

Author: jeromekelleher

c/dev-tools/dev-cli.c

@@ -157,29 +157,6 @@ print_ld_matrix(tsk_treeseq_t *ts)
     }
 }
 
-static void
-print_stats(tsk_treeseq_t *ts)
-{
-    int ret = 0;
-    uint32_t j;
-    size_t num_samples = tsk_treeseq_get_num_samples(ts) / 2;
-    tsk_id_t *sample = malloc(num_samples * sizeof(tsk_id_t));
-    double pi;
-
-    if (sample == NULL) {
-        fatal_error("no memory");
-    }
-    for (j = 0; j < num_samples; j++) {
-        sample[j] = (tsk_id_t) j;
-    }
-    ret = tsk_treeseq_get_pairwise_diversity(ts, sample, num_samples, &pi);
-    if (ret != 0) {
-        fatal_library_error(ret, "get_pairwise_diversity");
-    }
-    printf("pi = %f\n", pi);
-    free(sample);
-}
-
 static void
 print_newick_trees(tsk_treeseq_t *ts)
 {
@@ -292,16 +269,6 @@ run_newick(const char *filename, int TSK_UNUSED(verbose))
     tsk_treeseq_free(&ts);
 }
 
-static void
-run_stats(const char *filename, int TSK_UNUSED(verbose))
-{
-    tsk_treeseq_t ts;
-
-    load_tree_sequence(&ts, filename);
-    print_stats(&ts);
-    tsk_treeseq_free(&ts);
-}
-
 static void
 run_simplify(const char *input_filename, const char *output_filename, size_t num_samples,
         bool filter_sites, int verbose)
@@ -399,14 +366,6 @@ main(int argc, char** argv)
     void* argtable6[] = {cmd6, verbose6, infiles6, end6};
     int nerrors6;
 
-    /* SYNTAX 7: stats [-v] <input-file> */
-    struct arg_rex *cmd7 = arg_rex1(NULL, NULL, "stats", NULL, REG_ICASE, NULL);
-    struct arg_lit *verbose7 = arg_lit0("v", "verbose", NULL);
-    struct arg_file *infiles7 = arg_file1(NULL, NULL, NULL, NULL);
-    struct arg_end *end7 = arg_end(20);
-    void* argtable7[] = {cmd7, verbose7, infiles7, end7};
-    int nerrors7;
-
     int exitcode = EXIT_SUCCESS;
     const char *progname = "main";
 
@@ -419,7 +378,6 @@ main(int argc, char** argv)
     nerrors4 = arg_parse(argc, argv, argtable4);
     nerrors5 = arg_parse(argc, argv, argtable5);
     nerrors6 = arg_parse(argc, argv, argtable6);
-    nerrors7 = arg_parse(argc, argv, argtable7);
 
     if (nerrors1 == 0) {
         run_simplify(infiles1->filename[0], outfiles1->filename[0],
@@ -432,11 +390,9 @@ main(int argc, char** argv)
     } else if (nerrors4 == 0) {
         run_variants(infiles4->filename[0], verbose4->count);
     } else if (nerrors5 == 0) {
-        run_stats(infiles5->filename[0], verbose5->count);
+        run_print(infiles5->filename[0], verbose5->count);
     } else if (nerrors6 == 0) {
-        run_print(infiles6->filename[0], verbose6->count);
-    } else if (nerrors7 == 0) {
-        run_newick(infiles7->filename[0], verbose7->count);
+        run_newick(infiles6->filename[0], verbose6->count);
     } else {
         /* We get here if the command line matched none of the possible syntaxes */
         if (cmd1->count > 0) {
@@ -463,10 +419,6 @@ main(int argc, char** argv)
             arg_print_errors(stdout, end6, progname);
             printf("usage: %s ", progname);
             arg_print_syntax(stdout, argtable6, "\n");
-        } else if (cmd7->count > 0) {
-            arg_print_errors(stdout, end7, progname);
-            printf("usage: %s ", progname);
-            arg_print_syntax(stdout, argtable7, "\n");
         } else {
             /* no correct cmd literals were given, so we cant presume which syntax was intended */
             printf("%s: missing command.\n",progname);
@@ -476,7 +428,6 @@ main(int argc, char** argv)
             printf("usage 4: %s ", progname);  arg_print_syntax(stdout, argtable4, "\n");
             printf("usage 5: %s ", progname);  arg_print_syntax(stdout, argtable5, "\n");
             printf("usage 6: %s ", progname);  arg_print_syntax(stdout, argtable6, "\n");
-            printf("usage 7: %s ", progname);  arg_print_syntax(stdout, argtable7, "\n");
         }
     }
 
@@ -486,7 +437,6 @@ main(int argc, char** argv)
     arg_freetable(argtable4, sizeof(argtable4) / sizeof(argtable4[0]));
     arg_freetable(argtable5, sizeof(argtable5) / sizeof(argtable5[0]));
     arg_freetable(argtable6, sizeof(argtable6) / sizeof(argtable6[0]));
-    arg_freetable(argtable7, sizeof(argtable7) / sizeof(argtable7[0]));
 
     return exitcode;
 }

c/tests/test_stats.c

@@ -29,22 +29,6 @@
 #include <stdlib.h>
 #include <float.h>
 
-static tsk_size_t
-get_max_site_mutations(tsk_treeseq_t *ts)
-{
-    int ret;
-    tsk_id_t j;
-    tsk_size_t max_mutations = 0;
-    tsk_site_t site;
-
-    for (j = 0; j < (tsk_id_t) tsk_treeseq_get_num_sites(ts); j++) {
-        ret = tsk_treeseq_get_site(ts, j, &site);
-        CU_ASSERT_EQUAL_FATAL(ret, 0);
-        max_mutations = TSK_MAX(max_mutations, site.mutations_length);
-    }
-    return max_mutations;
-}
-
 static bool
 multi_mutations_exist(tsk_treeseq_t *ts, tsk_id_t start, tsk_id_t end)
 {
@@ -227,36 +211,6 @@ verify_ld(tsk_treeseq_t *ts)
     free(num_site_mutations);
 }
 
-static void
-verify_pairwise_diversity(tsk_treeseq_t *ts)
-{
-    int ret;
-    size_t num_samples = tsk_treeseq_get_num_samples(ts);
-    tsk_id_t *samples;
-    uint32_t j;
-    double pi;
-    tsk_size_t max_site_mutations = get_max_site_mutations(ts);
-
-    ret = tsk_treeseq_get_pairwise_diversity(ts, NULL, 0, &pi);
-    CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_BAD_PARAM_VALUE);
-    ret = tsk_treeseq_get_pairwise_diversity(ts, NULL, 1, &pi);
-    CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_BAD_PARAM_VALUE);
-    ret = tsk_treeseq_get_pairwise_diversity(ts, NULL, num_samples + 1, &pi);
-    CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_BAD_PARAM_VALUE);
-
-    samples = tsk_treeseq_get_samples(ts);
-
-    for (j = 2; j < num_samples; j++) {
-        ret = tsk_treeseq_get_pairwise_diversity(ts, samples, j, &pi);
-        if (max_site_mutations <= 1) {
-            CU_ASSERT_EQUAL_FATAL(ret, 0);
-            CU_ASSERT_TRUE_FATAL(pi >= 0);
-        } else {
-            CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_ONLY_INFINITE_SITES);
-        }
-    }
-}
-
 /* FIXME: this test is weak and should check the return value somehow.
  * We should also have simplest and single tree tests along with separate
  * tests for the error conditions. This should be done as part of the general
@@ -1023,17 +977,6 @@ test_single_tree_ld(void)
     tsk_treeseq_free(&ts);
 }
 
-static void
-test_single_tree_pairwise_diversity(void)
-{
-    tsk_treeseq_t ts;
-
-    tsk_treeseq_from_text(&ts, 1, single_tree_ex_nodes, single_tree_ex_edges,
-            NULL, single_tree_ex_sites, single_tree_ex_mutations, NULL, NULL);
-    verify_pairwise_diversity(&ts);
-    tsk_treeseq_free(&ts);
-}
-
 static void
 test_single_tree_mean_descendants(void)
 {
@@ -1095,17 +1038,6 @@ test_paper_ex_ld(void)
     tsk_treeseq_free(&ts);
 }
 
-static void
-test_paper_ex_pairwise_diversity(void)
-{
-    tsk_treeseq_t ts;
-
-    tsk_treeseq_from_text(&ts, 10, paper_ex_nodes, paper_ex_edges,
-            NULL, paper_ex_sites, paper_ex_mutations, paper_ex_individuals, NULL);
-    verify_pairwise_diversity(&ts);
-    tsk_treeseq_free(&ts);
-}
-
 static void
 test_paper_ex_mean_descendants(void)
 {
@@ -1173,26 +1105,23 @@ test_paper_ex_diversity(void)
     tsk_treeseq_t ts;
     tsk_id_t samples[] = {0, 1, 2, 3};
     tsk_size_t sample_set_sizes = 4;
-    double pi1, pi2;
+    double pi;
     int ret;
 
     tsk_treeseq_from_text(&ts, 10, paper_ex_nodes, paper_ex_edges,
             NULL, paper_ex_sites, paper_ex_mutations, paper_ex_individuals, NULL);
 
     ret = tsk_treeseq_diversity(&ts, 1, &sample_set_sizes, samples, 0, NULL,
-            &pi1, TSK_STAT_SITE);
+            &pi, TSK_STAT_SITE);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
-    /* This function will probably be removed, but it's a handy test for now */
-    ret = tsk_treeseq_get_pairwise_diversity(&ts, samples, 4, &pi2);
-    CU_ASSERT_EQUAL_FATAL(ret, 0);
-    CU_ASSERT_DOUBLE_EQUAL_FATAL(pi1, pi2, 1e-6);
+    CU_ASSERT_DOUBLE_EQUAL_FATAL(pi, 1.5, 1e-6);
 
     /* A sample set size of 1 leads to NaN */
     sample_set_sizes = 1;
     ret = tsk_treeseq_diversity(&ts, 1, &sample_set_sizes, samples, 0, NULL,
-            &pi1, TSK_STAT_SITE);
+            &pi, TSK_STAT_SITE);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
-    CU_ASSERT(isnan((float) pi1));
+    CU_ASSERT(isnan((float) pi));
 
     tsk_treeseq_free(&ts);
 }
@@ -1681,17 +1610,6 @@ test_nonbinary_ex_ld(void)
     tsk_treeseq_free(&ts);
 }
 
-static void
-test_nonbinary_ex_pairwise_diversity(void)
-{
-    tsk_treeseq_t ts;
-
-    tsk_treeseq_from_text(&ts, 100, nonbinary_ex_nodes, nonbinary_ex_edges,
-            NULL, nonbinary_ex_sites, nonbinary_ex_mutations, NULL, NULL);
-    verify_pairwise_diversity(&ts);
-    tsk_treeseq_free(&ts);
-}
-
 static void
 test_nonbinary_ex_mean_descendants(void)
 {
@@ -1756,15 +1674,13 @@ main(int argc, char **argv)
         {"test_empty_ts_afs", test_empty_ts_afs},
 
         {"test_single_tree_ld", test_single_tree_ld},
-        {"test_single_tree_pairwise_diversity", test_single_tree_pairwise_diversity},
         {"test_single_tree_mean_descendants", test_single_tree_mean_descendants},
         {"test_single_tree_genealogical_nearest_neighbours",
             test_single_tree_genealogical_nearest_neighbours},
         {"test_single_tree_general_stat", test_single_tree_general_stat},
         {"test_single_tree_general_stat_errors", test_single_tree_general_stat_errors},
 
         {"test_paper_ex_ld", test_paper_ex_ld},
-        {"test_paper_ex_pairwise_diversity", test_paper_ex_pairwise_diversity},
         {"test_paper_ex_mean_descendants", test_paper_ex_mean_descendants},
         {"test_paper_ex_genealogical_nearest_neighbours",
             test_paper_ex_genealogical_nearest_neighbours},
@@ -1798,7 +1714,6 @@ main(int argc, char **argv)
         {"test_paper_ex_afs", test_paper_ex_afs},
 
         {"test_nonbinary_ex_ld", test_nonbinary_ex_ld},
-        {"test_nonbinary_ex_pairwise_diversity", test_nonbinary_ex_pairwise_diversity},
         {"test_nonbinary_ex_mean_descendants", test_nonbinary_ex_mean_descendants},
         {"test_nonbinary_ex_genealogical_nearest_neighbours",
             test_nonbinary_ex_genealogical_nearest_neighbours},

c/tskit/trees.c

@@ -566,66 +566,8 @@ tsk_treeseq_is_sample(tsk_treeseq_t *self, tsk_id_t u)
     return ret;
 }
 
-/* Accessors for records */
+/* Stats functions */
 
-int TSK_WARN_UNUSED
-tsk_treeseq_get_pairwise_diversity(tsk_treeseq_t *self,
-    tsk_id_t *samples, size_t num_samples, double *pi)
-{
-    int ret = 0;
-    tsk_tree_t *tree = NULL;
-    double result, denom, n, count;
-    tsk_site_t *sites;
-    tsk_size_t j, k, num_sites;
-
-    if (num_samples < 2 || num_samples > self->num_samples) {
-        ret = TSK_ERR_BAD_PARAM_VALUE;
-        goto out;
-    }
-    n = (double) num_samples;
-    tree = malloc(sizeof(tsk_tree_t));
-    if (tree == NULL) {
-        ret = TSK_ERR_NO_MEMORY;
-        goto out;
-    }
-    ret = tsk_tree_init(tree, self, TSK_SAMPLE_COUNTS);
-    if (ret != 0) {
-        goto out;
-    }
-    ret = tsk_tree_set_tracked_samples(tree, num_samples, samples);
-    if (ret != 0) {
-        goto out;
-    }
-    /* Allocation done; move onto main algorithm. */
-    result = 0.0;
-    for (ret = tsk_tree_first(tree); ret == 1; ret = tsk_tree_next(tree)) {
-        ret = tsk_tree_get_sites(tree, &sites, &num_sites);
-        if (ret != 0) {
-            goto out;
-        }
-        for (j = 0; j < num_sites; j++) {
-            if (sites[j].mutations_length != 1) {
-                ret = TSK_ERR_ONLY_INFINITE_SITES;
-                goto out;
-            }
-            for (k = 0; k < sites[j].mutations_length; k++) {
-                count = (double) tree->num_tracked_samples[sites[j].mutations[k].node];
-                result += count * (n - count);
-            }
-        }
-    }
-    if (ret != 0) {
-        goto out;
-    }
-    denom = (n * (n - 1)) / 2.0;
-    *pi = result / denom;
-out:
-    if (tree != NULL) {
-        tsk_tree_free(tree);
-        free(tree);
-    }
-    return ret;
-}
 
 #define GET_2D_ROW(array, row_len, row) (array + (((size_t) (row_len)) * (size_t) row))
 

c/tskit/trees.h

@@ -196,8 +196,6 @@ int tsk_treeseq_simplify(tsk_treeseq_t *self, tsk_id_t *samples,
 /* TODO do these belong in trees or stats? They should probably be in stats.
  * Keep them here for now until we figure out the correct interface.
  */
-int tsk_treeseq_get_pairwise_diversity(tsk_treeseq_t *self,
-    tsk_id_t *samples, size_t num_samples, double *pi);
 int tsk_treeseq_genealogical_nearest_neighbours(tsk_treeseq_t *self,
         tsk_id_t *focal, size_t num_focal,
         tsk_id_t **reference_sets, size_t *reference_set_size, size_t num_reference_sets,