website/archive/binaries/mac/src/muscle/tracebacksw.cpp

   1 #include "muscle.h"\r
   2 #include "profile.h"\r
   3 #include "pwpath.h"\r
   4 #include <math.h>\r
   5 \r
   6 #define TRACE   0\r
   7 \r
   8 #define EQ(a, b)        (fabs(a-b) < 0.1)\r
   9 \r
  10 void TraceBackSW(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
  11   unsigned uLengthB, const SCORE *DPM_, const SCORE *DPD_, const SCORE *DPI_,\r
  12   unsigned uPrefixLengthAMax, unsigned uPrefixLengthBMax, PWPath &Path)\r
  13         {\r
  14 #if     TRACE\r
  15         Log("\n");\r
  16         Log("TraceBackSW LengthA=%u LengthB=%u PLAMax=%u PLBMax=%u\n",\r
  17           uLengthA, uLengthB, uPrefixLengthAMax, uPrefixLengthBMax);\r
  18 #endif\r
  19         assert(uLengthB > 0 && uLengthA > 0);\r
  20 \r
  21         const unsigned uPrefixCountA = uLengthA + 1;\r
  22         const unsigned uPrefixCountB = uLengthB + 1;\r
  23 \r
  24         Path.Clear();\r
  25 \r
  26         unsigned uPrefixLengthA = uPrefixLengthAMax;\r
  27         unsigned uPrefixLengthB = uPrefixLengthBMax;\r
  28 \r
  29         SCORE scoreMax = DPM(uPrefixLengthA, uPrefixLengthB);\r
  30         char cEdgeType = 'M';\r
  31 \r
  32         for (;;)\r
  33                 {\r
  34                 if ('S' == cEdgeType)\r
  35                         break;\r
  36 \r
  37                 PWEdge Edge;\r
  38                 Edge.cType = cEdgeType;\r
  39                 Edge.uPrefixLengthA = uPrefixLengthA;\r
  40                 Edge.uPrefixLengthB = uPrefixLengthB;\r
  41                 Path.PrependEdge(Edge);\r
  42 \r
  43                 char cPrevEdgeType;\r
  44                 unsigned uPrevPrefixLengthA = uPrefixLengthA;\r
  45                 unsigned uPrevPrefixLengthB = uPrefixLengthB;\r
  46 \r
  47                 switch (cEdgeType)\r
  48                         {\r
  49                 case 'M':\r
  50                         {\r
  51                         assert(uPrefixLengthA > 0);\r
  52                         assert(uPrefixLengthB > 0);\r
  53                         const ProfPos &PPA = PA[uPrefixLengthA - 1];\r
  54                         const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
  55 \r
  56                         const SCORE Score = DPM(uPrefixLengthA, uPrefixLengthB);\r
  57                         const SCORE scoreMatch = ScoreProfPos2(PPA, PPB);\r
  58 \r
  59                         SCORE scoreSM;\r
  60                         if (1 == uPrefixLengthA && 1 == uPrefixLengthB)\r
  61                                 scoreSM = scoreMatch;\r
  62                         else\r
  63                                 scoreSM = MINUS_INFINITY;\r
  64 \r
  65                         SCORE scoreMM = MINUS_INFINITY;\r
  66                         SCORE scoreDM = MINUS_INFINITY;\r
  67                         SCORE scoreIM = MINUS_INFINITY;\r
  68                         if (uPrefixLengthA > 1 && uPrefixLengthB > 1)\r
  69                                 {\r
  70                                 SCORE scoreTrans = DPM(uPrefixLengthA-1, uPrefixLengthB-1);\r
  71                                 scoreMM = scoreTrans + scoreMatch;\r
  72                                 }\r
  73                         if (uPrefixLengthA > 1)\r
  74                                 {\r
  75                                 SCORE scoreTransDM = PA[uPrefixLengthA-2].m_scoreGapClose;\r
  76                                 scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreTransDM + scoreMatch;\r
  77                                 }\r
  78                         if (uPrefixLengthB > 1)\r
  79                                 {\r
  80                                 SCORE scoreTransIM = PB[uPrefixLengthB-2].m_scoreGapClose;\r
  81                                 scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreTransIM + scoreMatch;\r
  82                                 }\r
  83 \r
  84                         if (EQ(scoreMM, Score))\r
  85                                 cPrevEdgeType = 'M';\r
  86                         else if (EQ(scoreDM, Score))\r
  87                                 cPrevEdgeType = 'D';\r
  88                         else if (EQ(scoreIM, Score))\r
  89                                 cPrevEdgeType = 'I';\r
  90                         else if (EQ(scoreSM, Score))\r
  91                                 cPrevEdgeType = 'S';\r
  92                         else if (EQ(scoreMatch, Score))\r
  93                                 cPrevEdgeType = 'S';\r
  94                         else\r
  95                                 Quit("TraceBack2: failed to match M score=%g M=%g D=%g I=%g S=%g",\r
  96                                   Score, scoreMM, scoreDM, scoreIM, scoreSM);\r
  97 \r
  98                         --uPrevPrefixLengthA;\r
  99                         --uPrevPrefixLengthB;\r
 100                         break;\r
 101                         }\r
 102 \r
 103                 case 'D':\r
 104                         {\r
 105                         assert(uPrefixLengthA > 0);\r
 106                         const SCORE Score = DPD(uPrefixLengthA, uPrefixLengthB);\r
 107 \r
 108                         SCORE scoreMD = MINUS_INFINITY;\r
 109                         SCORE scoreDD = MINUS_INFINITY;\r
 110                         SCORE scoreSD = MINUS_INFINITY;\r
 111                         if (uPrefixLengthB == 0)\r
 112                                 {\r
 113                                 if (uPrefixLengthA == 1)\r
 114                                         scoreSD = PA[0].m_scoreGapOpen;\r
 115                                 else\r
 116                                         scoreSD = DPD(uPrefixLengthA - 1, 0);\r
 117                                 }\r
 118                         if (uPrefixLengthA > 1)\r
 119                                 {\r
 120                                 const ProfPos &PPA = PA[uPrefixLengthA - 1];\r
 121                                 SCORE scoreTransMD = PPA.m_scoreGapOpen;\r
 122                                 scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) + scoreTransMD;\r
 123                                 scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB);\r
 124                                 }\r
 125 \r
 126                         if (EQ(Score, scoreMD))\r
 127                                 cPrevEdgeType = 'M';\r
 128                         else if (EQ(Score, scoreDD))\r
 129                                 cPrevEdgeType = 'D';\r
 130                         else if (EQ(Score, scoreSD))\r
 131                                 cPrevEdgeType = 'S';\r
 132                         else\r
 133                                 Quit("TraceBack2: failed to match D");\r
 134 \r
 135                         --uPrevPrefixLengthA;\r
 136                         break;\r
 137                         }\r
 138 \r
 139                 case 'I':\r
 140                         {\r
 141                         assert(uPrefixLengthB > 0);\r
 142                         const SCORE Score = DPI(uPrefixLengthA, uPrefixLengthB);\r
 143 \r
 144                         SCORE scoreMI = MINUS_INFINITY;\r
 145                         SCORE scoreII = MINUS_INFINITY;\r
 146                         SCORE scoreSI = MINUS_INFINITY;\r
 147                         if (uPrefixLengthA == 0)\r
 148                                 {\r
 149                                 if (uPrefixLengthB == 1)\r
 150                                         scoreSI = PB[0].m_scoreGapOpen;\r
 151                                 else\r
 152                                         scoreSI = DPI(0, uPrefixLengthB - 1);\r
 153                                 }\r
 154                         if (uPrefixLengthB > 1)\r
 155                                 {\r
 156                                 const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
 157                                 SCORE scoreTransMI = PPB.m_scoreGapOpen;\r
 158                                 scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) + scoreTransMI;\r
 159                                 scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1);\r
 160                                 }\r
 161 \r
 162                         if (EQ(Score, scoreMI))\r
 163                                 cPrevEdgeType = 'M';\r
 164                         else if (EQ(Score, scoreII))\r
 165                                 cPrevEdgeType = 'I';\r
 166                         else if (EQ(Score, scoreSI))\r
 167                                 cPrevEdgeType = 'S';\r
 168                         else\r
 169                                 Quit("TraceBack2: failed to match I");\r
 170 \r
 171                         --uPrevPrefixLengthB;\r
 172                         break;\r
 173                         }\r
 174 \r
 175                 default:\r
 176                         assert(false);\r
 177                         }\r
 178 #if     TRACE\r
 179                 Log("Edge %c%c%u.%u", cPrevEdgeType, cEdgeType, uPrefixLengthA, uPrefixLengthB);\r
 180                 Log("\n");\r
 181 #endif\r
 182                 cEdgeType = cPrevEdgeType;\r
 183                 uPrefixLengthA = uPrevPrefixLengthA;\r
 184                 uPrefixLengthB = uPrevPrefixLengthB;\r
 185                 }\r
 186         }\r