Logo Search packages:      
Sourcecode: ugene version File versions

nwdasmall.cpp

/*****************************************************************
* Unipro UGENE - Integrated Bioinformatics Suite
* Copyright (C) 2008 Unipro, Russia (http://ugene.unipro.ru)
* All Rights Reserved
* 
*     This source code is distributed under the terms of the
*     GNU General Public License. See the files COPYING and LICENSE
*     for details.
*****************************************************************/

#include "muscle.h"
#include <math.h>
#include "pwpath.h"
#include "profile.h"
#include <stdio.h>

#if   DOUBLE_AFFINE

// NW double affine small memory, term gaps fully penalized
// (so up to caller to adjust in profile if desired).

#define     TRACE 0

#define MIN(x, y) ((x) < (y) ? (x) : (y))

#if   TRACE
extern bool g_bKeepSimpleDP;
extern SCORE *g_DPM;
extern SCORE *g_DPD;
extern SCORE *g_DPE;
extern SCORE *g_DPI;
extern SCORE *g_DPJ;
extern char *g_TBM;
extern char *g_TBD;
extern char *g_TBE;
extern char *g_TBI;
extern char *g_TBJ;
#endif

#if   TRACE
#define ALLOC_TRACE()                                             \
      const SCORE UNINIT = MINUS_INFINITY;                  \
      const size_t LM = uPrefixCountA*uPrefixCountB;  \
                                                                              \
      SCORE *DPM_ = new SCORE[LM];                          \
      SCORE *DPD_ = new SCORE[LM];                          \
      SCORE *DPE_ = new SCORE[LM];                          \
      SCORE *DPI_ = new SCORE[LM];                          \
      SCORE *DPJ_ = new SCORE[LM];                          \
                                                                              \
      char *TBM_ = new char[LM];                                  \
      char *TBD_ = new char[LM];                                  \
      char *TBE_ = new char[LM];                                  \
      char *TBI_ = new char[LM];                                  \
      char *TBJ_ = new char[LM];                                  \
                                                                              \
      memset(TBM_, '?', LM);                                      \
      memset(TBD_, '?', LM);                                      \
      memset(TBE_, '?', LM);                                      \
      memset(TBI_, '?', LM);                                      \
      memset(TBJ_, '?', LM);                                      \
                                                                              \
      for (unsigned i = 0; i <= uLengthA; ++i)        \
            for (unsigned j = 0; j <= uLengthB; ++j)  \
                  {                                                           \
                  DPM(i, j) = UNINIT;                                   \
                  DPD(i, j) = UNINIT;                                   \
                  DPE(i, j) = UNINIT;                                   \
                  DPI(i, j) = UNINIT;                                   \
                  DPJ(i, j) = UNINIT;                                   \
                  }
#else
#define ALLOC_TRACE()
#endif

#if   TRACE
#define SetDPM(i, j, x)       DPM(i, j) = x
#define SetDPD(i, j, x)       DPD(i, j) = x
#define SetDPE(i, j, x)       DPE(i, j) = x
#define SetDPI(i, j, x)       DPI(i, j) = x
#define SetDPJ(i, j, x)       DPJ(i, j) = x
#define SetTBM(i, j, x)       TBM(i, j) = x
#define SetTBD(i, j, x)       TBD(i, j) = x
#define SetTBE(i, j, x)       TBE(i, j) = x
#define SetTBI(i, j, x)       TBI(i, j) = x
#define SetTBJ(i, j, x)       TBJ(i, j) = x
#else
#define SetDPM(i, j, x)       /* empty  */
#define SetDPD(i, j, x)       /* empty  */
#define SetDPE(i, j, x)       /* empty  */
#define SetDPI(i, j, x)       /* empty  */
#define SetDPJ(i, j, x)       /* empty  */
#define SetTBM(i, j, x)       /* empty  */
#define SetTBD(i, j, x)       /* empty  */
#define SetTBE(i, j, x)       /* empty  */
#define SetTBI(i, j, x)       /* empty  */
#define SetTBJ(i, j, x)       /* empty  */
#endif

