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

#include "muscle.h"
#include "diaglist.h"
#include "pwpath.h"

#define MAX(x, y)       ((x) > (y) ? (x) : (y))
#define MIN(x, y)       ((x) < (y) ? (x) : (y))

void DiagList::Add(const Diag &d)
        {
        if (m_uCount == MAX_DIAGS)
                Quit("DiagList::Add, overflow %u", m_uCount);
        m_Diags[m_uCount] = d;
        ++m_uCount;
        }

void DiagList::Add(unsigned uStartPosA, unsigned uStartPosB, unsigned uLength)
        {
        Diag d;
        d.m_uStartPosA = uStartPosA;
        d.m_uStartPosB = uStartPosB;
        d.m_uLength = uLength;
        Add(d);
        }

const Diag &DiagList::Get(unsigned uIndex) const
        {
        if (uIndex >= m_uCount)
                Quit("DiagList::Get(%u), count=%u", uIndex, m_uCount);
        return m_Diags[uIndex];
        }

void DiagList::LogMe() const
        {
        Log("DiagList::LogMe, count=%u\n", m_uCount);
        Log("  n  StartA  StartB  Length\n");
        Log("---  ------  ------  ------\n");
        for (unsigned n = 0; n < m_uCount; ++n)
                {
                const Diag &d = m_Diags[n];
                Log("%3u  %6u  %6u  %6u\n",
                  n, d.m_uStartPosA, d.m_uStartPosB, d.m_uLength);
                }
        }

