Fix bug in pair coal counts when window breakpoint is in missing interval

Author: nspope

c/tskit/trees.c

@@ -9586,8 +9586,9 @@ tsk_treeseq_pair_coalescence_stat(const tsk_treeseq_t *self, tsk_size_t num_samp
                 window_span = windows[w + 1] - windows[w] - missing_span;
                 missing_span = 0.0;
                 if (num_edges == 0) {
+                    /* missing interval, so remove overcounted missing span */
                     remaining_span = right - windows[w + 1];
-                    window_span -= remaining_span;
+                    window_span += remaining_span;
                     missing_span += remaining_span;
                 }
                 for (i = 0; i < (tsk_id_t) num_set_indexes; i++) {

python/CHANGELOG.rst

@@ -7,6 +7,12 @@
 - ``TreeSequence.map_to_vcf_model`` now also returns the transformed positions and
   contig length. (:user:`benjeffery`, :pr:`XXXX`, :issue:`3173`)
 
+**Bugfixes**
+
+- Fix bug in ``TreeSequence.pair_coalescence_counts`` when ``span_normalise=True``
+  and a window breakpoint falls within an internal missing interval.
+  (:user:`nspope`, :pr:`3176`, :issue:`3175`)
+
 --------------------
 [0.6.4] - 2025-05-21
 --------------------

python/tests/test_coalrate.py

@@ -1294,9 +1294,9 @@ def test_span_normalise(self):
         proto = proto_pair_coalescence_counts(ts, windows=windows, span_normalise=False)
         np.testing.assert_allclose(proto, check)
 
-    def test_span_normalise_with_missing(self):
+    def test_span_normalise_with_missing_flanks(self):
         """
-        test case where span is normalised and there are intervals without trees
+        test case where span is normalised and there are flanking intervals without trees
         """
         ts = self.example_ts()
         missing = np.array([[0.0, 0.1], [0.8, 1.0]]) * ts.sequence_length
@@ -1313,6 +1313,30 @@ def test_span_normalise_with_missing(self):
         proto = proto_pair_coalescence_counts(ts, windows=windows, span_normalise=True)
         np.testing.assert_allclose(proto, check)
 
+    def test_span_normalise_with_missing_interior(self):
+        """
+        test that span normalisation correctly calculates internal missing data
+        """
+        ts = msprime.sim_ancestry(samples=1, discrete_genome=False)
+        missing_interval = np.array([[0.3, 0.6]]) * ts.sequence_length
+        windows = np.array([0.0, 0.31, 1.0]) * ts.sequence_length
+        time_windows = np.array([0.0, np.inf])
+        ts = ts.delete_intervals(missing_interval)
+        check = np.ones(windows.size - 1)
+        implm = ts.pair_coalescence_counts(
+            windows=windows,
+            time_windows=time_windows,
+            span_normalise=True,
+        ).flatten()
+        np.testing.assert_array_almost_equal(implm, check)
+        proto = proto_pair_coalescence_counts(
+            ts,
+            windows=windows,
+            time_windows=time_windows,
+            span_normalise=True,
+        ).flatten()
+        np.testing.assert_array_almost_equal(proto, check)
+
     def test_empty_windows(self):
         """
         test that windows without nodes contain zeros