Skip to content

Commit

Permalink
Fix issues with C implm; refactor lowlevel tests; remove nodes mode
Browse files Browse the repository at this point in the history
  • Loading branch information
nspope committed Aug 1, 2024
1 parent 5d00273 commit 6318fad
Show file tree
Hide file tree
Showing 6 changed files with 321 additions and 252 deletions.
68 changes: 54 additions & 14 deletions c/tests/test_stats.c
Original file line number Diff line number Diff line change
Expand Up @@ -306,8 +306,7 @@ verify_divergence_matrix(tsk_treeseq_t *ts, tsk_flags_t options)

/* Check coalescence counts against naive implementation */
static void
verify_pair_coalescence_counts(tsk_treeseq_t *ts, tsk_size_t num_time_windows,
double *time_windows, tsk_flags_t options)
verify_pair_coalescence_counts(tsk_treeseq_t *ts, tsk_flags_t options)
{
int ret;
const tsk_size_t n = tsk_treeseq_get_num_samples(ts);
Expand All @@ -317,12 +316,17 @@ verify_pair_coalescence_counts(tsk_treeseq_t *ts, tsk_size_t num_time_windows,
const double *breakpoints = tsk_treeseq_get_breakpoints(ts);
const tsk_size_t P = 2;
const tsk_size_t I = P * (P + 1) / 2;
tsk_id_t sample_sets[n];
tsk_size_t sample_set_sizes[P];
tsk_id_t index_tuples[2 * I];
tsk_size_t dim = T * N * I;
double C1[dim]; //, C2[dim];
tsk_size_t i, j, k;

for (i = 0; i < n; i++) {
sample_sets[i] = samples[i];
}

for (i = 0; i < P; i++) {
sample_set_sizes[i] = 0;
}
Expand All @@ -338,11 +342,52 @@ verify_pair_coalescence_counts(tsk_treeseq_t *ts, tsk_size_t num_time_windows,
}
}

ret = tsk_treeseq_pair_coalescence_stat(ts, P, sample_set_sizes, samples, I,
index_tuples, T, breakpoints, num_time_windows, time_windows, options, C1);
double max_time = tsk_treeseq_get_max_time(ts);
double time_windows[3] = { 0.0, max_time / 2, INFINITY };

ret = tsk_treeseq_pair_coalescence_stat(ts, P, sample_set_sizes, sample_sets, I,
index_tuples, T, breakpoints, 2, time_windows, options, C1);
CU_ASSERT_EQUAL_FATAL(ret, 0);

double trunc_time_windows[3] = { max_time * 0.2, max_time * 0.5, max_time * 0.8 };
ret = tsk_treeseq_pair_coalescence_stat(ts, P, sample_set_sizes, sample_sets, I,
index_tuples, T, breakpoints, 2, trunc_time_windows, options, C1);
CU_ASSERT_EQUAL_FATAL(ret, 0);

/* TODO: check against tree by tree implementation here */
/* cover errors */
double nil_time_windows[1] = { 0.0 };
ret = tsk_treeseq_pair_coalescence_stat(ts, P, sample_set_sizes, sample_sets, I,
index_tuples, T, breakpoints, 0, nil_time_windows, options, C1);
CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_BAD_NUM_TIME_WINDOWS);

double bad_time_windows[2] = { 10.0, 0.0 };
ret = tsk_treeseq_pair_coalescence_stat(ts, P, sample_set_sizes, sample_sets, I,
index_tuples, T, breakpoints, 1, bad_time_windows, options, C1);
CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_BAD_TIME_WINDOWS);

double bad_breakpoints[2] = { breakpoints[1], 0.0 };
ret = tsk_treeseq_pair_coalescence_stat(ts, P, sample_set_sizes, sample_sets, I,
index_tuples, 1, bad_breakpoints, 1, time_windows, options, C1);
CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_BAD_WINDOWS);

index_tuples[0] = (tsk_id_t) P;
ret = tsk_treeseq_pair_coalescence_stat(ts, P, sample_set_sizes, sample_sets, I,
index_tuples, 1, breakpoints, 1, time_windows, options, C1);
CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_BAD_SAMPLE_SET_INDEX);
index_tuples[0] = 0;

tsk_size_t tmp = sample_set_sizes[0];
sample_set_sizes[0] = 0;
ret = tsk_treeseq_pair_coalescence_stat(ts, P, sample_set_sizes, sample_sets, I,
index_tuples, 1, breakpoints, 1, time_windows, options, C1);
CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_EMPTY_SAMPLE_SET);
sample_set_sizes[0] = tmp;

sample_sets[1] = 0;
ret = tsk_treeseq_pair_coalescence_stat(ts, P, sample_set_sizes, sample_sets, I,
index_tuples, 1, breakpoints, 1, time_windows, options, C1);
CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_DUPLICATE_SAMPLE);
sample_sets[1] = 1;
}

typedef struct {
Expand Down Expand Up @@ -2833,15 +2878,10 @@ test_pair_coalescence_counts(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, 0);
double max_time = tsk_treeseq_get_max_time(&ts);
double time_windows[3] = { 0.0, max_time / 2, INFINITY };
double trunc_windows[3] = { max_time * 0.2, max_time * 0.5, max_time * 0.8 };
verify_pair_coalescence_counts(&ts, 0, NULL, 0);
verify_pair_coalescence_counts(&ts, 0, NULL, TSK_STAT_SPAN_NORMALISE);
verify_pair_coalescence_counts(&ts, 2, time_windows, 0);
verify_pair_coalescence_counts(&ts, 2, time_windows, TSK_STAT_SPAN_NORMALISE);
verify_pair_coalescence_counts(&ts, 2, trunc_windows, 0);
verify_pair_coalescence_counts(&ts, 2, trunc_windows, TSK_STAT_SPAN_NORMALISE);
verify_pair_coalescence_counts(&ts, TSK_STAT_NODE);
verify_pair_coalescence_counts(&ts, TSK_STAT_NODE | TSK_STAT_SPAN_NORMALISE);
verify_pair_coalescence_counts(&ts, 0);
verify_pair_coalescence_counts(&ts, TSK_STAT_SPAN_NORMALISE);
tsk_treeseq_free(&ts);
}

Expand Down
Loading

0 comments on commit 6318fad

Please sign in to comment.