void DiagList::FromPath(const PWPath &Path)
        {
        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 DiagList::NonZeroIntersection(const Diag &d) const
        {
        for (unsigned n = 0; n < m_uCount; ++n)
                {
                const Diag &d2 = m_Diags[n];
                if (DiagOverlap(d, d2) > 0)
                        return true;
                }
        return false;
        }

// DialogOverlap returns the length of the overlapping
// section of the two diagonals along the diagonals
// themselves; in other words, the length of
// the intersection of the two sets of cells in
// the matrix.
unsigned DiagOverlap(const Diag &d1, const Diag &d2)
        {
// Determine where the diagonals intersect the A
// axis (extending them if required). If they
// intersect at different points, they do not
// overlap. Coordinates on a diagonal are
// given by B = A + c where c is the value of
// A at the intersection with the A axis.
// Hence, c = B - A for any point on the diagonal.
        int c1 = (int) d1.m_uStartPosB - (int) d1.m_uStartPosA;
        int c2 = (int) d2.m_uStartPosB - (int) d2.m_uStartPosA;
        if (c1 != c2)
                return 0;

        assert(DiagOverlapA(d1, d2) == DiagOverlapB(d1, d2));
        return DiagOverlapA(d1, d2);
        }

// DialogOverlapA returns the length of the overlapping
// section of the projection of the two diagonals onto
// the A axis.
unsigned DiagOverlapA(const Diag &d1, const Diag &d2)
        {
        unsigned uMaxStart = MAX(d1.m_uStartPosA, d2.m_uStartPosA);
        unsigned uMinEnd = MIN(d1.m_uStartPosA + d1.m_uLength - 1,
          d2.m_uStartPosA + d2.m_uLength - 1);

        int iLength = (int) uMinEnd - (int) uMaxStart + 1;
        if (iLength < 0)
                return 0;
        return (unsigned) iLength;
        }

// DialogOverlapB returns the length of the overlapping
// section of the projection of the two diagonals onto
// the B axis.
unsigned DiagOverlapB(const Diag &d1, const Diag &d2)
        {
        unsigned uMaxStart = MAX(d1.m_uStartPosB, d2.m_uStartPosB);
        unsigned uMinEnd = MIN(d1.m_uStartPosB + d1.m_uLength - 1,
          d2.m_uStartPosB + d2.m_uLength - 1);

        int iLength = (int) uMinEnd - (int) uMaxStart + 1;
        if (iLength < 0)
                return 0;
        return (unsigned) iLength;
        }

// Returns true if the two diagonals can be on the
// same path through the DP matrix. If DiagCompatible
// returns false, they cannot be in the same path
// and hence "contradict" each other.
bool DiagCompatible(const Diag &d1, const Diag &d2)
        {
        if (DiagOverlap(d1, d2) > 0)
                return true;
        return 0 == DiagOverlapA(d1, d2) && 0 == DiagOverlapB(d1, d2);
        }

// Returns the length of the "break" between two diagonals.
unsigned DiagBreak(const Diag &d1, const Diag &d2)
        {
        int c1 = (int) d1.m_uStartPosB - (int) d1.m_uStartPosA;
        int c2 = (int) d2.m_uStartPosB - (int) d2.m_uStartPosA;
        if (c1 != c2)
                return 0;

        int iMaxStart = MAX(d1.m_uStartPosA, d2.m_uStartPosA);
        int iMinEnd = MIN(d1.m_uStartPosA + d1.m_uLength - 1,
          d2.m_uStartPosA + d1.m_uLength - 1);
        int iBreak = iMaxStart - iMinEnd - 1;
        if (iBreak < 0)
                return 0;
        return (unsigned) iBreak;
        }

// Merge diagonals that are continuations of each other with
// short breaks of up to length g_uMaxDiagBreak.
// In a sorted list of diagonals, we only have to check
// consecutive entries.
void MergeDiags(DiagList &DL)
        {
        return;
#if     DEBUG
        if (!DL.IsSorted())
                Quit("MergeDiags: !IsSorted");
#endif

// TODO: Fix this!
// Breaks must be with no offset (no gaps)
        const unsigned uCount = DL.GetCount();
        if (uCount <= 1)
                return;

        DiagList NewList;

        Diag MergedDiag;
        const Diag *ptrPrev = &DL.Get(0);
        for (unsigned i = 1; i < uCount; ++i)
                {
                const Diag *ptrDiag = &DL.Get(i);
                unsigned uBreakLength = DiagBreak(*ptrPrev, *ptrDiag);
                if (uBreakLength <= g_uMaxDiagBreak)
                        {
                        MergedDiag.m_uStartPosA = ptrPrev->m_uStartPosA;
                        MergedDiag.m_uStartPosB = ptrPrev->m_uStartPosB;
                        MergedDiag.m_uLength = ptrPrev->m_uLength + ptrDiag->m_uLength
                          + uBreakLength;
                        ptrPrev = &MergedDiag;
                        }
                else
                        {
                        NewList.Add(*ptrPrev);
                        ptrPrev = ptrDiag;
                        }
                }
        NewList.Add(*ptrPrev);
        DL.Copy(NewList);
        }

void DiagList::DeleteIncompatible()
        {
        assert(IsSorted());

        if (m_uCount < 2)
                return;

        bool *bFlagForDeletion = new bool[m_uCount];
        for (unsigned i = 0; i < m_uCount; ++i)
                bFlagForDeletion[i] = false;

        for (unsigned i = 0; i < m_uCount; ++i)
                {
                const Diag &di = m_Diags[i];
                for (unsigned j = i + 1; j < m_uCount; ++j)
                        {
                        const Diag &dj = m_Diags[j];

                // Verify sorted correctly
                        assert(di.m_uStartPosA <= dj.m_uStartPosA);

                // If two diagonals are incompatible and
                // one is is much longer than the other,
                // keep the longer one.
                        if (!DiagCompatible(di, dj))
                                {
                                if (di.m_uLength > dj.m_uLength*4)
                                        bFlagForDeletion[j] = true;
                                else if (dj.m_uLength > di.m_uLength*4)
                                        bFlagForDeletion[i] = true;
                                else
                                        {
                                        bFlagForDeletion[i] = true;
                                        bFlagForDeletion[j] = true;
                                        }
                                }
                        }
                }

        for (unsigned i = 0; i < m_uCount; ++i)
                {
                const Diag &di = m_Diags[i];
                if (bFlagForDeletion[i])
                        continue;

                for (unsigned j = i + 1; j < m_uCount; ++j)
                        {
                        const Diag &dj = m_Diags[j];
                        if (bFlagForDeletion[j])
                                continue;

                // Verify sorted correctly
                        assert(di.m_uStartPosA <= dj.m_uStartPosA);

                // If sort order in B different from sorted order in A,
                // either diags are incompatible or we detected a repeat
                // or permutation.
                        if (di.m_uStartPosB >= dj.m_uStartPosB || !DiagCompatible(di, dj))
                                {
                                bFlagForDeletion[i] = true;
                                bFlagForDeletion[j] = true;
                                }
                        }
                }

        unsigned uNewCount = 0;
        Diag *NewDiags = new Diag[m_uCount];
        for (unsigned i = 0; i < m_uCount; ++i)
                {
                if (bFlagForDeletion[i])
                        continue;

                const Diag &d = m_Diags[i];
                NewDiags[uNewCount] = d;
                ++uNewCount;
                }
        memcpy(m_Diags, NewDiags, uNewCount*sizeof(Diag));
        m_uCount = uNewCount;
        delete[] NewDiags;
        }

void DiagList::Copy(const DiagList &DL)
        {
        Clear();
        unsigned uCount = DL.GetCount();
        for (unsigned i = 0; i < uCount; ++i)
                Add(DL.Get(i));
        }

// Check if sorted in increasing order of m_uStartPosA
bool DiagList::IsSorted() const
        {
        return true;
        unsigned uCount = GetCount();
        for (unsigned i = 1; i < uCount; ++i)
                if (m_Diags[i-1].m_uStartPosA > m_Diags[i].m_uStartPosA)
                        return false;
        return true;
        }

// Sort in increasing order of m_uStartPosA
// Dumb bubble sort, but don't care about speed
// because don't get long lists.
void DiagList::Sort()
        {
        if (m_uCount < 2)
                return;

        bool bContinue = true;
        while (bContinue)
                {
                bContinue = false;
                for (unsigned i = 0; i < m_uCount - 1; ++i)
                        {
                        if (m_Diags[i].m_uStartPosA > m_Diags[i+1].m_uStartPosA)
                                {
                                Diag Tmp = m_Diags[i];
                                m_Diags[i] = m_Diags[i+1];
                                m_Diags[i+1] = Tmp;
                                bContinue = true;
                                }
                        }
                }
        }

//void TestDiag()
//      {
//      Diag d1;
//      Diag d2;
//      Diag d3;
//
//      d1.m_uStartPosA = 0;
//      d1.m_uStartPosB = 1;
//      d1.m_uLength = 32;
//
//      d2.m_uStartPosA = 55;
//      d2.m_uStartPosB = 70;
//      d2.m_uLength = 36;
//
//      d3.m_uStartPosA = 102;
//      d3.m_uStartPosB = 122;
//      d3.m_uLength = 50;
//
//      DiagList DL;
//      DL.Add(d1);
//      DL.Add(d2);
//      DL.Add(d3);
//
//      Log("Before DeleteIncompatible:\n");
//      DL.LogMe();
//      DL.DeleteIncompatible();
//
//      Log("After DeleteIncompatible:\n");
//      DL.LogMe();
//
//      MergeDiags(DL);
//      Log("After Merge:\n");
//      DL.LogMe();
//
//      DPRegionList RL;
//      DiagListToDPRegionList(DL, RL, 200, 200);
//      RL.LogMe();
//      }