C++ PWPath类代码示例

SCORE GlobalAlign(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
  unsigned uLengthB, PWPath &Path)
	{
#if	TIMING
	TICKS t1 = GetClockTicks();
#endif
	g_bKeepSimpleDP = true;
	PWPath SimplePath;
	GlobalAlignSimple(PA, uLengthA, PB, uLengthB, SimplePath);

	SCORE Score = NWSmall(PA, uLengthA, PB, uLengthB, Path);

	if (!Path.Equal(SimplePath))
		{
		Log("Simple:\n");
		SimplePath.LogMe();
		Log("Small:\n");
		Path.LogMe();
		Quit("Paths differ");
		}

#if	TIMING
	TICKS t2 = GetClockTicks();
	g_ticksDP += (t2 - t1);
#endif
	return Score;
	}

static void AppendRegPath(PWPath &Path, const PWPath &RegPath)
	{
	const unsigned uRegEdgeCount = RegPath.GetEdgeCount();
	for (unsigned uRegEdgeIndex = 0; uRegEdgeIndex < uRegEdgeCount; ++uRegEdgeIndex)
		{
		const PWEdge &RegEdge = RegPath.GetEdge(uRegEdgeIndex);
		Path.AppendEdge(RegEdge);
		}
	}

static void AllDeletes(PWPath &Path, unsigned uLengthA)
	{
	Path.Clear();
	PWEdge Edge;
	Edge.cType = 'D';
	Edge.uPrefixLengthB = 0;
	for (unsigned uPrefixLengthA = 1; uPrefixLengthA <= uLengthA; ++uPrefixLengthA)
		{
		Edge.uPrefixLengthA = uPrefixLengthA;
		Path.AppendEdge(Edge);
		}
	}

static void AllInserts(PWPath &Path, unsigned uLengthB)
	{
	Path.Clear();
	PWEdge Edge;
	Edge.cType = 'I';
	Edge.uPrefixLengthA = 0;
	for (unsigned uPrefixLengthB = 1; uPrefixLengthB <= uLengthB; ++uPrefixLengthB)
		{
		Edge.uPrefixLengthB = uPrefixLengthB;
		Path.AppendEdge(Edge);
		}
	}

static void DiagToPath(const Diag &d, PWPath &Path)
	{
	Path.Clear();
	const unsigned uLength = d.m_uLength;
	for (unsigned i = 0; i < uLength; ++i)
		{
		PWEdge Edge;
		Edge.cType = 'M';
		Edge.uPrefixLengthA = d.m_uStartPosA + i + 1;
		Edge.uPrefixLengthB = d.m_uStartPosB + i + 1;
		Path.AppendEdge(Edge);
		}
	}

static void OffsetPath(PWPath &Path, unsigned uOffsetA, unsigned uOffsetB)
	{
	const unsigned uEdgeCount = Path.GetEdgeCount();
	for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)
		{
		const PWEdge &Edge = Path.GetEdge(uEdgeIndex);

	// Nasty hack -- poke new values back into path, circumventing class
		PWEdge &NonConstEdge = (PWEdge &) Edge;
		NonConstEdge.uPrefixLengthA += uOffsetA;
		NonConstEdge.uPrefixLengthB += uOffsetB;
		}
	}

static void PathSeq(const Seq &s, const PWPath &Path, bool bRight, Seq &sOut)
	{
	short *esA;
	short *esB;
	PathToEstrings(Path, &esA, &esB);

	const unsigned uSeqLength = s.Length();
	const unsigned uEdgeCount = Path.GetEdgeCount();

	sOut.Clear();
	sOut.SetName(s.GetName());
	unsigned uPos = 0;
	for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)
		{
		const PWEdge &Edge = Path.GetEdge(uEdgeIndex);
		char cType = Edge.cType;
		if (bRight)
			{
			if (cType == 'I')
				cType = 'D';
			else if (cType == 'D')
				cType = 'I';
			}
		switch (cType)
			{
		case 'M':
			sOut.AppendChar(s[uPos++]);
			break;
		case 'D':
			sOut.AppendChar('-');
			break;
		case 'I':
			sOut.AppendChar(s[uPos++]);
			break;
		default:
			Quit("PathSeq, invalid edge type %c", cType);
			}
		}
	}

