Skip to content

Commit 461cec4

Browse files
nspopenspope
committed
Debugging
1 parent a7be352 commit 461cec4

1 file changed

Lines changed: 0 additions & 249 deletions

File tree

‎c/tskit/trees.c‎

Lines changed: 0 additions & 249 deletions
Original file line numberDiff line numberDiff line change
@@ -8157,254 +8157,5 @@ tsk_treeseq_pair_coalescence_stat(const tsk_treeseq_t *self, tsk_size_t num_samp
81578157
const double *windows, tsk_flags_t options, double *result)
81588158
{
81598159
int ret = 0;
8160-
double left, right, remaining_span, window_span;
8161-
tsk_id_t e, p, c, u, v, w, i, j, k;
8162-
double x;
8163-
tsk_size_t total_samples, num_outputs;
8164-
tsk_tree_position_t tree_pos;
8165-
const tsk_table_collection_t *tables = self->tables;
8166-
const tsk_size_t num_nodes = tables->nodes.num_rows;
8167-
const double sequence_length = tables->sequence_length;
8168-
tsk_id_t *nodes_sample_set = NULL;
8169-
tsk_id_t *nodes_parent = NULL;
8170-
double *nodes_sample = NULL;
8171-
double *sample_count = NULL;
8172-
double *coalescing_pairs = NULL;
8173-
double *nodes_weight = NULL;
8174-
double *outside = NULL;
8175-
8176-
/* used as row pointers */
8177-
double *inside = NULL;
8178-
double *weight = NULL;
8179-
double *output = NULL;
8180-
double *above = NULL;
8181-
double *below = NULL;
8182-
double *state = NULL;
8183-
double *pairs = NULL;
8184-
8185-
tsk_memset(&tree_pos, 0, sizeof(tree_pos));
8186-
8187-
/* check inputs */
8188-
ret = tsk_treeseq_check_windows(self, num_windows, windows, options);
8189-
if (ret != 0) {
8190-
goto out;
8191-
}
8192-
ret = check_set_indexes(num_sample_sets, 2 * num_set_indexes, set_indexes);
8193-
if (ret != 0) {
8194-
goto out;
8195-
}
8196-
ret = tsk_treeseq_check_sample_sets(
8197-
self, num_sample_sets, sample_set_sizes, sample_sets);
8198-
if (ret != 0) {
8199-
goto out;
8200-
}
8201-
8202-
/* TODO: set defaults? */
8203-
8204-
/* map nodes to sample sets and output bins */
8205-
nodes_sample_set = tsk_malloc(num_nodes * sizeof(*nodes_sample_set));
8206-
if (nodes_sample_set == NULL) {
8207-
ret = TSK_ERR_NO_MEMORY;
8208-
goto out;
8209-
}
8210-
ret = get_sample_set_index_map(self, num_sample_sets, sample_set_sizes, sample_sets,
8211-
&total_samples, nodes_sample_set);
8212-
if (ret != 0) {
8213-
goto out;
8214-
}
8215-
num_outputs = num_nodes;
8216-
8217-
/* initialize internal state */
8218-
nodes_parent = tsk_malloc(num_nodes * sizeof(*nodes_parent));
8219-
nodes_sample = tsk_calloc(num_nodes * num_sample_sets, sizeof(*nodes_sample));
8220-
sample_count = tsk_malloc(num_nodes * num_sample_sets * sizeof(*sample_count));
8221-
coalescing_pairs
8222-
= tsk_calloc(num_outputs * num_set_indexes, sizeof(*coalescing_pairs));
8223-
nodes_weight = tsk_malloc(num_outputs * num_set_indexes * sizeof(*nodes_weight));
8224-
outside = tsk_malloc(num_set_indexes * sizeof(*outside));
8225-
if (nodes_parent == NULL || nodes_sample == NULL || sample_count == NULL
8226-
|| coalescing_pairs == NULL || nodes_weight == NULL || outside == NULL) {
8227-
ret = TSK_ERR_NO_MEMORY;
8228-
goto out;
8229-
}
8230-
8231-
for (c = 0; c < (tsk_id_t) num_nodes; c++) {
8232-
i = nodes_sample_set[c];
8233-
if (i != TSK_NULL) {
8234-
state = GET_2D_ROW(nodes_sample, num_sample_sets, c);
8235-
state[i] = 1.0;
8236-
}
8237-
nodes_parent[c] = TSK_NULL;
8238-
}
8239-
tsk_memcpy(
8240-
sample_count, nodes_sample, num_nodes * num_sample_sets * sizeof(*sample_count));
8241-
8242-
ret = tsk_tree_position_init(&tree_pos, self, 0);
8243-
if (ret != 0) {
8244-
goto out;
8245-
}
8246-
8247-
w = 0;
8248-
while (true) {
8249-
tsk_tree_position_next(&tree_pos);
8250-
if (tree_pos.index == TSK_NULL) {
8251-
break;
8252-
}
8253-
8254-
left = tree_pos.interval.left;
8255-
right = tree_pos.interval.right;
8256-
remaining_span = sequence_length - left;
8257-
8258-
for (u = tree_pos.out.start; u != tree_pos.out.stop; u++) {
8259-
e = tree_pos.out.order[u];
8260-
p = tables->edges.parent[e];
8261-
c = tables->edges.child[e];
8262-
nodes_parent[c] = TSK_NULL;
8263-
inside = GET_2D_ROW(sample_count, num_sample_sets, c);
8264-
while (p != TSK_NULL) {
8265-
v = p;
8266-
if (v != TSK_NULL) {
8267-
above = GET_2D_ROW(sample_count, num_sample_sets, p);
8268-
below = GET_2D_ROW(sample_count, num_sample_sets, c);
8269-
state = GET_2D_ROW(nodes_sample, num_sample_sets, p);
8270-
for (i = 0; i < (tsk_id_t) num_sample_sets; i++) {
8271-
outside[i] = above[i] - below[i] - state[i];
8272-
}
8273-
pairs = GET_2D_ROW(coalescing_pairs, num_set_indexes, v);
8274-
for (i = 0; i < (tsk_id_t) num_set_indexes; i++) {
8275-
j = set_indexes[2 * i];
8276-
k = set_indexes[2 * i + 1];
8277-
x = outside[j] * inside[k];
8278-
if (j != k) {
8279-
x += outside[k] * inside[j];
8280-
}
8281-
pairs[i] -= x * remaining_span;
8282-
}
8283-
}
8284-
c = p;
8285-
p = nodes_parent[c];
8286-
}
8287-
p = tables->edges.parent[e];
8288-
while (p != TSK_NULL) {
8289-
above = GET_2D_ROW(sample_count, num_sample_sets, p);
8290-
for (i = 0; i < (tsk_id_t) num_sample_sets; i++) {
8291-
above[i] -= inside[i];
8292-
}
8293-
p = nodes_parent[p];
8294-
}
8295-
}
8296-
8297-
for (u = tree_pos.in.start; u != tree_pos.in.stop; u++) {
8298-
e = tree_pos.in.order[u];
8299-
p = tables->edges.parent[e];
8300-
c = tables->edges.child[e];
8301-
nodes_parent[c] = p;
8302-
inside = GET_2D_ROW(sample_count, num_sample_sets, c);
8303-
while (p != TSK_NULL) {
8304-
above = GET_2D_ROW(sample_count, num_sample_sets, p);
8305-
for (i = 0; i < (tsk_id_t) num_sample_sets; i++) {
8306-
above[i] += inside[i];
8307-
}
8308-
p = nodes_parent[p];
8309-
}
8310-
p = tables->edges.parent[e];
8311-
while (p != TSK_NULL) {
8312-
v = p;
8313-
if (v != TSK_NULL) {
8314-
above = GET_2D_ROW(sample_count, num_sample_sets, p);
8315-
below = GET_2D_ROW(sample_count, num_sample_sets, c);
8316-
state = GET_2D_ROW(nodes_sample, num_sample_sets, p);
8317-
for (i = 0; i < (tsk_id_t) num_sample_sets; i++) {
8318-
outside[i] = above[i] - below[i] - state[i];
8319-
}
8320-
pairs = GET_2D_ROW(coalescing_pairs, num_set_indexes, v);
8321-
for (i = 0; i < (tsk_id_t) num_set_indexes; i++) {
8322-
j = set_indexes[2 * i];
8323-
k = set_indexes[2 * i + 1];
8324-
x = outside[j] * inside[k];
8325-
if (j != k) {
8326-
x += outside[k] * inside[j];
8327-
}
8328-
pairs[i] += x * remaining_span;
8329-
}
8330-
}
8331-
c = p;
8332-
p = nodes_parent[c];
8333-
}
8334-
}
8335-
8336-
/* flush windows */
8337-
while (w < (tsk_id_t) num_windows && windows[w + 1] <= right) {
8338-
remaining_span = sequence_length - windows[w + 1];
8339-
tsk_memcpy(nodes_weight, coalescing_pairs,
8340-
num_outputs * num_set_indexes * sizeof(*nodes_weight));
8341-
tsk_memset(coalescing_pairs, 0,
8342-
num_outputs * num_set_indexes * sizeof(*coalescing_pairs));
8343-
for (c = 0; c < (tsk_id_t) num_nodes; c++) {
8344-
/* TODO: better to loop over only those nodes in tree, by
8345-
* following nodes_parent up from samples; this should always
8346-
* be fine, I think, as disconnected nodes will have zero
8347-
* values anyway? */
8348-
p = nodes_parent[c];
8349-
if (p == TSK_NULL) {
8350-
continue;
8351-
}
8352-
v = p;
8353-
if (v == TSK_NULL) {
8354-
continue;
8355-
}
8356-
above = GET_2D_ROW(sample_count, num_sample_sets, p);
8357-
below = GET_2D_ROW(sample_count, num_sample_sets, c);
8358-
state = GET_2D_ROW(nodes_sample, num_sample_sets, p);
8359-
for (i = 0; i < (tsk_id_t) num_sample_sets; i++) {
8360-
outside[i] = above[i] - below[i] - state[i];
8361-
}
8362-
inside = GET_2D_ROW(sample_count, num_sample_sets, c);
8363-
weight = GET_2D_ROW(nodes_weight, num_set_indexes, v);
8364-
pairs = GET_2D_ROW(coalescing_pairs, num_set_indexes, v);
8365-
for (i = 0; i < (tsk_id_t) num_set_indexes; i++) {
8366-
j = set_indexes[2 * i];
8367-
k = set_indexes[2 * i + 1];
8368-
x = inside[j] * outside[k];
8369-
if (j != k) {
8370-
x += inside[k] * outside[j];
8371-
}
8372-
weight[i] -= x * remaining_span / 2;
8373-
pairs[i] += x * remaining_span / 2;
8374-
}
8375-
}
8376-
if (options & TSK_STAT_SPAN_NORMALISE) {
8377-
window_span = windows[w + 1] - windows[w];
8378-
for (v = 0; v < (tsk_id_t) num_outputs; v++) {
8379-
weight = GET_2D_ROW(nodes_weight, num_set_indexes, v);
8380-
for (i = 0; i < (tsk_id_t) num_set_indexes; i++) {
8381-
weight[i] /= window_span;
8382-
}
8383-
}
8384-
}
8385-
// TODO:
8386-
// for (i = 0; i < (tsk_id_t) num_set_indexes; i++) {
8387-
// reduce(i, col, nodes_weight, &result[w * row_dim * col_dim])
8388-
// };
8389-
for (v = 0; v < (tsk_id_t) num_outputs; v++) {
8390-
output = GET_3D_ROW(
8391-
result, num_outputs, num_set_indexes, (tsk_size_t) w, v);
8392-
weight = GET_2D_ROW(nodes_weight, num_set_indexes, v);
8393-
for (i = 0; i < (tsk_id_t) num_set_indexes; i++) {
8394-
output[i] = weight[i];
8395-
}
8396-
}
8397-
w += 1;
8398-
}
8399-
}
8400-
out:
8401-
tsk_tree_position_free(&tree_pos);
8402-
tsk_safe_free(nodes_sample_set);
8403-
tsk_safe_free(nodes_parent);
8404-
tsk_safe_free(nodes_sample);
8405-
tsk_safe_free(sample_count);
8406-
tsk_safe_free(coalescing_pairs);
8407-
tsk_safe_free(nodes_weight);
8408-
tsk_safe_free(outside);
84098160
return ret;
84108161
}

0 commit comments

Comments
 (0)