--- /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