Efficient metadata

QtSLiM/recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py

@@ -7,16 +7,18 @@
 # Load the .trees file
 ts = tskit.load("decap.trees")    # no simplify!
 
-# Calculate tree heights, giving uncoalesced sites the maximum time
+# Calculate tree heights
 def tree_heights(ts):
     heights = np.zeros(ts.num_trees + 1)
     for tree in ts.trees():
-        if tree.num_roots > 1:  # not fully coalesced
-            heights[tree.index] = ts.metadata['SLiM']['tick']
-        else:
-            children = tree.children(tree.root)
-            real_root = tree.root if len(children) > 1 else children[0]
-            heights[tree.index] = tree.time(real_root)
+        roots = tree.roots
+        if len(roots) == 1:
+            children = tree.children(roots[0])
+            if len(children) == 1:
+                roots = children
+        root_heights = [tree.time(r) for r in roots]
+        assert len(set(root_heights)) == 1
+        heights[tree.index] = root_heights[0]
     heights[-1] = heights[-2]  # repeat the last entry for plotting with step
     return heights
 
@@ -35,5 +37,3 @@ def tree_heights(ts):
 heights = tree_heights(recap)
 plt.step(breakpoints, heights, where='post')
 plt.show()
-
-

SLiMgui/Recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py

@@ -7,16 +7,18 @@
 # Load the .trees file
 ts = tskit.load("decap.trees")    # no simplify!
 
-# Calculate tree heights, giving uncoalesced sites the maximum time
+# Calculate tree heights
 def tree_heights(ts):
     heights = np.zeros(ts.num_trees + 1)
     for tree in ts.trees():
-        if tree.num_roots > 1:  # not fully coalesced
-            heights[tree.index] = ts.metadata['SLiM']['tick']
-        else:
-            children = tree.children(tree.root)
-            real_root = tree.root if len(children) > 1 else children[0]
-            heights[tree.index] = tree.time(real_root)
+        roots = tree.roots
+        if len(roots) == 1:
+            children = tree.children(roots[0])
+            if len(children) == 1:
+                roots = children
+        root_heights = [tree.time(r) for r in roots]
+        assert len(set(root_heights)) == 1
+        heights[tree.index] = root_heights[0]
     heights[-1] = heights[-2]  # repeat the last entry for plotting with step
     return heights
 
@@ -35,5 +37,3 @@ def tree_heights(ts):
 heights = tree_heights(recap)
 plt.step(breakpoints, heights, where='post')
 plt.show()
-
-

VERSIONS

@@ -204,6 +204,7 @@ multitrait branch:
 	add a built-in "Multispecies Multitrait Phenotype ~ Time" plot for phenotype ~ time, with one line per trait across all species
 	add a built-in "Phenotype ~ Time" plot for phenotype ~ time, focusing on a single user-selected trait across subpopulations
 	clean up the way pedigree IDs get recorded, for better efficiency and parallelization
+	fix recipe 18.10 to reduce metadata accesses (see discussion at https://github.com/tskit-dev/tskit/issues/3472#issuecomment-5210734749)
 
 
 version 5.2 (Eidos version 4.2):