--- /dev/null
+#include "muscle.h"\r
+#include "msa.h"\r
+#include "profile.h"\r
+#include "objscore.h"\r
+\r
+#if DOUBLE_AFFINE\r
+\r
+#define TRACE 0\r
+#define TEST_SPFAST 0\r
+\r
+static SCORE GapPenalty(unsigned uLength, bool Term, SCORE g, SCORE e)\r
+ {\r
+ //if (Term)\r
+ // {\r
+ // switch (g_TermGap)\r
+ // {\r
+ // case TERMGAP_Full:\r
+ // return g + (uLength - 1)*e;\r
+\r
+ // case TERMGAP_Half:\r
+ // return g/2 + (uLength - 1)*e;\r
+\r
+ // case TERMGAP_Ext:\r
+ // return uLength*e;\r
+ // }\r
+ // Quit("Bad termgap");\r
+ // }\r
+ //else\r
+ // return g + (uLength - 1)*e;\r
+ //return MINUS_INFINITY;\r
+ return g + (uLength - 1)*e;\r
+ }\r
+\r
+static SCORE GapPenalty(unsigned uLength, bool Term)\r
+ {\r
+ SCORE s1 = GapPenalty(uLength, Term, g_scoreGapOpen, g_scoreGapExtend);\r
+#if DOUBLE_AFFINE\r
+ SCORE s2 = GapPenalty(uLength, Term, g_scoreGapOpen2, g_scoreGapExtend2);\r
+ if (s1 > s2)\r
+ return s1;\r
+ return s2;\r
+#else\r
+ return s1;\r
+#endif\r
+ }\r
+\r
+static const MSA *g_ptrMSA1;\r
+static const MSA *g_ptrMSA2;\r
+static unsigned g_uSeqIndex1;\r
+static unsigned g_uSeqIndex2;\r
+\r
+static void LogGap(unsigned uStart, unsigned uEnd, unsigned uGapLength,\r
+ bool bNTerm, bool bCTerm)\r
+ {\r
+ Log("%16.16s ", "");\r
+ for (unsigned i = 0; i < uStart; ++i)\r
+ Log(" ");\r
+ unsigned uMyLength = 0;\r
+ for (unsigned i = uStart; i <= uEnd; ++i)\r
+ {\r
+ bool bGap1 = g_ptrMSA1->IsGap(g_uSeqIndex1, i);\r
+ bool bGap2 = g_ptrMSA2->IsGap(g_uSeqIndex2, i);\r
+ if (!bGap1 && !bGap2)\r
+ Quit("Error -- neither gapping");\r
+ if (bGap1 && bGap2)\r
+ Log(".");\r
+ else\r
+ {\r
+ ++uMyLength;\r
+ Log("-");\r
+ }\r
+ }\r
+ SCORE s = GapPenalty(uGapLength, bNTerm || bCTerm);\r
+ Log(" L=%d N%d C%d s=%.3g", uGapLength, bNTerm, bCTerm, s);\r
+ Log("\n");\r
+ if (uMyLength != uGapLength)\r
+ Quit("Lengths differ");\r
+\r
+ }\r
+\r
+static SCORE ScoreSeqPair(const MSA &msa1, unsigned uSeqIndex1,\r
+ const MSA &msa2, unsigned uSeqIndex2, SCORE *ptrLetters, SCORE *ptrGaps)\r
+ {\r
+ g_ptrMSA1 = &msa1;\r
+ g_ptrMSA2 = &msa2;\r
+ g_uSeqIndex1 = uSeqIndex1;\r
+ g_uSeqIndex2 = uSeqIndex2;\r
+\r
+ const unsigned uColCount = msa1.GetColCount();\r
+ const unsigned uColCount2 = msa2.GetColCount();\r
+ if (uColCount != uColCount2)\r
+ Quit("ScoreSeqPair, different lengths");\r
+\r
+#if TRACE\r
+ Log("ScoreSeqPair\n");\r
+ Log("%16.16s ", msa1.GetSeqName(uSeqIndex1));\r
+ for (unsigned i = 0; i < uColCount; ++i)\r
+ Log("%c", msa1.GetChar(uSeqIndex1, i));\r
+ Log("\n");\r
+ Log("%16.16s ", msa2.GetSeqName(uSeqIndex2));\r
+ for (unsigned i = 0; i < uColCount; ++i)\r
+ Log("%c", msa1.GetChar(uSeqIndex2, i));\r
+ Log("\n");\r
+#endif\r
+\r
+ SCORE scoreTotal = 0;\r
+\r
+// Substitution scores\r
+ unsigned uFirstLetter1 = uInsane;\r
+ unsigned uFirstLetter2 = uInsane;\r
+ unsigned uLastLetter1 = uInsane;\r
+ unsigned uLastLetter2 = uInsane;\r
+ for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
+ {\r
+ bool bGap1 = msa1.IsGap(uSeqIndex1, uColIndex);\r
+ bool bGap2 = msa2.IsGap(uSeqIndex2, uColIndex);\r
+ bool bWildcard1 = msa1.IsWildcard(uSeqIndex1, uColIndex);\r
+ bool bWildcard2 = msa2.IsWildcard(uSeqIndex2, uColIndex);\r
+\r
+ if (!bGap1)\r
+ {\r
+ if (uInsane == uFirstLetter1)\r
+ uFirstLetter1 = uColIndex;\r
+ uLastLetter1 = uColIndex;\r
+ }\r
+ if (!bGap2)\r
+ {\r
+ if (uInsane == uFirstLetter2)\r
+ uFirstLetter2 = uColIndex;\r
+ uLastLetter2 = uColIndex;\r
+ }\r
+\r
+ if (bGap1 || bGap2 || bWildcard1 || bWildcard2)\r
+ continue;\r
+\r
+ unsigned uLetter1 = msa1.GetLetter(uSeqIndex1, uColIndex);\r
+ unsigned uLetter2 = msa2.GetLetter(uSeqIndex2, uColIndex);\r
+\r
+ SCORE scoreMatch = (*g_ptrScoreMatrix)[uLetter1][uLetter2];\r
+ scoreTotal += scoreMatch;\r
+#if TRACE\r
+ Log("%c <-> %c = %7.1f %10.1f\n",\r
+ msa1.GetChar(uSeqIndex1, uColIndex),\r
+ msa2.GetChar(uSeqIndex2, uColIndex),\r
+ scoreMatch,\r
+ scoreTotal);\r
+#endif\r
+ }\r
+ \r
+ *ptrLetters = scoreTotal;\r
+\r
+// Gap penalties\r
+ unsigned uGapLength = uInsane;\r
+ unsigned uGapStartCol = uInsane;\r
+ bool bGapping1 = false;\r
+ bool bGapping2 = false;\r
+\r
+ for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
+ {\r
+ bool bGap1 = msa1.IsGap(uSeqIndex1, uColIndex);\r
+ bool bGap2 = msa2.IsGap(uSeqIndex2, uColIndex);\r
+\r
+ if (bGap1 && bGap2)\r
+ continue;\r
+\r
+ if (bGapping1)\r
+ {\r
+ if (bGap1)\r
+ ++uGapLength;\r
+ else\r
+ {\r
+ bGapping1 = false;\r
+ bool bNTerm = (uFirstLetter2 == uGapStartCol);\r
+ bool bCTerm = (uLastLetter2 + 1 == uColIndex);\r
+ SCORE scoreGap = GapPenalty(uGapLength, bNTerm || bCTerm);\r
+ scoreTotal += scoreGap;\r
+#if TRACE\r
+ LogGap(uGapStartCol, uColIndex - 1, uGapLength, bNTerm, bCTerm);\r
+ Log("GAP %7.1f %10.1f\n",\r
+ scoreGap,\r
+ scoreTotal);\r
+#endif\r
+ }\r
+ continue;\r
+ }\r
+ else\r
+ {\r
+ if (bGap1)\r
+ {\r
+ uGapStartCol = uColIndex;\r
+ bGapping1 = true;\r
+ uGapLength = 1;\r
+ continue;\r
+ }\r
+ }\r
+\r
+ if (bGapping2)\r
+ {\r
+ if (bGap2)\r
+ ++uGapLength;\r
+ else\r
+ {\r
+ bGapping2 = false;\r
+ bool bNTerm = (uFirstLetter1 == uGapStartCol);\r
+ bool bCTerm = (uLastLetter1 + 1 == uColIndex);\r
+ SCORE scoreGap = GapPenalty(uGapLength, bNTerm || bCTerm);\r
+ scoreTotal += scoreGap;\r
+#if TRACE\r
+ LogGap(uGapStartCol, uColIndex - 1, uGapLength, bNTerm, bCTerm);\r
+ Log("GAP %7.1f %10.1f\n",\r
+ scoreGap,\r
+ scoreTotal);\r
+#endif\r
+ }\r
+ }\r
+ else\r
+ {\r
+ if (bGap2)\r
+ {\r
+ uGapStartCol = uColIndex;\r
+ bGapping2 = true;\r
+ uGapLength = 1;\r
+ }\r
+ }\r
+ }\r
+\r
+ if (bGapping1 || bGapping2)\r
+ {\r
+ SCORE scoreGap = GapPenalty(uGapLength, true);\r
+ scoreTotal += scoreGap;\r
+#if TRACE\r
+ LogGap(uGapStartCol, uColCount - 1, uGapLength, false, true);\r
+ Log("GAP %7.1f %10.1f\n",\r
+ scoreGap,\r
+ scoreTotal);\r
+#endif\r
+ }\r
+ *ptrGaps = scoreTotal - *ptrLetters;\r
+ return scoreTotal;\r
+ }\r
+\r
+// The usual sum-of-pairs objective score: sum the score\r
+// of the alignment of each pair of sequences.\r
+SCORE ObjScoreDA(const MSA &msa, SCORE *ptrLetters, SCORE *ptrGaps)\r
+ {\r
+ const unsigned uSeqCount = msa.GetSeqCount();\r
+ SCORE scoreTotal = 0;\r
+ unsigned uPairCount = 0;\r
+#if TRACE\r
+ msa.LogMe();\r
+ Log(" Score Weight Weight Total\n");\r
+ Log("---------- ------ ------ ----------\n");\r
+#endif\r
+ SCORE TotalLetters = 0;\r
+ SCORE TotalGaps = 0;\r
+ for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount; ++uSeqIndex1)\r
+ {\r
+ const WEIGHT w1 = msa.GetSeqWeight(uSeqIndex1);\r
+ for (unsigned uSeqIndex2 = uSeqIndex1 + 1; uSeqIndex2 < uSeqCount; ++uSeqIndex2)\r
+ {\r
+ const WEIGHT w2 = msa.GetSeqWeight(uSeqIndex2);\r
+ const WEIGHT w = w1*w2;\r
+ SCORE Letters;\r
+ SCORE Gaps;\r
+ SCORE scorePair = ScoreSeqPair(msa, uSeqIndex1, msa, uSeqIndex2,\r
+ &Letters, &Gaps);\r
+ scoreTotal += w1*w2*scorePair;\r
+ TotalLetters += w1*w2*Letters;\r
+ TotalGaps += w1*w2*Gaps;\r
+ ++uPairCount;\r
+#if TRACE\r
+ Log("%10.2f %6.3f %6.3f %10.2f %d=%s %d=%s\n",\r
+ scorePair,\r
+ w1,\r
+ w2,\r
+ scorePair*w1*w2,\r
+ uSeqIndex1,\r
+ msa.GetSeqName(uSeqIndex1),\r
+ uSeqIndex2,\r
+ msa.GetSeqName(uSeqIndex2));\r
+#endif\r
+ }\r
+ }\r
+ *ptrLetters = TotalLetters;\r
+ *ptrGaps = TotalGaps;\r
+ return scoreTotal;\r
+ }\r
+\r
+#endif // DOUBLE_AFFINE\r
git://source.jalview.org / jabaws.git / blobdiff