phy.cpp

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

	{
	assert(IsRooted());

// Descend via left branches until we hit a leaf
	unsigned uNodeIndex = m_uRootNodeIndex;
	while (!IsLeaf(uNodeIndex))
		uNodeIndex = GetLeft(uNodeIndex);
	return uNodeIndex;
	}

unsigned Tree::FirstDepthFirstNodeR() const
	{
	assert(IsRooted());

// Descend via left branches until we hit a leaf
	unsigned uNodeIndex = m_uRootNodeIndex;
	while (!IsLeaf(uNodeIndex))
		uNodeIndex = GetRight(uNodeIndex);
	return uNodeIndex;
	}

unsigned Tree::NextDepthFirstNode(unsigned uNodeIndex) const
	{
#if	TRACE
	Log("NextDepthFirstNode(%3u) ", uNodeIndex);
#endif

	assert(IsRooted());
	assert(uNodeIndex < m_uNodeCount);

	if (IsRoot(uNodeIndex))
		{
#if	TRACE
		Log(">> Node %u is root, end of traversal\n", uNodeIndex);
#endif
		return NULL_NEIGHBOR;
		}

	unsigned uParent = GetParent(uNodeIndex);
	if (GetRight(uParent) == uNodeIndex)
		{
#if	TRACE
		Log(">> Is right branch, return parent=%u\n", uParent);
#endif
		return uParent;
		}

	uNodeIndex = GetRight(uParent);
#if	TRACE
		Log(">> Descend left from right sibling=%u ... ", uNodeIndex);
#endif
	while (!IsLeaf(uNodeIndex))
		uNodeIndex = GetLeft(uNodeIndex);

#if	TRACE
	Log("bottom out at leaf=%u\n", uNodeIndex);
#endif
	return uNodeIndex;
	}

unsigned Tree::NextDepthFirstNodeR(unsigned uNodeIndex) const
	{
#if	TRACE
	Log("NextDepthFirstNode(%3u) ", uNodeIndex);
#endif

	assert(IsRooted());
	assert(uNodeIndex < m_uNodeCount);

	if (IsRoot(uNodeIndex))
		{
#if	TRACE
		Log(">> Node %u is root, end of traversal\n", uNodeIndex);
#endif
		return NULL_NEIGHBOR;
		}

	unsigned uParent = GetParent(uNodeIndex);
	if (GetLeft(uParent) == uNodeIndex)
		{
#if	TRACE
		Log(">> Is left branch, return parent=%u\n", uParent);
#endif
		return uParent;
		}

	uNodeIndex = GetLeft(uParent);
#if	TRACE
		Log(">> Descend right from left sibling=%u ... ", uNodeIndex);
#endif
	while (!IsLeaf(uNodeIndex))
		uNodeIndex = GetRight(uNodeIndex);

#if	TRACE
	Log("bottom out at leaf=%u\n", uNodeIndex);
#endif
	return uNodeIndex;
	}

void Tree::UnrootByDeletingRoot()
	{
	assert(IsRooted());
	assert(m_uNodeCount >= 3);

	const unsigned uLeft = GetLeft(m_uRootNodeIndex);
	const unsigned uRight = GetRight(m_uRootNodeIndex);

	m_uNeighbor1[uLeft] = uRight;
	m_uNeighbor1[uRight] = uLeft;

	bool bHasEdgeLength = HasEdgeLength(m_uRootNodeIndex, uLeft) &&
	  HasEdgeLength(m_uRootNodeIndex, uRight);
	if (bHasEdgeLength)
		{
		double dEdgeLength = GetEdgeLength(m_uRootNodeIndex, uLeft) +
		  GetEdgeLength(m_uRootNodeIndex, uRight);
		m_dEdgeLength1[uLeft] = dEdgeLength;
		m_dEdgeLength1[uRight] = dEdgeLength;
		}

// Remove root node entry from arrays
	const unsigned uMoveCount = m_uNodeCount - m_uRootNodeIndex;
	const unsigned uUnsBytes = uMoveCount*sizeof(unsigned);
	memmove(m_uNeighbor1 + m_uRootNodeIndex, m_uNeighbor1 + m_uRootNodeIndex + 1,
	  uUnsBytes);
	memmove(m_uNeighbor2 + m_uRootNodeIndex, m_uNeighbor2 + m_uRootNodeIndex + 1,
	  uUnsBytes);
	memmove(m_uNeighbor3 + m_uRootNodeIndex, m_uNeighbor3 + m_uRootNodeIndex + 1,
	  uUnsBytes);

	const unsigned uDoubleBytes = uMoveCount*sizeof(double);
	memmove(m_dEdgeLength1 + m_uRootNodeIndex, m_dEdgeLength1 + m_uRootNodeIndex + 1,
	  uDoubleBytes);
	memmove(m_dEdgeLength2 + m_uRootNodeIndex, m_dEdgeLength2 + m_uRootNodeIndex + 1,
	  uDoubleBytes);
	memmove(m_dEdgeLength3 + m_uRootNodeIndex, m_dEdgeLength3 + m_uRootNodeIndex + 1,
	  uDoubleBytes);

	const unsigned uBoolBytes = uMoveCount*sizeof(bool);
	memmove(m_bHasEdgeLength1 + m_uRootNodeIndex, m_bHasEdgeLength1 + m_uRootNodeIndex + 1,
	  uBoolBytes);
	memmove(m_bHasEdgeLength2 + m_uRootNodeIndex, m_bHasEdgeLength2 + m_uRootNodeIndex + 1,
	  uBoolBytes);
	memmove(m_bHasEdgeLength3 + m_uRootNodeIndex, m_bHasEdgeLength3 + m_uRootNodeIndex + 1,
	  uBoolBytes);

	const unsigned uPtrBytes = uMoveCount*sizeof(char *);
	memmove(m_ptrName + m_uRootNodeIndex, m_ptrName + m_uRootNodeIndex + 1, uPtrBytes);

	--m_uNodeCount;
	m_bRooted = false;

// Fix up table entries
	for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)
		{
#define DEC(x)	if (x != NULL_NEIGHBOR && x > m_uRootNodeIndex) --x;
		DEC(m_uNeighbor1[uNodeIndex])
		DEC(m_uNeighbor2[uNodeIndex])
		DEC(m_uNeighbor3[uNodeIndex])
#undef	DEC
		}

	Validate();
	}