void DiagList::FromPath(const PWPath &Path)
	{
    unsigned &g_uMinDiagLength = getMuscleContext()->params.g_uMinDiagLength;

	Clear();

	const unsigned uEdgeCount = Path.GetEdgeCount();
	unsigned uLength = 0;
	unsigned uStartPosA;
	unsigned uStartPosB;
	for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)
		{
		const PWEdge &Edge = Path.GetEdge(uEdgeIndex);

	// Typical cases
		if (Edge.cType == 'M')
			{
			if (0 == uLength)
				{
				uStartPosA = Edge.uPrefixLengthA - 1;
				uStartPosB = Edge.uPrefixLengthB - 1;
				}
			++uLength;
			}
		else
			{
			if (uLength >= g_uMinDiagLength)
				Add(uStartPosA, uStartPosB, uLength);
			uLength = 0;
			}
		}

// Special case for last edge
	if (uLength >= g_uMinDiagLength)
		Add(uStartPosA, uStartPosB, uLength);
	}

    bool TryRealign(MSA &msaIn, const Tree &tree, const unsigned Leaves1[],
        unsigned uCount1, const unsigned Leaves2[], unsigned uCount2,
        SCORE *ptrscoreBefore, SCORE *ptrscoreAfter,
        bool bLockLeft, bool bLockRight)
    {
#if	TRACE
        Log("TryRealign, msaIn=\n");
#endif

        MuscleContext *ctx = getMuscleContext();

        const unsigned uSeqCount = msaIn.GetSeqCount();

        unsigned *Ids1 = new unsigned[uSeqCount];
        unsigned *Ids2 = new unsigned[uSeqCount];

        LeafIndexesToIds(tree, Leaves1, uCount1, Ids1);
        LeafIndexesToIds(tree, Leaves2, uCount2, Ids2);

        MSA msa1;
        MSA msa2;

        MSASubsetByIds(msaIn, Ids1, uCount1, msa1);
        MSASubsetByIds(msaIn, Ids2, uCount2, msa2);

#if	DEBUG
        ValidateMuscleIds(msa1);
        ValidateMuscleIds(msa2);
#endif

        // Computing the objective score may be expensive for
        // large numbers of sequences. As a speed optimization,
        // we check whether the alignment changes. If it does
        // not change, there is no need to compute the objective
        // score. We test for the alignment changing by comparing
        // the Viterbi paths before and after re-aligning.
        PWPath pathBefore;
        pathBefore.FromMSAPair(msa1, msa2);

        DeleteGappedCols(msa1);
        DeleteGappedCols(msa2);

        if (0 == msa1.GetColCount() || 0 == msa2.GetColCount()) {
            delete[] Ids1;
            delete[] Ids2;
            return false;
        }

        MSA msaRealigned;
        PWPath pathAfter;

        AlignTwoMSAs(msa1, msa2, msaRealigned, pathAfter, bLockLeft, bLockRight);

        bool bAnyChanges = !pathAfter.Equal(pathBefore);
        unsigned uDiffCount1;
        unsigned uDiffCount2;
        unsigned* Edges1 = ctx->refinehoriz.Edges1;
        unsigned* Edges2 = ctx->refinehoriz.Edges2;
        DiffPaths(pathBefore, pathAfter, Edges1, &uDiffCount1, Edges2, &uDiffCount2);

#if	TRACE
        Log("TryRealign, msa1=\n");
        Log("\nmsa2=\n");
        Log("\nRealigned (changes %s)=\n", bAnyChanges ? "TRUE" : "FALSE");
#endif

        if (!bAnyChanges)
        {
            *ptrscoreBefore = 0;
            *ptrscoreAfter = 0;
            delete[] Ids1;
            delete[] Ids2;
            return false;
        }

        SetMSAWeightsMuscle(msaIn);
        SetMSAWeightsMuscle(msaRealigned);

#if	DIFFOBJSCORE
        const SCORE scoreDiff = DiffObjScore(msaIn, pathBefore, Edges1, uDiffCount1,
            msaRealigned, pathAfter, Edges2, uDiffCount2);
        bool bAccept = (scoreDiff > 0);
        *ptrscoreBefore = 0;
        *ptrscoreAfter = scoreDiff;
        //const SCORE scoreBefore = ObjScoreIds(msaIn, Ids1, uCount1, Ids2, uCount2);
        //const SCORE scoreAfter = ObjScoreIds(msaRealigned, Ids1, uCount1, Ids2, uCount2);
        //Log("Diff = %.3g %.3g\n", scoreDiff, scoreAfter - scoreBefore);
#else
        const SCORE scoreBefore = ObjScoreIds(msaIn, Ids1, uCount1, Ids2, uCount2);
        const SCORE scoreAfter = ObjScoreIds(msaRealigned, Ids1, uCount1, Ids2, uCount2);

        bool bAccept = (scoreAfter > scoreBefore);

#if	TRACE
        Log("Score %g -> %g Accept %s\n", scoreBefore, scoreAfter, bAccept ? "TRUE" : "FALSE");
#endif

        *ptrscoreBefore = scoreBefore;
        *ptrscoreAfter = scoreAfter;
#endif
//.........这里部分代码省略.........

ProgNode *ProgressiveAlignE(const SeqVect &v, const Tree &GuideTree, MSA &a)
	{
	assert(GuideTree.IsRooted());

#if	TRACE
	Log("GuideTree:\n");
	GuideTree.LogMe();
#endif

	const unsigned uSeqCount = v.Length();
	const unsigned uNodeCount = 2*uSeqCount - 1;
	const unsigned uIterCount = uSeqCount - 1;

	WEIGHT *Weights = new WEIGHT[uSeqCount];
	CalcClustalWWeights(GuideTree, Weights);

	ProgNode *ProgNodes = new ProgNode[uNodeCount];

	unsigned uJoin = 0;
	unsigned uTreeNodeIndex = GuideTree.FirstDepthFirstNode();
	SetProgressDesc("Align node");
	do
		{
		if (GuideTree.IsLeaf(uTreeNodeIndex))
			{
			if (uTreeNodeIndex >= uNodeCount)
				Quit("TreeNodeIndex=%u NodeCount=%u\n", uTreeNodeIndex, uNodeCount);
			ProgNode &Node = ProgNodes[uTreeNodeIndex];
			unsigned uId = GuideTree.GetLeafId(uTreeNodeIndex);
			if (uId >= uSeqCount)
				Quit("Seq index out of range");
			const Seq &s = *(v[uId]);
			Node.m_MSA.FromSeq(s);
			Node.m_MSA.SetSeqId(0, uId);
			Node.m_uLength = Node.m_MSA.GetColCount();
			Node.m_Weight = Weights[uId];
		// TODO: Term gaps settable
			Node.m_Prof = ProfileFromMSA(Node.m_MSA);
			Node.m_EstringL = 0;
			Node.m_EstringR = 0;
#if	TRACE
			Log("Leaf id=%u\n", uId);
			Log("MSA=\n");
			Node.m_MSA.LogMe();
			Log("Profile (from MSA)=\n");
			ListProfile(Node.m_Prof, Node.m_uLength, &Node.m_MSA);
#endif
			}
		else
			{
			Progress(uJoin, uSeqCount - 1);
			++uJoin;

			const unsigned uMergeNodeIndex = uTreeNodeIndex;
			ProgNode &Parent = ProgNodes[uMergeNodeIndex];

			const unsigned uLeft = GuideTree.GetLeft(uTreeNodeIndex);
			const unsigned uRight = GuideTree.GetRight(uTreeNodeIndex);

			if (g_bVerbose)
				{
				Log("Align: (");
				LogLeafNames(GuideTree, uLeft);
				Log(") (");
				LogLeafNames(GuideTree, uRight);
				Log(")\n");
				}

			ProgNode &Node1 = ProgNodes[uLeft];
			ProgNode &Node2 = ProgNodes[uRight];

#if	TRACE
			Log("AlignTwoMSAs:\n");
#endif
			AlignTwoProfs(
			  Node1.m_Prof, Node1.m_uLength, Node1.m_Weight,
			  Node2.m_Prof, Node2.m_uLength, Node2.m_Weight,
			  Parent.m_Path,
			  &Parent.m_Prof, &Parent.m_uLength);
#if	TRACE_LENGTH_DELTA
			{
			unsigned L = Node1.m_uLength;
			unsigned R = Node2.m_uLength;
			unsigned P = Parent.m_Path.GetEdgeCount();
			unsigned Max = L > R ? L : R;
			unsigned d = P - Max;
			Log("LD%u;%u;%u;%u\n", L, R, P, d);
			}
#endif
			PathToEstrings(Parent.m_Path, &Parent.m_EstringL, &Parent.m_EstringR);

			Parent.m_Weight = Node1.m_Weight + Node2.m_Weight;

#if	VALIDATE
			{
#if	TRACE
			Log("AlignTwoMSAs:\n");
#endif
			PWPath TmpPath;
			AlignTwoMSAs(Node1.m_MSA, Node2.m_MSA, Parent.m_MSA, TmpPath);
//.........这里部分代码省略.........

void PathToEstrings(const PWPath &Path, short **ptresA, short **ptresB)
	{
// First pass to determine size of estrings esA and esB
	const unsigned uEdgeCount = Path.GetEdgeCount();
	if (0 == uEdgeCount)
		{
		short *esA = new short[1];
		short *esB = new short[1];
		esA[0] = 0;
		esB[0] = 0;
		*ptresA = esA;
		*ptresB = esB;
		return;
		}

	unsigned iLengthA = 1;
	unsigned iLengthB = 1;
	const char cFirstEdgeType = Path.GetEdge(0).cType;
	char cPrevEdgeType = cFirstEdgeType;
	for (unsigned uEdgeIndex = 1; uEdgeIndex < uEdgeCount; ++uEdgeIndex)
		{
		const PWEdge &Edge = Path.GetEdge(uEdgeIndex);
		char cEdgeType = Edge.cType;

		switch (c2(cPrevEdgeType, cEdgeType))
			{
		case c2('M', 'M'):
		case c2('D', 'D'):
		case c2('I', 'I'):
			break;

		case c2('D', 'M'):
		case c2('M', 'D'):
			++iLengthB;
			break;

		case c2('I', 'M'):
		case c2('M', 'I'):
			++iLengthA;
			break;

		case c2('I', 'D'):
		case c2('D', 'I'):
			++iLengthB;
			++iLengthA;
			break;

		default:
			assert(false);
			}
		cPrevEdgeType = cEdgeType;
		}

// Pass2 for seq A
	{
	short *esA = new short[iLengthA+1];
	unsigned iA = 0;
	switch (Path.GetEdge(0).cType)
		{
	case 'M':
	case 'D':
		esA[0] = 1;
		break;

	case 'I':
		esA[0] = -1;
		break;

	default:
		assert(false);
		}

	char cPrevEdgeType = cFirstEdgeType;
	for (unsigned uEdgeIndex = 1; uEdgeIndex < uEdgeCount; ++uEdgeIndex)
		{
		const PWEdge &Edge = Path.GetEdge(uEdgeIndex);
		char cEdgeType = Edge.cType;

		switch (c2(cPrevEdgeType, cEdgeType))
			{
		case c2('M', 'M'):
		case c2('D', 'D'):
		case c2('D', 'M'):
		case c2('M', 'D'):
			++(esA[iA]);
			break;

		case c2('I', 'D'):
		case c2('I', 'M'):
			++iA;
			esA[iA] = 1;
			break;

		case c2('M', 'I'):
		case c2('D', 'I'):
			++iA;
			esA[iA] = -1;
			break;

		case c2('I', 'I'):
//.........这里部分代码省略.........

void AlignTwoMSAsGivenPathSW(const PWPath &Path, const MSA &msaA, const MSA &msaB,
  MSA &msaCombined)
	{
	msaCombined.Clear();

#if	TRACE
	Log("AlignTwoMSAsGivenPathSW\n");
	Log("Template A:\n");
	msaA.LogMe();
	Log("Template B:\n");
	msaB.LogMe();
#endif

	const unsigned uColCountA = msaA.GetColCount();
	const unsigned uColCountB = msaB.GetColCount();

	const unsigned uSeqCountA = msaA.GetSeqCount();
	const unsigned uSeqCountB = msaB.GetSeqCount();

	msaCombined.SetSeqCount(uSeqCountA + uSeqCountB);

// Copy sequence names into combined MSA
	for (unsigned uSeqIndexA = 0; uSeqIndexA < uSeqCountA; ++uSeqIndexA)
		{
		msaCombined.SetSeqName(uSeqIndexA, msaA.GetSeqName(uSeqIndexA));
		msaCombined.SetSeqId(uSeqIndexA, msaA.GetSeqId(uSeqIndexA));
		}

	for (unsigned uSeqIndexB = 0; uSeqIndexB < uSeqCountB; ++uSeqIndexB)
		{
		msaCombined.SetSeqName(uSeqCountA + uSeqIndexB, msaB.GetSeqName(uSeqIndexB));
		msaCombined.SetSeqId(uSeqCountA + uSeqIndexB, msaB.GetSeqId(uSeqIndexB));
		}

	unsigned uColIndexA = 0;
	unsigned uColIndexB = 0;
	unsigned uColIndexCombined = 0;
	const unsigned uEdgeCount = Path.GetEdgeCount();
	for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)
		{
		const PWEdge &Edge = Path.GetEdge(uEdgeIndex);
#if	TRACE
		Log("\nEdge %u %c%u.%u\n",
		  uEdgeIndex,
		  Edge.cType,
		  Edge.uPrefixLengthA,
		  Edge.uPrefixLengthB);
#endif
		const char cType = Edge.cType;
		const unsigned uPrefixLengthA = Edge.uPrefixLengthA;
		unsigned uColCountA = 0;
		if (uPrefixLengthA > 0)
			{
			const unsigned uNodeIndexA = uPrefixLengthA - 1;
			const unsigned uTplColIndexA = uNodeIndexA;
			if (uTplColIndexA > uColIndexA)
				uColCountA = uTplColIndexA - uColIndexA;
			}

		const unsigned uPrefixLengthB = Edge.uPrefixLengthB;
		unsigned uColCountB = 0;
		if (uPrefixLengthB > 0)
			{
			const unsigned uNodeIndexB = uPrefixLengthB - 1;
			const unsigned uTplColIndexB = uNodeIndexB;
			if (uTplColIndexB > uColIndexB)
				uColCountB = uTplColIndexB - uColIndexB;
			}

		AppendUnalignedTerminals(msaA, uColIndexA, uColCountA, msaB, uColIndexB, uColCountB,
		  uSeqCountA, uSeqCountB, msaCombined, uColIndexCombined);

		switch (cType)
			{
		case 'M':
			{
			assert(uPrefixLengthA > 0);
			assert(uPrefixLengthB > 0);
			const unsigned uColA = uPrefixLengthA - 1;
			const unsigned uColB = uPrefixLengthB - 1;
			assert(uColIndexA == uColA);
			assert(uColIndexB == uColB);
			AppendMatch(msaA, uColIndexA, msaB, uColIndexB, uSeqCountA, uSeqCountB,
			  msaCombined, uColIndexCombined);
			break;
			}
		case 'D':
			{
			assert(uPrefixLengthA > 0);
			const unsigned uColA = uPrefixLengthA - 1;
			assert(uColIndexA == uColA);
			AppendDelete(msaA, uColIndexA, uSeqCountA, uSeqCountB, msaCombined, uColIndexCombined);
			break;
			}
		case 'I':
			{
			assert(uPrefixLengthB > 0);
			const unsigned uColB = uPrefixLengthB - 1;
			assert(uColIndexB == uColB);
			AppendInsert(msaB, uColIndexB, uSeqCountA, uSeqCountB, msaCombined, uColIndexCombined);
//.........这里部分代码省略.........

SCORE SW(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
  unsigned uLengthB, PWPath &Path)
	{
	assert(uLengthB > 0 && uLengthA > 0);

	const unsigned uPrefixCountA = uLengthA + 1;
	const unsigned uPrefixCountB = uLengthB + 1;

// Allocate DP matrices
	const size_t LM = uPrefixCountA*uPrefixCountB;
	SCORE *DPM_ = new SCORE[LM];
	SCORE *DPD_ = new SCORE[LM];
	SCORE *DPI_ = new SCORE[LM];

	DPM(0, 0) = 0;
	DPD(0, 0) = MINUS_INFINITY;
	DPI(0, 0) = MINUS_INFINITY;

	DPM(1, 0) = MINUS_INFINITY;
	DPD(1, 0) = MINUS_INFINITY;
	DPI(1, 0) = MINUS_INFINITY;

	DPM(0, 1) = MINUS_INFINITY;
	DPD(0, 1) = MINUS_INFINITY;
	DPI(0, 1) = MINUS_INFINITY;

// Empty prefix of B is special case
	for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
		{
	// M=LetterA+LetterB, impossible with empty prefix
		DPM(uPrefixLengthA, 0) = MINUS_INFINITY;

	// D=LetterA+GapB, never optimal in local alignment with gap penalties
		DPD(uPrefixLengthA, 0) = MINUS_INFINITY;

	// I=GapA+LetterB, impossible with empty prefix
		DPI(uPrefixLengthA, 0) = MINUS_INFINITY;
		}

// Empty prefix of A is special case
	for (unsigned uPrefixLengthB = 2; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
		{
	// M=LetterA+LetterB, impossible with empty prefix
		DPM(0, uPrefixLengthB) = MINUS_INFINITY;

	// D=LetterA+GapB, impossible with empty prefix
		DPD(0, uPrefixLengthB) = MINUS_INFINITY;

	// I=GapA+LetterB, never optimal in local alignment with gap penalties
		DPI(0, uPrefixLengthB) = MINUS_INFINITY;
		}

	SCORE scoreMax = MINUS_INFINITY;
	unsigned uPrefixLengthAMax = uInsane;
	unsigned uPrefixLengthBMax = uInsane;

// ============
// Main DP loop
// ============
	SCORE scoreGapCloseB = MINUS_INFINITY;
	for (unsigned uPrefixLengthB = 1; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
		{
		const ProfPos &PPB = PB[uPrefixLengthB - 1];

		SCORE scoreGapCloseA = MINUS_INFINITY;
		for (unsigned uPrefixLengthA = 1; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
			{
			const ProfPos &PPA = PA[uPrefixLengthA - 1];

			{
		// Match M=LetterA+LetterB
			SCORE scoreLL = ScoreProfPos2(PPA, PPB);

			SCORE scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1);
			SCORE scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseA;
			SCORE scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseB;

			SCORE scoreBest;
			if (scoreMM >= scoreDM && scoreMM >= scoreIM)
				scoreBest = scoreMM;
			else if (scoreDM >= scoreMM && scoreDM >= scoreIM)
				scoreBest = scoreDM;
			else 
				{
				assert(scoreIM >= scoreMM && scoreIM >= scoreDM);
				scoreBest = scoreIM;
				}
			if (scoreBest < 0)
				scoreBest = 0;
			scoreBest += scoreLL;
			if (scoreBest > scoreMax)
				{
				scoreMax = scoreBest;
				uPrefixLengthAMax = uPrefixLengthA;
				uPrefixLengthBMax = uPrefixLengthB;
				}
			DPM(uPrefixLengthA, uPrefixLengthB) = scoreBest;
			}

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

static SCORE TraceBackDimer(  const SCORE *DPM_, const SCORE *DPD_, const SCORE *DPI_,
                              const char *TBM_, const char *TBD_, const char *TBI_,
                              unsigned uLengthA, unsigned uLengthB, PWPath &Path)
{
    const unsigned uPrefixCountA = uLengthA + 1;

    unsigned uPrefixLengthA = uLengthA;
    unsigned uPrefixLengthB = uLengthB;

    char cEdge = 'M';
    SCORE scoreMax = DPM(uLengthA, uLengthB);
    if (DPD(uLengthA, uLengthB) > scoreMax)
    {
        scoreMax = DPD(uLengthA, uLengthB);
        cEdge = 'D';
    }
    if (DPI(uLengthA, uLengthB) > scoreMax)
    {
        scoreMax = DPI(uLengthA, uLengthB);
        cEdge = 'I';
    }

    for (;;)
    {
        if (0 == uPrefixLengthA && 0 == uPrefixLengthB)
            break;

        PWEdge Edge;
        Edge.cType = cEdge;
        Edge.uPrefixLengthA = uPrefixLengthA;
        Edge.uPrefixLengthB = uPrefixLengthB;
        Path.PrependEdge(Edge);

#if TRACE
        Log("PLA=%u PLB=%u Edge=%c\n", uPrefixLengthA, uPrefixLengthB, cEdge);
#endif
        switch (cEdge)
        {
        case 'M':
            assert(uPrefixLengthA > 0 && uPrefixLengthB > 0);
            cEdge = TBM(uPrefixLengthA, uPrefixLengthB);
            --uPrefixLengthA;
            --uPrefixLengthB;
            break;
        case 'D':
            assert(uPrefixLengthA > 0);
            cEdge = TBD(uPrefixLengthA, uPrefixLengthB);
            --uPrefixLengthA;
            break;
        case 'I':
            assert(uPrefixLengthB > 0);
            cEdge = TBI(uPrefixLengthA, uPrefixLengthB);
            --uPrefixLengthB;
            break;
        default:
            Quit("Invalid edge PLA=%u PLB=%u %c", uPrefixLengthA, uPrefixLengthB, cEdge);
        }
    }
#if	TRACE
    Path.LogMe();
#endif
    return scoreMax;
}

SCORE GlobalAlignDiags(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
  unsigned uLengthB, PWPath &Path)
	{
#if	LIST_DIAGS
	TICKS t1 = GetClockTicks();
#endif

	DiagList DL;

	if (ALPHA_Amino == g_Alpha)
		FindDiags(PA, uLengthA, PB, uLengthB, DL);
	else if (ALPHA_DNA == g_Alpha || ALPHA_RNA == g_Alpha)
		FindDiagsNuc(PA, uLengthA, PB, uLengthB, DL);
	else
		Quit("GlobalAlignDiags: bad alpha");

#if	TRACE
	Log("GlobalAlignDiags, diag list:\n");
	DL.LogMe();
#endif

	DL.Sort();
	DL.DeleteIncompatible();

#if	TRACE
	Log("After DeleteIncompatible:\n");
	DL.LogMe();
#endif

	MergeDiags(DL);

#if	TRACE
	Log("After MergeDiags:\n");
	DL.LogMe();
#endif

	DPRegionList RL;
	DiagListToDPRegionList(DL, RL, uLengthA, uLengthB);

#if	TRACE
	Log("RegionList:\n");
	RL.LogMe();
#endif

#if	LIST_DIAGS
	{
	TICKS t2 = GetClockTicks();
	unsigned uArea = RL.GetDPArea();
	Log("ticks=%ld\n", (long) (t2 - t1));
	Log("area=%u\n", uArea);
	}
#endif

	g_dDPAreaWithoutDiags += uLengthA*uLengthB;

	double dDPAreaWithDiags = 0.0;
	const unsigned uRegionCount = RL.GetCount();
	for (unsigned uRegionIndex = 0; uRegionIndex < uRegionCount; ++uRegionIndex)
		{
		const DPRegion &r = RL.Get(uRegionIndex);

		PWPath RegPath;
		if (DPREGIONTYPE_Diag == r.m_Type)
			{
			DiagToPath(r.m_Diag, RegPath);
#if	TRACE_PATH
			Log("DiagToPath, path=\n");
			RegPath.LogMe();
#endif
			}
		else if (DPREGIONTYPE_Rect == r.m_Type)
			{
			const unsigned uRegStartPosA = r.m_Rect.m_uStartPosA;
			const unsigned uRegStartPosB = r.m_Rect.m_uStartPosB;
			const unsigned uRegLengthA = r.m_Rect.m_uLengthA;
			const unsigned uRegLengthB = r.m_Rect.m_uLengthB;
			const ProfPos *RegPA = PA + uRegStartPosA;
			const ProfPos *RegPB = PB + uRegStartPosB;

			dDPAreaWithDiags += uRegLengthA*uRegLengthB;
			GlobalAlignNoDiags(RegPA, uRegLengthA, RegPB, uRegLengthB, RegPath);
#if	TRACE_PATH
			Log("GlobalAlignNoDiags RegPath=\n");
			RegPath.LogMe();
#endif
			OffsetPath(RegPath, uRegStartPosA, uRegStartPosB);
#if	TRACE_PATH
			Log("After offset path, RegPath=\n");
			RegPath.LogMe();
#endif
			}
		else
			Quit("GlobalAlignDiags, Invalid region type %u", r.m_Type);

		AppendRegPath(Path, RegPath);
#if	TRACE_PATH
		Log("After AppendPath, path=");
		Path.LogMe();
#endif
		}
//.........这里部分代码省略.........