fix: reset small negative branch lengths to 0
@@ -465,19 +465,22 @@ class ClockTree(TreeAnc):
# NOTE the Lx distribution is the likelihood, given the position of the parent
# (Lx.x = parent position, Lx.y = LH of the node_pos given Lx.x,
# the length of the branch corresponding to the most likely
- # subtree is node.Cx(node.time_before_present))
+ # subtree is node.Cx(node.up.time_before_present))
# subtree_LH = node.joint_pos_Lx(node.up.time_before_present)
node.branch_length = node.joint_pos_Cx(max(node.joint_pos_Cx.xmin,
- node.up.time_before_present)+ttconf.TINY_NUMBER)
- node.time_before_present = node.up.time_before_present - node.branch_length
- node.clock_length = node.branch_length
+ node.up.time_before_present))
- # just sanity check, should never happen:
- if node.branch_length < 0 or node.time_before_present < 0:
- if node.branch_length<0 and node.branch_length>-ttconf.TINY_NUMBER:
- self.logger("ClockTree - Joint reconstruction: correcting rounding error of %s"%node.name, 4)
+ # clean up tiny negative branch length, warn against bigger ones.
+ if node.branch_length<0:
+ if node.branch_length>-2*ttconf.TINY_NUMBER:
+ self.logger(f"ClockTree - Joint reconstruction: correcting rounding error of {node.name} bl={node.branch_length:1.2e}", 4)
+ else:
+ self.logger(f"ClockTree - Joint reconstruction: NEGATIVE BRANCH LENGTH {node.name} bl={node.branch_length:1.2e}", 2, warn=True)
node.branch_length = 0
+ node.time_before_present = node.up.time_before_present - node.branch_length
+ node.clock_length = node.branch_length
+
self.tree.positional_joint_LH = self.timetree_likelihood()
# cleanup, if required
if not self.debug:
