phy.cpp

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

#include "muscle.h"
#include "tree.h"
#include <math.h>

#define TRACE 0

/***
Node has 0 to 3 neighbors:
	0 neighbors:	singleton root
	1 neighbor:		leaf, neighbor is parent
	2 neigbors:		non-singleton root
	3 neighbors:	internal node (other than root)

Minimal rooted tree is single node.
Minimal unrooted tree is single edge.
Leaf node always has nulls in neighbors 2 and 3, neighbor 1 is parent.
When tree is rooted, neighbor 1=parent, 2=left, 3=right.
***/

void Tree::AssertAreNeighbors(unsigned uNodeIndex1, unsigned uNodeIndex2) const
	{
	if (uNodeIndex1 >= m_uNodeCount || uNodeIndex2 >= m_uNodeCount)
		Quit("AssertAreNeighbors(%u,%u), are %u nodes",
		  uNodeIndex1, uNodeIndex2, m_uNodeCount);

	if (m_uNeighbor1[uNodeIndex1] != uNodeIndex2 &&
	  m_uNeighbor2[uNodeIndex1] != uNodeIndex2 &&
	  m_uNeighbor3[uNodeIndex1] != uNodeIndex2)
		{
		LogMe();
		Quit("AssertAreNeighbors(%u,%u) failed", uNodeIndex1, uNodeIndex2);
		}

	if (m_uNeighbor1[uNodeIndex2] != uNodeIndex1 &&
	  m_uNeighbor2[uNodeIndex2] != uNodeIndex1 &&
	  m_uNeighbor3[uNodeIndex2] != uNodeIndex1)
		{
		LogMe();
		Quit("AssertAreNeighbors(%u,%u) failed", uNodeIndex1, uNodeIndex2);
		}

	bool Has12 = HasEdgeLength(uNodeIndex1, uNodeIndex2);
	bool Has21 = HasEdgeLength(uNodeIndex2, uNodeIndex1);
	if (Has12 != Has21)
		{
		HasEdgeLength(uNodeIndex1, uNodeIndex2);
		HasEdgeLength(uNodeIndex2, uNodeIndex1);
		LogMe();
		Log("HasEdgeLength(%u, %u)=%c HasEdgeLength(%u, %u)=%c\n",
		  uNodeIndex1,
		  uNodeIndex2,
		  Has12 ? 'T' : 'F',
		  uNodeIndex2,
		  uNodeIndex1,
		  Has21 ? 'T' : 'F');

		Quit("Tree::AssertAreNeighbors, HasEdgeLength not symmetric");
		}

	if (Has12)
		{
		double d12 = GetEdgeLength(uNodeIndex1, uNodeIndex2);
		double d21 = GetEdgeLength(uNodeIndex2, uNodeIndex1);
		if (d12 != d21)
			{
			LogMe();
			Quit("Tree::AssertAreNeighbors, Edge length disagrees %u-%u=%.3g, %u-%u=%.3g",
			  uNodeIndex1, uNodeIndex2, d12,
			  uNodeIndex2, uNodeIndex1, d21);
			}
		}
	}