#define RECURSE_D(i, j)                   \
      {                                               \
      SCORE DD = DRow[j] + e;             \
      SCORE MD = MPrev[j] + PA[i-1].m_scoreGapOpen;\
      if (DD > MD)                              \
            {                                         \
            DRow[j] = DD;                       \
            SetTBD(i, j, 'D');                  \
            }                                         \
      else                                      \
            {                                         \
            DRow[j] = MD;                       \
            SetBitTBD(TB, i, j, 'M');     \
            SetTBD(i, j, 'M');                  \
            }                                         \
      SetDPD(i, j, DRow[j]);              \
      }

#define RECURSE_E(i, j)                   \
      {                                               \
      SCORE EE = ERow[j] + e2;            \
      SCORE ME = MPrev[j] + PA[i-1].m_scoreGapOpen2;\
      if (EE > ME)                              \
            {                                         \
            ERow[j] = EE;                       \
            SetTBE(i, j, 'E');                  \
            }                                         \
      else                                      \
            {                                         \
            ERow[j] = ME;                       \
            SetBitTBE(TB, i, j, 'M');     \
            SetTBE(i, j, 'M');                  \
            }                                         \
      SetDPE(i, j, ERow[j]);              \
      }

#define RECURSE_D_ATerm(j)    RECURSE_D(uLengthA, j)
#define RECURSE_E_ATerm(j)    RECURSE_E(uLengthA, j)

#define RECURSE_D_BTerm(j)    RECURSE_D(i, uLengthB)
#define RECURSE_E_BTerm(j)    RECURSE_E(i, uLengthB)

#define RECURSE_I(i, j)                   \
      {                                               \
      Iij += e;                                 \
      SCORE MI = MCurr[j-1] + PB[j-1].m_scoreGapOpen;\
      if (MI >= Iij)                            \
            {                                         \
            Iij = MI;                           \
            SetBitTBI(TB, i, j, 'M');     \
            SetTBI(i, j, 'M');                  \
            }                                         \
      else                                      \
            SetTBI(i, j, 'I');                  \
      SetDPI(i, j, Iij);                        \
      }

#define RECURSE_J(i, j)                   \
      {                                               \
      Jij += e2;                                \
      SCORE MJ = MCurr[j-1] + PB[j-1].m_scoreGapOpen2;\
      if (MJ >= Jij)                            \
            {                                         \
            Jij = MJ;                           \
            SetBitTBJ(TB, i, j, 'M');     \
            SetTBJ(i, j, 'M');                  \
            }                                         \
      else                                      \
            SetTBJ(i, j, 'I');                  \
      SetDPJ(i, j, Jij);                        \
      }

#define RECURSE_I_ATerm(j)    RECURSE_I(uLengthA, j)
#define RECURSE_J_ATerm(j)    RECURSE_J(uLengthA, j)

#define RECURSE_I_BTerm(j)    RECURSE_I(i, uLengthB)
#define RECURSE_J_BTerm(j)    RECURSE_J(i, uLengthB)

#define RECURSE_M(i, j)                                                 \
      {                                                                             \
      SCORE Best = MCurr[j];  /*  MM  */                          \
      SetTBM(i+1, j+1, 'M');                                            \
      SetBitTBM(TB, i+1, j+1, 'M');                               \
                                                                                    \
      SCORE DM = DRow[j] + PA[i-1].m_scoreGapClose;         \
      if (DM > Best)                                                          \
            {                                                                       \
            Best = DM;                                                        \
            SetTBM(i+1, j+1, 'D');                                      \
            SetBitTBM(TB, i+1, j+1, 'D');                         \
            }                                                                       \
                                                                                    \
      SCORE EM = ERow[j] + PA[i-1].m_scoreGapClose2;        \
      if (EM > Best)                                                          \
            {                                                                       \
            Best = EM;                                                        \
            SetTBM(i+1, j+1, 'E');                                      \
            SetBitTBM(TB, i+1, j+1, 'E');                         \
            }                                                                       \
                                                                                    \
      SCORE IM = Iij + PB[j-1].m_scoreGapClose;             \
      if (IM > Best)                                                          \
            {                                                                       \
            Best = IM;                                                        \
            SetTBM(i+1, j+1, 'I');                                      \
            SetBitTBM(TB, i+1, j+1, 'I');                         \
            }                                                                       \
                                                                                    \
      SCORE JM = Jij + PB[j-1].m_scoreGapClose2;                  \
      if (JM > Best)                                                          \
            {                                                                       \
            Best = JM;                                                        \
            SetTBM(i+1, j+1, 'J');                                      \
            SetBitTBM(TB, i+1, j+1, 'J');                         \
            }                                                                       \
      MNext[j+1] += Best;                                                     \
      SetDPM(i+1, j+1, MNext[j+1]);                               \
      }

