+++ /dev/null
-/* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
-
-/* This a mix of tree functions and data-structures from
- * Bob Edgar's Muscle (version 3.7) ported to pure C, topped up with
- * some of our own stuff.
- *
- * Used files: phy.cpp, tree.h, phytofile.cpp and phyfromclust.cpp
- *
- * Muscle's code is public domain and so is this code here.
-
- * From http://www.drive5.com/muscle/license.htm:
- * """
- * MUSCLE is public domain software
- *
- * The MUSCLE software, including object and source code and
- * documentation, is hereby donated to the public domain.
- *
- * Disclaimer of warranty
- * THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
- * EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION IMPLIED
- * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * """
- *
- */
-
-/*
- * RCS $Id: muscle_tree.c 230 2011-04-09 15:37:50Z andreas $
- */
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <limits.h>
-#include <assert.h>
-#include <ctype.h>
-
-#include "util.h"
-#include "log.h"
-#include "muscle_tree.h"
-
-
-static const double VERY_NEGATIVE_DOUBLE = -9e29;
-/*const double dInsane = VERY_NEGATIVE_DOUBLE;*/
-static const double dInsane = -9e29;
-static const unsigned uInsane = 8888888;
-
-typedef enum
-{
- NTT_Unknown,
-
- /* Returned from Tree::GetToken: */
- NTT_Lparen,
- NTT_Rparen,
- NTT_Colon,
- NTT_Comma,
- NTT_Semicolon,
- NTT_String,
-
- /* Following are never returned from Tree::GetToken: */
- NTT_SingleQuotedString,
- NTT_DoubleQuotedString,
- NTT_Comment
-} NEWICK_TOKEN_TYPE;
-
-
-static void
-InitCache(uint uCacheCount, tree_t *tree);
-static void
-TreeZero(tree_t *tree);
-static uint
-GetNeighborCount(unsigned uNodeIndex, tree_t *tree);
-static bool
-IsEdge(unsigned uNodeIndex1, unsigned uNodeIndex2, tree_t *tree);
-static bool
-HasEdgeLength(uint uNodeIndex1, uint uNodeIndex2, tree_t *tree);
-static void
-TreeToFileNodeRooted(tree_t *tree, uint m_uRootNodeIndex, FILE *fp);
-static void
-ValidateNode(uint uNodeIndex, tree_t *tree);
-static void
-AssertAreNeighbors(unsigned uNodeIndex1, unsigned uNodeIndex2, tree_t *tree);
-static void
-ExpandCache(tree_t *tree);
-static void
-TreeCreateRooted(tree_t *tree);
-static bool
-GetGroupFromFile(FILE *fp, uint uNodeIndex, double *ptrdEdgeLength, tree_t *tree);
-static NEWICK_TOKEN_TYPE
-GetToken(FILE *fp, char szToken[], uint uBytes);
-/* stuff from textfile.cpp */
-static void
-FileSkipWhite(FILE *fp);
-static bool
-FileSkipWhiteX(FILE *fp);
-static void
-SetLeafName(uint uNodeIndex, const char *ptrName, tree_t *tree);
-uint
-AppendBranch(tree_t *tree, uint uExistingLeafIndex);
-static void
-SetEdgeLength(uint uNodeIndex1, uint uNodeIndex2,
- double dLength, tree_t *tree);
-static uint
-UnrootFromFile(tree_t *tree);
-uint
-GetNeighbor(uint uNodeIndex, uint uNeighborSubscript, tree_t *tree);
-static void
-InitNode(tree_t *prTree, uint uNodeIndex);
-
-
-/**
- * @param[in] uNodeIndex
- * node to examine
- * @param[in] tree
- * tree to examine
- * @return id of left node
- * @note called GetRight in Muscle3.7
- */
-uint
-GetLeft(uint uNodeIndex, tree_t *tree)
-{
- assert(NULL != tree);
- assert(tree->m_bRooted && uNodeIndex < tree->m_uNodeCount);
- return tree->m_uNeighbor2[uNodeIndex];
-}
-/*** end: GetLeft ***/
-
-
-
-/**
- * @param[in] uNodeIndex
- * node to examine
- * @param[in] tree
- * tree to examine
- * @return id of right node
- * @note called GetRight in Muscle3.7
- */
-uint
-GetRight(uint uNodeIndex, tree_t *tree)
-{
- assert(NULL != tree);
- assert(tree->m_bRooted && uNodeIndex < tree->m_uNodeCount);
- return tree->m_uNeighbor3[uNodeIndex];
-}
-/*** end: GetRight ***/
-
-
-
-/**
- * @param[in] uNodeIndex
- * node to examine
- * @param[in] tree
- * tree to examine
- * @return leaf id of current node
- */
-uint
-GetLeafId(uint uNodeIndex, tree_t *tree)
-{
- assert(NULL != tree);
- assert(uNodeIndex < tree->m_uNodeCount);
- assert(IsLeaf(uNodeIndex, tree));
-
- return tree->m_Ids[uNodeIndex];
-}
-/*** end: GetLeafId ***/
-
-
-
-/**
- * @note originally called GetLeafName
- *
- */
-char *
-GetLeafName(unsigned uNodeIndex, tree_t *tree)
-{
- assert(NULL != tree);
- assert(uNodeIndex < tree->m_uNodeCount);
- assert(IsLeaf(uNodeIndex, tree));
- return tree->m_ptrName[uNodeIndex];
-}
-/*** end: GetLeafName ***/
-
-
-
-/**
- * @brief returns first leaf node for a depth-first traversal of tree
- *
- * @param[in] tree
- * tree to traverse
- *
- * @return node index of first leaf node for depth-first traversal
- *
- * @note called FirstDepthFirstNode in Muscle3.7
- *
- */
-uint
-FirstDepthFirstNode(tree_t *tree)
-{
- uint uNodeIndex;
-
- assert(NULL != tree);
- assert(IsRooted(tree));
-
- /* Descend via left branches until we hit a leaf */
- uNodeIndex = tree->m_uRootNodeIndex;
- while (!IsLeaf(uNodeIndex, tree))
- uNodeIndex = GetLeft(uNodeIndex, tree);
- return uNodeIndex;
-}
-/*** end: FirstDepthFirstNode ***/
-
-
-/**
- * @brief returns next leaf node index for depth-first traversal of
- * tree
- *
- * @param[in] tree
- * tree to traverse
- * @param[in] uNodeIndex
- * Current node index
- *
- * @return node index of next leaf node during depth-first traversal
- *
- * @note called NextDepthFirstNode in Muscle3.7
- */
-uint
-NextDepthFirstNode(uint uNodeIndex, tree_t *tree)
-{
- uint uParent;
-
- assert(NULL != tree);
- assert(IsRooted(tree));
- assert(uNodeIndex < tree->m_uNodeCount);
-
- if (IsRoot(uNodeIndex, tree))
- {
- /* Node uNodeIndex is root, end of traversal */
- return NULL_NEIGHBOR;
- }
-
- uParent = GetParent(uNodeIndex, tree);
- if (GetRight(uParent, tree) == uNodeIndex) {
- /* Is right branch, return parent=uParent */
- return uParent;
- }
-
- uNodeIndex = GetRight(uParent, tree);
- /* Descend left from right sibling uNodeIndex */
- while (!IsLeaf(uNodeIndex, tree)) {
- uNodeIndex = GetLeft(uNodeIndex, tree);
- }
-
- /* bottom out at leaf uNodeIndex */
- return uNodeIndex;
-}
-/*** end: NextDepthFirstNode ***/
-
-
-
-/**
- * @brief check if tree is a rooted tree
- * @param[in] tree
- * tree to check
- * @return TRUE if given tree is rooted, FALSE otherwise
- *
- */
-bool
-IsRooted(tree_t *tree)
-{
- assert(NULL != tree);
- return tree->m_bRooted;
-}
-/*** end: IsRooted ***/
-
-
-/**
- *
- */
-void
-FreeMuscleTree(tree_t *tree)
-{
- uint i;
-
- assert(tree!=NULL);
-
-
- /* FIXME use GetLeafNodeIndex? or
- for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)
- {
- if (tree.IsLeaf(uNodeIndex))
- {
- const char *ptrName =
- tree.GetLeafName(uNodeIndex);
- */
- /* IsLeaf needs m_uNodeCount and all m_uNeighbor's
- * so free first
- */
- for (i=0; i<tree->m_uNodeCount; i++) {
- /* IsLeaf needs neighbour count, so free those guys later */
- if (IsLeaf(i, tree)) {
- CKFREE(tree->m_ptrName[i]);
- }
- }
- CKFREE(tree->m_ptrName);
-
- CKFREE(tree->m_uNeighbor1);
- CKFREE(tree->m_uNeighbor2);
- CKFREE(tree->m_uNeighbor3);
-
- CKFREE(tree->m_Ids);
-
- CKFREE(tree->m_dEdgeLength1);
- CKFREE(tree->m_dEdgeLength2);
- CKFREE(tree->m_dEdgeLength3);
-#if USE_HEIGHT
- CKFREE(tree->m_dHeight);
- CKFREE(tree->m_bHasHeight);
-#endif
- CKFREE(tree->m_bHasEdgeLength1);
- CKFREE(tree->m_bHasEdgeLength2);
- CKFREE(tree->m_bHasEdgeLength3);
-
- TreeZero(tree);
-
- free(tree);
-}
-/*** end: FreeMuscleTree ***/
-
-
-
-/**
- * @brief create a muscle tree
- *
- * @note Original comment in Muscle code: "Create rooted tree from a
- * vector description. Node indexes are 0..N-1 for leaves, N..2N-2 for
- * internal nodes. Vector subscripts are i-N and have values for
- * internal nodes only, but those values are node indexes 0..2N-2. So
- * e.g. if N=6 and Left[2]=1, this means that the third internal node
- * (node index 8) has the second leaf (node index 1) as its left
- * child. uRoot gives the vector subscript of the root, so add N to
- * get the node index."
- *
- * @param[out] tree
- * newly created tree
- * @param[in] uLeafCount
- * number of leaf nodes
- * @param[in] uRoot
- * Internal node index of root node
- * @param[in] Left
- * Node index of left sibling of an internal node.
- * Index range: 0--uLeafCount-1
- * @param[in] Right
- * Node index of right sibling of an internal node.
- * Index range: 0--uLeafCount-1
- * @param[in] LeftLength
- * Branch length of left branch of an internal node.
- * Index range: 0--uLeafCount-1
- * @param[in] RightLength
- * Branch length of right branch of an internal node.
- * Index range: 0--uLeafCount-1
- * @param[in] LeafIds
- * Leaf id. Index range: 0--uLeafCount
- * @param[in] LeafNames
- * Leaf label. Index range: 0--uLeafCount
- *
- */
-void
-MuscleTreeCreate(tree_t *tree,
- uint uLeafCount, uint uRoot,
- const uint *Left, const uint *Right,
- const float *LeftLength, const float* RightLength,
- const uint *LeafIds, char **LeafNames)
-{
- uint uNodeIndex;
-
- TreeZero(tree);
- tree->m_uNodeCount = 2*uLeafCount - 1;
- InitCache(tree->m_uNodeCount, tree);
-
- for (uNodeIndex = 0; uNodeIndex < uLeafCount; ++uNodeIndex) {
- tree->m_Ids[uNodeIndex] = LeafIds[uNodeIndex];
- tree->m_ptrName[uNodeIndex] = CkStrdup(LeafNames[uNodeIndex]);
- }
-
- for (uNodeIndex = uLeafCount; uNodeIndex < tree->m_uNodeCount; ++uNodeIndex) {
- uint v = uNodeIndex - uLeafCount;
- uint uLeft = Left[v];
- uint uRight = Right[v];
- float fLeft = LeftLength[v];
- float fRight = RightLength[v];
-
- tree->m_uNeighbor2[uNodeIndex] = uLeft;
- tree->m_uNeighbor3[uNodeIndex] = uRight;
-
- tree->m_bHasEdgeLength2[uNodeIndex] = TRUE;
- tree->m_bHasEdgeLength3[uNodeIndex] = TRUE;
-
- tree->m_dEdgeLength2[uNodeIndex] = fLeft;
- tree->m_dEdgeLength3[uNodeIndex] = fRight;
-
- tree->m_uNeighbor1[uLeft] = uNodeIndex;
- tree->m_uNeighbor1[uRight] = uNodeIndex;
-
- tree->m_dEdgeLength1[uLeft] = fLeft;
- tree->m_dEdgeLength1[uRight] = fRight;
-
- tree->m_bHasEdgeLength1[uLeft] = TRUE;
- tree->m_bHasEdgeLength1[uRight] = TRUE;
- }
-
- tree->m_bRooted = TRUE;
- tree->m_uRootNodeIndex = uRoot + uLeafCount;
-#ifndef NDEBUG
- TreeValidate(tree);
-#endif
-}
-/*** end: MuscleTreeCreate ***/
-
-
-
-/**
- * @param[in] tree
- * tree to write
- * @param[out] fp
- * filepointer to write to2
- *
- * @brief write a muscle tree to a file in newick format (distances
- * and all names)
- *
- *
- */
-void
-MuscleTreeToFile(FILE *fp, tree_t *tree)
-{
- assert(NULL != tree);
- if (IsRooted(tree)) {
- TreeToFileNodeRooted(tree, tree->m_uRootNodeIndex, fp);
- fprintf(fp, ";\n");
- return;
- } else {
- Log(&rLog, LOG_FATAL, "FIXME: output of unrooted muscle trees not implemented");
- }
-}
-/*** end: MuscleTreeToFile ***/
-
-
-
-/**
- * @brief check if given node is a leaf node
- *
- * @param[in] uNodeIndex
- * the node index
- * @param tree
- * the tree
- *
- * @return TRUE if given node is a leaf, FALSE otherwise
- *
- * @note called IsLeaf in Muscle3.7. See tree.h in original code
- */
-bool
-IsLeaf(uint uNodeIndex, tree_t *tree)
-{
- assert(NULL != tree);
- assert(uNodeIndex < tree->m_uNodeCount);
- if (1 == tree->m_uNodeCount)
- return TRUE;
- return 1 == GetNeighborCount(uNodeIndex, tree);
-}
-/*** end: IsLeaf ***/
-
-
-
-/**
- */
-bool
-IsRoot(uint uNodeIndex, tree_t *tree)
-{
- assert(NULL != tree);
- return IsRooted(tree) && tree->m_uRootNodeIndex == uNodeIndex;
-}
-/*** end: IsRoot ***/
-
-
-
-/**
- */
-uint
-GetNeighborCount(uint uNodeIndex, tree_t *tree)
-{
- uint n1, n2, n3;
- assert(uNodeIndex < tree->m_uNodeCount);
- assert(NULL != tree);
- assert(NULL != tree->m_uNeighbor1);
- assert(NULL != tree->m_uNeighbor2);
- assert(NULL != tree->m_uNeighbor3);
- n1 = tree->m_uNeighbor1[uNodeIndex];
- n2 = tree->m_uNeighbor2[uNodeIndex];
- n3 = tree->m_uNeighbor3[uNodeIndex];
- return (NULL_NEIGHBOR != n1) + (NULL_NEIGHBOR != n2) + (NULL_NEIGHBOR != n3);
-}
-/*** end: GetNeighborCount ***/
-
-
-/**
- */
-uint
-GetParent(unsigned uNodeIndex, tree_t *tree)
-{
- assert(NULL != tree);
- assert(tree->m_bRooted && uNodeIndex < tree->m_uNodeCount);
- return tree->m_uNeighbor1[uNodeIndex];
-}
-/*** end: GetParent ***/
-
-
-
-/**
- */
-bool
-IsEdge(unsigned uNodeIndex1, unsigned uNodeIndex2, tree_t *tree)
-{
- assert(uNodeIndex1 < tree->m_uNodeCount && uNodeIndex2 < tree->m_uNodeCount);
- assert(NULL != tree);
-
- return tree->m_uNeighbor1[uNodeIndex1] == uNodeIndex2 ||
- tree->m_uNeighbor2[uNodeIndex1] == uNodeIndex2 ||
- tree->m_uNeighbor3[uNodeIndex1] == uNodeIndex2;
-}
-/*** end: IsEdge ***/
-
-
-
-/**
- */
-bool
-HasEdgeLength(uint uNodeIndex1, uint uNodeIndex2, tree_t *tree)
-{
- assert(NULL != tree);
- assert(uNodeIndex1 < tree->m_uNodeCount);
- assert(uNodeIndex2 < tree->m_uNodeCount);
- assert(IsEdge(uNodeIndex1, uNodeIndex2, tree));
-
- if (tree->m_uNeighbor1[uNodeIndex1] == uNodeIndex2)
- return tree->m_bHasEdgeLength1[uNodeIndex1];
- else if (tree->m_uNeighbor2[uNodeIndex1] == uNodeIndex2)
- return tree->m_bHasEdgeLength2[uNodeIndex1];
- assert(tree->m_uNeighbor3[uNodeIndex1] == uNodeIndex2);
- return tree->m_bHasEdgeLength3[uNodeIndex1];
-}
-/*** end: ***/
-
-
-/**
- */
-double
-GetEdgeLength(uint uNodeIndex1, uint uNodeIndex2, tree_t *tree)
-{
- assert(NULL != tree);
- assert(uNodeIndex1 < tree->m_uNodeCount && uNodeIndex2 < tree->m_uNodeCount);
- if (!HasEdgeLength(uNodeIndex1, uNodeIndex2, tree))
- {
- Log(&rLog, LOG_FATAL, "Missing edge length in tree %u-%u", uNodeIndex1, uNodeIndex2);
- }
-
- if (tree->m_uNeighbor1[uNodeIndex1] == uNodeIndex2)
- return tree->m_dEdgeLength1[uNodeIndex1];
- else if (tree->m_uNeighbor2[uNodeIndex1] == uNodeIndex2)
- return tree->m_dEdgeLength2[uNodeIndex1];
- assert(tree->m_uNeighbor3[uNodeIndex1] == uNodeIndex2);
- return tree->m_dEdgeLength3[uNodeIndex1];
-}
-/*** end: GetEdgeLength ***/
-
-
-/**
- *
- */
-void
-InitNode(tree_t *prTree, uint uNodeIndex)
-{
- prTree->m_uNeighbor1[uNodeIndex] = NULL_NEIGHBOR;
- prTree->m_uNeighbor2[uNodeIndex] = NULL_NEIGHBOR;
- prTree->m_uNeighbor3[uNodeIndex] = NULL_NEIGHBOR;
- prTree->m_bHasEdgeLength1[uNodeIndex] = FALSE;
- prTree->m_bHasEdgeLength2[uNodeIndex] = FALSE;
- prTree->m_bHasEdgeLength3[uNodeIndex] = FALSE;
-#if USE_HEIGHT
- prTree->m_bHasHeight[uNodeIndex] = FALSE;
- prTree->m_dHeight[uNodeIndex] = dInsane;
-#endif
- prTree->m_dEdgeLength1[uNodeIndex] = dInsane;
- prTree->m_dEdgeLength2[uNodeIndex] = dInsane;
- prTree->m_dEdgeLength3[uNodeIndex] = dInsane;
-
- prTree->m_ptrName[uNodeIndex] = NULL;
- prTree->m_Ids[uNodeIndex] = uInsane;
-}
-/*** end: InitNode ***/
-
-
-/**
- *
- */
-void
-InitCache(uint uCacheCount, tree_t *tree)
-{
- uint uNodeIndex;
-
- tree->m_uCacheCount = uCacheCount;
-
- tree->m_uNeighbor1 = (uint *) CKMALLOC(sizeof(uint) * tree->m_uCacheCount);
- tree->m_uNeighbor2 = (uint *) CKMALLOC(sizeof(uint) * tree->m_uCacheCount);
- tree->m_uNeighbor3 = (uint *) CKMALLOC(sizeof(uint) * tree->m_uCacheCount);
-
- tree->m_Ids = (uint *) CKMALLOC(sizeof(uint) * tree->m_uCacheCount);
-
- tree->m_dEdgeLength1 = (double *) CKMALLOC(sizeof(double) * tree->m_uCacheCount);
- tree->m_dEdgeLength2 = (double *) CKMALLOC(sizeof(double) * tree->m_uCacheCount);
- tree->m_dEdgeLength3 = (double *) CKMALLOC(sizeof(double) * tree->m_uCacheCount);
-#if USE_HEIGHT
- tree->m_dHeight = (double *) CKMALLOC(sizeof(double) * tree->m_uCacheCount);
- tree->m_bHasHeight = (bool *) CKMALLOC(sizeof(bool) * tree->m_uCacheCount);
-#endif
- tree->m_bHasEdgeLength1 = (bool *) CKMALLOC(sizeof(bool) * tree->m_uCacheCount);
- tree->m_bHasEdgeLength2 = (bool *) CKMALLOC(sizeof(bool) * tree->m_uCacheCount);
- tree->m_bHasEdgeLength3 = (bool *) CKMALLOC(sizeof(bool) * tree->m_uCacheCount);
-
- tree->m_ptrName = (char **) CKMALLOC(sizeof(char *) * tree->m_uCacheCount);
-
- for (uNodeIndex = 0; uNodeIndex < tree->m_uNodeCount; ++uNodeIndex) {
- InitNode(tree, uNodeIndex);
- }
-}
-/*** end: InitCache ***/
-
-
-
-
-/**
- *
- * @note Replacing Tree::Clear but no freeing of memory! Just setting
- * everything to 0/NULL
- */
-void
-TreeZero(tree_t *tree)
-{
- assert(NULL != tree);
- tree->m_uNodeCount = 0;
- tree->m_uCacheCount = 0;
-
- tree->m_uNeighbor1 = NULL;
- tree->m_uNeighbor2 = NULL;
- tree->m_uNeighbor3 = NULL;
-
- tree->m_dEdgeLength1 = NULL;
- tree->m_dEdgeLength2 = NULL;
- tree->m_dEdgeLength3 = NULL;
-
-#if USE_HEIGHT
- tree->m_dHeight = NULL;
- tree->m_bHasHeight = NULL;
-#endif
- tree->m_bHasEdgeLength1 = NULL;
- tree->m_bHasEdgeLength2 = NULL;
- tree->m_bHasEdgeLength3 = NULL;
-
- tree->m_ptrName = NULL;
- tree->m_Ids = NULL;
-
- tree->m_bRooted = FALSE;
- tree->m_uRootNodeIndex = 0;
-}
-/* end: TreeZero */
-
-
-
-/**
- *
- */
-void
-TreeToFileNodeRooted(tree_t *tree, uint uNodeIndex, FILE *fp)
-{
- bool bGroup;
-
- assert(IsRooted(tree));
- assert(NULL != tree);
- bGroup = !IsLeaf(uNodeIndex, tree) || IsRoot(uNodeIndex, tree);
-
- if (bGroup)
- fprintf(fp, "(\n");
-
-
- if (IsLeaf(uNodeIndex, tree)) {
- /* File.PutString(GetName(uNodeIndex)); */
- fprintf(fp, "%s", tree->m_ptrName[uNodeIndex]);
- } else {
- TreeToFileNodeRooted(tree, GetLeft(uNodeIndex, tree), fp);
- fprintf(fp, ",\n");
- TreeToFileNodeRooted(tree, GetRight(uNodeIndex, tree), fp);
- }
-
- if (bGroup)
- fprintf(fp, ")");
-
- if (!IsRoot(uNodeIndex, tree)) {
- uint uParent = GetParent(uNodeIndex, tree);
- if (HasEdgeLength(uNodeIndex, uParent, tree))
- fprintf(fp, ":%g", GetEdgeLength(uNodeIndex, uParent, tree));
- }
- fprintf(fp, "\n");
-}
-/*** end: TreeToFileNodeRooted ***/
-
-
-/**
- *
- *
- */
-void
-TreeValidate(tree_t *tree)
-{
- uint uNodeIndex;
- assert(NULL != tree);
- for (uNodeIndex = 0; uNodeIndex < tree->m_uNodeCount; ++uNodeIndex) {
- ValidateNode(uNodeIndex, tree);
- }
- /* FIXME: maybe set negative length to zero? What impact would
- * that have? */
-}
-/*** end TreeValidate ***/
-
-
-
-/**
- *
- *
- */
-void
-ValidateNode(uint uNodeIndex, tree_t *tree)
-{
- uint uNeighborCount;
- uint n1, n2, n3;
- assert(NULL != tree);
-
- if (uNodeIndex >= tree->m_uNodeCount)
- Log(&rLog, LOG_FATAL, "ValidateNode(%u), %u nodes", uNodeIndex, tree->m_uNodeCount);
-
- uNeighborCount = GetNeighborCount(uNodeIndex, tree);
-
-
- if (2 == uNeighborCount) {
- if (!tree->m_bRooted) {
- Log(&rLog, LOG_FATAL, "Tree::ValidateNode: Node %u has two neighbors, tree is not rooted",
- uNodeIndex);
- }
- if (uNodeIndex != tree->m_uRootNodeIndex) {
- Log(&rLog, LOG_FATAL, "Tree::ValidateNode: Node %u has two neighbors, but not root node=%u",
- uNodeIndex, tree->m_uRootNodeIndex);
- }
- }
-
- n1 = tree->m_uNeighbor1[uNodeIndex];
- n2 = tree->m_uNeighbor2[uNodeIndex];
- n3 = tree->m_uNeighbor3[uNodeIndex];
-
- if (NULL_NEIGHBOR == n2 && NULL_NEIGHBOR != n3) {
- Log(&rLog, LOG_FATAL, "Tree::ValidateNode, n2=null, n3!=null", uNodeIndex);
- }
- if (NULL_NEIGHBOR == n3 && NULL_NEIGHBOR != n2) {
- Log(&rLog, LOG_FATAL, "Tree::ValidateNode, n3=null, n2!=null", uNodeIndex);
- }
-
- if (n1 != NULL_NEIGHBOR)
- AssertAreNeighbors(uNodeIndex, n1, tree);
- if (n2 != NULL_NEIGHBOR)
- AssertAreNeighbors(uNodeIndex, n2, tree);
- if (n3 != NULL_NEIGHBOR)
- AssertAreNeighbors(uNodeIndex, n3, tree);
-
-
-
- if (n1 != NULL_NEIGHBOR && (n1 == n2 || n1 == n3)) {
- Log(&rLog, LOG_FATAL, "Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);
- }
- if (n2 != NULL_NEIGHBOR && (n2 == n1 || n2 == n3)) {
- Log(&rLog, LOG_FATAL, "Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);
- }
- if (n3 != NULL_NEIGHBOR && (n3 == n1 || n3 == n2)) {
- Log(&rLog, LOG_FATAL, "Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);
- }
-
-
- if (IsRooted(tree)) {
- if (NULL_NEIGHBOR == GetParent(uNodeIndex, tree)) {
- if (uNodeIndex != tree->m_uRootNodeIndex) {
- Log(&rLog, LOG_FATAL, "Tree::ValiateNode(%u), no parent", uNodeIndex);
- }
- } else if (GetLeft(GetParent(uNodeIndex, tree), tree) != uNodeIndex &&
- GetRight(GetParent(uNodeIndex, tree), tree) != uNodeIndex) {
- Log(&rLog, LOG_FATAL, "Tree::ValidateNode(%u), parent / child mismatch", uNodeIndex);
- }
- }
-}
-/*** end: ValidateNode ***/
-
-
-/**
- *
- *
- */
-void
-AssertAreNeighbors(unsigned uNodeIndex1, unsigned uNodeIndex2, tree_t *tree)
-{
- bool Has12, Has21;
- assert(NULL != tree);
-
- if (uNodeIndex1 >= tree->m_uNodeCount || uNodeIndex2 >= tree->m_uNodeCount)
- Log(&rLog, LOG_FATAL, "AssertAreNeighbors(%u,%u), are %u nodes",
- uNodeIndex1, uNodeIndex2, tree->m_uNodeCount);
-
- if (tree->m_uNeighbor1[uNodeIndex1] != uNodeIndex2 &&
- tree->m_uNeighbor2[uNodeIndex1] != uNodeIndex2 &&
- tree->m_uNeighbor3[uNodeIndex1] != uNodeIndex2) {
- Log(&rLog, LOG_FATAL, "AssertAreNeighbors(%u,%u) failed", uNodeIndex1, uNodeIndex2);
- }
-
-
- if (tree->m_uNeighbor1[uNodeIndex2] != uNodeIndex1 &&
- tree->m_uNeighbor2[uNodeIndex2] != uNodeIndex1 &&
- tree->m_uNeighbor3[uNodeIndex2] != uNodeIndex1) {
- Log(&rLog, LOG_FATAL, "AssertAreNeighbors(%u,%u) failed", uNodeIndex1, uNodeIndex2);
- }
-
-
- Has12 = HasEdgeLength(uNodeIndex1, uNodeIndex2, tree);
- Has21 = HasEdgeLength(uNodeIndex2, uNodeIndex1, tree);
- if (Has12 != Has21) {
- HasEdgeLength(uNodeIndex1, uNodeIndex2, tree);
- HasEdgeLength(uNodeIndex2, uNodeIndex1, tree);
- Log(&rLog, LOG_ERROR, "HasEdgeLength(%u, %u)=%c HasEdgeLength(%u, %u)=%c\n",
- uNodeIndex1,
- uNodeIndex2,
- Has12 ? 'T' : 'F',
- uNodeIndex2,
- uNodeIndex1,
- Has21 ? 'T' : 'F');
- Log(&rLog, LOG_FATAL, "Tree::AssertAreNeighbors, HasEdgeLength not symmetric");
- }
-
-
- if (Has12) {
- double d12 = GetEdgeLength(uNodeIndex1, uNodeIndex2, tree);
- double d21 = GetEdgeLength(uNodeIndex2, uNodeIndex1, tree);
- if (d12 != d21) {
- Log(&rLog, LOG_FATAL, "Tree::AssertAreNeighbors, Edge length disagrees %u-%u=%.3g, %u-%u=%.3g",
- uNodeIndex1, uNodeIndex2, d12,
- uNodeIndex2, uNodeIndex1, d21);
- }
- }
-}
-/*** end: AssertAreNeighbors ***/
-
-
-
-/**
- *
- * @note see phyfromfile.cpp in Muscle3.7. Tree has to be a pointer to
- * an already allocated tree structure.
- *
- * return non-Zero on failure
- *
- * leafids will not be set here (FIXME:CHECK if true)
- */
-int
-MuscleTreeFromFile(tree_t *tree, char *ftree)
-{
- double dEdgeLength;
- bool bEdgeLength;
- char szToken[16];
- NEWICK_TOKEN_TYPE NTT;
- unsigned uThirdNode;
- FILE *fp = NULL;
-
- assert(tree!=NULL);
- assert(ftree!=NULL);
-
- if (NULL == (fp=fopen(ftree, "r"))) {
- Log(&rLog, LOG_ERROR, "Couldn't open tree-file '%s' for reading. Skipping", ftree);
- return -1;
- }
-
- /* Assume rooted.
- * If we discover that it is unrooted, will convert on the fly.
- */
- TreeCreateRooted(tree);
-
- bEdgeLength = GetGroupFromFile(fp, 0, &dEdgeLength, tree);
-
-
- /* Next token should be either ';' for rooted tree or ',' for
- * unrooted.
- */
- NTT = GetToken(fp, szToken, sizeof(szToken));
-
- /* If rooted, all done. */
- if (NTT_Semicolon == NTT) {
- if (bEdgeLength)
- Log(&rLog, LOG_WARN, " *** Warning *** edge length on root group in Newick file %s\n", ftree);
- TreeValidate(tree);
- fclose(fp);
- return 0;
- }
-
- if (NTT_Comma != NTT)
- Log(&rLog, LOG_FATAL, "Tree::FromFile, expected ';' or ',', got '%s'", szToken);
-
- uThirdNode = UnrootFromFile(tree);
- bEdgeLength = GetGroupFromFile(fp, uThirdNode, &dEdgeLength, tree);
- if (bEdgeLength)
- SetEdgeLength(0, uThirdNode, dEdgeLength, tree);
- TreeValidate(tree);
-
- fclose(fp);
- return 0;
-}
-/*** end MuscleTreeFromFile ***/
-
-
-
-/**
- *
- */
-void
-ExpandCache(tree_t *tree)
-{
- const uint uNodeCount = 100;
- uint uNewCacheCount;
- uint *uNewNeighbor1, *uNewNeighbor2, *uNewNeighbor3;
- uint *uNewIds;
- double *dNewEdgeLength1, *dNewEdgeLength2, *dNewEdgeLength3;
-#if USE_HEIGHT
- double *dNewHeight;
- bool *bNewHasHeight;
-#endif
- bool *bNewHasEdgeLength1, *bNewHasEdgeLength2, *bNewHasEdgeLength3;
- char **ptrNewName;
-
- assert(NULL != tree);
- uNewCacheCount = tree->m_uCacheCount + uNodeCount;
- uNewNeighbor1 = (uint *) CKMALLOC(
- uNewCacheCount * sizeof(uint));
- uNewNeighbor2 = (uint *) CKMALLOC(
- uNewCacheCount * sizeof(uint));
- uNewNeighbor3 = (uint *) CKMALLOC(
- uNewCacheCount * sizeof(uint));
-
- uNewIds = (uint *) CKCALLOC(
- uNewCacheCount, sizeof(uint));
-
- dNewEdgeLength1 = (double *) CKMALLOC(
- uNewCacheCount * sizeof(double));
- dNewEdgeLength2 = (double *) CKMALLOC(
- uNewCacheCount * sizeof(double));
- dNewEdgeLength3 = (double *) CKMALLOC(
- uNewCacheCount * sizeof(double));
-#if USE_HEIGHT
- dNewHeight = (double *) CKMALLOC(
- uNewCacheCount * sizeof(double));
- bNewHasHeight = (bool *) CKMALLOC(
- uNewCacheCount * sizeof(bool));
-#endif
-
- bNewHasEdgeLength1 = (bool *) CKMALLOC(
- uNewCacheCount * sizeof(bool));
- bNewHasEdgeLength2 = (bool *) CKMALLOC(
- uNewCacheCount * sizeof(bool));
- bNewHasEdgeLength3 = (bool *) CKMALLOC(
- uNewCacheCount * sizeof(bool));
- ptrNewName = (char **) CKCALLOC(uNewCacheCount, sizeof(char*));
-
- if (tree->m_uCacheCount > 0) {
- uint uUnsignedBytes, uEdgeBytes;
- uint uBoolBytes, uNameBytes;
-
- uUnsignedBytes = tree->m_uCacheCount*sizeof(uint);
-
- memcpy(uNewNeighbor1, tree->m_uNeighbor1, uUnsignedBytes);
- memcpy(uNewNeighbor2, tree->m_uNeighbor2, uUnsignedBytes);
- memcpy(uNewNeighbor3, tree->m_uNeighbor3, uUnsignedBytes);
-
- memcpy(uNewIds, tree->m_Ids, uUnsignedBytes);
-
- uEdgeBytes = tree->m_uCacheCount*sizeof(double);
- memcpy(dNewEdgeLength1, tree->m_dEdgeLength1, uEdgeBytes);
- memcpy(dNewEdgeLength2, tree->m_dEdgeLength2, uEdgeBytes);
- memcpy(dNewEdgeLength3, tree->m_dEdgeLength3, uEdgeBytes);
-#if USE_HEIGHT
- memcpy(dNewHeight, tree->m_dHeight, uEdgeBytes);
-#endif
-
- uBoolBytes = tree->m_uCacheCount*sizeof(bool);
- memcpy(bNewHasEdgeLength1, tree->m_bHasEdgeLength1, uBoolBytes);
- memcpy(bNewHasEdgeLength2, tree->m_bHasEdgeLength2, uBoolBytes);
- memcpy(bNewHasEdgeLength3, tree->m_bHasEdgeLength3, uBoolBytes);
-#if USE_HEIGHT
- memcpy(bNewHasHeight, tree->m_bHasHeight, uBoolBytes);
-#endif
- uNameBytes = tree->m_uCacheCount*sizeof(char *);
- memcpy(ptrNewName, tree->m_ptrName, uNameBytes);
-
- /* similiar to FreeMuscleTree
- */
-
- /* IsLeaf needs m_uNodeCount and all m_uNeighbor's
- * so free first
- */
-#if 0
- for (i=0; i<tree->m_uNodeCount; i++) {
- if (IsLeaf(i, tree)) {
-#ifndef NDEBUG
- if (NULL==tree->m_ptrName[i]) {
- Log(&rLog, LOG_WARN, "FIXME tree->m_ptrName[%d] is already NULL", i);
- }
-#endif
- CKFREE(tree->m_ptrName[i]);
- }
- }
-#endif
- CKFREE(tree->m_ptrName);
-
- CKFREE(tree->m_uNeighbor1);
- CKFREE(tree->m_uNeighbor2);
- CKFREE(tree->m_uNeighbor3);
-
- CKFREE(tree->m_Ids);
-
- CKFREE(tree->m_dEdgeLength1);
- CKFREE(tree->m_dEdgeLength2);
- CKFREE(tree->m_dEdgeLength3);
-
- CKFREE(tree->m_bHasEdgeLength1);
- CKFREE(tree->m_bHasEdgeLength2);
- CKFREE(tree->m_bHasEdgeLength3);
-#if USE_HEIGHT
- CKFREE(tree->m_bHasHeight);
- CKFREE(tree->m_dHeight);
-#endif
- }
-
- tree->m_uCacheCount = uNewCacheCount;
- tree->m_uNeighbor1 = uNewNeighbor1;
- tree->m_uNeighbor2 = uNewNeighbor2;
- tree->m_uNeighbor3 = uNewNeighbor3;
- tree->m_Ids = uNewIds;
- tree->m_dEdgeLength1 = dNewEdgeLength1;
- tree->m_dEdgeLength2 = dNewEdgeLength2;
- tree->m_dEdgeLength3 = dNewEdgeLength3;
-
-#ifdef USE_HEIGHT
- tree->m_dHeight = dNewHeight;
- tree->m_bHasHeight = bNewHasHeight;
-#endif
- tree->m_bHasEdgeLength1 = bNewHasEdgeLength1;
- tree->m_bHasEdgeLength2 = bNewHasEdgeLength2;
- tree->m_bHasEdgeLength3 = bNewHasEdgeLength3;
-
- tree->m_ptrName = ptrNewName;
-
-}
-/*** end: ExpandCache ***/
-
-
-
-/**
- *
- * Tree must be pointer to an already allocated tree structure
- *
- */
-void
-TreeCreateRooted(tree_t *tree)
-{
- TreeZero(tree);
- ExpandCache(tree);
- tree->m_uNodeCount = 1;
-
- tree->m_uNeighbor1[0] = NULL_NEIGHBOR;
- tree->m_uNeighbor2[0] = NULL_NEIGHBOR;
- tree->m_uNeighbor3[0] = NULL_NEIGHBOR;
-
- tree->m_bHasEdgeLength1[0] = FALSE;
- tree->m_bHasEdgeLength2[0] = FALSE;
- tree->m_bHasEdgeLength3[0] = FALSE;
-#if USE_HEIGHT
- tree->m_bHasHeight[0] = FALSE;
-#endif
-
- tree->m_uRootNodeIndex = 0;
- tree->m_bRooted = TRUE;
-
-#ifndef NDEBUG
- TreeValidate(tree);
-#endif
-}
-/*** end: TreeCreateRooted ***/
-
-
-/**
- *
- */
-uint
-UnrootFromFile(tree_t *tree)
-{
- uint uThirdNode;
-#if TRACE
- Log("Before unroot:\n");
- LogMe();
-#endif
-
- if (!tree->m_bRooted)
- Log(&rLog, LOG_FATAL, "Tree::Unroot, not rooted");
-
- /* Convention: root node is always node zero */
- assert(IsRoot(0, tree));
- assert(NULL_NEIGHBOR == tree->m_uNeighbor1[0]);
-
- uThirdNode = tree->m_uNodeCount++;
-
- tree->m_uNeighbor1[0] = uThirdNode;
- tree->m_uNeighbor1[uThirdNode] = 0;
-
- tree->m_uNeighbor2[uThirdNode] = NULL_NEIGHBOR;
- tree->m_uNeighbor3[uThirdNode] = NULL_NEIGHBOR;
-
- tree->m_dEdgeLength1[0] = 0;
- tree->m_dEdgeLength1[uThirdNode] = 0;
- tree->m_bHasEdgeLength1[uThirdNode] = TRUE;
-
- tree->m_bRooted = FALSE;
-
-#if TRACE
- Log("After unroot:\n");
- LogMe();
-#endif
-
- return uThirdNode;
-}
-/*** end: UnrootFromFile ***/
-
-
-
-/**
- *
- *
- */
-bool
-GetGroupFromFile(FILE *fp, uint uNodeIndex, double *ptrdEdgeLength, tree_t *tree)
-{
- char szToken[1024];
- NEWICK_TOKEN_TYPE NTT = GetToken(fp, szToken, sizeof(szToken));
- bool bRightLength;
- bool bEof;
- char c;
-
- /* Group is either leaf name or (left, right). */
- if (NTT_String == NTT) {
- SetLeafName(uNodeIndex, szToken, tree);
-#if TRACE
- Log(&rLog, LOG_FORCED_DEBUG, "Group is leaf '%s'", szToken);
-#endif
- } else if (NTT_Lparen == NTT) {
- const unsigned uLeft = AppendBranch(tree, uNodeIndex);
- const unsigned uRight = uLeft + 1;
- double dEdgeLength;
- bool bLeftLength = GetGroupFromFile(fp, uLeft, &dEdgeLength, tree);
-
- /* Left sub-group...
- */
-#if TRACE
- Log(&rLog, LOG_FORCED_DEBUG, "%s", "Got '(', group is compound, expect left sub-group");
-
- if (bLeftLength) {
- Log(&rLog, LOG_FORCED_DEBUG, "Edge length for left sub-group: %.3g", dEdgeLength);
- } else {
- Log(&rLog, LOG_FORCED_DEBUG, "%s", "No edge length for left sub-group");
- }
-#endif
- if (bLeftLength)
- SetEdgeLength(uNodeIndex, uLeft, dEdgeLength, tree);
-
- /* ... then comma ...
- */
-#if TRACE
- Log(&rLog, LOG_FORCED_DEBUG, "%s", "Expect comma");
-#endif
- NTT = GetToken(fp, szToken, sizeof(szToken));
- if (NTT_Comma != NTT)
- Log(&rLog, LOG_FATAL, "Tree::GetGroupFromFile, expected ',', got '%s'", szToken);
-
- /* ...then right sub-group...
- */
-#if TRACE
- Log(&rLog, LOG_FORCED_DEBUG, "%s", "Expect right sub-group");
-#endif
- bRightLength = GetGroupFromFile(fp, uRight, &dEdgeLength, tree);
- if (bRightLength)
- SetEdgeLength(uNodeIndex, uRight, dEdgeLength, tree);
-
-#if TRACE
- if (bRightLength)
- Log(&rLog, LOG_FORCED_DEBUG, "Edge length for right sub-group: %.3g", dEdgeLength);
- else
- Log(&rLog, LOG_FORCED_DEBUG, "%s", "No edge length for right sub-group");
-#endif
-
- /* ... then closing parenthesis.
- */
-#if TRACE
- Log(&rLog, LOG_FORCED_DEBUG, "%s", "Expect closing parenthesis (or comma if > 2-ary)");
-#endif
- NTT = GetToken(fp, szToken, sizeof(szToken));
- if (NTT_Rparen == NTT)
- ;
- else if (NTT_Comma == NTT) {
- if (ungetc(',', fp)==EOF)
- Log(&rLog, LOG_FATAL, "%s" "ungetc failed");
- return FALSE;
- } else
- Log(&rLog, LOG_FATAL, "Tree::GetGroupFromFile, expected ')' or ',', got '%s'", szToken);
- } else {
- Log(&rLog, LOG_FATAL, "Tree::GetGroupFromFile, expected '(' or leaf name, got '%s'",
- szToken);
- }
-
- /* Group may optionally be followed by edge length.
- */
- bEof = FileSkipWhiteX(fp);
- if (bEof)
- return FALSE;
- if ((c = fgetc(fp))==EOF) /* GetCharX */
- Log(&rLog, LOG_FATAL, "%s", "fgetc reached end of file");
-#if TRACE
- Log(&rLog, LOG_FORCED_DEBUG, "Character following group, could be colon, is '%c'", c);
-#endif
- if (':' == c) {
- NTT = GetToken(fp, szToken, sizeof(szToken));
- if (NTT_String != NTT)
- Log(&rLog, LOG_FATAL, "Tree::GetGroupFromFile, expected edge length, got '%s'", szToken);
- *ptrdEdgeLength = atof(szToken);
- return TRUE;
- }
- if (ungetc(c, fp)==EOF)
- Log(&rLog, LOG_FATAL, "%s" "ungetc failed");
-
- return FALSE;
-}
-/*** end: GetGroupFromFile ***/
-
-
-
-
-/**
- *
- */
-void
-FileSkipWhite(FILE *fp)
-{
- bool bEof = FileSkipWhiteX(fp);
- if (bEof)
- Log(&rLog, LOG_FATAL, "%s", "End-of-file skipping white space");
-}
-/*** end: FileSkipWhite ***/
-
-
-
-
-/**
- *
- */
-bool
-FileSkipWhiteX(FILE *fp)
-{
- for (;;) {
- int c;
- bool bEof;
-
- /* GetChar */
- if ((c = fgetc(fp))==EOF) {
- bEof = TRUE;
- } else {
- bEof = FALSE;
- }
-
- if (bEof)
- return TRUE;
- if (!isspace(c)) {
- if (ungetc(c, fp)==EOF)
- Log(&rLog, LOG_FATAL, "%s" "ungetc failed");
- break;
- }
- }
- return FALSE;
-}
-/*** end: FileSkipWhiteX ***/
-
-
-
-
-/**
- *
- */
-NEWICK_TOKEN_TYPE
-GetToken(FILE *fp, char szToken[], uint uBytes)
-{
- char c;
- unsigned uBytesCopied = 0;
- NEWICK_TOKEN_TYPE TT;
-
- /* Skip leading white space */
- FileSkipWhite(fp);
-
- if ((c = fgetc(fp))==EOF) /* GetCharX */
- Log(&rLog, LOG_FATAL, "%s", "fgetc reached end of file");
-
- /* In case a single-character token */
- szToken[0] = c;
- szToken[1] = 0;
-
- switch (c) {
-
- case '(':
- return NTT_Lparen;
-
- case ')':
- return NTT_Rparen;
-
- case ':':
- return NTT_Colon;
-
- case ';':
- return NTT_Semicolon;
-
- case ',':
- return NTT_Comma;
-
- case '\'':
- TT = NTT_SingleQuotedString;
- if ((c = fgetc(fp))==EOF) /* GetCharX */
- Log(&rLog, LOG_FATAL, "%s", "fgetc reached end of file");
- break;
-
- case '"':
- TT = NTT_DoubleQuotedString;
- if ((c = fgetc(fp))==EOF) /* GetCharX */
- Log(&rLog, LOG_FATAL, "%s", "fgetc reached end of file");
- break;
-
- case '[':
- TT = NTT_Comment;
- break;
-
- default:
- TT = NTT_String;
- break;
- }
-
- for (;;)
- {
- bool bEof;
- if (TT != NTT_Comment) {
- if (uBytesCopied < uBytes - 2) {
- szToken[uBytesCopied++] = c;
- szToken[uBytesCopied] = 0;
- } else {
- Log(&rLog, LOG_FATAL, "Tree::GetToken: input buffer too small, token so far='%s'", szToken);
- }
- }
- c = fgetc(fp); /* GetChar */
- bEof = (c==EOF ? TRUE : FALSE);
- if (bEof)
- return TT;
-
- switch (TT) {
-
- case NTT_String:
- if (0 != strchr("():;,", c)) {
- if (ungetc(c, fp)==EOF)
- Log(&rLog, LOG_FATAL, "%s" "ungetc failed");
- return NTT_String;
- }
- if (isspace(c))
- return NTT_String;
- break;
-
- case NTT_SingleQuotedString:
- if ('\'' == c)
- return NTT_String;
- break;
-
- case NTT_DoubleQuotedString:
- if ('"' == c)
- return NTT_String;
- break;
-
- case NTT_Comment:
- if (']' == c)
- return GetToken(fp, szToken, uBytes);
- break;
-
- default:
- Log(&rLog, LOG_FATAL, "Tree::GetToken, invalid TT=%u", TT);
- }
- }
-}
-/*** end: GetToken ***/
-
-
-
-/*** SetLeafName
- *
- */
-void
-SetLeafName(unsigned uNodeIndex, const char *ptrName, tree_t *tree)
-{
- assert(uNodeIndex < tree->m_uNodeCount);
- assert(IsLeaf(uNodeIndex, tree));
- free(tree->m_ptrName[uNodeIndex]);
- /*LOG_DEBUG("Setting tree->m_ptrName[uNodeIndex=%d] to %s", uNodeIndex, ptrName);*/
- tree->m_ptrName[uNodeIndex] = CkStrdup(ptrName);
-}
-/*** end: SetLeafName ***/
-
-
-
-
-/**
- *
- * Append a new branch. This adds two new nodes and joins them to an
- * existing leaf node. Return value is k, new nodes have indexes k and
- * k+1 respectively.
- *
- */
-uint
-AppendBranch(tree_t *tree, uint uExistingLeafIndex)
-{
- uint uNewLeaf1;
- uint uNewLeaf2;
-
- assert(tree!=NULL);
- if (0 == tree->m_uNodeCount) {
- Log(&rLog, LOG_FATAL, "%s(): %s", __FUNCTION__, "tree has not been created");
- }
- assert(NULL_NEIGHBOR == tree->m_uNeighbor2[uExistingLeafIndex]);
- assert(NULL_NEIGHBOR == tree->m_uNeighbor3[uExistingLeafIndex]);
- assert(uExistingLeafIndex < tree->m_uNodeCount);
-#ifndef NDEBUG
- if (!IsLeaf(uExistingLeafIndex, tree)) {
- Log(&rLog, LOG_FATAL, "AppendBranch(%u): not leaf", uExistingLeafIndex);
- }
-#endif
-
- if (tree->m_uNodeCount >= tree->m_uCacheCount - 2) {
- ExpandCache(tree);
- }
- uNewLeaf1 = tree->m_uNodeCount;
- uNewLeaf2 = tree->m_uNodeCount + 1;
-
- tree->m_uNodeCount += 2;
-
- tree->m_uNeighbor2[uExistingLeafIndex] = uNewLeaf1;
- tree->m_uNeighbor3[uExistingLeafIndex] = uNewLeaf2;
-
- tree->m_uNeighbor1[uNewLeaf1] = uExistingLeafIndex;
- tree->m_uNeighbor1[uNewLeaf2] = uExistingLeafIndex;
-
- tree->m_uNeighbor2[uNewLeaf1] = NULL_NEIGHBOR;
- tree->m_uNeighbor2[uNewLeaf2] = NULL_NEIGHBOR;
-
- tree->m_uNeighbor3[uNewLeaf1] = NULL_NEIGHBOR;
- tree->m_uNeighbor3[uNewLeaf2] = NULL_NEIGHBOR;
-
- tree->m_dEdgeLength2[uExistingLeafIndex] = 0;
- tree->m_dEdgeLength3[uExistingLeafIndex] = 0;
-
- tree->m_dEdgeLength1[uNewLeaf1] = 0;
- tree->m_dEdgeLength2[uNewLeaf1] = 0;
- tree->m_dEdgeLength3[uNewLeaf1] = 0;
-
- tree->m_dEdgeLength1[uNewLeaf2] = 0;
- tree->m_dEdgeLength2[uNewLeaf2] = 0;
- tree->m_dEdgeLength3[uNewLeaf2] = 0;
-
- tree->m_bHasEdgeLength1[uNewLeaf1] = FALSE;
- tree->m_bHasEdgeLength2[uNewLeaf1] = FALSE;
- tree->m_bHasEdgeLength3[uNewLeaf1] = FALSE;
-
- tree->m_bHasEdgeLength1[uNewLeaf2] = FALSE;
- tree->m_bHasEdgeLength2[uNewLeaf2] = FALSE;
- tree->m_bHasEdgeLength3[uNewLeaf2] = FALSE;
-
-#if USE_HEIGHT
- tree->m_bHasHeight[uNewLeaf1] = FALSE;
- tree->m_bHasHeight[uNewLeaf2] = FALSE;
-#endif
- tree->m_Ids[uNewLeaf1] = uInsane;
- tree->m_Ids[uNewLeaf2] = uInsane;
-
- return uNewLeaf1;
-}
-/*** end: AppendBranch ***/
-
-
-/**
- *
- *
- */
-void
-SetEdgeLength(uint uNodeIndex1, uint uNodeIndex2,
- double dLength, tree_t *tree)
-{
- assert(uNodeIndex1 < tree->m_uNodeCount && uNodeIndex2 < tree->m_uNodeCount);
- assert(IsEdge(uNodeIndex1, uNodeIndex2, tree));
-
- if (tree->m_uNeighbor1[uNodeIndex1] == uNodeIndex2) {
- tree->m_dEdgeLength1[uNodeIndex1] = dLength;
- tree->m_bHasEdgeLength1[uNodeIndex1] = TRUE;
- } else if (tree->m_uNeighbor2[uNodeIndex1] == uNodeIndex2) {
- tree->m_dEdgeLength2[uNodeIndex1] = dLength;
- tree->m_bHasEdgeLength2[uNodeIndex1] = TRUE;
-
- } else {
- assert(tree->m_uNeighbor3[uNodeIndex1] == uNodeIndex2);
- tree->m_dEdgeLength3[uNodeIndex1] = dLength;
- tree->m_bHasEdgeLength3[uNodeIndex1] = TRUE;
- }
-
- if (tree->m_uNeighbor1[uNodeIndex2] == uNodeIndex1) {
- tree->m_dEdgeLength1[uNodeIndex2] = dLength;
- tree->m_bHasEdgeLength1[uNodeIndex2] = TRUE;
- } else if (tree->m_uNeighbor2[uNodeIndex2] == uNodeIndex1) {
- tree->m_dEdgeLength2[uNodeIndex2] = dLength;
- tree->m_bHasEdgeLength2[uNodeIndex2] = TRUE;
- } else {
- assert(tree->m_uNeighbor3[uNodeIndex2] == uNodeIndex1);
- tree->m_dEdgeLength3[uNodeIndex2] = dLength;
- tree->m_bHasEdgeLength3[uNodeIndex2] = TRUE;
- }
-}
-/*** end: SetEdgeLength ***/
-
-
-
-/**
- *
- * Debug output
- *
- * LogMe in phy.cpp
- *
- */
-void
-LogTree(tree_t *tree, FILE *fp)
-{
- uint uNodeIndex;
- uint n1, n2, n3;
- char *ptrName;
-
- fprintf(fp, "This is a tree with %u nodes, which is ", tree->m_uNodeCount);
-
- if (IsRooted(tree)) {
- fprintf(fp, "rooted:\n");
- fprintf(fp, "Index Parnt LengthP Left LengthL Right LengthR Id Name\n");
- fprintf(fp, "----- ----- ------- ---- ------- ----- ------- ----- ----\n");
-
- } else {
- fprintf(fp, "unrooted;\n");
- fprintf(fp, "Index Nbr_1 Length1 Nbr_2 Length2 Nbr_3 Length3 Id Name\n");
- fprintf(fp, "----- ----- ------- ----- ------- ----- ------- ----- ----\n");
- }
-
- for (uNodeIndex = 0; uNodeIndex < tree->m_uNodeCount; ++uNodeIndex) {
- fprintf(fp, "%5u ", uNodeIndex);
- n1 = tree->m_uNeighbor1[uNodeIndex];
- n2 = tree->m_uNeighbor2[uNodeIndex];
- n3 = tree->m_uNeighbor3[uNodeIndex];
-
- if (NULL_NEIGHBOR != n1) {
- fprintf(fp, "%5u ", n1);
- if (tree->m_bHasEdgeLength1[uNodeIndex])
- fprintf(fp, "%7.3g ", tree->m_dEdgeLength1[uNodeIndex]);
- else
- fprintf(fp, " * ");
- } else {
- fprintf(fp, " ");
- }
-
- if (NULL_NEIGHBOR != n2) {
- fprintf(fp, "%5u ", n2);
- if (tree->m_bHasEdgeLength2[uNodeIndex])
- fprintf(fp, "%7.3g ", tree->m_dEdgeLength2[uNodeIndex]);
- else
- fprintf(fp, " * ");
- } else {
- fprintf(fp, " ");
- }
-
- if (NULL_NEIGHBOR != n3) {
- fprintf(fp, "%5u ", n3);
- if (tree->m_bHasEdgeLength3[uNodeIndex])
- fprintf(fp, "%7.3g ", tree->m_dEdgeLength3[uNodeIndex]);
- else
- fprintf(fp, " * ");
- } else {
- fprintf(fp, " ");
- }
-
- if (tree->m_Ids != 0 && IsLeaf(uNodeIndex, tree)) {
- unsigned uId = tree->m_Ids[uNodeIndex];
- if (uId == uInsane)
- fprintf(fp, " *");
- else
- fprintf(fp, "%5u", uId);
- } else {
- fprintf(fp, " ");
- }
-
- if (tree->m_bRooted && uNodeIndex == tree->m_uRootNodeIndex)
- fprintf(fp, " [ROOT] ");
- ptrName = tree->m_ptrName[uNodeIndex];
- if (ptrName != 0)
- fprintf(fp, " %s", ptrName);
-
- fprintf(fp, "\n");
- }
-}
-/*** end: LogTree ***/
-
-
-
-/**
- *
- * replaces m_uLeafCount
- */
-uint
-GetLeafCount(tree_t *tree)
-{
- assert(tree!=NULL);
-
- return (tree->m_uNodeCount+1)/2;
-}
-/*** GetLeafCount ***/
-
-
-
-/**
- *
- */
-uint
-GetNodeCount(tree_t *tree)
-{
- assert(tree!=NULL);
-
- return 2*(GetLeafCount(tree)) - 1;
-}
-/*** end: GetNodeCount ***/
-
-
-/**
- *
- */
-uint
-GetNeighbor(uint uNodeIndex, uint uNeighborSubscript, tree_t *prTree)
-{
- assert(uNodeIndex < prTree->m_uNodeCount);
- switch (uNeighborSubscript)
- {
- case 0:
- return prTree->m_uNeighbor1[uNodeIndex];
- case 1:
- return prTree->m_uNeighbor2[uNodeIndex];
- case 2:
- return prTree->m_uNeighbor3[uNodeIndex];
- }
- Log(&rLog, LOG_FATAL, "Internal error in %s: sub=%u", __FUNCTION__, uNeighborSubscript);
- return NULL_NEIGHBOR;
-}
-/*** end: GetNeighbor ***/
-
-
-
-
-
-/**
- *
- */
-void
-SetLeafId(tree_t *tree, uint uNodeIndex, uint uId)
-{
- assert(uNodeIndex < tree->m_uNodeCount);
- assert(IsLeaf(uNodeIndex, tree));
- tree->m_Ids[uNodeIndex] = uId;
-}
-/*** end: SetLeafId ***/
-
-
-/**
- *
- */
-uint
-GetRootNodeIndex(tree_t *tree)
-{
- assert(NULL!=tree);
- return tree->m_uRootNodeIndex;
-}
-/*** end: GetRootNodeIndex ***/
-
-
-
-/**
- * @note avoid usage if you want to iterate over all indices, because
- * this will be slow
- *
- */
-uint
-LeafIndexToNodeIndex(uint uLeafIndex, tree_t *prTree) {
- uint uLeafCount = 0;
- unsigned uNodeCount = GetNodeCount(prTree);
- uint uNodeIndex;
-
- for (uNodeIndex = 0; uNodeIndex < uNodeCount; uNodeIndex++) {
- if (IsLeaf(uNodeIndex, prTree)) {
- if (uLeafCount == uLeafIndex) {
- return uNodeIndex;
- } else {
- uLeafCount++;
- }
- }
- }
- Log(&rLog, LOG_FATAL, "Internal error: node index out of range");
- return 0;
-}
-/*** end: LeafIndexToNodeIndex ***/
-
-
-
-
-/**
- * @brief Append a (source) tree to a (dest) tree to a given node
- * which will be replaced. All other nodes in that tree will stay the
- * same.
- *
- * @param[out] prDstTree
- * The tree to append to
- * @param[in] uDstTreeReplaceNodeIndex
- * Dest tree node to which source tree will be appended
- * @param[in] prSrcTree
- * The tree to append
- *
- * @note No nodes inside prDstTree will change except
- * uDstTreeReplaceNodeIndex
- *
- *
- * @warning: Function won't check or touch the m_Ids/leaf-indices!
- * That means if you want to join two trees with leaf indices 1-10 and
- * 1-10 your m_Ids/leaf-indices won't be unique anymore and the
- * association between your sequences and the tree are broken. Make
- * sure m_Ids are unique before calling me.
- *
- * The easiest would have been to do this by recursively calling
- * AppendBranch() (after adding uSrcTreeNodeIndex as extra argument to
- * this function). But recursion is evil. Yet another version would be
- * to setup all the data and call MuscleTreeCreate() to create a third
- * tree, which seems complicated and wasteful. The approach taken here
- * is the following: increase dest tree memory, copy over each src
- * tree node data and update the indices and counters. This is tricky
- * and has a lot of potential for bugs if tree interface is changed
- * (and even if it isn't).
- *
- */
-void
-AppendTree(tree_t *prDstTree, uint uDstTreeReplaceNodeIndex, tree_t *prSrcTree)
-{
- uint uSrcTreeNodeIndex;
- uint uOrgDstTreeNodeCount;
-
- assert(NULL!=prDstTree);
- assert(NULL!=prSrcTree);
- assert(uDstTreeReplaceNodeIndex < prDstTree->m_uNodeCount);
- assert(IsLeaf(uDstTreeReplaceNodeIndex, prDstTree));
- assert(IsRooted(prDstTree) && IsRooted(prSrcTree));
-
-
- uOrgDstTreeNodeCount = prDstTree->m_uNodeCount;
-
-
- /* increase dest tree memory
- */
- while (prDstTree->m_uCacheCount
- <
- (GetNodeCount(prDstTree) + GetNodeCount(prSrcTree))) {
- ExpandCache(prDstTree);
- }
-
-
- /* copy all src tree nodes
- *
- */
- for (uSrcTreeNodeIndex=0;
- uSrcTreeNodeIndex<GetNodeCount(prSrcTree); uSrcTreeNodeIndex++) {
- uint uNewDstNodeIndex = prDstTree->m_uNodeCount;
-
- /* distinguish 4 cases for src nodes to copy:
- *
- * 1. src node is the only node, i.e. root & leaf
- *
- * 2. src node is root: set only left & right, but not parent
- * and just replace the given dest index
- *
- * 3. src node is leaf: set only parent
- *
- * 4. src node is internal node: update all three neighbours
- *
- * FIXME: this is messy. Is there a cleaner way to do this by
- * merging all cases?
- *
- */
- if (IsRoot(uSrcTreeNodeIndex, prSrcTree) && IsLeaf(uSrcTreeNodeIndex, prSrcTree)) {
- /* special case: if this is the only member in
- * tree, i.e. it's root and leaf. Copy leaf name and leaf
- * id. No neighbours to update
- */
-
- /* free dst node name if set */
- if (NULL != prDstTree->m_ptrName[uDstTreeReplaceNodeIndex]) {
- CKFREE(prDstTree->m_ptrName[uDstTreeReplaceNodeIndex]);
- }
-
- prDstTree->m_ptrName[uDstTreeReplaceNodeIndex] =
- CkStrdup(GetLeafName(uSrcTreeNodeIndex, prSrcTree));
-
- prDstTree->m_Ids[uDstTreeReplaceNodeIndex] =
- prSrcTree->m_Ids[uSrcTreeNodeIndex];
-
- /* no updated of uNodeCount, because we used the replace node */
-
-#if TRACE
- Log(&rLog, LOG_FORCED_DEBUG, "Updated dst rpl node %d with the only src leaf node: parent=%d (%f)",
- uDstTreeReplaceNodeIndex,
- prDstTree->m_uNeighbor1[uDstTreeReplaceNodeIndex], prDstTree->m_dEdgeLength1[uDstTreeReplaceNodeIndex]);
-#endif
-
- } else if (IsRoot(uSrcTreeNodeIndex, prSrcTree)) {
- /* src node is root: replace uDstTreeReplaceNodeIndex
- * (not uNewDstNodeIndex) with src root, i.e. the
- * uDstTreeReplaceNodeIndex becomes an internal node now.
- *
- * We only have two neighbours 2 & 3 (no parent). Keep old
- * parent info (neighbor 1).
- */
-
- /* free dst node name if set */
- if (NULL != prDstTree->m_ptrName[uDstTreeReplaceNodeIndex]) {
- CKFREE(prDstTree->m_ptrName[uDstTreeReplaceNodeIndex]);
- }
-
- prDstTree->m_uNeighbor2[uDstTreeReplaceNodeIndex] =
- prSrcTree->m_uNeighbor2[uSrcTreeNodeIndex] + uOrgDstTreeNodeCount;
- prDstTree->m_uNeighbor3[uDstTreeReplaceNodeIndex] =
- prSrcTree->m_uNeighbor3[uSrcTreeNodeIndex] + uOrgDstTreeNodeCount;
-
- prDstTree->m_bHasEdgeLength2[uDstTreeReplaceNodeIndex] =
- prSrcTree->m_bHasEdgeLength2[uSrcTreeNodeIndex];
- prDstTree->m_bHasEdgeLength3[uDstTreeReplaceNodeIndex] =
- prSrcTree->m_bHasEdgeLength3[uSrcTreeNodeIndex];
-
- prDstTree->m_dEdgeLength2[uDstTreeReplaceNodeIndex] =
- prSrcTree->m_dEdgeLength2[uSrcTreeNodeIndex];
- prDstTree->m_dEdgeLength3[uDstTreeReplaceNodeIndex] =
- prSrcTree->m_dEdgeLength3[uSrcTreeNodeIndex];
-
- /* make Id invalid */
- prDstTree->m_Ids[uDstTreeReplaceNodeIndex] = uInsane;
-
- /* no updated of uNodeCount, because we used the replace node */
-
-#if TRACE
- Log(&rLog, LOG_FORCED_DEBUG, "Updated dst rpl node %d with the src root node: (untouched) parent=%d (%f) left=%d (%f) right=%d (%f)",
- uDstTreeReplaceNodeIndex,
- prDstTree->m_uNeighbor1[uDstTreeReplaceNodeIndex], prDstTree->m_dEdgeLength1[uDstTreeReplaceNodeIndex],
- prDstTree->m_uNeighbor2[uDstTreeReplaceNodeIndex], prDstTree->m_dEdgeLength2[uDstTreeReplaceNodeIndex],
- prDstTree->m_uNeighbor3[uDstTreeReplaceNodeIndex], prDstTree->m_dEdgeLength3[uDstTreeReplaceNodeIndex]);
-#endif
-
- } else if (IsLeaf(uSrcTreeNodeIndex, prSrcTree)) {
- /* src node is a leaf, which means we only have one
- * neighbour, and that is its parent, i.e. n1
- *
- */
-
- /* initialise/zero new node to default values
- */
- InitNode(prDstTree, uNewDstNodeIndex);
-
-
- /* update m_ptrName/leaf name
- */
- prDstTree->m_ptrName[uNewDstNodeIndex] =
- CkStrdup(GetLeafName(uSrcTreeNodeIndex, prSrcTree));
-
- /* update parent node (beware of special case: parent was
- src tree root */
- if (IsRoot(prSrcTree->m_uNeighbor1[uSrcTreeNodeIndex], prSrcTree)) {
- prDstTree->m_uNeighbor1[uNewDstNodeIndex] =
- uDstTreeReplaceNodeIndex;
- } else {
- prDstTree->m_uNeighbor1[uNewDstNodeIndex] =
- prSrcTree->m_uNeighbor1[uSrcTreeNodeIndex] + uOrgDstTreeNodeCount;
- }
-
- /* update edge length info to parent
- */
- prDstTree->m_bHasEdgeLength1[uNewDstNodeIndex] =
- prSrcTree->m_bHasEdgeLength1[uSrcTreeNodeIndex];
- prDstTree->m_dEdgeLength1[uNewDstNodeIndex] =
- prSrcTree->m_dEdgeLength1[uSrcTreeNodeIndex];
-
- /* update sequence/object id
- */
- prDstTree->m_Ids[uNewDstNodeIndex] =
- prSrcTree->m_Ids[uSrcTreeNodeIndex];
-
- /* we used a new node so increase their count */
- prDstTree->m_uNodeCount += 1;
-
-#if TRACE
- Log(&rLog, LOG_FORCED_DEBUG, "Updated dst node %d with a src leaf node: parent=%d (%f)",
- uNewDstNodeIndex,
- prDstTree->m_uNeighbor1[uNewDstNodeIndex], prDstTree->m_dEdgeLength1[uNewDstNodeIndex]);
-#endif
-
- } else {
- /* src node is not root neither leaf, means we have an
- * internal node. Update all neighbour info
- *
- */
-
- /* initialise/zero node values to default values
- */
- InitNode(prDstTree, uNewDstNodeIndex);
-
- /* update neigbours
- */
- /* parent: special case if parent was src tree root */
- if (IsRoot(prSrcTree->m_uNeighbor1[uSrcTreeNodeIndex], prSrcTree)) {
- prDstTree->m_uNeighbor1[uNewDstNodeIndex] =
- uDstTreeReplaceNodeIndex;
- } else {
- prDstTree->m_uNeighbor1[uNewDstNodeIndex] =
- prSrcTree->m_uNeighbor1[uSrcTreeNodeIndex] + uOrgDstTreeNodeCount;
- }
- /* left */
- prDstTree->m_uNeighbor2[uNewDstNodeIndex] =
- prSrcTree->m_uNeighbor2[uSrcTreeNodeIndex] + uOrgDstTreeNodeCount;
- /* right */
- prDstTree->m_uNeighbor3[uNewDstNodeIndex] =
- prSrcTree->m_uNeighbor3[uSrcTreeNodeIndex] + uOrgDstTreeNodeCount;
-
- /* update edge length info
- */
- /* parent */
- prDstTree->m_bHasEdgeLength1[uNewDstNodeIndex] =
- prSrcTree->m_bHasEdgeLength1[uSrcTreeNodeIndex];
- prDstTree->m_dEdgeLength1[uNewDstNodeIndex] =
- prSrcTree->m_dEdgeLength1[uSrcTreeNodeIndex];
- /* left */
- prDstTree->m_bHasEdgeLength2[uNewDstNodeIndex] =
- prSrcTree->m_bHasEdgeLength2[uSrcTreeNodeIndex];
- prDstTree->m_dEdgeLength2[uNewDstNodeIndex] =
- prSrcTree->m_dEdgeLength2[uSrcTreeNodeIndex];
- /* right */
- prDstTree->m_bHasEdgeLength3[uNewDstNodeIndex] =
- prSrcTree->m_bHasEdgeLength3[uSrcTreeNodeIndex];
- prDstTree->m_dEdgeLength3[uNewDstNodeIndex] =
- prSrcTree->m_dEdgeLength3[uSrcTreeNodeIndex];
-
- /* we used a new node so increase their count */
- prDstTree->m_uNodeCount += 1;
-
-#if TRACE
- Log(&rLog, LOG_FORCED_DEBUG, "Updated dst node %d with an internal src node: parent=%d (%f) left=%d (%f) right=%d (%f)",
- uNewDstNodeIndex,
- prDstTree->m_uNeighbor1[uNewDstNodeIndex], prDstTree->m_dEdgeLength1[uNewDstNodeIndex],
- prDstTree->m_uNeighbor2[uNewDstNodeIndex], prDstTree->m_dEdgeLength2[uNewDstNodeIndex],
- prDstTree->m_uNeighbor3[uNewDstNodeIndex], prDstTree->m_dEdgeLength3[uNewDstNodeIndex]);
-#endif
- }
-
- }
- /* end for each src tree node */
-
-
- /*
- * m_uRootNodeIndex stays the same.
- *
- * No need to touch m_uCacheCount.
- *
- */
-#if USE_HEIGHT
- Log(&rLog, LOG_FATAL, "Internal error: Height usage not implemented in %s", __FUNCTION__);
-#endif
-
-
-#ifndef NDEBUG
- TreeValidate(prDstTree);
-#endif
-
- return;
-}
-/*** end: AppendTree() ***/
-
git://source.jalview.org / jabaws.git / blobdiff