Set root time values to something reasonable

Author: hyanwong

CHANGELOG.rst

@@ -1,3 +1,22 @@
+********************
+[0.2.2] - 2021-0X-XX
+********************
+
+**Bugfixes**:
+
+- Mutations at non-inference sites are now guaranteed to be fully parsimonious.
+  Previous versions required a mutation above the root when the input ancestral state
+  disagreed with the ancestral state produced by the parsimony algorithm. Now fixed by
+  using the new map_mutations code from tskit 0.3.7 (:pr:`557`, :user:`hyanwong`)
+
+**New Features**:
+
+**Breaking changes**:
+
+- Oldest nodes in a standard inferred tree sequence are no longer set to frequencies ~2
+  and ~3 (i.e. 2 or 3 times as old as all the other nodes), but are spaced above the
+  others by the mean time between unique ancestor ages (:pr:`485`, :user:`hyanwong`)
+
 ********************
 [0.2.1] - 2021-05-26
 ********************

tests/test_inference.py

@@ -1515,6 +1515,26 @@ def test_random_data(self):
         ancestor_data = tsinfer.generate_ancestors(sample_data)
         self.verify_ancestors(sample_data, ancestor_data)
 
+    def test_oldest_ancestors(self):
+        m = 50
+        G, positions = get_random_data_example(20, m)
+        sample_data = tsinfer.SampleData(sequence_length=m)
+        for genotypes, position in zip(G, positions):
+            sample_data.add_site(position, genotypes)
+        sample_data.finalise()
+        ancestor_data = tsinfer.generate_ancestors(sample_data)
+        assert ancestor_data.num_ancestors > 2
+        times = ancestor_data.ancestors_time[:]
+        unique_times = np.unique(times)
+        ultimate_root_time = unique_times[-1]
+        root_time = unique_times[-2]
+        oldest_non_root = unique_times[-3]
+        assert np.sum(times == root_time) == 1  # root ancestor at unique time
+        assert np.sum(times == ultimate_root_time) == 1  # ultimate anc at unique time
+        expected_time_diff = oldest_non_root / len(unique_times[:-2])
+        assert np.isclose(ultimate_root_time - root_time, expected_time_diff)
+        assert np.isclose(root_time - oldest_non_root, expected_time_diff)
+
 
 class TestAncestorsTreeSequence:
     """

tsinfer/inference.py

@@ -946,8 +946,10 @@ def run(self):
             logger.info(f"Starting build for {self.num_ancestors} ancestors")
             progress = self.progress_monitor.get("ga_generate", self.num_ancestors)
             a = np.zeros(self.num_sites, dtype=np.int8)
-            root_time = max(self.timepoint_to_epoch.keys()) + 1
-            ultimate_ancestor_time = root_time + 1
+            root_time = max(self.timepoint_to_epoch.keys())
+            av_timestep = root_time / len(self.timepoint_to_epoch)
+            root_time += av_timestep  # Add a root a bit older than the oldest ancestor
+            ultimate_ancestor_time = root_time + av_timestep
             # Add the ultimate ancestor. This is an awkward hack really; we don't
             # ever insert this ancestor. The only reason to add it here is that
             # it makes sure that the ancestor IDs we have in the ancestor file are