code: def felsensteins(binaryTree, subMatrices, ancestorProbs, leaves, alphabetSize):
      
      l = {}
      def upPass(binaryTree):
          if binaryTree.internal: #is internal binaryTree
              i = branchUp(binaryTree.left)
              j = branchUp(binaryTree.right)
              k = multiplyWV(i, j, alphabetSize)
              l[binaryTree.traversalID.mid] = (k, i, j)
              return k
          l[binaryTree.traversalID.mid] = leaves[binaryTree.traversalID.leafNo]
          return leaves[binaryTree.traversalID.leafNo]
      def downPass(binaryTree, ancestorProbs):
          if binaryTree.internal: #is internal binaryTree
              i = l[binaryTree.traversalID.mid]
              l[binaryTree.traversalID.mid] = multiplyWV(ancestorProbs, i[0], alphabetSize)
              branchDown(binaryTree.left, multiplyWV(ancestorProbs, i[2], alphabetSize))
              branchDown(binaryTree.right, multiplyWV(ancestorProbs, i[1], alphabetSize))
      def branchUp(binaryTree):
          return transformByDistance(upPass(binaryTree), subMatrices[binaryTree.traversalID.mid], alphabetSize)
      def branchDown(binaryTree, ancestorProbs):
          downPass(binaryTree, transformByDistance(ancestorProbs, subMatrices[binaryTree.traversalID.mid], alphabetSize))
      upPass(binaryTree)
      downPass(binaryTree, ancestorProbs)
      return l
docs: calculates the un-normalised probabilties of each non-gap residue position