C++ Tree::GetEdgeLength方法代码示例

void ApplyMinEdgeLength(Tree &tree, double dMinEdgeLength)
	{
	const unsigned uNodeCount = tree.GetNodeCount();
	for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
		{
		const unsigned uNeighborCount = tree.GetNeighborCount(uNodeIndex);
		for (unsigned n = 0; n < uNeighborCount; ++n)
			{
			const unsigned uNeighborNodeIndex = tree.GetNeighbor(uNodeIndex, n);
			if (!tree.HasEdgeLength(uNodeIndex, uNeighborNodeIndex))
				continue;
			if (tree.GetEdgeLength(uNodeIndex, uNeighborNodeIndex) < dMinEdgeLength)
				tree.SetEdgeLength(uNodeIndex, uNeighborNodeIndex, dMinEdgeLength);
			}
		}
	}

void CalcClustalWWeights(const Tree &tree, WEIGHT Weights[])
	{
#if	TRACE
	Log("CalcClustalWWeights\n");
	tree.LogMe();
#endif

	const unsigned uLeafCount = tree.GetLeafCount();
	if (0 == uLeafCount)
		return;
	else if (1 == uLeafCount)
		{
		Weights[0] = (WEIGHT) 1.0;
		return;
		}
	else if (2 == uLeafCount)
		{
		Weights[0] = (WEIGHT) 0.5;
		Weights[1] = (WEIGHT) 0.5;
		return;
		}

	if (!tree.IsRooted())
		Quit("CalcClustalWWeights requires rooted tree");

	const unsigned uNodeCount = tree.GetNodeCount();
	unsigned *LeavesUnderNode = new unsigned[uNodeCount];
	memset(LeavesUnderNode, 0, uNodeCount*sizeof(unsigned));

	const unsigned uRootNodeIndex = tree.GetRootNodeIndex();
	unsigned uLeavesUnderRoot = CountLeaves(tree, uRootNodeIndex, LeavesUnderNode);
	if (uLeavesUnderRoot != uLeafCount)
		Quit("WeightsFromTreee: Internal error, root count %u %u",
		  uLeavesUnderRoot, uLeafCount);

#if	TRACE
	Log("Node  Leaves    Length  Strength\n");
	Log("----  ------  --------  --------\n");
	//    1234  123456  12345678  12345678
#endif

	double *Strengths = new double[uNodeCount];
	for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
		{
		if (tree.IsRoot(uNodeIndex))
			{
			Strengths[uNodeIndex] = 0.0;
			continue;
			}
		const unsigned uParent = tree.GetParent(uNodeIndex);
		const double dLength = tree.GetEdgeLength(uNodeIndex, uParent);
		const unsigned uLeaves = LeavesUnderNode[uNodeIndex];
		const double dStrength = dLength / (double) uLeaves;
		Strengths[uNodeIndex] = dStrength;
#if	TRACE
		Log("%4u  %6u  %8g  %8g\n", uNodeIndex, uLeaves, dLength, dStrength);
#endif
		}

#if	TRACE
	Log("\n");
	Log("                 Seq  Path..Weight\n");
	Log("--------------------  ------------\n");
#endif
	for (unsigned n = 0; n < uLeafCount; ++n)
		{
		const unsigned uLeafNodeIndex = tree.LeafIndexToNodeIndex(n);
#if	TRACE
		Log("%20.20s  %4u ", tree.GetLeafName(uLeafNodeIndex), uLeafNodeIndex);
#endif
		if (!tree.IsLeaf(uLeafNodeIndex))
			Quit("CalcClustalWWeights: leaf");

		double dWeight = 0;
		unsigned uNode = uLeafNodeIndex;
		while (!tree.IsRoot(uNode))
			{
			dWeight += Strengths[uNode];
			uNode = tree.GetParent(uNode);
#if	TRACE
			Log("->%u(%g)", uNode, Strengths[uNode]);
#endif
			}
		if (dWeight < 0.0001)
			{
#if	TRACE
			Log("zero->one");
#endif
			dWeight = 1.0;
			}
		Weights[n] = (WEIGHT) dWeight;
#if	TRACE
		Log(" = %g\n", dWeight);
#endif
		}

	delete[] Strengths;
	delete[] LeavesUnderNode;

	Normalize(Weights, uLeafCount);
//.........这里部分代码省略.........

void Tree::PruneTree(const Tree &tree, unsigned Subfams[],
                     unsigned uSubfamCount)
{
    if (!tree.IsRooted())
        Quit("Tree::PruneTree: requires rooted tree");

    Clear();

    m_uNodeCount = 2*uSubfamCount - 1;
    InitCache(m_uNodeCount);

    const unsigned uUnprunedNodeCount = tree.GetNodeCount();

    unsigned *uUnprunedToPrunedIndex = new unsigned[uUnprunedNodeCount];
    unsigned *uPrunedToUnprunedIndex = new unsigned[m_uNodeCount];

    for (unsigned n = 0; n < uUnprunedNodeCount; ++n)
        uUnprunedToPrunedIndex[n] = NULL_NEIGHBOR;

    for (unsigned n = 0; n < m_uNodeCount; ++n)
        uPrunedToUnprunedIndex[n] = NULL_NEIGHBOR;

// Create mapping between unpruned and pruned node indexes
    unsigned uInternalNodeIndex = uSubfamCount;
    for (unsigned uSubfamIndex = 0; uSubfamIndex < uSubfamCount; ++uSubfamIndex)
    {
        unsigned uUnprunedNodeIndex = Subfams[uSubfamIndex];
        uUnprunedToPrunedIndex[uUnprunedNodeIndex] = uSubfamIndex;
        uPrunedToUnprunedIndex[uSubfamIndex] = uUnprunedNodeIndex;
        for (;;)
        {
            uUnprunedNodeIndex = tree.GetParent(uUnprunedNodeIndex);
            if (tree.IsRoot(uUnprunedNodeIndex))
                break;

            // Already visited this node?
            if (NULL_NEIGHBOR != uUnprunedToPrunedIndex[uUnprunedNodeIndex])
                break;

            uUnprunedToPrunedIndex[uUnprunedNodeIndex] = uInternalNodeIndex;
            uPrunedToUnprunedIndex[uInternalNodeIndex] = uUnprunedNodeIndex;

            ++uInternalNodeIndex;
        }
    }

    const unsigned uUnprunedRootIndex = tree.GetRootNodeIndex();
    uUnprunedToPrunedIndex[uUnprunedRootIndex] = uInternalNodeIndex;
    uPrunedToUnprunedIndex[uInternalNodeIndex] = uUnprunedRootIndex;

#if	TRACE
    {
        Log("Pruned to unpruned:\n");
        for (unsigned i = 0; i < m_uNodeCount; ++i)
            Log(" [%u]=%u", i, uPrunedToUnprunedIndex[i]);
        Log("\n");
        Log("Unpruned to pruned:\n");
        for (unsigned i = 0; i < uUnprunedNodeCount; ++i)
        {
            unsigned n = uUnprunedToPrunedIndex[i];
            if (n != NULL_NEIGHBOR)
                Log(" [%u]=%u", i, n);
        }
        Log("\n");
    }
#endif

    if (uInternalNodeIndex != m_uNodeCount - 1)
        Quit("Tree::PruneTree, Internal error");

// Nodes 0, 1 ... are the leaves
    for (unsigned uSubfamIndex = 0; uSubfamIndex < uSubfamCount; ++uSubfamIndex)
    {
        char szName[32];
        sprintf(szName, "Subfam_%u", uSubfamIndex + 1);
        m_ptrName[uSubfamIndex] = strsave(szName);
    }

    for (unsigned uPrunedNodeIndex = uSubfamCount; uPrunedNodeIndex < m_uNodeCount;
            ++uPrunedNodeIndex)
    {
        unsigned uUnprunedNodeIndex = uPrunedToUnprunedIndex[uPrunedNodeIndex];

        const unsigned uUnprunedLeft = tree.GetLeft(uUnprunedNodeIndex);
        const unsigned uUnprunedRight = tree.GetRight(uUnprunedNodeIndex);

        const unsigned uPrunedLeft = uUnprunedToPrunedIndex[uUnprunedLeft];
        const unsigned uPrunedRight = uUnprunedToPrunedIndex[uUnprunedRight];

        const double dLeftLength =
            tree.GetEdgeLength(uUnprunedNodeIndex, uUnprunedLeft);
        const double dRightLength =
            tree.GetEdgeLength(uUnprunedNodeIndex, uUnprunedRight);

        m_uNeighbor2[uPrunedNodeIndex] = uPrunedLeft;
        m_uNeighbor3[uPrunedNodeIndex] = uPrunedRight;

        m_dEdgeLength1[uPrunedLeft] = dLeftLength;
        m_dEdgeLength1[uPrunedRight] = dRightLength;

//.........这里部分代码省略.........

static void CalcInfo(const Tree &tree, unsigned uNode1, unsigned uNode2, EdgeInfo **EIs)
	{
	const unsigned uNeighborIndex = tree.GetNeighborSubscript(uNode1, uNode2);
	EdgeInfo &EI = EIs[uNode1][uNeighborIndex];
	EI.m_uNode1 = uNode1;
	EI.m_uNode2 = uNode2;

	if (tree.IsLeaf(uNode2))
		{
		EI.m_dMaxDistToLeaf = 0;
		EI.m_dTotalDistToLeaves = 0;
		EI.m_uMaxStep = NULL_NEIGHBOR;
		EI.m_uMostDistantLeaf = uNode2;
		EI.m_uLeafCount = 1;
		EI.m_bSet = true;
		return;
		}

	double dMaxDistToLeaf = -1e29;
	double dTotalDistToLeaves = 0.0;
	unsigned uLeafCount = 0;
	unsigned uMostDistantLeaf = NULL_NEIGHBOR;
	unsigned uMaxStep = NULL_NEIGHBOR;

	const unsigned uNeighborCount = tree.GetNeighborCount(uNode2);
	for (unsigned uSub = 0; uSub < uNeighborCount; ++uSub)
		{
		const unsigned uNode3 = tree.GetNeighbor(uNode2, uSub);
		if (uNode3 == uNode1)
			continue;
		const EdgeInfo &EINext = EIs[uNode2][uSub];
		if (!EINext.m_bSet)
			Quit("CalcInfo: internal error, dist %u->%u not known",
				uNode2, uNode3);


		uLeafCount += EINext.m_uLeafCount;

		const double dEdgeLength = tree.GetEdgeLength(uNode2, uNode3);
		const double dTotalDist = EINext.m_dTotalDistToLeaves +
		  EINext.m_uLeafCount*dEdgeLength;
		dTotalDistToLeaves += dTotalDist;

		const double dDist = EINext.m_dMaxDistToLeaf + dEdgeLength;
		if (dDist > dMaxDistToLeaf)
			{
			dMaxDistToLeaf = dDist;
			uMostDistantLeaf = EINext.m_uMostDistantLeaf;
			uMaxStep = uNode3;
			}
		}
	if (NULL_NEIGHBOR == uMaxStep || NULL_NEIGHBOR == uMostDistantLeaf ||
	  0 == uLeafCount)
		Quit("CalcInfo: internal error 2");

	const double dThisDist = tree.GetEdgeLength(uNode1, uNode2);
	EI.m_dMaxDistToLeaf = dMaxDistToLeaf;
	EI.m_dTotalDistToLeaves = dTotalDistToLeaves;
	EI.m_uMaxStep = uMaxStep;
	EI.m_uMostDistantLeaf = uMostDistantLeaf;
	EI.m_uLeafCount = uLeafCount;
	EI.m_bSet = true;
	}

static void RootByMinAvgLeafDist(const Tree &tree, EdgeInfo **EIs,
  unsigned *ptruNode1, unsigned *ptruNode2,
  double *ptrdLength1, double *ptrdLength2)
	{
	const unsigned uNodeCount = tree.GetNodeCount();
	const unsigned uLeafCount = tree.GetLeafCount();
	unsigned uNode1 = NULL_NEIGHBOR;
	unsigned uNode2 = NULL_NEIGHBOR;
	double dMinHeight = VERY_LARGE_DOUBLE;
	double dBestLength1 = VERY_LARGE_DOUBLE;
	double dBestLength2 = VERY_LARGE_DOUBLE;

	for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
		{
		const unsigned uNeighborCount = tree.GetNeighborCount(uNodeIndex);
		for (unsigned uSub = 0; uSub < uNeighborCount; ++uSub)
			{
			const unsigned uNeighborIndex = tree.GetNeighbor(uNodeIndex, uSub);

		// Avoid visiting same edge a second time in reversed order.
			if (uNeighborIndex < uNodeIndex)
				continue;

			const unsigned uSubRev = tree.GetNeighborSubscript(uNeighborIndex, uNodeIndex);
			if (NULL_NEIGHBOR == uSubRev)
				Quit("RootByMinAvgLeafDist, internal error 1");

		// Get info for edges Node1->Node2 and Node2->Node1 (reversed)
			const EdgeInfo &EI = EIs[uNodeIndex][uSub];
			const EdgeInfo &EIRev = EIs[uNeighborIndex][uSubRev];

			if (EI.m_uNode1 != uNodeIndex || EI.m_uNode2 != uNeighborIndex ||
			  EIRev.m_uNode1 != uNeighborIndex || EIRev.m_uNode2 != uNodeIndex)
				Quit("RootByMinAvgLeafDist, internal error 2");
			if (!EI.m_bSet)
				Quit("RootByMinAvgLeafDist, internal error 3");
			if (uLeafCount != EI.m_uLeafCount + EIRev.m_uLeafCount)
				Quit("RootByMinAvgLeafDist, internal error 4");

			const double dEdgeLength = tree.GetEdgeLength(uNodeIndex, uNeighborIndex);
			if (dEdgeLength != tree.GetEdgeLength(uNeighborIndex, uNodeIndex))
				Quit("RootByMinAvgLeafDist, internal error 5");

		// Consider point p on edge 12 in tree (1=Node, 2=Neighbor).
		//
        //	-----         ----
        //	     |       |
        //	     1----p--2
        //	     |       |
        //	-----         ----
		//
		// Define:
		//    ADLp = average distance to leaves to left of point p.
		//	  ADRp = average distance to leaves to right of point p.
		//	  L = edge length = distance 12
		//    x = distance 1p
		// So distance p2 = L - x.
		// Average distance from p to leaves on left of p is:
		//		ADLp = ADL1 + x
		// Average distance from p to leaves on right of p is:
		//		ADRp = ADR2 + (L - x)
		// To be a root, we require these two distances to be equal,
		//		ADLp = ADRp
		//		ADL1 + x = ADR2 + (L - x)
		// Solving for x,
		//		x = (ADR2 - ADL1 + L)/2
		// If 0 <= x <= L, we can place the root on edge 12.

			const double ADL1 = EI.m_dTotalDistToLeaves / EI.m_uLeafCount;
			const double ADR2 = EIRev.m_dTotalDistToLeaves / EIRev.m_uLeafCount;

			const double x = (ADR2 - ADL1 + dEdgeLength)/2.0;
			if (x >= 0 && x <= dEdgeLength)
				{
				const double dLength1 = x;
				const double dLength2 = dEdgeLength - x;
				const double dHeight1 = EI.m_dMaxDistToLeaf + dLength1;
				const double dHeight2 = EIRev.m_dMaxDistToLeaf + dLength2;
				const double dHeight = dHeight1 >= dHeight2 ? dHeight1 : dHeight2;
#if	TRACE
				Log("Candidate root Node1=%u Node2=%u Height=%g\n",
				  uNodeIndex, uNeighborIndex, dHeight);
#endif
				if (dHeight < dMinHeight)
					{
					uNode1 = uNodeIndex;
					uNode2 = uNeighborIndex;
					dBestLength1 = dLength1;
					dBestLength2 = dLength2;
					dMinHeight = dHeight;
					}
				}
			}
		}

	if (NULL_NEIGHBOR == uNode1 || NULL_NEIGHBOR == uNode2)
		Quit("RootByMinAvgLeafDist, internal error 6");

#if	TRACE
	Log("Best root Node1=%u Node2=%u Length1=%g Length2=%g Height=%g\n",
//.........这里部分代码省略.........

static void RootByMidLongestSpan(const Tree &tree, EdgeInfo **EIs,
  unsigned *ptruNode1, unsigned *ptruNode2,
  double *ptrdLength1, double *ptrdLength2)
	{
	const unsigned uNodeCount = tree.GetNodeCount();

	unsigned uLeaf1 = NULL_NEIGHBOR;
	unsigned uMostDistantLeaf = NULL_NEIGHBOR;
	double dMaxDist = -VERY_LARGE_DOUBLE;
	for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
		{
		if (!tree.IsLeaf(uNodeIndex))
			continue;

		const unsigned uNode2 = tree.GetNeighbor1(uNodeIndex);
		if (NULL_NEIGHBOR == uNode2)
			Quit("RootByMidLongestSpan: internal error 0");
		const double dEdgeLength = tree.GetEdgeLength(uNodeIndex, uNode2);
		const EdgeInfo &EI = EIs[uNodeIndex][0];
		if (!EI.m_bSet)
			Quit("RootByMidLongestSpan: internal error 1");
		if (EI.m_uNode1 != uNodeIndex || EI.m_uNode2 != uNode2)
			Quit("RootByMidLongestSpan: internal error 2");
		const double dSpanLength = dEdgeLength + EI.m_dMaxDistToLeaf;
		if (dSpanLength > dMaxDist)
			{
			dMaxDist = dSpanLength;
			uLeaf1 = uNodeIndex;
			uMostDistantLeaf = EI.m_uMostDistantLeaf;
			}
		}
	
	if (NULL_NEIGHBOR == uLeaf1)
		Quit("RootByMidLongestSpan: internal error 3");

	const double dTreeHeight = dMaxDist/2.0;
	unsigned uNode1 = uLeaf1;
	unsigned uNode2 = tree.GetNeighbor1(uLeaf1);
	double dAccumSpanLength = 0;

#if	TRACE
	Log("RootByMidLongestSpan: span=%u", uLeaf1);
#endif

	for (;;)
		{
		const double dEdgeLength = tree.GetEdgeLength(uNode1, uNode2);
#if	TRACE
		Log("->%u(%g;%g)", uNode2, dEdgeLength, dAccumSpanLength);
#endif
		if (dAccumSpanLength + dEdgeLength >= dTreeHeight)
			{
			*ptruNode1 = uNode1;
			*ptruNode2 = uNode2;
			*ptrdLength1 = dTreeHeight - dAccumSpanLength;
			*ptrdLength2 = dEdgeLength - *ptrdLength1;
#if	TRACE
			{
			const EdgeInfo &EI = EIs[uLeaf1][0];
			Log("...\n");
			Log("Midpoint: Leaf1=%u Leaf2=%u Node1=%u Node2=%u Length1=%g Length2=%g\n",
			  uLeaf1, EI.m_uMostDistantLeaf, *ptruNode1, *ptruNode2, *ptrdLength1, *ptrdLength2);
			}
#endif
			return;
			}

		if (tree.IsLeaf(uNode2))
			Quit("RootByMidLongestSpan: internal error 4");

		dAccumSpanLength += dEdgeLength;
		const unsigned uSub = tree.GetNeighborSubscript(uNode1, uNode2);
		const EdgeInfo &EI = EIs[uNode1][uSub];
		if (!EI.m_bSet)
			Quit("RootByMidLongestSpan: internal error 5");

		uNode1 = uNode2;
		uNode2 = EI.m_uMaxStep;
		}
	}