void Tree::ValidateNode(unsigned uNodeIndex) const
	{
	if (uNodeIndex >= m_uNodeCount)
		Quit("ValidateNode(%u), %u nodes", uNodeIndex, m_uNodeCount);

	const unsigned uNeighborCount = GetNeighborCount(uNodeIndex);

	if (2 == uNeighborCount)
		{
		if (!m_bRooted)
			{
			LogMe();
			Quit("Tree::ValidateNode: Node %u has two neighbors, tree is not rooted",
			 uNodeIndex);
			}
		if (uNodeIndex != m_uRootNodeIndex)
			{
			LogMe();
			Quit("Tree::ValidateNode: Node %u has two neighbors, but not root node=%u",
			 uNodeIndex, m_uRootNodeIndex);
			}
		}

	const unsigned n1 = m_uNeighbor1[uNodeIndex];
	const unsigned n2 = m_uNeighbor2[uNodeIndex];
	const unsigned n3 = m_uNeighbor3[uNodeIndex];

	if (NULL_NEIGHBOR == n2 && NULL_NEIGHBOR != n3)
		{
		LogMe();
		Quit("Tree::ValidateNode, n2=null, n3!=null", uNodeIndex);
		}
	if (NULL_NEIGHBOR == n3 && NULL_NEIGHBOR != n2)
		{
		LogMe();
		Quit("Tree::ValidateNode, n3=null, n2!=null", uNodeIndex);
		}

	if (n1 != NULL_NEIGHBOR)
		AssertAreNeighbors(uNodeIndex, n1);
	if (n2 != NULL_NEIGHBOR)
		AssertAreNeighbors(uNodeIndex, n2);
	if (n3 != NULL_NEIGHBOR)
		AssertAreNeighbors(uNodeIndex, n3);

	if (n1 != NULL_NEIGHBOR && (n1 == n2 || n1 == n3))
		{
		LogMe();
		Quit("Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);
		}
	if (n2 != NULL_NEIGHBOR && (n2 == n1 || n2 == n3))
		{
		LogMe();
		Quit("Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);
		}
	if (n3 != NULL_NEIGHBOR && (n3 == n1 || n3 == n2))
		{
		LogMe();
		Quit("Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);
		}

	if (IsRooted())
		{
		if (NULL_NEIGHBOR == GetParent(uNodeIndex))
			{
			if (uNodeIndex != m_uRootNodeIndex)
				{
				LogMe();
				Quit("Tree::ValiateNode(%u), no parent", uNodeIndex);
				}
			}
		else if (GetLeft(GetParent(uNodeIndex)) != uNodeIndex &&
		  GetRight(GetParent(uNodeIndex)) != uNodeIndex)
			{
			LogMe();
			Quit("Tree::ValidateNode(%u), parent / child mismatch", uNodeIndex);
			}
		}
	}

void Tree::Validate() const
	{
	for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)
		ValidateNode(uNodeIndex);
	}

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

	return m_uNeighbor1[uNodeIndex1] == uNodeIndex2 ||
	  m_uNeighbor2[uNodeIndex1] == uNodeIndex2 ||
	  m_uNeighbor3[uNodeIndex1] == uNodeIndex2;
	}

double Tree::GetEdgeLength(unsigned uNodeIndex1, unsigned uNodeIndex2) const
	{
	assert(uNodeIndex1 < m_uNodeCount && uNodeIndex2 < m_uNodeCount);
	if (!HasEdgeLength(uNodeIndex1, uNodeIndex2))
		{
		LogMe();
		Quit("Missing edge length in tree %u-%u", uNodeIndex1, uNodeIndex2);
		}

	if (m_uNeighbor1[uNodeIndex1] == uNodeIndex2)
		return m_dEdgeLength1[uNodeIndex1];
	else if (m_uNeighbor2[uNodeIndex1] == uNodeIndex2)
		return m_dEdgeLength2[uNodeIndex1];
	assert(m_uNeighbor3[uNodeIndex1] == uNodeIndex2);
	return m_dEdgeLength3[uNodeIndex1];
	}

void Tree::ExpandCache()
	{
	const unsigned uNodeCount = 100;
	unsigned uNewCacheCount = m_uCacheCount + uNodeCount;
	unsigned *uNewNeighbor1 = new unsigned[uNewCacheCount];
	unsigned *uNewNeighbor2 = new unsigned[uNewCacheCount];
	unsigned *uNewNeighbor3 = new unsigned[uNewCacheCount];

	unsigned *uNewIds = new unsigned[uNewCacheCount];
	memset(uNewIds, 0xff, uNewCacheCount*sizeof(unsigned));

	double *dNewEdgeLength1 = new double[uNewCacheCount];
	double *dNewEdgeLength2 = new double[uNewCacheCount];
	double *dNewEdgeLength3 = new double[uNewCacheCount];
	double *dNewHeight = new double[uNewCacheCount];

	bool *bNewHasEdgeLength1 = new bool[uNewCacheCount];
	bool *bNewHasEdgeLength2 = new bool[uNewCacheCount];
	bool *bNewHasEdgeLength3 = new bool[uNewCacheCount];
	bool *bNewHasHeight = new bool[uNewCacheCount];

	char **ptrNewName = new char *[uNewCacheCount];
	memset(ptrNewName, 0, uNewCacheCount*sizeof(char *));

	if (m_uCacheCount > 0)
		{
		const unsigned uUnsignedBytes = m_uCacheCount*sizeof(unsigned);
		memcpy(uNewNeighbor1, m_uNeighbor1, uUnsignedBytes);
		memcpy(uNewNeighbor2, m_uNeighbor2, uUnsignedBytes);
		memcpy(uNewNeighbor3, m_uNeighbor3, uUnsignedBytes);

		memcpy(uNewIds, m_Ids, uUnsignedBytes);

		const unsigned uEdgeBytes = m_uCacheCount*sizeof(double);
		memcpy(dNewEdgeLength1, m_dEdgeLength1, uEdgeBytes);
		memcpy(dNewEdgeLength2, m_dEdgeLength2, uEdgeBytes);
		memcpy(dNewEdgeLength3, m_dEdgeLength3, uEdgeBytes);
		memcpy(dNewHeight, m_dHeight, uEdgeBytes);

		const unsigned uBoolBytes = m_uCacheCount*sizeof(bool);
		memcpy(bNewHasEdgeLength1, m_bHasEdgeLength1, uBoolBytes);
		memcpy(bNewHasEdgeLength2, m_bHasEdgeLength2, uBoolBytes);
		memcpy(bNewHasEdgeLength3, m_bHasEdgeLength3, uBoolBytes);
		memcpy(bNewHasHeight, m_bHasHeight, uBoolBytes);

		const unsigned uNameBytes = m_uCacheCount*sizeof(char *);
		memcpy(ptrNewName, m_ptrName, uNameBytes);

		delete[] m_uNeighbor1;
		delete[] m_uNeighbor2;
		delete[] m_uNeighbor3;

		delete[] m_Ids;

		delete[] m_dEdgeLength1;
		delete[] m_dEdgeLength2;
		delete[] m_dEdgeLength3;

		delete[] m_bHasEdgeLength1;
		delete[] m_bHasEdgeLength2;
		delete[] m_bHasEdgeLength3;
		delete[] m_bHasHeight;

		delete[] m_ptrName;
		}
	m_uCacheCount = uNewCacheCount;
	m_uNeighbor1 = uNewNeighbor1;
	m_uNeighbor2 = uNewNeighbor2;
	m_uNeighbor3 = uNewNeighbor3;
	m_Ids = uNewIds;
	m_dEdgeLength1 = dNewEdgeLength1;
	m_dEdgeLength2 = dNewEdgeLength2;
	m_dEdgeLength3 = dNewEdgeLength3;
	m_dHeight = dNewHeight;
	m_bHasEdgeLength1 = bNewHasEdgeLength1;
	m_bHasEdgeLength2 = bNewHasEdgeLength2;
	m_bHasEdgeLength3 = bNewHasEdgeLength3;
	m_bHasHeight = bNewHasHeight;
	m_ptrName = ptrNewName;
	}

// Creates tree with single node, no edges.
// Root node always has index 0.
void Tree::CreateRooted()
	{
	Clear();
	ExpandCache();
	m_uNodeCount = 1;

	m_uNeighbor1[0] = NULL_NEIGHBOR;
	m_uNeighbor2[0] = NULL_NEIGHBOR;
	m_uNeighbor3[0] = NULL_NEIGHBOR;

	m_bHasEdgeLength1[0] = false;
	m_bHasEdgeLength2[0] = false;
	m_bHasEdgeLength3[0] = false;
	m_bHasHeight[0] = false;

	m_uRootNodeIndex = 0;
	m_bRooted = true;

#if	DEBUG
	Validate();
#endif
	}

// Creates unrooted tree with single edge.
// Nodes for that edge are always 0 and 1.
void Tree::CreateUnrooted(double dEdgeLength)
	{
	Clear();
	ExpandCache();

	m_uNeighbor1[0] = 1;
	m_uNeighbor2[0] = NULL_NEIGHBOR;
	m_uNeighbor3[0] = NULL_NEIGHBOR;

	m_uNeighbor1[1] = 0;
	m_uNeighbor2[1] = NULL_NEIGHBOR;
	m_uNeighbor3[1] = NULL_NEIGHBOR;

	m_dEdgeLength1[0] = dEdgeLength;
	m_dEdgeLength1[1] = dEdgeLength;

	m_bHasEdgeLength1[0] = true;
	m_bHasEdgeLength1[1] = true;

	m_bRooted = false;

#if	DEBUG
	Validate();
#endif
	}

void Tree::SetLeafName(unsigned uNodeIndex, const char *ptrName)
	{
	assert(uNodeIndex < m_uNodeCount);
	assert(IsLeaf(uNodeIndex));
	free(m_ptrName[uNodeIndex]);
	m_ptrName[uNodeIndex] = strsave(ptrName);
	}

void Tree::SetLeafId(unsigned uNodeIndex, unsigned uId)
	{
	assert(uNodeIndex < m_uNodeCount);
	assert(IsLeaf(uNodeIndex));
	m_Ids[uNodeIndex] = uId;
	}

const char *Tree::GetLeafName(unsigned uNodeIndex) const
	{
	assert(uNodeIndex < m_uNodeCount);
	assert(IsLeaf(uNodeIndex));
	return m_ptrName[uNodeIndex];
	}

unsigned Tree::GetLeafId(unsigned uNodeIndex) const
	{
	assert(uNodeIndex < m_uNodeCount);
	assert(IsLeaf(uNodeIndex));
	return m_Ids[uNodeIndex];
	}

// 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.
unsigned Tree::AppendBranch(unsigned uExistingLeafIndex)
	{
	if (0 == m_uNodeCount)
		Quit("Tree::AppendBranch: tree has not been created");

#if	DEBUG
	assert(uExistingLeafIndex < m_uNodeCount);
	if (!IsLeaf(uExistingLeafIndex))
		{
		LogMe();
		Quit("AppendBranch(%u): not leaf", uExistingLeafIndex);
		}
#endif

	if (m_uNodeCount >= m_uCacheCount - 2)
		ExpandCache();

	const unsigned uNewLeaf1 = m_uNodeCount;
	const unsigned uNewLeaf2 = m_uNodeCount + 1;

	m_uNodeCount += 2;

	assert(m_uNeighbor2[uExistingLeafIndex] == NULL_NEIGHBOR);
	assert(m_uNeighbor3[uExistingLeafIndex] == NULL_NEIGHBOR);

	m_uNeighbor2[uExistingLeafIndex] = uNewLeaf1;