From 1a1214afa866a090fdc9114730ab682981e6810d Mon Sep 17 00:00:00 2001 From: peter Date: Fri, 14 Aug 2026 13:04:48 -0700 Subject: [PATCH 1/2] fixup of python recipe --- ...rn-in after simulation with recapitation II.py | 15 ++++++++------- ...rn-in after simulation with recapitation II.py | 15 ++++++++------- 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/QtSLiM/recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py b/QtSLiM/recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py index a56d3c0a..ea43ff2c 100644 --- a/QtSLiM/recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py +++ b/QtSLiM/recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py @@ -9,15 +9,16 @@ # Calculate tree heights, giving uncoalesced sites the maximum time def tree_heights(ts): - uncoalesced_height = ts.metadata['SLiM']['tick'] heights = np.zeros(ts.num_trees + 1) for tree in ts.trees(): - if tree.num_roots > 1: # not fully coalesced - heights[tree.index] = uncoalesced_height - 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 tree.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 diff --git a/SLiMgui/Recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py b/SLiMgui/Recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py index a56d3c0a..ea43ff2c 100644 --- a/SLiMgui/Recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py +++ b/SLiMgui/Recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py @@ -9,15 +9,16 @@ # Calculate tree heights, giving uncoalesced sites the maximum time def tree_heights(ts): - uncoalesced_height = ts.metadata['SLiM']['tick'] heights = np.zeros(ts.num_trees + 1) for tree in ts.trees(): - if tree.num_roots > 1: # not fully coalesced - heights[tree.index] = uncoalesced_height - 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 tree.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 From d6def0541bc3957313c25cc4f662c20bb69046a1 Mon Sep 17 00:00:00 2001 From: peter Date: Sat, 15 Aug 2026 09:02:18 -0700 Subject: [PATCH 2/2] bug fixup --- ...l burn-in after simulation with recapitation II.py | 11 ++++++----- ...l burn-in after simulation with recapitation II.py | 11 ++++++----- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/QtSLiM/recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py b/QtSLiM/recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py index ea43ff2c..51fd497f 100644 --- a/QtSLiM/recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py +++ b/QtSLiM/recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py @@ -16,7 +16,7 @@ def tree_heights(ts): children = tree.children(roots[0]) if len(children) == 1: roots = children - root_heights = [tree.time(r) for r in tree.roots] + 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 @@ -25,8 +25,7 @@ def tree_heights(ts): # Plot tree heights before recapitation breakpoints = list(ts.breakpoints()) heights = tree_heights(ts) -plt.step(breakpoints, heights, where='post') -plt.show() +plt.step(breakpoints, heights, where='post', label="SLiM") # Recapitate! recap = pyslim.recapitate(ts, ancestral_Ne=1e5, recombination_rate=3e-10, random_seed=1) @@ -35,7 +34,9 @@ def tree_heights(ts): # Plot the tree heights after recapitation breakpoints = list(recap.breakpoints()) heights = tree_heights(recap) -plt.step(breakpoints, heights, where='post') -plt.show() +plt.step(breakpoints, heights, where='post', label="recapitated") +plt.xlabel("genome"); plt.ylabel("time ago (generations)") +plt.savefig("tree_heights.png") + diff --git a/SLiMgui/Recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py b/SLiMgui/Recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py index ea43ff2c..51fd497f 100644 --- a/SLiMgui/Recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py +++ b/SLiMgui/Recipes/Recipe 18.10 - Adding a neutral burn-in after simulation with recapitation II.py @@ -16,7 +16,7 @@ def tree_heights(ts): children = tree.children(roots[0]) if len(children) == 1: roots = children - root_heights = [tree.time(r) for r in tree.roots] + 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 @@ -25,8 +25,7 @@ def tree_heights(ts): # Plot tree heights before recapitation breakpoints = list(ts.breakpoints()) heights = tree_heights(ts) -plt.step(breakpoints, heights, where='post') -plt.show() +plt.step(breakpoints, heights, where='post', label="SLiM") # Recapitate! recap = pyslim.recapitate(ts, ancestral_Ne=1e5, recombination_rate=3e-10, random_seed=1) @@ -35,7 +34,9 @@ def tree_heights(ts): # Plot the tree heights after recapitation breakpoints = list(recap.breakpoints()) heights = tree_heights(recap) -plt.step(breakpoints, heights, where='post') -plt.show() +plt.step(breakpoints, heights, where='post', label="recapitated") +plt.xlabel("genome"); plt.ylabel("time ago (generations)") +plt.savefig("tree_heights.png") +