jeetsukumaran / DendroPy

A Python library for phylogenetic scripting, simulation, data processing and manipulation.
https://pypi.org/project/DendroPy/.
BSD 3-Clause "New" or "Revised" License
210 stars 61 forks source link

treesim.birth_death_tree() raises division by zero error when used to extend an existing tree #188

Closed Androstane closed 6 months ago

Androstane commented 6 months ago

Hello,

Thanks for developing the tool! When I was using Treesim to simulate a birth-death tree with multiple stages, an error had occurred. Specifically, when treesim.birth_death_tree() raises division by zero error when the tree argument is supplied. Here's a minimal example:

import random
import dendropy
from dendropy.simulate import treesim

tree = dendropy.Tree()
tree = treesim.birth_death_tree(1.0, 0.5, max_time=random.randint(1, 8), tree=tree, repeat_until_success=True)

When running the code, the following error would arise:

ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-2-a488d6f12df8> in <module>
      3 from dendropy.simulate import treesim
      4 tree = dendropy.Tree()
----> 5 tree = treesim.birth_death_tree(1.0,0.5,max_time=random.randint(1,8),tree=tree,repeat_until_success=True)

[~/anaconda3/lib/python3.8/site-packages/dendropy/model/birthdeath.py](https://vscode-remote+ssh-002dremote-002b7b22686f73744e616d65223a226b696d7572612e63732e726963652e656475222c2275736572223a22796c313833227d.vscode-resource.vscode-cdn.net/home/yl183/sampled_ancestor/preprocessing/~/anaconda3/lib/python3.8/site-packages/dendropy/model/birthdeath.py) in birth_death_tree(birth_rate, death_rate, birth_rate_sd, death_rate_sd, **kwargs)
    329 
    330         # waiting time based on above probability
--> 331         waiting_time = rng.expovariate(rate_of_any_event)
    332 
    333         if ( (gsa_ntax is not None)

[~/anaconda3/lib/python3.8/random.py](https://vscode-remote+ssh-002dremote-002b7b22686f73744e616d65223a226b696d7572612e63732e726963652e656475222c2275736572223a22796c313833227d.vscode-resource.vscode-cdn.net/home/yl183/sampled_ancestor/preprocessing/~/anaconda3/lib/python3.8/random.py) in expovariate(self, lambd)
    493         # we use 1-random() instead of random() to preclude the
    494         # possibility of taking the log of zero.
--> 495         return -_log(1.0 - self.random())[/](https://vscode-remote+ssh-002dremote-002b7b22686f73744e616d65223a226b696d7572612e63732e726963652e656475222c2275736572223a22796c313833227d.vscode-resource.vscode-cdn.net/)lambd
    496 
    497 ## -------------------- von Mises distribution --------------------

ZeroDivisionError: float division by zero
mmore500 commented 6 months ago

Thanks for the report! Do you happen to know if this error occurs when working with a non-empty tree?

As a pointer for myself in returning to this issue, here's a link to this method's documentation https://jeetsukumaran.github.io/DendroPy/library/birthdeath.html#dendropy.model.birthdeath.birth_death_tree

Androstane commented 6 months ago

It seems to me that the error would occur when working with empty and non-empty trees. A quick fix I have found is to remove the condition check at:

https://github.com/jeetsukumaran/DendroPy/blob/a8ee43179bd11822e4368d56371e7d6dd85b0e07/src/dendropy/model/birthdeath.py#L703

as with it, the list of the extant tips would always be empty even if a non-empty tree is supplied, leading to a division by zero issues later when sampling the event node from the list of extant tips. It is probably because, by default, the is_extinct attribute is not added to the node. I am not sure if this fix would lead to other unexpected performance in the code though.

mmore500 commented 6 months ago

Hi @Androstane --- I believe I found the root cause of the issue and fixed it. The fix will appear in DendroPy 4.6.2 which will release shortly!

Thank you for the report, it is much appreciated :+1:

mmore500 commented 6 months ago

Quick update, looks like we need to get the maintainers to add verified email addresses before publishing due to a new PyPI policy. Hope to get the release out soon.