#if   TRACE
static bool LocalEq(BASETYPE b1, BASETYPE b2)
      {
      if (b1 < -100000 && b2 < -100000)
            return true;
      double diff = fabs(b1 - b2);
      if (diff < 0.0001)
            return true;
      double sum = fabs(b1) + fabs(b2);
      return diff/sum < 0.005;
      }

static char Get_M_Char(char Bits)
      {
      switch (Bits & BIT_xM)
            {
      case BIT_MM:
            return 'M';
      case BIT_DM:
            return 'D';
      case BIT_EM:
            return 'E';
      case BIT_IM:
            return 'I';
      case BIT_JM:
            return 'J';
            }
      Quit("Huh?");
      return '?';
      }

static char Get_D_Char(char Bits)
      {
      return (Bits & BIT_xD) ? 'M' : 'D';
      }

static char Get_E_Char(char Bits)
      {
      return (Bits & BIT_xE) ? 'M' : 'E';
      }

static char Get_I_Char(char Bits)
      {
      return (Bits & BIT_xI) ? 'M' : 'I';
      }

static char Get_J_Char(char Bits)
      {
      return (Bits & BIT_xJ) ? 'M' : 'J';
      }

static bool DPEq(char c, SCORE *g_DP, SCORE *DPD_,
  unsigned uPrefixCountA, unsigned uPrefixCountB)
      {
      if (0 == g_DP)
            {
            Log("***DPDIFF*** DP%c=NULL\n", c);
            return true;
            }

      SCORE *DPM_ = g_DP;
      for (unsigned i = 0; i < uPrefixCountA; ++i)
            for (unsigned j = 0; j < uPrefixCountB; ++j)
                  if (!LocalEq(DPM(i, j), DPD(i, j)))
                        {
                        Log("***DPDIFF*** DP%c(%d, %d) Simple = %.2g, Small = %.2g\n",
                          c, i, j, DPM(i, j), DPD(i, j));
                        return false;
                        }
      return true;
      }

