nwsmall.cpp

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

		Log("\n");
		}
	}

static void ListTB(char *TBM_, const ProfPos *PA, const ProfPos *PB,
  unsigned uPrefixCountA, unsigned uPrefixCountB)
	{
	Log("        ");
	for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
		{
		char c = ' ';
		if (uPrefixLengthB > 0)
			c = ConsensusChar(PB[uPrefixLengthB - 1]);
		Log(" %4u:%c", uPrefixLengthB, c);
		}
	Log("\n");
	for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
		{
		char c = ' ';
		if (uPrefixLengthA > 0)
			c = ConsensusChar(PA[uPrefixLengthA - 1]);
		Log("%4u:%c  ", uPrefixLengthA, c);
		for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
			{
			char c = TBM(uPrefixLengthA, uPrefixLengthB);
			Log(" %6c", c);
			}
		Log("\n");
		}
	}

static const char *BitsToStr(char Bits)
	{
	static char Str[9];

	sprintf(Str, "%cM %cD %cI",
	  Get_M_Char(Bits),
	  Get_D_Char(Bits),
	  Get_I_Char(Bits));
	}
#endif	// TRACE

static inline void SetBitTBM(char **TB, unsigned i, unsigned j, char c)
	{
	char Bit;
	switch (c)
		{
	case 'M':
		Bit = BIT_MM;
		break;
	case 'D':
		Bit = BIT_DM;
		break;
	case 'I':
		Bit = BIT_IM;
		break;
	default:
		Quit("Huh?!");
		}
	TB[i][j] &= ~BIT_xM;
	TB[i][j] |= Bit;
	}

static inline void SetBitTBD(char **TB, unsigned i, unsigned j, char c)
	{
	char Bit;
	switch (c)
		{
	case 'M':
		Bit = BIT_MD;
		break;
	case 'D':
		Bit = BIT_DD;
		break;
	default:
		Quit("Huh?!");
		}
	TB[i][j] &= ~BIT_xD;
	TB[i][j] |= Bit;
	}

static inline void SetBitTBI(char **TB, unsigned i, unsigned j, char c)
	{
	char Bit;
	switch (c)
		{
	case 'M':
		Bit = BIT_MI;
		break;
	case 'I':
		Bit = BIT_II;
		break;
	default:
		Quit("Huh?!");
		}
	TB[i][j] &= ~BIT_xI;
	TB[i][j] |= Bit;
	}

#if	TRACE
#define LogMatrices()											\
	{															\
	Log("Bit DPM:\n");											\
	LogDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);			\
	Log("Bit DPD:\n");											\
	LogDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);			\
	Log("Bit DPI:\n");											\
	LogDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);			\
	Log("Bit TB:\n");											\
	LogBitTB(TB, PA, PB, uPrefixCountA, uPrefixCountB);			\
	bool Same;													\
	Same = DPEq('M', g_DPM, DPM_, uPrefixCountA, uPrefixCountB);\
	if (Same)													\
		Log("DPM success\n");									\
	Same = DPEq('D', g_DPD, DPD_, uPrefixCountA, uPrefixCountB);\
	if (Same)													\
		Log("DPD success\n");									\
	Same = DPEq('I', g_DPI, DPI_, uPrefixCountA, uPrefixCountB);\
	if (Same)													\
		Log("DPI success\n");									\
	CompareTB(TB, g_TBM, g_TBD, g_TBI, uPrefixCountA, uPrefixCountB);\
	}
#else
#define LogMatrices()	/* empty */
#endif

static unsigned uCachePrefixCountB;
static unsigned uCachePrefixCountA;
static SCORE *CacheMCurr;
static SCORE *CacheMNext;
static SCORE *CacheMPrev;
static SCORE *CacheDRow;
static char **CacheTB;

static void AllocCache(unsigned uPrefixCountA, unsigned uPrefixCountB)
	{
	if (uPrefixCountA <= uCachePrefixCountA && uPrefixCountB <= uCachePrefixCountB)
		return;

	delete[] CacheMCurr;
	delete[] CacheMNext;
	delete[] CacheMPrev;
	delete[] CacheDRow;
	for (unsigned i = 0; i < uCachePrefixCountA; ++i)
		delete[] CacheTB[i];
	delete[] CacheTB;

	uCachePrefixCountA = uPrefixCountA + 1024;
	uCachePrefixCountB = uPrefixCountB + 1024;

	CacheMCurr = new SCORE[uCachePrefixCountB];
	CacheMNext = new SCORE[uCachePrefixCountB];
	CacheMPrev = new SCORE[uCachePrefixCountB];
	CacheDRow = new SCORE[uCachePrefixCountB];

	CacheTB = new char *[uCachePrefixCountA];
	for (unsigned i = 0; i < uCachePrefixCountA; ++i)
		CacheTB[i] = new char [uCachePrefixCountB];
	}

