6 #include "scorehistory.h"
\r
7 #include "objscore.h"
\r
9 unsigned g_uRefineHeightSubtree;
\r
10 unsigned g_uRefineHeightSubtreeTotal;
\r
13 #define DIFFOBJSCORE 0
\r
15 static bool TryRealign(MSA &msaIn, const Tree &tree, const unsigned Leaves1[],
\r
16 unsigned uCount1, const unsigned Leaves2[], unsigned uCount2,
\r
17 SCORE *ptrscoreBefore, SCORE *ptrscoreAfter,
\r
18 bool bLockLeft, bool bLockRight)
\r
21 Log("TryRealign, msaIn=\n");
\r
25 const unsigned uSeqCount = msaIn.GetSeqCount();
\r
27 unsigned *Ids1 = new unsigned[uSeqCount];
\r
28 unsigned *Ids2 = new unsigned[uSeqCount];
\r
30 LeafIndexesToIds(tree, Leaves1, uCount1, Ids1);
\r
31 LeafIndexesToIds(tree, Leaves2, uCount2, Ids2);
\r
36 MSASubsetByIds(msaIn, Ids1, uCount1, msa1);
\r
37 MSASubsetByIds(msaIn, Ids2, uCount2, msa2);
\r
40 ValidateMuscleIds(msa1);
\r
41 ValidateMuscleIds(msa2);
\r
44 // Computing the objective score may be expensive for
\r
45 // large numbers of sequences. As a speed optimization,
\r
46 // we check whether the alignment changes. If it does
\r
47 // not change, there is no need to compute the objective
\r
48 // score. We test for the alignment changing by comparing
\r
49 // the Viterbi paths before and after re-aligning.
\r
51 pathBefore.FromMSAPair(msa1, msa2);
\r
53 DeleteGappedCols(msa1);
\r
54 DeleteGappedCols(msa2);
\r
56 if (0 == msa1.GetColCount() || 0 == msa2.GetColCount())
\r
62 AlignTwoMSAs(msa1, msa2, msaRealigned, pathAfter, bLockLeft, bLockRight);
\r
64 bool bAnyChanges = !pathAfter.Equal(pathBefore);
\r
65 unsigned uDiffCount1;
\r
66 unsigned uDiffCount2;
\r
67 static unsigned Edges1[10000];
\r
68 static unsigned Edges2[10000];
\r
69 DiffPaths(pathBefore, pathAfter, Edges1, &uDiffCount1, Edges2, &uDiffCount2);
\r
72 Log("TryRealign, msa1=\n");
\r
76 Log("\nRealigned (changes %s)=\n", bAnyChanges ? "TRUE" : "FALSE");
\r
77 msaRealigned.LogMe();
\r
82 *ptrscoreBefore = 0;
\r
87 SetMSAWeightsMuscle(msaIn);
\r
88 SetMSAWeightsMuscle(msaRealigned);
\r
91 const SCORE scoreDiff = DiffObjScore(msaIn, pathBefore, Edges1, uDiffCount1,
\r
92 msaRealigned, pathAfter, Edges2, uDiffCount2);
\r
93 bool bAccept = (scoreDiff > 0);
\r
94 *ptrscoreBefore = 0;
\r
95 *ptrscoreAfter = scoreDiff;
\r
96 //const SCORE scoreBefore = ObjScoreIds(msaIn, Ids1, uCount1, Ids2, uCount2);
\r
97 //const SCORE scoreAfter = ObjScoreIds(msaRealigned, Ids1, uCount1, Ids2, uCount2);
\r
98 //Log("Diff = %.3g %.3g\n", scoreDiff, scoreAfter - scoreBefore);
\r
100 const SCORE scoreBefore = ObjScoreIds(msaIn, Ids1, uCount1, Ids2, uCount2);
\r
101 const SCORE scoreAfter = ObjScoreIds(msaRealigned, Ids1, uCount1, Ids2, uCount2);
\r
103 bool bAccept = (scoreAfter > scoreBefore);
\r
106 Log("Score %g -> %g Accept %s\n", scoreBefore, scoreAfter, bAccept ? "TRUE" : "FALSE");
\r
109 *ptrscoreBefore = scoreBefore;
\r
110 *ptrscoreAfter = scoreAfter;
\r
114 msaIn.Copy(msaRealigned);
\r
120 static void RefineHeightParts(MSA &msaIn, const Tree &tree,
\r
121 const unsigned InternalNodeIndexes[], bool bReversed, bool bRight,
\r
123 ScoreHistory &History,
\r
124 bool *ptrbAnyChanges, bool *ptrbOscillating, bool bLockLeft, bool bLockRight)
\r
126 *ptrbOscillating = false;
\r
128 const unsigned uSeqCount = msaIn.GetSeqCount();
\r
129 const unsigned uInternalNodeCount = uSeqCount - 1;
\r
131 unsigned *Leaves1 = new unsigned[uSeqCount];
\r
132 unsigned *Leaves2 = new unsigned[uSeqCount];
\r
134 const unsigned uRootNodeIndex = tree.GetRootNodeIndex();
\r
135 bool bAnyAccepted = false;
\r
136 for (unsigned i = 0; i < uInternalNodeCount; ++i)
\r
138 const unsigned uInternalNodeIndex = InternalNodeIndexes[i];
\r
139 unsigned uNeighborNodeIndex;
\r
140 if (tree.IsRoot(uInternalNodeIndex) && !bRight)
\r
143 uNeighborNodeIndex = tree.GetRight(uInternalNodeIndex);
\r
145 uNeighborNodeIndex = tree.GetLeft(uInternalNodeIndex);
\r
147 g_uTreeSplitNode1 = uInternalNodeIndex;
\r
148 g_uTreeSplitNode2 = uNeighborNodeIndex;
\r
153 GetLeaves(tree, uNeighborNodeIndex, Leaves1, &uCount1);
\r
154 GetLeavesExcluding(tree, uRootNodeIndex, uNeighborNodeIndex,
\r
155 Leaves2, &uCount2);
\r
158 Log("\nRefineHeightParts node %u\n", uInternalNodeIndex);
\r
160 for (unsigned n = 0; n < uCount1; ++n)
\r
161 Log(" %u(%s)", Leaves1[n], tree.GetName(Leaves1[n]));
\r
164 for (unsigned n = 0; n < uCount2; ++n)
\r
165 Log(" %u(%s)", Leaves2[n], tree.GetName(Leaves2[n]));
\r
171 bool bAccepted = TryRealign(msaIn, tree, Leaves1, uCount1, Leaves2, uCount2,
\r
172 &scoreBefore, &scoreAfter, bLockLeft, bLockRight);
\r
173 SetCurrentAlignment(msaIn);
\r
175 ++g_uRefineHeightSubtree;
\r
176 Progress(g_uRefineHeightSubtree, g_uRefineHeightSubtreeTotal);
\r
180 Log("Before %g %g\n", scoreBefore,
\r
181 History.GetScore(uIter - 1, uInternalNodeIndex, bReversed, bRight));
\r
183 SCORE scoreMax = scoreAfter > scoreBefore? scoreAfter : scoreBefore;
\r
184 bool bRepeated = History.SetScore(uIter, uInternalNodeIndex, bRight, scoreMax);
\r
187 *ptrbOscillating = true;
\r
192 bAnyAccepted = true;
\r
198 *ptrbAnyChanges = bAnyAccepted;
\r
201 // Return true if any changes made
\r
202 bool RefineHoriz(MSA &msaIn, const Tree &tree, unsigned uIters, bool bLockLeft,
\r
209 if (!tree.IsRooted())
\r
210 Quit("RefineHeight: requires rooted tree");
\r
212 const unsigned uSeqCount = msaIn.GetSeqCount();
\r
216 const unsigned uInternalNodeCount = uSeqCount - 1;
\r
217 unsigned *InternalNodeIndexes = new unsigned[uInternalNodeCount];
\r
218 unsigned *InternalNodeIndexesR = new unsigned[uInternalNodeCount];
\r
220 GetInternalNodesInHeightOrder(tree, InternalNodeIndexes);
\r
222 ScoreHistory History(uIters, 2*uSeqCount - 1);
\r
224 bool bAnyChangesAnyIter = false;
\r
225 for (unsigned n = 0; n < uInternalNodeCount; ++n)
\r
226 InternalNodeIndexesR[uInternalNodeCount - 1 - n] = InternalNodeIndexes[n];
\r
228 for (unsigned uIter = 0; uIter < uIters; ++uIter)
\r
230 bool bAnyChangesThisIter = false;
\r
232 SetProgressDesc("Refine biparts");
\r
233 g_uRefineHeightSubtree = 0;
\r
234 g_uRefineHeightSubtreeTotal = uInternalNodeCount*2 - 1;
\r
236 bool bReverse = (uIter%2 != 0);
\r
237 unsigned *Internals;
\r
239 Internals = InternalNodeIndexesR;
\r
241 Internals = InternalNodeIndexes;
\r
244 for (unsigned i = 0; i < 2; ++i)
\r
246 bool bAnyChanges = false;
\r
257 Quit("RefineHeight default case");
\r
259 RefineHeightParts(msaIn, tree, Internals, bReverse, bRight,
\r
262 &bAnyChanges, &bOscillating, bLockLeft, bLockRight);
\r
265 ProgressStepsDone();
\r
270 bAnyChangesThisIter = true;
\r
271 bAnyChangesAnyIter = true;
\r
275 ProgressStepsDone();
\r
279 if (!bAnyChangesThisIter)
\r
284 delete[] InternalNodeIndexes;
\r
285 delete[] InternalNodeIndexesR;
\r
287 return bAnyChangesAnyIter;
\r