source: trunk/GDE/MUSCLE/src/realigndiffse.cpp

#include "muscle.h"
#include "msa.h"
#include "tree.h"
#include "profile.h"
#include "pwpath.h"
#include "seqvect.h"
#include "estring.h"

#define TRACE           0

void DeleteProgNode(ProgNode &Node)
        {
        delete[] Node.m_Prof;
        delete[] Node.m_EstringL;
        delete[] Node.m_EstringR;

        Node.m_Prof = 0;
        Node.m_EstringL = 0;
        Node.m_EstringR = 0;
        }

static void MakeNode(ProgNode &OldNode, ProgNode &NewNode, bool bSwapLR)
        {
        if (bSwapLR)
                {
                NewNode.m_EstringL = OldNode.m_EstringR;
                NewNode.m_EstringR = OldNode.m_EstringL;
                }
        else
                {
                NewNode.m_EstringL = OldNode.m_EstringL;
                NewNode.m_EstringR = OldNode.m_EstringR;
                }
        NewNode.m_Prof = OldNode.m_Prof;
        NewNode.m_uLength = OldNode.m_uLength;
        NewNode.m_Weight = OldNode.m_Weight;

        OldNode.m_Prof = 0;
        OldNode.m_EstringL = 0;
        OldNode.m_EstringR = 0;
        }

void RealignDiffsE(const MSA &msaIn, const SeqVect &v,
  const Tree &NewTree, const Tree &OldTree, 
  const unsigned uNewNodeIndexToOldNodeIndex[],
  MSA &msaOut, ProgNode *OldProgNodes)
        {
        assert(OldProgNodes != 0);

        const unsigned uNodeCount = NewTree.GetNodeCount();
        if (uNodeCount%2 == 0)
                Quit("RealignDiffs: Expected odd number of nodes");

        const unsigned uMergeCount = (uNodeCount - 1)/2;
        ProgNode *NewProgNodes = new ProgNode[uNodeCount];

        for (unsigned uNewNodeIndex = 0; uNewNodeIndex < uNodeCount; ++uNewNodeIndex)
                {
                if (NODE_CHANGED == uNewNodeIndexToOldNodeIndex[uNewNodeIndex])
                        continue;

                unsigned uOldNodeIndex = uNewNodeIndexToOldNodeIndex[uNewNodeIndex];
                assert(uNewNodeIndex < uNodeCount);
                assert(uOldNodeIndex < uNodeCount);

                ProgNode &NewNode = NewProgNodes[uNewNodeIndex];
                ProgNode &OldNode = OldProgNodes[uOldNodeIndex];
                bool bSwapLR = false;
                if (!NewTree.IsLeaf(uNewNodeIndex))
                        {
                        unsigned uNewLeft = NewTree.GetLeft(uNewNodeIndex);
                        unsigned uNewRight = NewTree.GetRight(uNewNodeIndex);
                        unsigned uOld = uNewNodeIndexToOldNodeIndex[uNewNodeIndex];
                        unsigned uOldLeft = OldTree.GetLeft(uOld);
                        unsigned uOldRight = OldTree.GetRight(uOld);
                        assert(uOldLeft < uNodeCount && uOldRight < uNodeCount);
                        if (uOldLeft != uNewNodeIndexToOldNodeIndex[uNewLeft])
                                {
                                assert(uOldLeft == uNewNodeIndexToOldNodeIndex[uNewRight]);
                                bSwapLR = true;
                                }
                        }
                MakeNode(OldNode, NewNode, bSwapLR);
#if     TRACE
                Log("MakeNode old=%u new=%u swap=%d length=%u weight=%.3g\n",
                  uOldNodeIndex, uNewNodeIndex, bSwapLR, NewNode.m_uLength, NewNode.m_Weight);
#endif
                }

        unsigned uJoin = 0;
        SetProgressDesc("Refine tree");
        for (unsigned uNewNodeIndex = NewTree.FirstDepthFirstNode();
          NULL_NEIGHBOR != uNewNodeIndex;
          uNewNodeIndex = NewTree.NextDepthFirstNode(uNewNodeIndex))
                {
                if (NODE_CHANGED != uNewNodeIndexToOldNodeIndex[uNewNodeIndex])
                        continue;

                Progress(uJoin, uMergeCount - 1);
                ++uJoin;

                const unsigned uMergeNodeIndex = uNewNodeIndex;
                ProgNode &Parent = NewProgNodes[uMergeNodeIndex];

                const unsigned uLeft = NewTree.GetLeft(uNewNodeIndex);
                const unsigned uRight = NewTree.GetRight(uNewNodeIndex);

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

                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);
                PathToEstrings(Parent.m_Path, &Parent.m_EstringL, &Parent.m_EstringR);

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

                delete[] Node1.m_Prof;
                delete[] Node2.m_Prof;

                Node1.m_Prof = 0;
                Node2.m_Prof = 0;
                }

        ProgressStepsDone();

        if (g_bBrenner)
                MakeRootMSABrenner((SeqVect &) v, NewTree, NewProgNodes, msaOut);
        else
                MakeRootMSA(v, NewTree, NewProgNodes, msaOut);

#if     DEBUG
        AssertMSAEqIgnoreCaseAndGaps(msaIn, msaOut);
#endif

        for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
                DeleteProgNode(NewProgNodes[uNodeIndex]);

        delete[] NewProgNodes;
        }