+++ /dev/null
-#include "muscle.h"\r
-#include "diaglist.h"\r
-#include "pwpath.h"\r
-\r
-#define MAX(x, y) ((x) > (y) ? (x) : (y))\r
-#define MIN(x, y) ((x) < (y) ? (x) : (y))\r
-\r
-void DiagList::Add(const Diag &d)\r
- {\r
- if (m_uCount == MAX_DIAGS)\r
- Quit("DiagList::Add, overflow %u", m_uCount);\r
- m_Diags[m_uCount] = d;\r
- ++m_uCount;\r
- }\r
-\r
-void DiagList::Add(unsigned uStartPosA, unsigned uStartPosB, unsigned uLength)\r
- {\r
- Diag d;\r
- d.m_uStartPosA = uStartPosA;\r
- d.m_uStartPosB = uStartPosB;\r
- d.m_uLength = uLength;\r
- Add(d);\r
- }\r
-\r
-const Diag &DiagList::Get(unsigned uIndex) const\r
- {\r
- if (uIndex >= m_uCount)\r
- Quit("DiagList::Get(%u), count=%u", uIndex, m_uCount);\r
- return m_Diags[uIndex];\r
- }\r
-\r
-void DiagList::LogMe() const\r
- {\r
- Log("DiagList::LogMe, count=%u\n", m_uCount);\r
- Log(" n StartA StartB Length\n");\r
- Log("--- ------ ------ ------\n");\r
- for (unsigned n = 0; n < m_uCount; ++n)\r
- {\r
- const Diag &d = m_Diags[n];\r
- Log("%3u %6u %6u %6u\n",\r
- n, d.m_uStartPosA, d.m_uStartPosB, d.m_uLength);\r
- }\r
- }\r
-\r
-void DiagList::FromPath(const PWPath &Path)\r
- {\r
- Clear();\r
-\r
- const unsigned uEdgeCount = Path.GetEdgeCount();\r
- unsigned uLength = 0;\r
- unsigned uStartPosA;\r
- unsigned uStartPosB;\r
- for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
- {\r
- const PWEdge &Edge = Path.GetEdge(uEdgeIndex);\r
-\r
- // Typical cases\r
- if (Edge.cType == 'M')\r
- {\r
- if (0 == uLength)\r
- {\r
- uStartPosA = Edge.uPrefixLengthA - 1;\r
- uStartPosB = Edge.uPrefixLengthB - 1;\r
- }\r
- ++uLength;\r
- }\r
- else\r
- {\r
- if (uLength >= g_uMinDiagLength)\r
- Add(uStartPosA, uStartPosB, uLength);\r
- uLength = 0;\r
- }\r
- }\r
-\r
-// Special case for last edge\r
- if (uLength >= g_uMinDiagLength)\r
- Add(uStartPosA, uStartPosB, uLength);\r
- }\r
-\r
-bool DiagList::NonZeroIntersection(const Diag &d) const\r
- {\r
- for (unsigned n = 0; n < m_uCount; ++n)\r
- {\r
- const Diag &d2 = m_Diags[n];\r
- if (DiagOverlap(d, d2) > 0)\r
- return true;\r
- }\r
- return false;\r
- }\r
-\r
-// DialogOverlap returns the length of the overlapping\r
-// section of the two diagonals along the diagonals\r
-// themselves; in other words, the length of\r
-// the intersection of the two sets of cells in\r
-// the matrix.\r
-unsigned DiagOverlap(const Diag &d1, const Diag &d2)\r
- {\r
-// Determine where the diagonals intersect the A\r
-// axis (extending them if required). If they\r
-// intersect at different points, they do not\r
-// overlap. Coordinates on a diagonal are\r
-// given by B = A + c where c is the value of\r
-// A at the intersection with the A axis.\r
-// Hence, c = B - A for any point on the diagonal.\r
- int c1 = (int) d1.m_uStartPosB - (int) d1.m_uStartPosA;\r
- int c2 = (int) d2.m_uStartPosB - (int) d2.m_uStartPosA;\r
- if (c1 != c2)\r
- return 0;\r
-\r
- assert(DiagOverlapA(d1, d2) == DiagOverlapB(d1, d2));\r
- return DiagOverlapA(d1, d2);\r
- }\r
-\r
-// DialogOverlapA returns the length of the overlapping\r
-// section of the projection of the two diagonals onto\r
-// the A axis.\r
-unsigned DiagOverlapA(const Diag &d1, const Diag &d2)\r
- {\r
- unsigned uMaxStart = MAX(d1.m_uStartPosA, d2.m_uStartPosA);\r
- unsigned uMinEnd = MIN(d1.m_uStartPosA + d1.m_uLength - 1,\r
- d2.m_uStartPosA + d2.m_uLength - 1);\r
-\r
- int iLength = (int) uMinEnd - (int) uMaxStart + 1;\r
- if (iLength < 0)\r
- return 0;\r
- return (unsigned) iLength;\r
- }\r
-\r
-// DialogOverlapB returns the length of the overlapping\r
-// section of the projection of the two diagonals onto\r
-// the B axis.\r
-unsigned DiagOverlapB(const Diag &d1, const Diag &d2)\r
- {\r
- unsigned uMaxStart = MAX(d1.m_uStartPosB, d2.m_uStartPosB);\r
- unsigned uMinEnd = MIN(d1.m_uStartPosB + d1.m_uLength - 1,\r
- d2.m_uStartPosB + d2.m_uLength - 1);\r
-\r
- int iLength = (int) uMinEnd - (int) uMaxStart + 1;\r
- if (iLength < 0)\r
- return 0;\r
- return (unsigned) iLength;\r
- }\r
-\r
-// Returns true if the two diagonals can be on the\r
-// same path through the DP matrix. If DiagCompatible\r
-// returns false, they cannot be in the same path\r
-// and hence "contradict" each other.\r
-bool DiagCompatible(const Diag &d1, const Diag &d2)\r
- {\r
- if (DiagOverlap(d1, d2) > 0)\r
- return true;\r
- return 0 == DiagOverlapA(d1, d2) && 0 == DiagOverlapB(d1, d2);\r
- }\r
-\r
-// Returns the length of the "break" between two diagonals.\r
-unsigned DiagBreak(const Diag &d1, const Diag &d2)\r
- {\r
- int c1 = (int) d1.m_uStartPosB - (int) d1.m_uStartPosA;\r
- int c2 = (int) d2.m_uStartPosB - (int) d2.m_uStartPosA;\r
- if (c1 != c2)\r
- return 0;\r
-\r
- int iMaxStart = MAX(d1.m_uStartPosA, d2.m_uStartPosA);\r
- int iMinEnd = MIN(d1.m_uStartPosA + d1.m_uLength - 1,\r
- d2.m_uStartPosA + d1.m_uLength - 1);\r
- int iBreak = iMaxStart - iMinEnd - 1;\r
- if (iBreak < 0)\r
- return 0;\r
- return (unsigned) iBreak;\r
- }\r
-\r
-// Merge diagonals that are continuations of each other with\r
-// short breaks of up to length g_uMaxDiagBreak.\r
-// In a sorted list of diagonals, we only have to check\r
-// consecutive entries.\r
-void MergeDiags(DiagList &DL)\r
- {\r
- return;\r
-#if DEBUG\r
- if (!DL.IsSorted())\r
- Quit("MergeDiags: !IsSorted");\r
-#endif\r
-\r
-// TODO: Fix this!\r
-// Breaks must be with no offset (no gaps)\r
- const unsigned uCount = DL.GetCount();\r
- if (uCount <= 1)\r
- return;\r
-\r
- DiagList NewList;\r
-\r
- Diag MergedDiag;\r
- const Diag *ptrPrev = &DL.Get(0);\r
- for (unsigned i = 1; i < uCount; ++i)\r
- {\r
- const Diag *ptrDiag = &DL.Get(i);\r
- unsigned uBreakLength = DiagBreak(*ptrPrev, *ptrDiag);\r
- if (uBreakLength <= g_uMaxDiagBreak)\r
- {\r
- MergedDiag.m_uStartPosA = ptrPrev->m_uStartPosA;\r
- MergedDiag.m_uStartPosB = ptrPrev->m_uStartPosB;\r
- MergedDiag.m_uLength = ptrPrev->m_uLength + ptrDiag->m_uLength\r
- + uBreakLength;\r
- ptrPrev = &MergedDiag;\r
- }\r
- else\r
- {\r
- NewList.Add(*ptrPrev);\r
- ptrPrev = ptrDiag;\r
- }\r
- }\r
- NewList.Add(*ptrPrev);\r
- DL.Copy(NewList);\r
- }\r
-\r
-void DiagList::DeleteIncompatible()\r
- {\r
- assert(IsSorted());\r
-\r
- if (m_uCount < 2)\r
- return;\r
-\r
- bool *bFlagForDeletion = new bool[m_uCount];\r
- for (unsigned i = 0; i < m_uCount; ++i)\r
- bFlagForDeletion[i] = false;\r
-\r
- for (unsigned i = 0; i < m_uCount; ++i)\r
- {\r
- const Diag &di = m_Diags[i];\r
- for (unsigned j = i + 1; j < m_uCount; ++j)\r
- {\r
- const Diag &dj = m_Diags[j];\r
-\r
- // Verify sorted correctly\r
- assert(di.m_uStartPosA <= dj.m_uStartPosA);\r
-\r
- // If two diagonals are incompatible and\r
- // one is is much longer than the other,\r
- // keep the longer one.\r
- if (!DiagCompatible(di, dj))\r
- {\r
- if (di.m_uLength > dj.m_uLength*4)\r
- bFlagForDeletion[j] = true;\r
- else if (dj.m_uLength > di.m_uLength*4)\r
- bFlagForDeletion[i] = true;\r
- else\r
- {\r
- bFlagForDeletion[i] = true;\r
- bFlagForDeletion[j] = true;\r
- }\r
- }\r
- }\r
- }\r
-\r
- for (unsigned i = 0; i < m_uCount; ++i)\r
- {\r
- const Diag &di = m_Diags[i];\r
- if (bFlagForDeletion[i])\r
- continue;\r
-\r
- for (unsigned j = i + 1; j < m_uCount; ++j)\r
- {\r
- const Diag &dj = m_Diags[j];\r
- if (bFlagForDeletion[j])\r
- continue;\r
-\r
- // Verify sorted correctly\r
- assert(di.m_uStartPosA <= dj.m_uStartPosA);\r
-\r
- // If sort order in B different from sorted order in A,\r
- // either diags are incompatible or we detected a repeat\r
- // or permutation.\r
- if (di.m_uStartPosB >= dj.m_uStartPosB || !DiagCompatible(di, dj))\r
- {\r
- bFlagForDeletion[i] = true;\r
- bFlagForDeletion[j] = true;\r
- }\r
- }\r
- }\r
-\r
- unsigned uNewCount = 0;\r
- Diag *NewDiags = new Diag[m_uCount];\r
- for (unsigned i = 0; i < m_uCount; ++i)\r
- {\r
- if (bFlagForDeletion[i])\r
- continue;\r
-\r
- const Diag &d = m_Diags[i];\r
- NewDiags[uNewCount] = d;\r
- ++uNewCount;\r
- }\r
- memcpy(m_Diags, NewDiags, uNewCount*sizeof(Diag));\r
- m_uCount = uNewCount;\r
- delete[] NewDiags;\r
- }\r
-\r
-void DiagList::Copy(const DiagList &DL)\r
- {\r
- Clear();\r
- unsigned uCount = DL.GetCount();\r
- for (unsigned i = 0; i < uCount; ++i)\r
- Add(DL.Get(i));\r
- }\r
-\r
-// Check if sorted in increasing order of m_uStartPosA\r
-bool DiagList::IsSorted() const\r
- {\r
- return true;\r
- unsigned uCount = GetCount();\r
- for (unsigned i = 1; i < uCount; ++i)\r
- if (m_Diags[i-1].m_uStartPosA > m_Diags[i].m_uStartPosA)\r
- return false;\r
- return true;\r
- }\r
-\r
-// Sort in increasing order of m_uStartPosA\r
-// Dumb bubble sort, but don't care about speed\r
-// because don't get long lists.\r
-void DiagList::Sort()\r
- {\r
- if (m_uCount < 2)\r
- return;\r
-\r
- bool bContinue = true;\r
- while (bContinue)\r
- {\r
- bContinue = false;\r
- for (unsigned i = 0; i < m_uCount - 1; ++i)\r
- {\r
- if (m_Diags[i].m_uStartPosA > m_Diags[i+1].m_uStartPosA)\r
- {\r
- Diag Tmp = m_Diags[i];\r
- m_Diags[i] = m_Diags[i+1];\r
- m_Diags[i+1] = Tmp;\r
- bContinue = true;\r
- }\r
- }\r
- }\r
- }\r
-\r
-//void TestDiag()\r
-// {\r
-// Diag d1;\r
-// Diag d2;\r
-// Diag d3;\r
-//\r
-// d1.m_uStartPosA = 0;\r
-// d1.m_uStartPosB = 1;\r
-// d1.m_uLength = 32;\r
-//\r
-// d2.m_uStartPosA = 55;\r
-// d2.m_uStartPosB = 70;\r
-// d2.m_uLength = 36;\r
-//\r
-// d3.m_uStartPosA = 102;\r
-// d3.m_uStartPosB = 122;\r
-// d3.m_uLength = 50;\r
-//\r
-// DiagList DL;\r
-// DL.Add(d1);\r
-// DL.Add(d2);\r
-// DL.Add(d3);\r
-//\r
-// Log("Before DeleteIncompatible:\n");\r
-// DL.LogMe();\r
-// DL.DeleteIncompatible();\r
-//\r
-// Log("After DeleteIncompatible:\n");\r
-// DL.LogMe();\r
-//\r
-// MergeDiags(DL);\r
-// Log("After Merge:\n");\r
-// DL.LogMe();\r
-//\r
-// DPRegionList RL;\r
-// DiagListToDPRegionList(DL, RL, 200, 200);\r
-// RL.LogMe();\r
-// }\r
git://source.jalview.org / jabaws.git / blobdiff