unsigned Tree::GetLeafParent(unsigned uNodeIndex) const
	{
	assert(IsLeaf(uNodeIndex));

	if (IsRooted())
		return GetParent(uNodeIndex);

	if (m_uNeighbor1[uNodeIndex] != NULL_NEIGHBOR)
		return m_uNeighbor1[uNodeIndex];
	if (m_uNeighbor2[uNodeIndex] != NULL_NEIGHBOR)
		return m_uNeighbor2[uNodeIndex];
	return m_uNeighbor3[uNodeIndex];
	}

// TODO: This is not efficient for large trees, should cache.
double Tree::GetNodeHeight(unsigned uNodeIndex) const
	{
	if (!IsRooted())
		Quit("Tree::GetNodeHeight: undefined unless rooted tree");
	
	if (IsLeaf(uNodeIndex))
		return 0.0;

	if (m_bHasHeight[uNodeIndex])
		return m_dHeight[uNodeIndex];

	const unsigned uLeft = GetLeft(uNodeIndex);
	const unsigned uRight = GetRight(uNodeIndex);
	double dLeftLength = GetEdgeLength(uNodeIndex, uLeft);
	double dRightLength = GetEdgeLength(uNodeIndex, uRight);

	if (dLeftLength < 0)
		dLeftLength = 0;
	if (dRightLength < 0)
		dRightLength = 0;

	const double dLeftHeight = dLeftLength + GetNodeHeight(uLeft);
	const double dRightHeight = dRightLength + GetNodeHeight(uRight);
	const double dHeight = (dLeftHeight + dRightHeight)/2;
	m_bHasHeight[uNodeIndex] = true;
	m_dHeight[uNodeIndex] = dHeight;
	return dHeight;
	}

unsigned Tree::GetNeighborSubscript(unsigned uNodeIndex, unsigned uNeighborIndex) const
	{
	assert(uNodeIndex < m_uNodeCount);
	assert(uNeighborIndex < m_uNodeCount);
	if (uNeighborIndex == m_uNeighbor1[uNodeIndex])
		return 0;
	if (uNeighborIndex == m_uNeighbor2[uNodeIndex])
		return 1;
	if (uNeighborIndex == m_uNeighbor3[uNodeIndex])
		return 2;
	return NULL_NEIGHBOR;
	}

unsigned Tree::GetNeighbor(unsigned uNodeIndex, unsigned uNeighborSubscript) const
	{
	switch (uNeighborSubscript)
		{
	case 0:
		return m_uNeighbor1[uNodeIndex];
	case 1:
		return m_uNeighbor2[uNodeIndex];
	case 2:
		return m_uNeighbor3[uNodeIndex];
		}
	Quit("Tree::GetNeighbor, sub=%u", uNeighborSubscript);
	return NULL_NEIGHBOR;
	}