static bool CompareTB(char **TB, char *TBM_, char *TBD_, char *TBE_, char *TBI_, char *TBJ_,
  unsigned uPrefixCountA, unsigned uPrefixCountB)
      {
      if (!g_bKeepSimpleDP)
            return true;
      SCORE *DPM_ = g_DPM;
      bool Eq = true;
      for (unsigned i = 0; i < uPrefixCountA; ++i)
            for (unsigned j = 0; j < uPrefixCountB; ++j)
                  {
                  char c1 = TBM(i, j);
                  char c2 = Get_M_Char(TB[i][j]);
                  if (c1 != '?' && c1 != c2 && DPM(i, j) > -100000)
                        {
                        Log("TBM(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
                        Eq = false;
                        goto D;
                        }
                  }

D:
      SCORE *DPD_ = g_DPD;
      for (unsigned i = 0; i < uPrefixCountA; ++i)
            for (unsigned j = 0; j < uPrefixCountB; ++j)
                  {
                  char c1 = TBD(i, j);
                  char c2 = Get_D_Char(TB[i][j]);
                  if (c1 != '?' && c1 != c2 && DPD(i, j) > -100000)
                        {
                        Log("TBD(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
                        Eq = false;
                        goto E;
                        }
                  }
E:
      SCORE *DPE_ = g_DPE;
      if (0 == TBE_)
            goto I;
      for (unsigned i = 0; i < uPrefixCountA; ++i)
            for (unsigned j = 0; j < uPrefixCountB; ++j)
                  {
                  char c1 = TBE(i, j);
                  char c2 = Get_E_Char(TB[i][j]);
                  if (c1 != '?' && c1 != c2 && DPE(i, j) > -100000)
                        {
                        Log("TBE(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
                        Eq = false;
                        goto I;
                        }
                  }
I:
      SCORE *DPI_ = g_DPI;
      for (unsigned i = 0; i < uPrefixCountA; ++i)
            for (unsigned j = 0; j < uPrefixCountB; ++j)
                  {
                  char c1 = TBI(i, j);
                  char c2 = Get_I_Char(TB[i][j]);
                  if (c1 != '?' && c1 != c2 && DPI(i, j) > -100000)
                        {
                        Log("TBI(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
                        Eq = false;
                        goto J;
                        }
                  }
J:
      SCORE *DPJ_ = g_DPJ;
      if (0 == DPJ_)
            goto Done;
      for (unsigned i = 0; i < uPrefixCountA; ++i)
            for (unsigned j = 0; j < uPrefixCountB; ++j)
                  {
                  char c1 = TBJ(i, j);
                  char c2 = Get_J_Char(TB[i][j]);
                  if (c1 != '?' && c1 != c2 && DPJ(i, j) > -100000)
                        {
                        Log("TBJ(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
                        Eq = false;
                        goto Done;
                        }
                  }
Done:
      if (Eq)
            Log("TB success\n");
      return Eq;
      }

static const char *LocalScoreToStr(SCORE s)
      {
      static char str[16];
      if (s < -100000)
            return "     *";
      sprintf(str, "%6.1f", s);
      return str;
      }

static void LogDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,
  unsigned uPrefixCountA, unsigned uPrefixCountB)
      {
      Log("        ");
      for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
            {
            char c = ' ';
            if (uPrefixLengthB > 0)
                  c = ConsensusChar(PB[uPrefixLengthB - 1]);
            Log(" %4u:%c", uPrefixLengthB, c);
            }
      Log("\n");
      for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
            {
            char c = ' ';
            if (uPrefixLengthA > 0)
                  c = ConsensusChar(PA[uPrefixLengthA - 1]);
            Log("%4u:%c  ", uPrefixLengthA, c);
            for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
                  Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));
            Log("\n");
            }
      }

static void LogBitTB(char **TB, const ProfPos *PA, const ProfPos *PB,
  unsigned uPrefixCountA, unsigned uPrefixCountB)
      {
      Log("        ");
      for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
            {
            char c = ' ';
            if (uPrefixLengthB > 0)
                  c = ConsensusChar(PB[uPrefixLengthB - 1]);
            Log(" %4u:%c", uPrefixLengthB, c);
            }
      Log("\n");
      Log("Bit TBM:\n");
      for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
            {
            char c = ' ';
            if (uPrefixLengthA > 0)
                  c = ConsensusChar(PA[uPrefixLengthA - 1]);
            Log("%4u:%c  ", uPrefixLengthA, c);
            for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
                  {
                  char c = Get_M_Char(TB[uPrefixLengthA][uPrefixLengthB]);
                  Log(" %6c", c);
                  }
            Log("\n");
            }

      Log("\n");
      Log("Bit TBD:\n");
      for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
            {
            char c = ' ';
            if (uPrefixLengthA > 0)
                  c = ConsensusChar(PA[uPrefixLengthA - 1]);
            Log("%4u:%c  ", uPrefixLengthA, c);
            for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
                  {
                  char c = Get_D_Char(TB[uPrefixLengthA][uPrefixLengthB]);
                  Log(" %6c", c);
                  }
            Log("\n");
            }

      Log("\n");
      Log("Bit TBE:\n");
      for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
            {
            char c = ' ';
            if (uPrefixLengthA > 0)
                  c = ConsensusChar(PA[uPrefixLengthA - 1]);
            Log("%4u:%c  ", uPrefixLengthA, c);
            for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
                  {
                  char c = Get_E_Char(TB[uPrefixLengthA][uPrefixLengthB]);
                  Log(" %6c", c);
                  }
            Log("\n");
            }

      Log("\n");
      Log("Bit TBI:\n");
      for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
            {
            char c = ' ';
            if (uPrefixLengthA > 0)
                  c = ConsensusChar(PA[uPrefixLengthA - 1]);
            Log("%4u:%c  ", uPrefixLengthA, c);
            for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
                  {
                  char c = Get_I_Char(TB[uPrefixLengthA][uPrefixLengthB]);
                  Log(" %6c", c);
                  }
            Log("\n");
            }

      Log("\n");
      Log("Bit TBJ:\n");
      for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
            {
            char c = ' ';
            if (uPrefixLengthA > 0)
                  c = ConsensusChar(PA[uPrefixLengthA - 1]);
            Log("%4u:%c  ", uPrefixLengthA, c);
            for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
                  {
                  char c = Get_J_Char(TB[uPrefixLengthA][uPrefixLengthB]);
                  Log(" %6c", c);
                  }
            Log("\n");
            }
      }

static void ListTB(char *TBM_, const ProfPos *PA, const ProfPos *PB,
  unsigned uPrefixCountA, unsigned uPrefixCountB)
      {
      Log("        ");
      for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
            {
            char c = ' ';
            if (uPrefixLengthB > 0)
                  c = ConsensusChar(PB[uPrefixLengthB - 1]);
            Log(" %4u:%c", uPrefixLengthB, c);
            }
      Log("\n");
      for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
            {
            char c = ' ';
            if (uPrefixLengthA > 0)
                  c = ConsensusChar(PA[uPrefixLengthA - 1]);
            Log("%4u:%c  ", uPrefixLengthA, c);
            for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
                  {
                  char c = TBM(uPrefixLengthA, uPrefixLengthB);
                  Log(" %6c", c);
                  }
            Log("\n");
            }
      }

static const char *BitsToStr(char Bits)
      {
      static char Str[32];

      sprintf(Str, "%cM %cD %cE %cI %cJ",
        Get_M_Char(Bits),
        Get_D_Char(Bits),
        Get_E_Char(Bits),
        Get_I_Char(Bits),
        Get_J_Char(Bits));
      }
#endif      // TRACE

static inline void SetBitTBM(char **TB, unsigned i, unsigned j, char c)
      {
      char Bit;
      switch (c)
            {
      case 'M':
            Bit = BIT_MM;
            break;
      case 'D':
            Bit = BIT_DM;
            break;
#if   DOUBLE_AFFINE
      case 'E':
            Bit = BIT_EM;
            break;
      case 'I':
            Bit = BIT_IM;
            break;
      case 'J':
            Bit = BIT_JM;
            break;
#endif
      default:
            Quit("Huh?!");
            }
      TB[i][j] &= ~BIT_xM;
      TB[i][j] |= Bit;
      }

static inline void SetBitTBD(char **TB, unsigned i, unsigned j, char c)
      {
      char Bit;
      switch (c)
            {
      case 'M':
            Bit = BIT_MD;
            break;
      case 'D':
            Bit = BIT_DD;
            break;
      default:
            Quit("Huh?!");
            }
      TB[i][j] &= ~BIT_xD;
      TB[i][j] |= Bit;
      }

static inline void SetBitTBI(char **TB, unsigned i, unsigned j, char c)
      {
      char Bit;
      switch (c)
            {
      case 'M':
            Bit = BIT_MI;
            break;
      case 'I':
            Bit = BIT_II;
            break;
      default:
            Quit("Huh?!");
            }
      TB[i][j] &= ~BIT_xI;
      TB[i][j] |= Bit;
      }

#if   DOUBLE_AFFINE
static inline void SetBitTBE(char **TB, unsigned i, unsigned j, char c)
      {
      char Bit;
      switch (c)
            {
      case 'M':
            Bit = BIT_ME;
            break;
      case 'E':
            Bit = BIT_EE;
            break;
      default:
            Quit("Huh?!");
            }
      TB[i][j] &= ~BIT_xE;
      TB[i][j] |= Bit;
      }

static inline void SetBitTBJ(char **TB, unsigned i, unsigned j, char c)
      {
      char Bit;
      switch (c)
            {
      case 'M':
            Bit = BIT_MJ;
            break;
      case 'J':
            Bit = BIT_JJ;
            break;
      default:
            Quit("Huh?!");
            }
      TB[i][j] &= ~BIT_xJ;
      TB[i][j] |= Bit;
      }
#endif

#if   TRACE
#define LogMatrices()                                                               \
      {                                                                                         \
      Log("Bit DPM:\n");                                                                  \
      LogDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);                \
      Log("Bit DPD:\n");                                                                  \
      LogDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);                \
      Log("Bit DPE:\n");                                                                  \
      LogDP(DPE_, PA, PB, uPrefixCountA, uPrefixCountB);                \
      Log("Bit DPI:\n");                                                                  \
      LogDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);                \
      Log("Bit DPJ:\n");                                                                  \
      LogDP(DPJ_, PA, PB, uPrefixCountA, uPrefixCountB);                \
      Log("Bit TB:\n");                                                             \
      LogBitTB(TB, PA, PB, uPrefixCountA, uPrefixCountB);               \
      bool Same;                                                                          \
      Same = DPEq('M', g_DPM, DPM_, uPrefixCountA, uPrefixCountB);\
      if (Same)                                                                           \
            Log("DPM success\n");                                                   \
      Same = DPEq('D', g_DPD, DPD_, uPrefixCountA, uPrefixCountB);\
      if (Same)                                                                           \
            Log("DPD success\n");                                                   \
      Same = DPEq('E', g_DPE, DPE_, uPrefixCountA, uPrefixCountB);\
      if (Same)                                                                           \
            Log("DPE success\n");                                                   \
      Same = DPEq('I', g_DPI, DPI_, uPrefixCountA, uPrefixCountB);\
      if (Same)                                                                           \
            Log("DPI success\n");                                                   \
      Same = DPEq('J', g_DPJ, DPJ_, uPrefixCountA, uPrefixCountB);\
      if (Same)                                                                           \
            Log("DPJ success\n");                                                   \
      CompareTB(TB, g_TBM, g_TBD, g_TBE, g_TBI, g_TBJ, uPrefixCountA, uPrefixCountB);\
      }
