docs: def calculate_top_down_likelihood(tree, character, frequencies, sf, kappa=None, model=F81):
      """
      Calculates the top-down likelihood for the given tree.
      The likelihood for each node is stored in the corresponding feature,
      given by get_personalised_feature_name(feature, TD_LH).
  
      To calculate the top-down likelihood of a node, we assume that the tree is rooted in this node
      and combine the likelihoods of the “up-subtrees”,
      e.g. to calculate the top-down likelihood of a node N1 being in a state i,
      given that its parent node is P and its brother node is N2, we imagine that the tree is re-rooted in N1,
      therefore P becoming the child of N1, and N2 its grandchild.
      We then calculate the bottom-up likelihood from the P subtree:
      L_top_down(N1, i) = \sum_j P(i -> j, dist(N1, P)) * L_top_down(P) * \sum_k P(j -> k, dist(N2, P)) * L_bottom_up (N2).
  
      For the root node we assume its top-down likelihood to be 1 for all the states.
  
      :param model: model of character evolution
      :type model: str
      :param sf: scaling factor
      :type sf: float
      :param character: character whose ancestral state likelihood is being calculated
      :type character: str
      :param tree: tree of interest (with bottom-up likelihood pre-calculated)
      :type tree: ete3.Tree
      :param frequencies: state frequencies
      :type frequencies: numpy.array
      :return: void, stores the node top-down likelihoods in the get_personalised_feature_name(feature, TD_LH) feature.
      """
  
      lh_feature = get_personalized_feature_name(character, TD_LH)
      lh_sf_feature = get_personalized_feature_name(character, TD_LH_SF)
      bu_lh_feature = get_personalized_feature_name(character, BU_LH)
      bu_lh_sf_feature = get_personalized_feature_name(character, BU_LH_SF)
  
      get_pij = get_pij_method(model, frequencies, kappa)
      for node in tree.traverse('preorder'):
          if node.is_root():
              node.add_feature(lh_feature, np.ones(len(frequencies), np.float64))
              node.add_feature(lh_sf_feature, 0)
              continue
  
          parent = node.up
          parent_bu_likelihood = getattr(parent, bu_lh_feature)
  
          node_pjis = np.transpose(get_pij(node.dist * sf))
          node_contribution = getattr(node, bu_lh_feature).dot(node_pjis)
  
          parent_likelihood = getattr(parent, lh_feature) * parent_bu_likelihood
          parent_likelihood[np.nonzero(parent_likelihood)] /= node_contribution[np.nonzero(parent_likelihood)]
          factors = getattr(parent, lh_sf_feature) + getattr(parent, bu_lh_sf_feature) - getattr(node, bu_lh_sf_feature)
  
          td_likelihood = parent_likelihood.dot(node_pjis)
          factors += rescale(td_likelihood, fraction_of_limit=len(node.children) if not node.is_leaf() else 1)
  
          node.add_feature(lh_feature, td_likelihood)
          node.add_feature(lh_sf_feature, factors)