// TODO: check if this is a performance issue, could cache a lookup table
unsigned Tree::LeafIndexToNodeIndex(unsigned uLeafIndex) const
	{
	const unsigned uNodeCount = GetNodeCount();
	unsigned uLeafCount = 0;
	for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
		{
		if (IsLeaf(uNodeIndex))
			{
			if (uLeafCount == uLeafIndex)
				return uNodeIndex;
			else
				++uLeafCount;
			}
		}
	Quit("LeafIndexToNodeIndex: out of range");
	return 0;
	}

unsigned Tree::GetLeafNodeIndex(const char *ptrName) const
	{
	const unsigned uNodeCount = GetNodeCount();
	for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
		{
		if (!IsLeaf(uNodeIndex))
			continue;
		const char *ptrLeafName = GetLeafName(uNodeIndex);
		if (0 == strcmp(ptrName, ptrLeafName))
			return uNodeIndex;
		}
	Quit("Tree::GetLeafNodeIndex, name not found");
	return 0;
	}

void Tree::Copy(const Tree &tree)
	{
	const unsigned uNodeCount = tree.GetNodeCount();
	InitCache(uNodeCount);

	m_uNodeCount = uNodeCount;

	const size_t UnsignedBytes = uNodeCount*sizeof(unsigned);
	const size_t DoubleBytes = uNodeCount*sizeof(double);
	const size_t BoolBytes = uNodeCount*sizeof(bool);

	memcpy(m_uNeighbor1, tree.m_uNeighbor1, UnsignedBytes);
	memcpy(m_uNeighbor2, tree.m_uNeighbor2, UnsignedBytes);
	memcpy(m_uNeighbor3, tree.m_uNeighbor3, UnsignedBytes);

	memcpy(m_Ids, tree.m_Ids, UnsignedBytes);

	memcpy(m_dEdgeLength1, tree.m_dEdgeLength1, DoubleBytes);
	memcpy(m_dEdgeLength2, tree.m_dEdgeLength2, DoubleBytes);
	memcpy(m_dEdgeLength3, tree.m_dEdgeLength3, DoubleBytes);

	memcpy(m_dHeight, tree.m_dHeight, DoubleBytes);

	memcpy(m_bHasEdgeLength1, tree.m_bHasEdgeLength1, BoolBytes);
	memcpy(m_bHasEdgeLength2, tree.m_bHasEdgeLength2, BoolBytes);
	memcpy(m_bHasEdgeLength3, tree.m_bHasEdgeLength3, BoolBytes);

	memcpy(m_bHasHeight, tree.m_bHasHeight, BoolBytes);

	m_uRootNodeIndex = tree.m_uRootNodeIndex;
	m_bRooted = tree.m_bRooted;

	for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)
		{
		if (tree.IsLeaf(uNodeIndex))
			{
			const char *ptrName = tree.GetLeafName(uNodeIndex);
			m_ptrName[uNodeIndex] = strsave(ptrName);
			}
		else
			m_ptrName[uNodeIndex] = 0;
		}

#if	DEBUG
	Validate();
#endif
	}

// 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.
void Tree::Create(unsigned uLeafCount, unsigned uRoot, const unsigned Left[],
  const unsigned Right[], const float LeftLength[], const float RightLength[],
  const unsigned LeafIds[], char **LeafNames)
	{
	Clear();

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

	for (unsigned uNodeIndex = 0; uNodeIndex < uLeafCount; ++uNodeIndex)
		{
		m_Ids[uNodeIndex] = LeafIds[uNodeIndex];
		m_ptrName[uNodeIndex] = strsave(LeafNames[uNodeIndex]);
		}

	for (unsigned uNodeIndex = uLeafCount; uNodeIndex < m_uNodeCount; ++uNodeIndex)
		{
		unsigned v = uNodeIndex - uLeafCount;
		unsigned uLeft = Left[v];
		unsigned uRight = Right[v];
		float fLeft = LeftLength[v];
		float fRight = RightLength[v];

		m_uNeighbor2[uNodeIndex] = uLeft;
		m_uNeighbor3[uNodeIndex] = uRight;

		m_bHasEdgeLength2[uNodeIndex] = true;
		m_bHasEdgeLength3[uNodeIndex] = true;

		m_dEdgeLength2[uNodeIndex] = fLeft;
		m_dEdgeLength3[uNodeIndex] = fRight;

		m_uNeighbor1[uLeft] = uNodeIndex;
		m_uNeighbor1[uRight] = uNodeIndex;

		m_dEdgeLength1[uLeft] = fLeft;
		m_dEdgeLength1[uRight] = fRight;

		m_bHasEdgeLength1[uLeft] = true;
		m_bHasEdgeLength1[uRight] = true;
		}

	m_bRooted = true;
	m_uRootNodeIndex = uRoot + uLeafCount;

	Validate();
	}