binaries/src/muscle/glbaligndiag.cpp

   1 #include "muscle.h"\r
   2 #include "dpreglist.h"\r
   3 #include "diaglist.h"\r
   4 #include "pwpath.h"\r
   5 #include "profile.h"\r
   6 #include "timing.h"\r
   7 \r
   8 #define TRACE           0\r
   9 #define TRACE_PATH      0\r
  10 #define LIST_DIAGS      0\r
  11 \r
  12 static double g_dDPAreaWithoutDiags = 0.0;\r
  13 static double g_dDPAreaWithDiags = 0.0;\r
  14 \r
  15 static void OffsetPath(PWPath &Path, unsigned uOffsetA, unsigned uOffsetB)\r
  16         {\r
  17         const unsigned uEdgeCount = Path.GetEdgeCount();\r
  18         for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
  19                 {\r
  20                 const PWEdge &Edge = Path.GetEdge(uEdgeIndex);\r
  21 \r
  22         // Nasty hack -- poke new values back into path, circumventing class\r
  23                 PWEdge &NonConstEdge = (PWEdge &) Edge;\r
  24                 NonConstEdge.uPrefixLengthA += uOffsetA;\r
  25                 NonConstEdge.uPrefixLengthB += uOffsetB;\r
  26                 }\r
  27         }\r
  28 \r
  29 static void DiagToPath(const Diag &d, PWPath &Path)\r
  30         {\r
  31         Path.Clear();\r
  32         const unsigned uLength = d.m_uLength;\r
  33         for (unsigned i = 0; i < uLength; ++i)\r
  34                 {\r
  35                 PWEdge Edge;\r
  36                 Edge.cType = 'M';\r
  37                 Edge.uPrefixLengthA = d.m_uStartPosA + i + 1;\r
  38                 Edge.uPrefixLengthB = d.m_uStartPosB + i + 1;\r
  39                 Path.AppendEdge(Edge);\r
  40                 }\r
  41         }\r
  42 \r
  43 static void AppendRegPath(PWPath &Path, const PWPath &RegPath)\r
  44         {\r
  45         const unsigned uRegEdgeCount = RegPath.GetEdgeCount();\r
  46         for (unsigned uRegEdgeIndex = 0; uRegEdgeIndex < uRegEdgeCount; ++uRegEdgeIndex)\r
  47                 {\r
  48                 const PWEdge &RegEdge = RegPath.GetEdge(uRegEdgeIndex);\r
  49                 Path.AppendEdge(RegEdge);\r
  50                 }\r
  51         }\r
  52 \r
  53 SCORE GlobalAlignDiags(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
  54   unsigned uLengthB, PWPath &Path)\r
  55         {\r
  56 #if     LIST_DIAGS\r
  57         TICKS t1 = GetClockTicks();\r
  58 #endif\r
  59 \r
  60         DiagList DL;\r
  61 \r
  62         if (ALPHA_Amino == g_Alpha)\r
  63                 FindDiags(PA, uLengthA, PB, uLengthB, DL);\r
  64         else if (ALPHA_DNA == g_Alpha || ALPHA_RNA == g_Alpha)\r
  65                 FindDiagsNuc(PA, uLengthA, PB, uLengthB, DL);\r
  66         else\r
  67                 Quit("GlobalAlignDiags: bad alpha");\r
  68 \r
  69 #if     TRACE\r
  70         Log("GlobalAlignDiags, diag list:\n");\r
  71         DL.LogMe();\r
  72 #endif\r
  73 \r
  74         DL.Sort();\r
  75         DL.DeleteIncompatible();\r
  76 \r
  77 #if     TRACE\r
  78         Log("After DeleteIncompatible:\n");\r
  79         DL.LogMe();\r
  80 #endif\r
  81 \r
  82         MergeDiags(DL);\r
  83 \r
  84 #if     TRACE\r
  85         Log("After MergeDiags:\n");\r
  86         DL.LogMe();\r
  87 #endif\r
  88 \r
  89         DPRegionList RL;\r
  90         DiagListToDPRegionList(DL, RL, uLengthA, uLengthB);\r
  91 \r
  92 #if     TRACE\r
  93         Log("RegionList:\n");\r
  94         RL.LogMe();\r
  95 #endif\r
  96 \r
  97 #if     LIST_DIAGS\r
  98         {\r
  99         TICKS t2 = GetClockTicks();\r
 100         unsigned uArea = RL.GetDPArea();\r
 101         Log("ticks=%ld\n", (long) (t2 - t1));\r
 102         Log("area=%u\n", uArea);\r
 103         }\r
 104 #endif\r
 105 \r
 106         g_dDPAreaWithoutDiags += uLengthA*uLengthB;\r
 107 \r
 108         double dDPAreaWithDiags = 0.0;\r
 109         const unsigned uRegionCount = RL.GetCount();\r
 110         for (unsigned uRegionIndex = 0; uRegionIndex < uRegionCount; ++uRegionIndex)\r
 111                 {\r
 112                 const DPRegion &r = RL.Get(uRegionIndex);\r
 113 \r
 114                 PWPath RegPath;\r
 115                 if (DPREGIONTYPE_Diag == r.m_Type)\r
 116                         {\r
 117                         DiagToPath(r.m_Diag, RegPath);\r
 118 #if     TRACE_PATH\r
 119                         Log("DiagToPath, path=\n");\r
 120                         RegPath.LogMe();\r
 121 #endif\r
 122                         }\r
 123                 else if (DPREGIONTYPE_Rect == r.m_Type)\r
 124                         {\r
 125                         const unsigned uRegStartPosA = r.m_Rect.m_uStartPosA;\r
 126                         const unsigned uRegStartPosB = r.m_Rect.m_uStartPosB;\r
 127                         const unsigned uRegLengthA = r.m_Rect.m_uLengthA;\r
 128                         const unsigned uRegLengthB = r.m_Rect.m_uLengthB;\r
 129                         const ProfPos *RegPA = PA + uRegStartPosA;\r
 130                         const ProfPos *RegPB = PB + uRegStartPosB;\r
 131 \r
 132                         dDPAreaWithDiags += uRegLengthA*uRegLengthB;\r
 133                         GlobalAlignNoDiags(RegPA, uRegLengthA, RegPB, uRegLengthB, RegPath);\r
 134 #if     TRACE_PATH\r
 135                         Log("GlobalAlignNoDiags RegPath=\n");\r
 136                         RegPath.LogMe();\r
 137 #endif\r
 138                         OffsetPath(RegPath, uRegStartPosA, uRegStartPosB);\r
 139 #if     TRACE_PATH\r
 140                         Log("After offset path, RegPath=\n");\r
 141                         RegPath.LogMe();\r
 142 #endif\r
 143                         }\r
 144                 else\r
 145                         Quit("GlobalAlignDiags, Invalid region type %u", r.m_Type);\r
 146 \r
 147                 AppendRegPath(Path, RegPath);\r
 148 #if     TRACE_PATH\r
 149                 Log("After AppendPath, path=");\r
 150                 Path.LogMe();\r
 151 #endif\r
 152                 }\r
 153 \r
 154 #if     TRACE\r
 155         {\r
 156         double dDPAreaWithoutDiags = uLengthA*uLengthB;\r
 157         Log("DP area with diags %.3g without %.3g pct saved %.3g %%\n",\r
 158           dDPAreaWithDiags, dDPAreaWithoutDiags, (1.0 - dDPAreaWithDiags/dDPAreaWithoutDiags)*100.0);\r
 159         }\r
 160 #endif\r
 161         g_dDPAreaWithDiags += dDPAreaWithDiags;\r
 162         return 0;\r
 163         }\r
 164 \r
 165 void ListDiagSavings()\r
 166         {\r
 167         if (!g_bVerbose || !g_bDiags)\r
 168                 return;\r
 169         double dAreaSaved = g_dDPAreaWithoutDiags - g_dDPAreaWithDiags;\r
 170         double dPct = dAreaSaved*100.0/g_dDPAreaWithoutDiags;\r
 171         Log("DP area saved by diagonals %-4.1f%%\n", dPct);\r
 172         }\r