binaries/src/muscle/diaglist.cpp

   1 #include "muscle.h"\r
   2 #include "diaglist.h"\r
   3 #include "pwpath.h"\r
   4 \r
   5 #define MAX(x, y)       ((x) > (y) ? (x) : (y))\r
   6 #define MIN(x, y)       ((x) < (y) ? (x) : (y))\r
   7 \r
   8 void DiagList::Add(const Diag &d)\r
   9         {\r
  10         if (m_uCount == MAX_DIAGS)\r
  11                 Quit("DiagList::Add, overflow %u", m_uCount);\r
  12         m_Diags[m_uCount] = d;\r
  13         ++m_uCount;\r
  14         }\r
  15 \r
  16 void DiagList::Add(unsigned uStartPosA, unsigned uStartPosB, unsigned uLength)\r
  17         {\r
  18         Diag d;\r
  19         d.m_uStartPosA = uStartPosA;\r
  20         d.m_uStartPosB = uStartPosB;\r
  21         d.m_uLength = uLength;\r
  22         Add(d);\r
  23         }\r
  24 \r
  25 const Diag &DiagList::Get(unsigned uIndex) const\r
  26         {\r
  27         if (uIndex >= m_uCount)\r
  28                 Quit("DiagList::Get(%u), count=%u", uIndex, m_uCount);\r
  29         return m_Diags[uIndex];\r
  30         }\r
  31 \r
  32 void DiagList::LogMe() const\r
  33         {\r
  34         Log("DiagList::LogMe, count=%u\n", m_uCount);\r
  35         Log("  n  StartA  StartB  Length\n");\r
  36         Log("---  ------  ------  ------\n");\r
  37         for (unsigned n = 0; n < m_uCount; ++n)\r
  38                 {\r
  39                 const Diag &d = m_Diags[n];\r
  40                 Log("%3u  %6u  %6u  %6u\n",\r
  41                   n, d.m_uStartPosA, d.m_uStartPosB, d.m_uLength);\r
  42                 }\r
  43         }\r
  44 \r
  45 void DiagList::FromPath(const PWPath &Path)\r
  46         {\r
  47         Clear();\r
  48 \r
  49         const unsigned uEdgeCount = Path.GetEdgeCount();\r
  50         unsigned uLength = 0;\r
  51         unsigned uStartPosA;\r
  52         unsigned uStartPosB;\r
  53         for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
  54                 {\r
  55                 const PWEdge &Edge = Path.GetEdge(uEdgeIndex);\r
  56 \r
  57         // Typical cases\r
  58                 if (Edge.cType == 'M')\r
  59                         {\r
  60                         if (0 == uLength)\r
  61                                 {\r
  62                                 uStartPosA = Edge.uPrefixLengthA - 1;\r
  63                                 uStartPosB = Edge.uPrefixLengthB - 1;\r
  64                                 }\r
  65                         ++uLength;\r
  66                         }\r
  67                 else\r
  68                         {\r
  69                         if (uLength >= g_uMinDiagLength)\r
  70                                 Add(uStartPosA, uStartPosB, uLength);\r
  71                         uLength = 0;\r
  72                         }\r
  73                 }\r
  74 \r
  75 // Special case for last edge\r
  76         if (uLength >= g_uMinDiagLength)\r
  77                 Add(uStartPosA, uStartPosB, uLength);\r
  78         }\r
  79 \r
  80 bool DiagList::NonZeroIntersection(const Diag &d) const\r
  81         {\r
  82         for (unsigned n = 0; n < m_uCount; ++n)\r
  83                 {\r
  84                 const Diag &d2 = m_Diags[n];\r
  85                 if (DiagOverlap(d, d2) > 0)\r
  86                         return true;\r
  87                 }\r
  88         return false;\r
  89         }\r
  90 \r
  91 // DialogOverlap returns the length of the overlapping\r
  92 // section of the two diagonals along the diagonals\r
  93 // themselves; in other words, the length of\r
  94 // the intersection of the two sets of cells in\r
  95 // the matrix.\r
  96 unsigned DiagOverlap(const Diag &d1, const Diag &d2)\r
  97         {\r
  98 // Determine where the diagonals intersect the A\r
  99 // axis (extending them if required). If they\r
 100 // intersect at different points, they do not\r
 101 // overlap. Coordinates on a diagonal are\r
 102 // given by B = A + c where c is the value of\r
 103 // A at the intersection with the A axis.\r
 104 // Hence, c = B - A for any point on the diagonal.\r
 105         int c1 = (int) d1.m_uStartPosB - (int) d1.m_uStartPosA;\r
 106         int c2 = (int) d2.m_uStartPosB - (int) d2.m_uStartPosA;\r
 107         if (c1 != c2)\r
 108                 return 0;\r
 109 \r
 110         assert(DiagOverlapA(d1, d2) == DiagOverlapB(d1, d2));\r
 111         return DiagOverlapA(d1, d2);\r
 112         }\r
 113 \r
 114 // DialogOverlapA returns the length of the overlapping\r
 115 // section of the projection of the two diagonals onto\r
 116 // the A axis.\r
 117 unsigned DiagOverlapA(const Diag &d1, const Diag &d2)\r
 118         {\r
 119         unsigned uMaxStart = MAX(d1.m_uStartPosA, d2.m_uStartPosA);\r
 120         unsigned uMinEnd = MIN(d1.m_uStartPosA + d1.m_uLength - 1,\r
 121           d2.m_uStartPosA + d2.m_uLength - 1);\r
 122 \r
 123         int iLength = (int) uMinEnd - (int) uMaxStart + 1;\r
 124         if (iLength < 0)\r
 125                 return 0;\r
 126         return (unsigned) iLength;\r
 127         }\r
 128 \r
 129 // DialogOverlapB returns the length of the overlapping\r
 130 // section of the projection of the two diagonals onto\r
 131 // the B axis.\r
 132 unsigned DiagOverlapB(const Diag &d1, const Diag &d2)\r
 133         {\r
 134         unsigned uMaxStart = MAX(d1.m_uStartPosB, d2.m_uStartPosB);\r
 135         unsigned uMinEnd = MIN(d1.m_uStartPosB + d1.m_uLength - 1,\r
 136           d2.m_uStartPosB + d2.m_uLength - 1);\r
 137 \r
 138         int iLength = (int) uMinEnd - (int) uMaxStart + 1;\r
 139         if (iLength < 0)\r
 140                 return 0;\r
 141         return (unsigned) iLength;\r
 142         }\r
 143 \r
 144 // Returns true if the two diagonals can be on the\r
 145 // same path through the DP matrix. If DiagCompatible\r
 146 // returns false, they cannot be in the same path\r
 147 // and hence "contradict" each other.\r
 148 bool DiagCompatible(const Diag &d1, const Diag &d2)\r
 149         {\r
 150         if (DiagOverlap(d1, d2) > 0)\r
 151                 return true;\r
 152         return 0 == DiagOverlapA(d1, d2) && 0 == DiagOverlapB(d1, d2);\r
 153         }\r
 154 \r
 155 // Returns the length of the "break" between two diagonals.\r
 156 unsigned DiagBreak(const Diag &d1, const Diag &d2)\r
 157         {\r
 158         int c1 = (int) d1.m_uStartPosB - (int) d1.m_uStartPosA;\r
 159         int c2 = (int) d2.m_uStartPosB - (int) d2.m_uStartPosA;\r
 160         if (c1 != c2)\r
 161                 return 0;\r
 162 \r
 163         int iMaxStart = MAX(d1.m_uStartPosA, d2.m_uStartPosA);\r
 164         int iMinEnd = MIN(d1.m_uStartPosA + d1.m_uLength - 1,\r
 165           d2.m_uStartPosA + d1.m_uLength - 1);\r
 166         int iBreak = iMaxStart - iMinEnd - 1;\r
 167         if (iBreak < 0)\r
 168                 return 0;\r
 169         return (unsigned) iBreak;\r
 170         }\r
 171 \r
 172 // Merge diagonals that are continuations of each other with\r
 173 // short breaks of up to length g_uMaxDiagBreak.\r
 174 // In a sorted list of diagonals, we only have to check\r
 175 // consecutive entries.\r
 176 void MergeDiags(DiagList &DL)\r
 177         {\r
 178         return;\r
 179 #if     DEBUG\r
 180         if (!DL.IsSorted())\r
 181                 Quit("MergeDiags: !IsSorted");\r
 182 #endif\r
 183 \r
 184 // TODO: Fix this!\r
 185 // Breaks must be with no offset (no gaps)\r
 186         const unsigned uCount = DL.GetCount();\r
 187         if (uCount <= 1)\r
 188                 return;\r
 189 \r
 190         DiagList NewList;\r
 191 \r
 192         Diag MergedDiag;\r
 193         const Diag *ptrPrev = &DL.Get(0);\r
 194         for (unsigned i = 1; i < uCount; ++i)\r
 195                 {\r
 196                 const Diag *ptrDiag = &DL.Get(i);\r
 197                 unsigned uBreakLength = DiagBreak(*ptrPrev, *ptrDiag);\r
 198                 if (uBreakLength <= g_uMaxDiagBreak)\r
 199                         {\r
 200                         MergedDiag.m_uStartPosA = ptrPrev->m_uStartPosA;\r
 201                         MergedDiag.m_uStartPosB = ptrPrev->m_uStartPosB;\r
 202                         MergedDiag.m_uLength = ptrPrev->m_uLength + ptrDiag->m_uLength\r
 203                           + uBreakLength;\r
 204                         ptrPrev = &MergedDiag;\r
 205                         }\r
 206                 else\r
 207                         {\r
 208                         NewList.Add(*ptrPrev);\r
 209                         ptrPrev = ptrDiag;\r
 210                         }\r
 211                 }\r
 212         NewList.Add(*ptrPrev);\r
 213         DL.Copy(NewList);\r
 214         }\r
 215 \r
 216 void DiagList::DeleteIncompatible()\r
 217         {\r
 218         assert(IsSorted());\r
 219 \r
 220         if (m_uCount < 2)\r
 221                 return;\r
 222 \r
 223         bool *bFlagForDeletion = new bool[m_uCount];\r
 224         for (unsigned i = 0; i < m_uCount; ++i)\r
 225                 bFlagForDeletion[i] = false;\r
 226 \r
 227         for (unsigned i = 0; i < m_uCount; ++i)\r
 228                 {\r
 229                 const Diag &di = m_Diags[i];\r
 230                 for (unsigned j = i + 1; j < m_uCount; ++j)\r
 231                         {\r
 232                         const Diag &dj = m_Diags[j];\r
 233 \r
 234                 // Verify sorted correctly\r
 235                         assert(di.m_uStartPosA <= dj.m_uStartPosA);\r
 236 \r
 237                 // If two diagonals are incompatible and\r
 238                 // one is is much longer than the other,\r
 239                 // keep the longer one.\r
 240                         if (!DiagCompatible(di, dj))\r
 241                                 {\r
 242                                 if (di.m_uLength > dj.m_uLength*4)\r
 243                                         bFlagForDeletion[j] = true;\r
 244                                 else if (dj.m_uLength > di.m_uLength*4)\r
 245                                         bFlagForDeletion[i] = true;\r
 246                                 else\r
 247                                         {\r
 248                                         bFlagForDeletion[i] = true;\r
 249                                         bFlagForDeletion[j] = true;\r
 250                                         }\r
 251                                 }\r
 252                         }\r
 253                 }\r
 254 \r
 255         for (unsigned i = 0; i < m_uCount; ++i)\r
 256                 {\r
 257                 const Diag &di = m_Diags[i];\r
 258                 if (bFlagForDeletion[i])\r
 259                         continue;\r
 260 \r
 261                 for (unsigned j = i + 1; j < m_uCount; ++j)\r
 262                         {\r
 263                         const Diag &dj = m_Diags[j];\r
 264                         if (bFlagForDeletion[j])\r
 265                                 continue;\r
 266 \r
 267                 // Verify sorted correctly\r
 268                         assert(di.m_uStartPosA <= dj.m_uStartPosA);\r
 269 \r
 270                 // If sort order in B different from sorted order in A,\r
 271                 // either diags are incompatible or we detected a repeat\r
 272                 // or permutation.\r
 273                         if (di.m_uStartPosB >= dj.m_uStartPosB || !DiagCompatible(di, dj))\r
 274                                 {\r
 275                                 bFlagForDeletion[i] = true;\r
 276                                 bFlagForDeletion[j] = true;\r
 277                                 }\r
 278                         }\r
 279                 }\r
 280 \r
 281         unsigned uNewCount = 0;\r
 282         Diag *NewDiags = new Diag[m_uCount];\r
 283         for (unsigned i = 0; i < m_uCount; ++i)\r
 284                 {\r
 285                 if (bFlagForDeletion[i])\r
 286                         continue;\r
 287 \r
 288                 const Diag &d = m_Diags[i];\r
 289                 NewDiags[uNewCount] = d;\r
 290                 ++uNewCount;\r
 291                 }\r
 292         memcpy(m_Diags, NewDiags, uNewCount*sizeof(Diag));\r
 293         m_uCount = uNewCount;\r
 294         delete[] NewDiags;\r
 295         }\r
 296 \r
 297 void DiagList::Copy(const DiagList &DL)\r
 298         {\r
 299         Clear();\r
 300         unsigned uCount = DL.GetCount();\r
 301         for (unsigned i = 0; i < uCount; ++i)\r
 302                 Add(DL.Get(i));\r
 303         }\r
 304 \r
 305 // Check if sorted in increasing order of m_uStartPosA\r
 306 bool DiagList::IsSorted() const\r
 307         {\r
 308         return true;\r
 309         unsigned uCount = GetCount();\r
 310         for (unsigned i = 1; i < uCount; ++i)\r
 311                 if (m_Diags[i-1].m_uStartPosA > m_Diags[i].m_uStartPosA)\r
 312                         return false;\r
 313         return true;\r
 314         }\r
 315 \r
 316 // Sort in increasing order of m_uStartPosA\r
 317 // Dumb bubble sort, but don't care about speed\r
 318 // because don't get long lists.\r
 319 void DiagList::Sort()\r
 320         {\r
 321         if (m_uCount < 2)\r
 322                 return;\r
 323 \r
 324         bool bContinue = true;\r
 325         while (bContinue)\r
 326                 {\r
 327                 bContinue = false;\r
 328                 for (unsigned i = 0; i < m_uCount - 1; ++i)\r
 329                         {\r
 330                         if (m_Diags[i].m_uStartPosA > m_Diags[i+1].m_uStartPosA)\r
 331                                 {\r
 332                                 Diag Tmp = m_Diags[i];\r
 333                                 m_Diags[i] = m_Diags[i+1];\r
 334                                 m_Diags[i+1] = Tmp;\r
 335                                 bContinue = true;\r
 336                                 }\r
 337                         }\r
 338                 }\r
 339         }\r
 340 \r
 341 //void TestDiag()\r
 342 //      {\r
 343 //      Diag d1;\r
 344 //      Diag d2;\r
 345 //      Diag d3;\r
 346 //\r
 347 //      d1.m_uStartPosA = 0;\r
 348 //      d1.m_uStartPosB = 1;\r
 349 //      d1.m_uLength = 32;\r
 350 //\r
 351 //      d2.m_uStartPosA = 55;\r
 352 //      d2.m_uStartPosB = 70;\r
 353 //      d2.m_uLength = 36;\r
 354 //\r
 355 //      d3.m_uStartPosA = 102;\r
 356 //      d3.m_uStartPosB = 122;\r
 357 //      d3.m_uLength = 50;\r
 358 //\r
 359 //      DiagList DL;\r
 360 //      DL.Add(d1);\r
 361 //      DL.Add(d2);\r
 362 //      DL.Add(d3);\r
 363 //\r
 364 //      Log("Before DeleteIncompatible:\n");\r
 365 //      DL.LogMe();\r
 366 //      DL.DeleteIncompatible();\r
 367 //\r
 368 //      Log("After DeleteIncompatible:\n");\r
 369 //      DL.LogMe();\r
 370 //\r
 371 //      MergeDiags(DL);\r
 372 //      Log("After Merge:\n");\r
 373 //      DL.LogMe();\r
 374 //\r
 375 //      DPRegionList RL;\r
 376 //      DiagListToDPRegionList(DL, RL, 200, 200);\r
 377 //      RL.LogMe();\r
 378 //      }\r