Felsensteins algorithms

felsenstein.py

def felsenstein(tree, states):
  # Initialize a dictionary to store the probabilities at each node
  probabilities = {}

  # Traverse the tree in postorder (i.e. children first, parent last)
  for node in postorder_traverse(tree):
    # If the node is a leaf, its probability is 1
    if is_leaf(node):
      probabilities[node] = 1
    else:
      # If the node is not a leaf, its probability is the product of
      # the probabilities of its children times the transition probabilities
      # between the states at the node and its children
      left_child = left_child(node)
      right_child = right_child(node)
      probabilities[node] = probabilities[left_child] * probabilities[right_child] * transition_probability(states[node], states[left_child], states[right_child])

  # Return the probability of the root node, which is the overall probability of the tree
  return probabilities[root(tree)]