phy.cpp

unix,linux下编译。用于蛋白质

	m_uNeighbor3[uExistingLeafIndex] = uNewLeaf2;

	m_uNeighbor1[uNewLeaf1] = uExistingLeafIndex;
	m_uNeighbor1[uNewLeaf2] = uExistingLeafIndex;

	m_uNeighbor2[uNewLeaf1] = NULL_NEIGHBOR;
	m_uNeighbor2[uNewLeaf2] = NULL_NEIGHBOR;

	m_uNeighbor3[uNewLeaf1] = NULL_NEIGHBOR;
	m_uNeighbor3[uNewLeaf2] = NULL_NEIGHBOR;

	m_dEdgeLength2[uExistingLeafIndex] = 0;
	m_dEdgeLength3[uExistingLeafIndex] = 0;

	m_dEdgeLength1[uNewLeaf1] = 0;
	m_dEdgeLength2[uNewLeaf1] = 0;
	m_dEdgeLength3[uNewLeaf1] = 0;

	m_dEdgeLength1[uNewLeaf2] = 0;
	m_dEdgeLength2[uNewLeaf2] = 0;
	m_dEdgeLength3[uNewLeaf2] = 0;

	m_bHasEdgeLength1[uNewLeaf1] = false;
	m_bHasEdgeLength2[uNewLeaf1] = false;
	m_bHasEdgeLength3[uNewLeaf1] = false;

	m_bHasEdgeLength1[uNewLeaf2] = false;
	m_bHasEdgeLength2[uNewLeaf2] = false;
	m_bHasEdgeLength3[uNewLeaf2] = false;

	m_bHasHeight[uNewLeaf1] = false;
	m_bHasHeight[uNewLeaf2] = false;

	m_Ids[uNewLeaf1] = uInsane;
	m_Ids[uNewLeaf2] = uInsane;
	return uNewLeaf1;
	}

void Tree::LogMe() const
	{
	Log("Tree::LogMe %u nodes, ", m_uNodeCount);

	if (IsRooted())
		{
		Log("rooted.\n");
		Log("\n");
		Log("Index  Parnt  LengthP  Left   LengthL  Right  LengthR     Id  Name\n");
		Log("-----  -----  -------  ----   -------  -----  -------  -----  ----\n");
		}
	else
		{
		Log("unrooted.\n");
		Log("\n");
		Log("Index  Nbr_1  Length1  Nbr_2  Length2  Nbr_3  Length3     Id  Name\n");
		Log("-----  -----  -------  -----  -------  -----  -------  -----  ----\n");
		}

	for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)
		{
		Log("%5u  ", uNodeIndex);
		const unsigned n1 = m_uNeighbor1[uNodeIndex];
		const unsigned n2 = m_uNeighbor2[uNodeIndex];
		const unsigned n3 = m_uNeighbor3[uNodeIndex];
		if (NULL_NEIGHBOR != n1)
			{
			Log("%5u  ", n1);
			if (m_bHasEdgeLength1[uNodeIndex])
				Log("%7.3g  ", m_dEdgeLength1[uNodeIndex]);
			else
				Log("      *  ");
			}
		else
			Log("                ");

		if (NULL_NEIGHBOR != n2)
			{
			Log("%5u  ", n2);
			if (m_bHasEdgeLength2[uNodeIndex])
				Log("%7.3g  ", m_dEdgeLength2[uNodeIndex]);
			else
				Log("      *  ");
			}
		else
			Log("                ");

		if (NULL_NEIGHBOR != n3)
			{
			Log("%5u  ", n3);
			if (m_bHasEdgeLength3[uNodeIndex])
				Log("%7.3g  ", m_dEdgeLength3[uNodeIndex]);
			else
				Log("      *  ");
			}
		else
			Log("                ");

		if (m_Ids != 0 && IsLeaf(uNodeIndex))
			{
			unsigned uId = m_Ids[uNodeIndex];
			if (uId == uInsane)
				Log("    *");
			else
				Log("%5u", uId);
			}
		else
			Log("     ");

		if (m_bRooted && uNodeIndex == m_uRootNodeIndex)
			Log("  [ROOT] ");
		const char *ptrName = m_ptrName[uNodeIndex];
		if (ptrName != 0)
			Log("  %s", ptrName);
		Log("\n");
		}
	}