SCORE NWSmall(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
  unsigned uLengthB, PWPath &Path)
	{
	if (0 == uLengthB || 0 == uLengthA )
		Quit("Internal error, NWSmall: length=0");

	SetTermGaps(PA, uLengthA);
	SetTermGaps(PB, uLengthB);

	const unsigned uPrefixCountA = uLengthA + 1;
	const unsigned uPrefixCountB = uLengthB + 1;
	const SCORE e = g_scoreGapExtend;

	ALLOC_TRACE()

	AllocCache(uPrefixCountA, uPrefixCountB);

	SCORE *MCurr = CacheMCurr;
	SCORE *MNext = CacheMNext;
	SCORE *MPrev = CacheMPrev;
	SCORE *DRow = CacheDRow;

	char **TB = CacheTB;
	for (unsigned i = 0; i < uPrefixCountA; ++i)
		memset(TB[i], 0, uPrefixCountB);

	SCORE Iij = MINUS_INFINITY;
	SetDPI(0, 0, Iij);

	Iij = PB[0].m_scoreGapOpen;
	SetDPI(0, 1, Iij);

	for (unsigned j = 2; j <= uLengthB; ++j)
		{
		Iij += e;
		SetDPI(0, j, Iij);
		SetTBI(0, j, 'I');
		}

	for (unsigned j = 0; j <= uLengthB; ++j)
		{
		DRow[j] = MINUS_INFINITY;
		SetDPD(0, j, DRow[j]);
		SetTBD(0, j, 'D');
		}

	MPrev[0] = 0;
	SetDPM(0, 0, MPrev[0]);
	for (unsigned j = 1; j <= uLengthB; ++j)
		{
		MPrev[j] = MINUS_INFINITY;
		SetDPM(0, j, MPrev[j]);
		}

	MCurr[0] = MINUS_INFINITY;
	SetDPM(1, 0, MCurr[0]);

	MCurr[1] = ScoreProfPos2(PA[0], PB[0]);
	SetDPM(1, 1, MCurr[1]);
	SetBitTBM(TB, 1, 1, 'M');
	SetTBM(1, 1, 'M');

	for (unsigned j = 2; j <= uLengthB; ++j)
		{
		MCurr[j] = ScoreProfPos2(PA[0], PB[j-1]) + PB[0].m_scoreGapOpen +
		  (j - 2)*e + PB[j-2].m_scoreGapClose;
		SetDPM(1, j, MCurr[j]);
		SetBitTBM(TB, 1, j, 'I');
		SetTBM(1, j, 'I');
		}

// Main DP loop
	for (unsigned i = 1; i < uLengthA; ++i)
		{
		char *TBRow = TB[i];

		Iij = MINUS_INFINITY;
		SetDPI(i, 0, Iij);

		DRow[0] = PA[0].m_scoreGapOpen + (i - 1)*e;
		SetDPD(i, 0, DRow[0]);

		MCurr[0] = MINUS_INFINITY; 
		if (i == 1)
			{
			MCurr[1] = ScoreProfPos2(PA[0], PB[0]);
			SetBitTBM(TB, i, 1, 'M');
			SetTBM(i, 1, 'M');
			}
		else
			{
			MCurr[1] = ScoreProfPos2(PA[i-1], PB[0]) + PA[0].m_scoreGapOpen +
			  (i - 2)*e + PA[i-2].m_scoreGapClose;
			SetBitTBM(TB, i, 1, 'D');
			SetTBM(i, 1, 'D');
			}
		SetDPM(i, 0, MCurr[0]);
		SetDPM(i, 1, MCurr[1]);

		for (unsigned j = 1; j < uLengthB; ++j)
			MNext[j+1] = ScoreProfPos2(PA[i], PB[j]);

		for (unsigned j = 1; j < uLengthB; ++j)
			{
			RECURSE_D(i, j)
			RECURSE_I(i, j)
			RECURSE_M(i, j)
			}
	// Special case for j=uLengthB
		RECURSE_D_BTerm(i)
		RECURSE_I_BTerm(i)

	// Prev := Curr, Curr := Next, Next := Prev
		Rotate(MPrev, MCurr, MNext);
		}

// Special case for i=uLengthA
	char *TBRow = TB[uLengthA];
	MCurr[0] = MINUS_INFINITY;
	if (uLengthA > 1)
		MCurr[1] = ScoreProfPos2(PA[uLengthA-1], PB[0]) + (uLengthA - 2)*e +
		  PA[0].m_scoreGapOpen + PA[uLengthA-2].m_scoreGapClose;
	else
		MCurr[1] = ScoreProfPos2(PA[uLengthA-1], PB[0]) + PA[0].m_scoreGapOpen +
		  PA[0].m_scoreGapClose;
	SetBitTBM(TB, uLengthA, 1, 'D');
	SetTBM(uLengthA, 1, 'D');
	SetDPM(uLengthA, 0, MCurr[0]);
	SetDPM(uLengthA, 1, MCurr[1]);

	DRow[0] = MINUS_INFINITY;
	SetDPD(uLengthA, 0, DRow[0]);
	for (unsigned j = 1; j <= uLengthB; ++j)
		RECURSE_D_ATerm(j);

	Iij = MINUS_INFINITY;
	for (unsigned j = 1; j <= uLengthB; ++j)
		RECURSE_I_ATerm(j)

	LogMatrices();

	SCORE MAB = MCurr[uLengthB];
	SCORE DAB = DRow[uLengthB];
	SCORE IAB = Iij;

	SCORE Score = MAB;
	char cEdgeType = 'M';
	if (DAB > Score)
		{
		Score = DAB;
		cEdgeType = 'D';
		}
	if (IAB > Score)
		{
		Score = IAB;
		cEdgeType = 'I';
		}

#if TRACE
	Log("    Fast: MAB=%.4g DAB=%.4g IAB=%.4g best=%c\n",
	  MAB, DAB, IAB, cEdgeType);
#endif

	BitTraceBack(TB, uLengthA, uLengthB, cEdgeType, Path);

#if	DBEUG
	Path.Validate();
#endif

	return 0;
	}