#else
#define LogMatrices()   /* empty */
#endif

SCORE NWDASmall(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
  unsigned uLengthB, PWPath &Path)
      {
      assert(uLengthB > 0 && uLengthA > 0);

      ProfPos *pa0 = (ProfPos *) PA;
      ProfPos *pb0 = (ProfPos *) PB;
      ProfPos *paa = (ProfPos *) (PA + uLengthA - 1);
      ProfPos *pbb = (ProfPos *) (PB + uLengthB - 1);

      pa0->m_scoreGapOpen *= -1;
      pb0->m_scoreGapOpen *= -1;

      paa->m_scoreGapClose *= -1;
      pbb->m_scoreGapClose *= -1;

      pa0->m_scoreGapOpen2 *= -1;
      pb0->m_scoreGapOpen2 *= -1;
      paa->m_scoreGapClose2 *= -1;
      pbb->m_scoreGapClose2 *= -1;

      const unsigned uPrefixCountA = uLengthA + 1;
      const unsigned uPrefixCountB = uLengthB + 1;
      const SCORE e = g_scoreGapExtend;

      const SCORE e2 = g_scoreGapExtend2;
      const SCORE min_e = MIN(g_scoreGapExtend, g_scoreGapExtend2);

      ALLOC_TRACE()

      SCORE *MCurr = new SCORE[uPrefixCountB];
      SCORE *MNext = new SCORE[uPrefixCountB];
      SCORE *MPrev = new SCORE[uPrefixCountB];
      SCORE *DRow = new SCORE[uPrefixCountB];
      SCORE *ERow = new SCORE[uPrefixCountB];

      char **TB = new char *[uPrefixCountA];
      for (unsigned i = 0; i < uPrefixCountA; ++i)
            {
            TB[i] = new char [uPrefixCountB];
            memset(TB[i], 0, uPrefixCountB);
            }

      SCORE Iij = MINUS_INFINITY;
      SetDPI(0, 0, Iij);

      SCORE Jij = MINUS_INFINITY;
      SetDPJ(0, 0, Jij);

      Iij = PB[0].m_scoreGapOpen;
      SetDPI(0, 1, Iij);

      Jij = PB[0].m_scoreGapOpen2;
      SetDPJ(0, 1, Jij);

      for (unsigned j = 2; j <= uLengthB; ++j)
            {
            Iij += e;
            Jij += e2;

            SetDPI(0, j, Iij);
            SetDPJ(0, j, Jij);

            SetTBI(0, j, 'I');
            SetTBJ(0, j, 'J');
            }

      for (unsigned j = 0; j <= uLengthB; ++j)
            {
            DRow[j] = MINUS_INFINITY;
            ERow[j] = MINUS_INFINITY;

            SetDPD(0, j, DRow[j]);
            SetDPE(0, j, ERow[j]);

            SetTBD(0, j, 'D');
            SetTBE(0, j, 'E');
            }

      MPrev[0] = 0;
      SetDPM(0, 0, MPrev[0]);
      for (unsigned j = 1; j <= uLengthB; ++j)
            {
            MPrev[j] = MINUS_INFINITY;
            SetDPM(0, j, MPrev[j]);
            }

      MCurr[0] = MINUS_INFINITY;
      SetDPM(1, 0, MCurr[0]);

      MCurr[1] = ScoreProfPos2(PA[0], PB[0]);
      SetDPM(1, 1, MCurr[1]);
      SetBitTBM(TB, 1, 1, 'M');
      SetTBM(1, 1, 'M');

      for (unsigned j = 2; j <= uLengthB; ++j)
            {
            SCORE M = ScoreProfPos2(PA[0], PB[j-1]) + PB[0].m_scoreGapOpen +
              (j - 2)*e + PB[j-2].m_scoreGapClose;
            SCORE M2 = ScoreProfPos2(PA[0], PB[j-1]) + PB[0].m_scoreGapOpen2 +
              (j - 2)*e2 + PB[j-2].m_scoreGapClose2;
            
            if (M >= M2)
                  {
                  MCurr[j] = M;
                  SetBitTBM(TB, 1, j, 'I');
                  SetTBM(1, j, 'I');
                  }
            else
                  {
                  MCurr[j] = M2;
                  SetBitTBM(TB, 1, j, 'J');
                  SetTBM(1, j, 'J');
                  }
            SetDPM(1, j, MCurr[j]);
            }

// Main DP loop
      for (unsigned i = 1; i < uLengthA; ++i)
            {
            Iij = MINUS_INFINITY;
            Jij = MINUS_INFINITY;
            SetDPI(i, 0, Iij);
            SetDPJ(i, 0, Jij);

            DRow[0] = PA[0].m_scoreGapOpen + (i - 1)*e;
            ERow[0] = PA[0].m_scoreGapOpen2 + (i - 1)*e2;
            SetDPD(i, 0, DRow[0]);
            SetDPE(i, 0, ERow[0]);

            MCurr[0] = MINUS_INFINITY; 
            if (i == 1)
                  {
                  MCurr[1] = ScoreProfPos2(PA[0], PB[0]);
                  SetBitTBM(TB, i, 1, 'M');
                  SetTBM(i, 1, 'M');
                  }
            else
                  {
                  SCORE M = ScoreProfPos2(PA[i-1], PB[0]) + PA[0].m_scoreGapOpen +
                    (i - 2)*e + PA[i-2].m_scoreGapClose;
                  SCORE M2 = ScoreProfPos2(PA[i-1], PB[0]) + PA[0].m_scoreGapOpen2 +
                    (i - 2)*e2 + PA[i-2].m_scoreGapClose2;
                  if (M >= M2)
                        {
                        MCurr[1] = M;
                        SetBitTBM(TB, i, 1, 'D');
                        SetTBM(i, 1, 'D');
                        }
                  else
                        {
                        MCurr[1] = M2;
                        SetBitTBM(TB, i, 1, 'E');
                        SetTBM(i, 1, 'E');
                        }
                  }
            SetDPM(i, 0, MCurr[0]);
            SetDPM(i, 1, MCurr[1]);

            for (unsigned j = 1; j < uLengthB; ++j)
                  MNext[j+1] = ScoreProfPos2(PA[i], PB[j]);

            for (unsigned j = 1; j < uLengthB; ++j)
                  {
                  RECURSE_D(i, j)
                  RECURSE_E(i, j)
                  RECURSE_I(i, j)
                  RECURSE_J(i, j)
                  RECURSE_M(i, j)
                  }
      // Special case for j=uLengthB
            RECURSE_D_BTerm(i)
            RECURSE_E_BTerm(i)
            RECURSE_I_BTerm(i)
            RECURSE_J_BTerm(i)

      // Prev := Curr, Curr := Next, Next := Prev
            Rotate(MPrev, MCurr, MNext);
            }

// Special case for i=uLengthA
      MCurr[0] = MINUS_INFINITY;
      SCORE M = ScoreProfPos2(PA[uLengthA-1], PB[0]) + (uLengthA - 2)*e +
        PA[0].m_scoreGapOpen + PA[uLengthA-2].m_scoreGapClose;
      SCORE M2 = ScoreProfPos2(PA[uLengthA-1], PB[0]) + (uLengthA - 2)*e +
        PA[0].m_scoreGapOpen + PA[uLengthA-2].m_scoreGapClose;
      if (M >= M2)
            {
            MCurr[1] = M;
            SetBitTBM(TB, uLengthA, 1, 'D');
            SetTBM(uLengthA, 1, 'D');
            }
      else
            {
            MCurr[1] = M2;
            SetBitTBM(TB, uLengthA, 1, 'E');
            SetTBM(uLengthA, 1, 'D');
            }
      SetDPM(uLengthA, 0, MCurr[0]);
      SetDPM(uLengthA, 1, MCurr[1]);

      DRow[0] = MINUS_INFINITY;
      ERow[0] = MINUS_INFINITY;

      SetDPD(uLengthA, 0, DRow[0]);
      SetDPE(uLengthA, 0, ERow[0]);

      for (unsigned j = 1; j <= uLengthB; ++j)
            {
            RECURSE_D_ATerm(j);
            RECURSE_E_ATerm(j);
            }

      Iij = MINUS_INFINITY;
      Jij = MINUS_INFINITY;

      for (unsigned j = 1; j <= uLengthB; ++j)
            {
            RECURSE_I_ATerm(j)
            RECURSE_J_ATerm(j)
            }

      LogMatrices();

      SCORE MAB = MCurr[uLengthB];
      SCORE DAB = DRow[uLengthB] + PA[uLengthA-1].m_scoreGapClose;
      SCORE EAB = ERow[uLengthB] + PA[uLengthA-1].m_scoreGapClose2;
      SCORE IAB = Iij + PB[uLengthB-1].m_scoreGapClose;
      SCORE JAB = Jij + PB[uLengthB-1].m_scoreGapClose2;

      SCORE Score = MAB;
      char cEdgeType = 'M';
      if (DAB > Score)
            {
            Score = DAB;
            cEdgeType = 'D';
            }
      if (EAB > Score)
            {
            Score = EAB;
            cEdgeType = 'E';
            }
      if (IAB > Score)
            {
            Score = IAB;
            cEdgeType = 'I';
            }
      if (JAB > Score)
            {
            Score = JAB;
            cEdgeType = 'J';
            }

#if TRACE
      Log("    Small: MAB=%.4g DAB=%.4g EAB=%.4g IAB=%.4g JAB=%.4g best=%c\n",
        MAB, DAB, EAB, IAB, JAB, cEdgeType);
#endif

      BitTraceBack(TB, uLengthA, uLengthB, cEdgeType, Path);

#if   DBEUG
      Path.Validate();
#endif

      delete[] MCurr;
      delete[] MNext;
      delete[] MPrev;
      delete[] DRow;
      delete[] ERow;
      for (unsigned i = 0; i < uPrefixCountA; ++i)
            delete[] TB[i];
      delete[] TB;

      return 0;
      }
#endif      // DOUBLE_AFFINE

Generated by  Doxygen 1.6.0   Back to index