void Tree::SetEdgeLength(unsigned uNodeIndex1, unsigned uNodeIndex2,
  double dLength)
	{
	assert(uNodeIndex1 < m_uNodeCount && uNodeIndex2 < m_uNodeCount);
	assert(IsEdge(uNodeIndex1, uNodeIndex2));

	if (m_uNeighbor1[uNodeIndex1] == uNodeIndex2)
		{
		m_dEdgeLength1[uNodeIndex1] = dLength;
		m_bHasEdgeLength1[uNodeIndex1] = true;
		}
	else if (m_uNeighbor2[uNodeIndex1] == uNodeIndex2)
		{
		m_dEdgeLength2[uNodeIndex1] = dLength;
		m_bHasEdgeLength2[uNodeIndex1] = true;

		}
	else
		{
		assert(m_uNeighbor3[uNodeIndex1] == uNodeIndex2);
		m_dEdgeLength3[uNodeIndex1] = dLength;
		m_bHasEdgeLength3[uNodeIndex1] = true;
		}

	if (m_uNeighbor1[uNodeIndex2] == uNodeIndex1)
		{
		m_dEdgeLength1[uNodeIndex2] = dLength;
		m_bHasEdgeLength1[uNodeIndex2] = true;
		}
	else if (m_uNeighbor2[uNodeIndex2] == uNodeIndex1)
		{
		m_dEdgeLength2[uNodeIndex2] = dLength;
		m_bHasEdgeLength2[uNodeIndex2] = true;
		}
	else
		{
		assert(m_uNeighbor3[uNodeIndex2] == uNodeIndex1);
		m_dEdgeLength3[uNodeIndex2] = dLength;
		m_bHasEdgeLength3[uNodeIndex2] = true;
		}
	}

unsigned Tree::UnrootFromFile()
	{
#if	TRACE
	Log("Before unroot:\n");
	LogMe();
#endif

	if (!m_bRooted)
		Quit("Tree::Unroot, not rooted");

// Convention: root node is always node zero
	assert(IsRoot(0));
	assert(NULL_NEIGHBOR == m_uNeighbor1[0]);

	const unsigned uThirdNode = m_uNodeCount++;

	m_uNeighbor1[0] = uThirdNode;
	m_uNeighbor1[uThirdNode] = 0;

	m_uNeighbor2[uThirdNode] = NULL_NEIGHBOR;
	m_uNeighbor3[uThirdNode] = NULL_NEIGHBOR;

	m_dEdgeLength1[0] = 0;
	m_dEdgeLength1[uThirdNode] = 0;
	m_bHasEdgeLength1[uThirdNode] = true;

	m_bRooted = false;

#if	TRACE
	Log("After unroot:\n");
	LogMe();
#endif

	return uThirdNode;
	}

// In an unrooted tree, equivalent of GetLeft/Right is 
// GetFirst/SecondNeighbor.
// uNeighborIndex must be a known neighbor of uNodeIndex.
// This is the way to find the other two neighbor nodes of
// an internal node.
// The labeling as "First" and "Second" neighbor is arbitrary.
// Calling these functions on a leaf returns NULL_NEIGHBOR, as
// for GetLeft/Right.
unsigned Tree::GetFirstNeighbor(unsigned uNodeIndex, unsigned uNeighborIndex) const
	{
	assert(uNodeIndex < m_uNodeCount);
	assert(uNeighborIndex < m_uNodeCount);
	assert(IsEdge(uNodeIndex, uNeighborIndex));

	for (unsigned n = 0; n < 3; ++n)
		{
		unsigned uNeighbor = GetNeighbor(uNodeIndex, n);
		if (NULL_NEIGHBOR != uNeighbor && uNeighborIndex != uNeighbor)
			return uNeighbor;
		}
	return NULL_NEIGHBOR;
	}

unsigned Tree::GetSecondNeighbor(unsigned uNodeIndex, unsigned uNeighborIndex) const
	{
	assert(uNodeIndex < m_uNodeCount);
	assert(uNeighborIndex < m_uNodeCount);
	assert(IsEdge(uNodeIndex, uNeighborIndex));

	bool bFoundOne = false;
	for (unsigned n = 0; n < 3; ++n)
		{
		unsigned uNeighbor = GetNeighbor(uNodeIndex, n);
		if (NULL_NEIGHBOR != uNeighbor && uNeighborIndex != uNeighbor)
			{
			if (bFoundOne)
				return uNeighbor;
			else
				bFoundOne = true;
			}
		}
	return NULL_NEIGHBOR;
	}

// Compute the number of leaves in the sub-tree defined by an edge
// in an unrooted tree. Conceptually, the tree is cut at this edge,
// and uNodeIndex2 considered the root of the sub-tree.
unsigned Tree::GetLeafCountUnrooted(unsigned uNodeIndex1, unsigned uNodeIndex2,
  double *ptrdTotalDistance) const
	{
	assert(!IsRooted());

	if (IsLeaf(uNodeIndex2))
		{
		*ptrdTotalDistance = GetEdgeLength(uNodeIndex1, uNodeIndex2);
		return 1;
		}

// Recurse down the rooted sub-tree defined by cutting the edge
// and considering uNodeIndex2 as the root.
	const unsigned uLeft = GetFirstNeighbor(uNodeIndex2, uNodeIndex1);
	const unsigned uRight = GetSecondNeighbor(uNodeIndex2, uNodeIndex1);

	double dLeftDistance;
	double dRightDistance;

	const unsigned uLeftCount = GetLeafCountUnrooted(uNodeIndex2, uLeft,
	  &dLeftDistance);
	const unsigned uRightCount = GetLeafCountUnrooted(uNodeIndex2, uRight,
	  &dRightDistance);

	*ptrdTotalDistance = dLeftDistance + dRightDistance;
	return uLeftCount + uRightCount;
	}

