Merge pull request #221 from tskit-dev/node_stats

Implement node statistics

Author: petrelharp

c/tskit/trees.c

@@ -1559,13 +1559,17 @@ tsk_treeseq_node_general_stat(tsk_treeseq_t *self,
 
     /* Set the initial conditions */
     for (j = 0; j < self->num_samples; j++) {
-        u = self->samples[j];
-        state_u = GET_2D_ROW(state, state_dim, u);
+        // total_weight is used in the next loop
         weight_u = GET_2D_ROW(sample_weights, state_dim, j);
-        memcpy(state_u, weight_u, state_dim * sizeof(*state_u));
         for (k = 0; k < state_dim; k++) {
             total_weight[k] += weight_u[k];
         }
+        u = self->samples[j];
+        state_u = GET_2D_ROW(state, state_dim, u);
+        weight_u = GET_2D_ROW(sample_weights, state_dim, j);
+        memcpy(state_u, weight_u, state_dim * sizeof(*state_u));
+    }
+    for (u = 0; u < (tsk_id_t) num_nodes; u++) {
         ret = update_node_summary(u, result_dim, node_summary, state, state_dim,
             f, f_params, total_weight, polarised);
         if (ret != 0) {

docs/stats.rst

@@ -129,7 +129,7 @@ then we'd just obtain the number of differing sites per window.
 
 And, a final note about "length": in tree sequences produced by ``msprime``
 coordinates along the sequence are **continuous**,
-so the the "lengths" used here may not correspond to distance along the genome in (say) base pairs.
+so the "lengths" used here may not correspond to distance along the genome in (say) base pairs.
 For instance, pairwise sequence divergence is usually a number between 0 and 1
 because it is the proportion of bases that differ;
 this will only be true if length is measured in base pairs
@@ -155,7 +155,7 @@ For instance, nucleotide divergence is defined for a *pair* of groups of samples
 so if you wanted to compute pairwise divergences between some groups of samples,
 you'd specify these as your ``sample_sets``.
 Then, if ``p[i]`` is the derived allele frequency in sample set ``i``,
-under the hood we (essentially) compute the divergency between sample sets ``i`` and ``j``
+under the hood we (essentially) compute the divergence between sample sets ``i`` and ``j``
 by averaging ``p[i] * (1 - p[j]) + (1 - p[i]) * p[j]`` across the genome.
 
 So, what if you
@@ -220,6 +220,20 @@ from ``windows[i]`` to ``windows[i + 1]`` (including the left but not the right
 
 The final dimension of the arrays in other cases is specified by the method.
 
+A note about **default values** and **division by zero**: 
+Under the hood, statistics computation fills in zeros everywhere, then updates these
+(since statistics are all **additive**, this makes sense).
+But now suppose that you've got a statistic that returns ``nan``
+("not a number") sometimes, like if you're taking the diversity of a sample set with only ``n=1`` sample,
+which involves dividing by ``n * (n - 1)``.
+Usually, you'll just get ``nan`` everywhere that the division by zero happens.
+But there's a couple of caveats.
+For ``site`` statistics, any windows without any sites in them never get touched,
+so they will have a value of 0.
+For ``branch`` statistics, any windows with **no branches** will similarly remain 0.
+That said, you should **not** rely on the specific behavior of whether ``0`` or ``nan`` is returned
+for "empty" cases like these: it is subject to change.
+
 
 .. Commenting these out for now as they are duplicates of the methods in the TreeSequence
    and sphinx is unhappy.