+++ /dev/null
-#include "muscle.h"\r
-#include "tree.h"\r
-#include "msa.h"\r
-#include "pwpath.h"\r
-#include "profile.h"\r
-#include "scorehistory.h"\r
-#include "objscore.h"\r
-\r
-unsigned g_uRefineHeightSubtree;\r
-unsigned g_uRefineHeightSubtreeTotal;\r
-\r
-#define TRACE 0\r
-#define DIFFOBJSCORE 0\r
-\r
-static bool TryRealign(MSA &msaIn, const Tree &tree, const unsigned Leaves1[],\r
- unsigned uCount1, const unsigned Leaves2[], unsigned uCount2,\r
- SCORE *ptrscoreBefore, SCORE *ptrscoreAfter,\r
- bool bLockLeft, bool bLockRight)\r
- {\r
-#if TRACE\r
- Log("TryRealign, msaIn=\n");\r
- msaIn.LogMe();\r
-#endif\r
-\r
- const unsigned uSeqCount = msaIn.GetSeqCount();\r
-\r
- unsigned *Ids1 = new unsigned[uSeqCount];\r
- unsigned *Ids2 = new unsigned[uSeqCount];\r
-\r
- LeafIndexesToIds(tree, Leaves1, uCount1, Ids1);\r
- LeafIndexesToIds(tree, Leaves2, uCount2, Ids2);\r
-\r
- MSA msa1;\r
- MSA msa2;\r
-\r
- MSASubsetByIds(msaIn, Ids1, uCount1, msa1);\r
- MSASubsetByIds(msaIn, Ids2, uCount2, msa2);\r
-\r
-#if DEBUG\r
- ValidateMuscleIds(msa1);\r
- ValidateMuscleIds(msa2);\r
-#endif\r
-\r
-// Computing the objective score may be expensive for\r
-// large numbers of sequences. As a speed optimization,\r
-// we check whether the alignment changes. If it does\r
-// not change, there is no need to compute the objective\r
-// score. We test for the alignment changing by comparing\r
-// the Viterbi paths before and after re-aligning.\r
- PWPath pathBefore;\r
- pathBefore.FromMSAPair(msa1, msa2);\r
-\r
- DeleteGappedCols(msa1);\r
- DeleteGappedCols(msa2);\r
-\r
- if (0 == msa1.GetColCount() || 0 == msa2.GetColCount())\r
- return false;\r
-\r
- MSA msaRealigned;\r
- PWPath pathAfter;\r
-\r
- AlignTwoMSAs(msa1, msa2, msaRealigned, pathAfter, bLockLeft, bLockRight);\r
-\r
- bool bAnyChanges = !pathAfter.Equal(pathBefore);\r
- unsigned uDiffCount1;\r
- unsigned uDiffCount2;\r
- static unsigned Edges1[10000];\r
- static unsigned Edges2[10000];\r
- DiffPaths(pathBefore, pathAfter, Edges1, &uDiffCount1, Edges2, &uDiffCount2);\r
-\r
-#if TRACE\r
- Log("TryRealign, msa1=\n");\r
- msa1.LogMe();\r
- Log("\nmsa2=\n");\r
- msa2.LogMe();\r
- Log("\nRealigned (changes %s)=\n", bAnyChanges ? "TRUE" : "FALSE");\r
- msaRealigned.LogMe();\r
-#endif\r
-\r
- if (!bAnyChanges)\r
- {\r
- *ptrscoreBefore = 0;\r
- *ptrscoreAfter = 0;\r
- return false;\r
- }\r
-\r
- SetMSAWeightsMuscle(msaIn);\r
- SetMSAWeightsMuscle(msaRealigned);\r
-\r
-#if DIFFOBJSCORE\r
- const SCORE scoreDiff = DiffObjScore(msaIn, pathBefore, Edges1, uDiffCount1,\r
- msaRealigned, pathAfter, Edges2, uDiffCount2);\r
- bool bAccept = (scoreDiff > 0);\r
- *ptrscoreBefore = 0;\r
- *ptrscoreAfter = scoreDiff;\r
- //const SCORE scoreBefore = ObjScoreIds(msaIn, Ids1, uCount1, Ids2, uCount2);\r
- //const SCORE scoreAfter = ObjScoreIds(msaRealigned, Ids1, uCount1, Ids2, uCount2);\r
- //Log("Diff = %.3g %.3g\n", scoreDiff, scoreAfter - scoreBefore);\r
-#else\r
- const SCORE scoreBefore = ObjScoreIds(msaIn, Ids1, uCount1, Ids2, uCount2);\r
- const SCORE scoreAfter = ObjScoreIds(msaRealigned, Ids1, uCount1, Ids2, uCount2);\r
-\r
- bool bAccept = (scoreAfter > scoreBefore);\r
-\r
-#if TRACE\r
- Log("Score %g -> %g Accept %s\n", scoreBefore, scoreAfter, bAccept ? "TRUE" : "FALSE");\r
-#endif\r
-\r
- *ptrscoreBefore = scoreBefore;\r
- *ptrscoreAfter = scoreAfter;\r
-#endif\r
-\r
- if (bAccept)\r
- msaIn.Copy(msaRealigned);\r
- delete[] Ids1;\r
- delete[] Ids2;\r
- return bAccept;\r
- }\r
-\r
-static void RefineHeightParts(MSA &msaIn, const Tree &tree,\r
- const unsigned InternalNodeIndexes[], bool bReversed, bool bRight,\r
- unsigned uIter, \r
- ScoreHistory &History,\r
- bool *ptrbAnyChanges, bool *ptrbOscillating, bool bLockLeft, bool bLockRight)\r
- {\r
- *ptrbOscillating = false;\r
-\r
- const unsigned uSeqCount = msaIn.GetSeqCount();\r
- const unsigned uInternalNodeCount = uSeqCount - 1;\r
-\r
- unsigned *Leaves1 = new unsigned[uSeqCount];\r
- unsigned *Leaves2 = new unsigned[uSeqCount];\r
-\r
- const unsigned uRootNodeIndex = tree.GetRootNodeIndex();\r
- bool bAnyAccepted = false;\r
- for (unsigned i = 0; i < uInternalNodeCount; ++i)\r
- {\r
- const unsigned uInternalNodeIndex = InternalNodeIndexes[i];\r
- unsigned uNeighborNodeIndex;\r
- if (tree.IsRoot(uInternalNodeIndex) && !bRight)\r
- continue;\r
- else if (bRight)\r
- uNeighborNodeIndex = tree.GetRight(uInternalNodeIndex);\r
- else\r
- uNeighborNodeIndex = tree.GetLeft(uInternalNodeIndex);\r
-\r
- g_uTreeSplitNode1 = uInternalNodeIndex;\r
- g_uTreeSplitNode2 = uNeighborNodeIndex;\r
-\r
- unsigned uCount1;\r
- unsigned uCount2;\r
-\r
- GetLeaves(tree, uNeighborNodeIndex, Leaves1, &uCount1);\r
- GetLeavesExcluding(tree, uRootNodeIndex, uNeighborNodeIndex,\r
- Leaves2, &uCount2);\r
-\r
-#if TRACE\r
- Log("\nRefineHeightParts node %u\n", uInternalNodeIndex);\r
- Log("Group1=");\r
- for (unsigned n = 0; n < uCount1; ++n)\r
- Log(" %u(%s)", Leaves1[n], tree.GetName(Leaves1[n]));\r
- Log("\n");\r
- Log("Group2=");\r
- for (unsigned n = 0; n < uCount2; ++n)\r
- Log(" %u(%s)", Leaves2[n], tree.GetName(Leaves2[n]));\r
- Log("\n");\r
-#endif\r
-\r
- SCORE scoreBefore;\r
- SCORE scoreAfter;\r
- bool bAccepted = TryRealign(msaIn, tree, Leaves1, uCount1, Leaves2, uCount2,\r
- &scoreBefore, &scoreAfter, bLockLeft, bLockRight);\r
- SetCurrentAlignment(msaIn);\r
-\r
- ++g_uRefineHeightSubtree;\r
- Progress(g_uRefineHeightSubtree, g_uRefineHeightSubtreeTotal);\r
-\r
-#if TRACE\r
- if (uIter > 0)\r
- Log("Before %g %g\n", scoreBefore,\r
- History.GetScore(uIter - 1, uInternalNodeIndex, bReversed, bRight));\r
-#endif\r
- SCORE scoreMax = scoreAfter > scoreBefore? scoreAfter : scoreBefore;\r
- bool bRepeated = History.SetScore(uIter, uInternalNodeIndex, bRight, scoreMax);\r
- if (bRepeated)\r
- {\r
- *ptrbOscillating = true;\r
- break;\r
- }\r
-\r
- if (bAccepted)\r
- bAnyAccepted = true;\r
- }\r
-\r
- delete[] Leaves1;\r
- delete[] Leaves2;\r
-\r
- *ptrbAnyChanges = bAnyAccepted;\r
- }\r
-\r
-// Return true if any changes made\r
-bool RefineHoriz(MSA &msaIn, const Tree &tree, unsigned uIters, bool bLockLeft,\r
- bool bLockRight)\r
- {\r
-#if TRACE\r
- tree.LogMe();\r
-#endif\r
-\r
- if (!tree.IsRooted())\r
- Quit("RefineHeight: requires rooted tree");\r
-\r
- const unsigned uSeqCount = msaIn.GetSeqCount();\r
- if (uSeqCount < 3)\r
- return false;\r
-\r
- const unsigned uInternalNodeCount = uSeqCount - 1;\r
- unsigned *InternalNodeIndexes = new unsigned[uInternalNodeCount];\r
- unsigned *InternalNodeIndexesR = new unsigned[uInternalNodeCount];\r
-\r
- GetInternalNodesInHeightOrder(tree, InternalNodeIndexes);\r
-\r
- ScoreHistory History(uIters, 2*uSeqCount - 1);\r
-\r
- bool bAnyChangesAnyIter = false;\r
- for (unsigned n = 0; n < uInternalNodeCount; ++n)\r
- InternalNodeIndexesR[uInternalNodeCount - 1 - n] = InternalNodeIndexes[n];\r
-\r
- for (unsigned uIter = 0; uIter < uIters; ++uIter)\r
- {\r
- bool bAnyChangesThisIter = false;\r
- IncIter();\r
- SetProgressDesc("Refine biparts");\r
- g_uRefineHeightSubtree = 0;\r
- g_uRefineHeightSubtreeTotal = uInternalNodeCount*2 - 1;\r
-\r
- bool bReverse = (uIter%2 != 0);\r
- unsigned *Internals;\r
- if (bReverse)\r
- Internals = InternalNodeIndexesR;\r
- else\r
- Internals = InternalNodeIndexes;\r
-\r
- bool bOscillating;\r
- for (unsigned i = 0; i < 2; ++i)\r
- {\r
- bool bAnyChanges = false;\r
- bool bRight;\r
- switch (i)\r
- {\r
- case 0:\r
- bRight = true;\r
- break;\r
- case 1:\r
- bRight = false;\r
- break;\r
- default:\r
- Quit("RefineHeight default case");\r
- }\r
- RefineHeightParts(msaIn, tree, Internals, bReverse, bRight,\r
- uIter, \r
- History, \r
- &bAnyChanges, &bOscillating, bLockLeft, bLockRight);\r
- if (bOscillating)\r
- {\r
- ProgressStepsDone();\r
- goto Osc;\r
- }\r
- if (bAnyChanges)\r
- {\r
- bAnyChangesThisIter = true;\r
- bAnyChangesAnyIter = true;\r
- }\r
- }\r
-\r
- ProgressStepsDone();\r
- if (bOscillating)\r
- break;\r
-\r
- if (!bAnyChangesThisIter)\r
- break;\r
- }\r
-\r
-Osc:\r
- delete[] InternalNodeIndexes;\r
- delete[] InternalNodeIndexesR;\r
-\r
- return bAnyChangesAnyIter;\r
- }\r
git://source.jalview.org / jabaws.git / blobdiff