void Tree::RootUnrootedTree(ROOT Method)
	{
	assert(!IsRooted());
#if TRACE
	Log("Tree::RootUnrootedTree, before:");
	LogMe();
#endif

	unsigned uNode1;
	unsigned uNode2;
	double dLength1;
	double dLength2;
	FindRoot(*this, &uNode1, &uNode2, &dLength1, &dLength2, Method);

	if (m_uNodeCount == m_uCacheCount)
		ExpandCache();
	m_uRootNodeIndex = m_uNodeCount++;

	double dEdgeLength = GetEdgeLength(uNode1, uNode2);

	m_uNeighbor1[m_uRootNodeIndex] = NULL_NEIGHBOR;
	m_uNeighbor2[m_uRootNodeIndex] = uNode1;
	m_uNeighbor3[m_uRootNodeIndex] = uNode2;

	if (m_uNeighbor1[uNode1] == uNode2)
		m_uNeighbor1[uNode1] = m_uRootNodeIndex;
	else if (m_uNeighbor2[uNode1] == uNode2)
		m_uNeighbor2[uNode1] = m_uRootNodeIndex;
	else
		{
		assert(m_uNeighbor3[uNode1] == uNode2);
		m_uNeighbor3[uNode1] = m_uRootNodeIndex;
		}

	if (m_uNeighbor1[uNode2] == uNode1)
		m_uNeighbor1[uNode2] = m_uRootNodeIndex;
	else if (m_uNeighbor2[uNode2] == uNode1)
		m_uNeighbor2[uNode2] = m_uRootNodeIndex;
	else
		{
		assert(m_uNeighbor3[uNode2] == uNode1);
		m_uNeighbor3[uNode2] = m_uRootNodeIndex;
		}

	OrientParent(uNode1, m_uRootNodeIndex);
	OrientParent(uNode2, m_uRootNodeIndex);

	SetEdgeLength(m_uRootNodeIndex, uNode1, dLength1);
	SetEdgeLength(m_uRootNodeIndex, uNode2, dLength2);

	m_bHasHeight[m_uRootNodeIndex] = false;

	m_ptrName[m_uRootNodeIndex] = 0;

	m_bRooted = true;

#if	TRACE
	Log("\nPhy::RootUnrootedTree, after:");
	LogMe();
#endif

	Validate();
	}

bool Tree::HasEdgeLength(unsigned uNodeIndex1, unsigned uNodeIndex2) const
	{
	assert(uNodeIndex1 < m_uNodeCount);
	assert(uNodeIndex2 < m_uNodeCount);
	assert(IsEdge(uNodeIndex1, uNodeIndex2));

	if (m_uNeighbor1[uNodeIndex1] == uNodeIndex2)
		return m_bHasEdgeLength1[uNodeIndex1];
	else if (m_uNeighbor2[uNodeIndex1] == uNodeIndex2)
		return m_bHasEdgeLength2[uNodeIndex1];
	assert(m_uNeighbor3[uNodeIndex1] == uNodeIndex2);
	return m_bHasEdgeLength3[uNodeIndex1];
	}

void Tree::OrientParent(unsigned uNodeIndex, unsigned uParentNodeIndex)
	{
	if (NULL_NEIGHBOR == uNodeIndex)
		return;

	if (m_uNeighbor1[uNodeIndex] == uParentNodeIndex)
		;
	else if (m_uNeighbor2[uNodeIndex] == uParentNodeIndex)
		{
		double dEdgeLength2 = m_dEdgeLength2[uNodeIndex];
		m_uNeighbor2[uNodeIndex] = m_uNeighbor1[uNodeIndex];
		m_dEdgeLength2[uNodeIndex] = m_dEdgeLength1[uNodeIndex];
		m_uNeighbor1[uNodeIndex] = uParentNodeIndex;
		m_dEdgeLength1[uNodeIndex] = dEdgeLength2;
		}
	else
		{
		assert(m_uNeighbor3[uNodeIndex] == uParentNodeIndex);
		double dEdgeLength3 = m_dEdgeLength3[uNodeIndex];
		m_uNeighbor3[uNodeIndex] = m_uNeighbor1[uNodeIndex];
		m_dEdgeLength3[uNodeIndex] = m_dEdgeLength1[uNodeIndex];
		m_uNeighbor1[uNodeIndex] = uParentNodeIndex;
		m_dEdgeLength1[uNodeIndex] = dEdgeLength3;
		}

	OrientParent(m_uNeighbor2[uNodeIndex], uNodeIndex);
	OrientParent(m_uNeighbor3[uNodeIndex], uNodeIndex);
	}

unsigned Tree::FirstDepthFirstNode() const