开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
mm_aligner.cpp
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：29k
源码类别：

生物技术

开发平台：

C/C++

mm_aligner.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: mm_aligner.cpp,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 18:04:46  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.20
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: mm_aligner.cpp,v 1000.1 2004/06/01 18:04:46 gouriano Exp $
  11.  * ===========================================================================
  12.  *
  13.  *                            PUBLIC DOMAIN NOTICE                          
  14.  *               National Center for Biotechnology Information
  15.  *                                                                          
  16.  *  This software/database is a "United States Government Work" under the   
  17.  *  terms of the United States Copyright Act.  It was written as part of    
  18.  *  the author's official duties as a United States Government employee and 
  19.  *  thus cannot be copyrighted.  This software/database is freely available 
  20.  *  to the public for use. The National Library of Medicine and the U.S.    
  21.  *  Government have not placed any restriction on its use or reproduction.  
  22.  *                                                                          
  23.  *  Although all reasonable efforts have been taken to ensure the accuracy  
  24.  *  and reliability of the software and data, the NLM and the U.S.           
  25.  *  Government do not and cannot warrant the performance or results that    
  26.  *  may be obtained by using this software or data. The NLM and the U.S.    
  27.  *  Government disclaim all warranties, express or implied, including       
  28.  *  warranties of performance, merchantability or fitness for any particular
  29.  *  purpose.                                                                
  30.  *                                                                          
  31.  *  Please cite the author in any work or product based on this material.   
  32.  *
  33.  * ===========================================================================
  34.  *
  35.  * Authors:  Yuri Kapustin
  36.  *
  37.  * File Description:  CMMAligner implementation
  38.  *                   
  39.  * ===========================================================================
  40.  *
  41.  */
  43. #include <ncbi_pch.hpp>
  44. #include "mm_aligner_threads.hpp"
  46. #include <corelib/ncbimtx.hpp>
  47. #include <corelib/ncbi_system.hpp>
  48. #include <algo/align/align_exception.hpp>
  50. BEGIN_NCBI_SCOPE
  52. CMMAligner::CMMAligner():
  53.     m_mt(false),
  54.     m_maxthreads(1)
  55. {
  56. }
  59. CMMAligner::CMMAligner( const char* seq1, size_t len1,
  60.                         const char* seq2, size_t len2,
  61.                         const SNCBIPackedScoreMatrix* scoremat )
  62.   : CNWAligner(seq1, len1, seq2, len2, scoremat),
  63.     m_mt(false),
  64.     m_maxthreads(1)
  65. {
  66. }
  69. void CMMAligner::EnableMultipleThreads(bool allow)
  70. {
  71.     m_maxthreads = (m_mt = allow)? GetCpuCount(): 1;
  72. }
  75. CNWAligner::TScore CMMAligner::x_Run()
  76. {
  77.     m_terminate = false;
  78.     if(m_prg_callback) {
  79.         m_prg_info.m_iter_total = 2*m_SeqLen1*m_SeqLen2;
  80.         m_prg_info.m_iter_done = 0;
  81.         m_terminate = m_prg_callback(&m_prg_info);
  82.     }
  84.     if(m_terminate) {
  85.         return m_score = 0;
  86.     }
  87.     
  88.     m_score = kMin_Int;
  89.     m_TransList.clear();
  90.     m_TransList.push_back(eTS_None);
  92.     SCoordRect m (0, 0, m_SeqLen1 - 1, m_SeqLen2 - 1);
  93.     x_DoSubmatrix(m, m_TransList.end(), false, false); // top-level call
  95.     if(m_terminate) {
  96.         return m_score = 0;
  97.     }
  99.     // reverse_copy not supported by some compilers
  100.     list<ETranscriptSymbol>::const_iterator ii = m_TransList.begin();
  101.     size_t nsize = m_TransList.size() - 1;
  102.     m_Transcript.clear();
  103.     m_Transcript.resize(nsize);
  104.     for(size_t k = 1; k <= nsize; ++k)
  105.         m_Transcript[nsize - k] = *++ii;
  107.     return m_score;
  108. }
  111. // trace bit coding (four bits per value): D E Ec Fc
  112. // D:  1 if diagonal; 0 - otherwise
  113. // E:  1 if space in 1st sequence; 0 if space in 2nd sequence
  114. // Ec: 1 if gap in 1st sequence was extended; 0 if it is was opened
  115. // Fc: 1 if gap in 2nd sequence was extended; 0 if it is was opened
  116. //
  118. const unsigned char kMaskFc  = 0x0001;
  119. const unsigned char kMaskEc  = 0x0002;
  120. const unsigned char kMaskE   = 0x0004;
  121. const unsigned char kMaskD   = 0x0008;
  123. DEFINE_STATIC_FAST_MUTEX (masterlist_mutex);
  125. void CMMAligner::x_DoSubmatrix( const SCoordRect& submatr,
  126.     list<ETranscriptSymbol>::iterator translist_pos,
  127.     bool left_top, bool right_bottom )
  128. {
  130.     if( m_terminate ) {
  131.         return;
  132.     }
  134.     const int dimI = submatr.i2 - submatr.i1 + 1;
  135.     const int dimJ = submatr.j2 - submatr.j1 + 1;
  136.     if(dimI < 1 || dimJ < 1) return;
  138.     bool top_level = submatr.i1 == 0 && submatr.j1 == 0 &&
  139.         submatr.i2 == m_SeqLen1-1 && submatr.j2 == m_SeqLen2-1;
  142.     if(dimI < 3 || dimJ < 3) { // terminal case
  143.         CFastMutexGuard guard (masterlist_mutex);
  144.         list<ETranscriptSymbol> lts;
  145.         TScore score = x_RunTerm(submatr, left_top, right_bottom, lts);
  146.         if(top_level) m_score = score;
  147.         m_TransList.splice(translist_pos, lts);
  148.         return;
  149.     }
  151.     const size_t I = submatr.i1 + dimI / 2;
  152.     const size_t dim = dimJ + 1;
  154.     // run Needleman-Wunsch without storing full traceback matrix
  155.     SCoordRect rtop (submatr.i1, submatr.j1, I, submatr.j2);
  156.     vector<TScore> vEtop (dim), vFtop (dim), vGtop (dim);
  157.     vector<unsigned char> trace_top (dim);
  159.     SCoordRect rbtm (I + 1, submatr.j1, submatr.i2 , submatr.j2);
  160.     vector<TScore> vEbtm (dim), vFbtm (dim), vGbtm (dim);
  161.     vector<unsigned char> trace_btm (dim);
  162.     
  163.     if( m_mt && m_maxthreads > 1 && MM_RequestNewThread(m_maxthreads) ) {
  164.         CThreadRunOnTop* thr2 = new CThreadRunOnTop ( this, &rtop,
  165.                                             &vEtop, &vFtop, &vGtop,
  166.                                             &trace_top, left_top );
  167.         thr2->Run();
  168.         x_RunBtm(rbtm, vEbtm, vFbtm, vGbtm, trace_btm, right_bottom);
  169.         thr2->Join(0);
  170.     }
  171.     else {
  172.         x_RunTop(rtop, vEtop, vFtop, vGtop, trace_top, left_top);
  173.         x_RunBtm(rbtm, vEbtm, vFbtm, vGbtm, trace_btm, right_bottom);
  174.     }
  176.     if( m_terminate ) {
  177.         return;
  178.     }
  180.     // locate the transition point
  181.     size_t trans_pos;
  182.     ETransitionType trans_type;
  183.     TScore score1 = x_FindBestJ ( vEtop, vFtop, vGtop, vEbtm, vFbtm, vGbtm,
  184.                                   trans_pos, trans_type );
  185.     if(top_level)
  186.         m_score = score1;
  188.     // modify traces at the transition point if applicable
  189.     if(trans_type == eII) {
  190.         unsigned char mask_top = kMaskE;
  191.         if(trace_top[trans_pos] & kMaskEc)
  192.             mask_top |= kMaskEc;
  193.         trace_top[trans_pos] = mask_top;
  194.         trace_btm[trans_pos] = kMaskE | kMaskEc;
  195.     }
  196.     if(trans_type == eDD) {
  197.         trace_top[trans_pos] = (trace_top[trans_pos] & kMaskFc)? kMaskFc: 0;
  198.         trace_btm[trans_pos] = kMaskFc;
  199.     }
  201.     // extend a subpath to both directions
  202.     vector<unsigned char>::const_iterator trace_it_top = 
  203.         trace_top.begin() + trans_pos;
  204.     list<ETranscriptSymbol> subpath_left;
  205.     size_t steps_left = x_ExtendSubpath(trace_it_top, false, subpath_left);
  207.     vector<unsigned char>::const_iterator trace_it_btm = 
  208.         trace_btm.begin() + trans_pos;
  209.     list<ETranscriptSymbol> subpath_right;
  210.     size_t steps_right = x_ExtendSubpath(trace_it_btm, true, subpath_right);
  212.     // check if we reached left or top line
  213.     bool bNoLT = false;
  214.     int nc0 = trans_pos - steps_left;
  215.     if(nc0 < 0) {
  216.         NCBI_THROW(CAlgoAlignException, eInternal,
  217.                    "Assertion: LT underflow");
  218.     }
  219.     bool bLeft = nc0 == 0;
  220.     bool bTop  = I == submatr.i1;
  221.     if(bLeft && !bTop) {
  222.         int jump = I - submatr.i1;
  223.         subpath_left.insert(subpath_left.begin(), jump, eTS_Delete);
  224.     }
  225.     if(!bLeft && bTop) {
  226.         subpath_left.insert(subpath_left.begin(), nc0, eTS_Insert);
  227.     }
  228.     if(bLeft || bTop) {
  229.         bNoLT = true;
  230.     }
  232.     // check if we reached right or bottom line
  233.     bool bNoRB = false;
  234.     int nc1 = trans_pos + steps_right;
  235.     if(nc1 > dimJ) {
  236.         NCBI_THROW(CAlgoAlignException, eInternal,
  237.                    "Assertion: RB overflow");
  238.     }
  239.     bool bRight = nc1 == dimJ;
  240.     bool bBottom  = I == submatr.i2 - 1;
  241.     if(bRight && !bBottom) {
  242.         int jump = submatr.i2 - I - 1;
  243.         subpath_right.insert(subpath_right.end(), jump, eTS_Delete);
  244.     }
  245.     if(!bRight && bBottom) {
  246.         int jump = dimJ - nc1;
  247.         subpath_right.insert(subpath_right.end(), jump, eTS_Insert);
  248.     }
  249.     if(bRight || bBottom) {
  250.         bNoRB = true;
  251.     }
  253.     // find out if we have special left-top and/or right-bottom
  254.     // on the new sub-matrices
  255.     bool rb = subpath_left.front() == eTS_Delete;
  256.     bool lt = subpath_right.back() == eTS_Delete;
  258.     // coalesce the pieces together and insert into the master list
  259.     subpath_left.splice( subpath_left.end(), subpath_right );
  261.     list<ETranscriptSymbol>::iterator ti0, ti1;
  262.     {{
  263.         CFastMutexGuard  guard (masterlist_mutex);
  264.         ti0 = translist_pos;
  265.         --ti0;
  266.         m_TransList.splice( translist_pos, subpath_left );
  267.         ++ti0;
  268.         ti1 = translist_pos;
  269.     }}
  271.     // Recurse submatrices:
  272.     SCoordRect rlt;
  273.     if(!bNoLT) {
  274.         // left top
  275.         size_t left_bnd = submatr.j1 + trans_pos - steps_left - 1;
  276.         if(left_bnd >= submatr.j1) {
  277.             rlt.i1 = submatr.i1;
  278.             rlt.j1 = submatr.j1;
  279.             rlt.i2 = I - 1;
  280.             rlt.j2 = left_bnd;
  281.         }
  282.         else {
  283.             NCBI_THROW(CAlgoAlignException, eInternal,
  284.                        "Assertion: Left boundary out of range");
  285.         }
  286.     }
  288.     SCoordRect rrb;
  289.     if(!bNoRB) {
  290.         // right bottom
  291.         size_t right_bnd = submatr.j1 + trans_pos + steps_right;
  292.         if(right_bnd <= submatr.j2) {
  293.             rrb.i1 = I + 2;
  294.             rrb.j1 = right_bnd;
  295.             rrb.i2 = submatr.i2;
  296.             rrb.j2 = submatr.j2;
  297.         }
  298.         else {
  299.             NCBI_THROW(CAlgoAlignException, eInternal,
  300.                        "Assertion: Right boundary out of range");
  301.         }
  302.     }
  304.     if(!bNoLT && !bNoRB) {
  305.         if( m_mt && m_maxthreads > 1 && MM_RequestNewThread(m_maxthreads) ) {
  306.             // find out which rect is larger
  307.             if(rlt.GetArea() > rrb.GetArea()) {
  308.                 // do rb in a second thread
  309.                 CThread* thr2 = new CThreadDoSM( this, &rrb, ti1,
  310.                                                  lt, right_bottom);
  311.                 thr2->Run();
  312.                 // do lt in this thread
  313.                 x_DoSubmatrix(rlt, ti0, left_top, rb);
  314.                 thr2->Join(0);
  315.             }
  316.             else {
  317.                 // do lt in a second thread
  318.                 CThread* thr2 = new CThreadDoSM( this, &rlt, ti0,
  319.                                                  left_top, rb);
  320.                 thr2->Run();
  321.                 // do rb in this thread
  322.                 x_DoSubmatrix(rrb, ti1, lt, right_bottom);
  323.                 thr2->Join(0);
  324.             }
  325.         }
  326.         else {  // just do both
  327.             x_DoSubmatrix(rlt, ti0, left_top, rb);
  328.             x_DoSubmatrix(rrb, ti1, lt, right_bottom);
  329.         }
  330.     }
  331.     else if(!bNoLT) {
  332.         x_DoSubmatrix(rlt, ti0, left_top, rb);
  333.     }
  334.     else if(!bNoRB) {
  335.         x_DoSubmatrix(rrb, ti1, lt, right_bottom);
  336.     }
  337. }
  340. CNWAligner::TScore CMMAligner::x_FindBestJ (
  341.     const vector<TScore>& vEtop, const vector<TScore>& vFtop,
  342.     const vector<TScore>& vGtop, const vector<TScore>& vEbtm,
  343.     const vector<TScore>& vFbtm, const vector<TScore>& vGbtm,
  344.     size_t& pos,
  345.     ETransitionType& trans_type ) const
  346. {
  347.     const size_t dim = vEtop.size();
  348.     TScore score = kMin_Int;
  349.     TScore trans_alts [9];
  351.     bool bFreeGapLeft2  = m_esf_L2 && dim == m_SeqLen2 + 1;
  352.     bool bFreeGapRight2 = m_esf_R2 && dim == m_SeqLen2 + 1;
  354.     for(size_t i = 0; i < dim ; ++i) {
  355.         trans_alts [0] = vEtop[i] + vEbtm[i] - m_Wg;   // II
  356.         trans_alts [1] = vFtop[i] + vEbtm[i];          // DI
  357.         trans_alts [2] = vGtop[i] + vEbtm[i];          // GI
  358.         trans_alts [3] = vEtop[i] + vFbtm[i];          // ID
  359.         TScore wg = (bFreeGapLeft2 && i == 0 || bFreeGapRight2 && i == dim -1)?
  360.                     0: m_Wg;
  361.         trans_alts [4] = vFtop[i] + vFbtm[i] - wg;     // DD
  362.         trans_alts [5] = vGtop[i] + vFbtm[i];          // GD
  363.         trans_alts [6] = vEtop[i] + vGbtm[i];          // IG
  364.         trans_alts [7] = vFtop[i] + vGbtm[i];          // DG
  365.         trans_alts [8] = vGtop[i] + vGbtm[i];          // GG
  367.         for(size_t k = 0; k < 9; ++k) {
  368.             if(trans_alts[k] > score) {
  369.                 score = trans_alts[k];
  370.                 pos = i;
  371.                 trans_type = (ETransitionType)k;
  372.             }
  373.         }
  374.     }
  376.     return score;
  377. }
  380. size_t CMMAligner::x_ExtendSubpath (
  381.     vector<unsigned char>::const_iterator trace_it,
  382.     bool direction,
  383.     list<ETranscriptSymbol>& subpath ) const
  384. {
  385.     subpath.clear();
  386.     size_t step_counter = 0;
  387.     if(direction) {
  388.         while(true) {
  389.             unsigned char Key = *trace_it;
  390.             if (Key & kMaskD) {
  391.                 subpath.push_back(eTS_Match);
  392.                 ++step_counter;
  393.                 break;
  394.             }
  395.             else if (Key & kMaskE) {
  396.                 subpath.push_back(eTS_Insert);
  397.                 ++step_counter;
  398.                 ++trace_it;
  399.                 while(Key & kMaskEc) {
  400.                     Key = *trace_it++;
  401.                     subpath.push_back(eTS_Insert);
  402.                     ++step_counter;
  403.                 }
  404.             }
  405.             else if ((Key & kMaskE) == 0) {
  406.                 subpath.push_back(eTS_Delete);
  407.                 break;
  408.             }
  409.             else {
  410.                 NCBI_THROW(CAlgoAlignException, eInternal,
  411.                            "Assertion: incorrect backtrace symbol "
  412.                            "(right expansion)");
  413.             }
  414.         }
  415.     }
  416.     else {
  417.         while(true) {
  418.             unsigned char Key = *trace_it;
  419.             if (Key & kMaskD) {
  420.                 subpath.push_front(eTS_Match);
  421.                 ++step_counter;
  422.                 break;
  423.             }
  424.             else if (Key & kMaskE) {
  425.                 subpath.push_front(eTS_Insert);
  426.                 ++step_counter;
  427.                 --trace_it;
  428.                 while(Key & kMaskEc) {
  429.                     Key = *trace_it--;
  430.                     subpath.push_front(eTS_Insert);
  431.                     ++step_counter;
  432.                 }
  433.             }
  434.             else if ((Key & kMaskE) == 0) {
  435.                 subpath.push_front(eTS_Delete);
  436.                 break;
  437.             }
  438.             else {
  439.                 NCBI_THROW(CAlgoAlignException, eInternal,
  440.                            "Assertion: incorrect backtrace symbol "
  441.                            "(left expansion)");
  442.             }
  443.         }
  444.     }
  445.     return step_counter;
  446. }
  449. #ifdef NCBI_THREADS
  450. DEFINE_STATIC_FAST_MUTEX (progress_mutex);
  451. #endif
  454. void CMMAligner::x_RunTop ( const SCoordRect& rect,
  455.              vector<TScore>& vE, vector<TScore>& vF, vector<TScore>& vG,
  456.              vector<unsigned char>& trace, bool lt ) const
  457. {
  458.     if( m_terminate ) {
  459.         return;
  460.     }
  462.     const size_t dim1 = rect.i2 - rect.i1 + 1;
  463.     const size_t dim2 = rect.j2 - rect.j1 + 1;
  464.     const size_t N1   = dim1 + 1;
  465.     const size_t N2   = dim2 + 1;
  467.     vector<TScore> stl_rowV (N2), stl_rowF (N2);
  468.     TScore* rowV    = &stl_rowV [0];
  469.     TScore* rowF    = &stl_rowF [0];
  471.     TScore* pV = rowV - 1;
  473.     const char* seq1 = m_Seq1 - 1 + rect.i1;
  474.     const char* seq2 = m_Seq2 - 1 + rect.j1;
  476.     const TNCBIScore (*sm) [NCBI_FSM_DIM] = m_ScoreMatrix.s;
  478.     bool bFreeGapLeft1  = m_esf_L1 && rect.i1 == 0;
  479.     bool bFreeGapLeft2  = m_esf_L2 && rect.j1 == 0;
  480.     bool bFreeGapRight2  = m_esf_R2 && rect.j2 == m_SeqLen2 - 1;
  482.     // progress reporting
  483.     const size_t prg_rep_rate = 100;
  484.     const size_t prg_rep_increment = prg_rep_rate*N2;
  486.     // first row
  488.     TScore wg = bFreeGapLeft1? 0: m_Wg;
  489.     TScore ws = bFreeGapLeft1? 0: m_Ws;
  491.     rowV[0] = wg;
  492.     size_t i, j;
  493.     for (j = 1; j < N2; ++j) {
  494.         rowV[j] = pV[j] + ws;
  495.         rowF[j] = kInfMinus;        
  496.     }
  497.     rowV[0] = 0;
  499.     // recurrences
  501.     wg = bFreeGapLeft2? 0: m_Wg;
  502.     ws = bFreeGapLeft2? 0: m_Ws;
  504.     TScore V  = 0;
  505.     TScore V0 = lt? 0: wg;
  506.     TScore E, G, n0;
  508.     for(i = 1;  i < N1 - 1;  ++i) {
  509.         
  510.         V = V0 += ws;
  511.         E = kInfMinus;
  512.         unsigned char ci = seq1[i];
  514.         TScore wg2 = m_Wg, ws2 = m_Ws;
  516.         for (j = 1; j < N2; ++j) {
  518.             G = pV[j] + sm[ci][(unsigned char)seq2[j]];
  519.             pV[j] = V;
  521.             n0 = V + m_Wg;
  522.             if(E >= n0)
  523.                 E += m_Ws;      // continue the gap
  524.             else
  525.                 E = n0 + m_Ws;  // open a new gap
  527.             if(j == N2 - 1 && bFreeGapRight2)
  528.                 wg2 = ws2 = 0;
  530.             n0 = rowV[j] + wg2;
  531.             if (rowF[j] > n0)
  532.                 rowF[j] += ws2;
  533.             else
  534.                 rowF[j] = n0 + ws2;
  536.             V = (E >= rowF[j])? (E >= G? E: G): (rowF[j] >= G? rowF[j]: G);
  537.         }
  538.         pV[j] = V;
  540.         if( m_prg_callback && i % prg_rep_rate == 0 ) {
  541. #ifdef NCBI_THREADS
  542.             CFastMutexGuard guard (progress_mutex);
  543. #endif
  544.             m_prg_info.m_iter_done += prg_rep_increment;
  545.             if( m_terminate = m_prg_callback(&m_prg_info) ) {
  546.                 break;
  547.             }
  548.         }
  549.     }
  551.     // the last row (i == N1 - 1)
  552.     if(!m_terminate) {
  554.         vG[0] = vE[0] = E = kInfMinus;
  555.         vF[0] = V = V0 += ws;
  556.         trace[0] = kMaskFc;
  557.         unsigned char ci = seq1[i];
  559.         TScore wg2 = m_Wg, ws2 = m_Ws;
  561.         unsigned char tracer;
  562.         for (j = 1; j < N2; ++j) {
  564.             vG[j] = G = pV[j] + sm[ci][(unsigned char)seq2[j]];
  565.             pV[j] = V;
  567.             n0 = V + m_Wg;
  568.             if(E >= n0) {
  569.                 E += m_Ws;      // continue the gap
  570.                 tracer = kMaskEc;
  571.             }
  572.             else {
  573.                 E = n0 + m_Ws;  // open a new gap
  574.                 tracer = 0;
  575.             }
  576.             vE[j] = E;
  578.             if(j == N2 - 1 && bFreeGapRight2) {
  579.                 wg2 = ws2 = 0;
  580.             }
  582.             n0 = rowV[j] + wg2;
  583.             if(rowF[j] >= n0) {
  584.                 rowF[j] += ws2;
  585.                 tracer |= kMaskFc;
  586.             }
  587.             else {
  588.                 rowF[j] = n0 + ws2;
  589.             }
  590.             vF[j] = rowF[j];
  592.             if (E >= rowF[j]) {
  593.                 if(E >= G) {
  594.                     V = E;
  595.                     tracer |= kMaskE;
  596.                 }
  597.                 else {
  598.                     V = G;
  599.                     tracer |= kMaskD;
  600.                 }
  601.             } else {
  602.                 if(rowF[j] >= G) {
  603.                     V = rowF[j];
  604.                 }
  605.                 else {
  606.                     V = G;
  607.                     tracer |= kMaskD;
  608.                 }
  609.             }
  610.             trace[j] = tracer;
  611.         }
  612.     }
  614.     if( m_prg_callback ) {
  615. #ifdef NCBI_THREADS
  616.         CFastMutexGuard guard (progress_mutex);
  617. #endif
  618.         m_prg_info.m_iter_done += (i % prg_rep_rate)*N2;
  619.         m_terminate = m_prg_callback(&m_prg_info);
  620.     }
  621. }
  623.  
  624. void CMMAligner::x_RunBtm(const SCoordRect& rect,
  625.              vector<TScore>& vE, vector<TScore>& vF, vector<TScore>& vG,
  626.              vector<unsigned char>& trace, bool rb) const
  627. {
  628.     if( m_terminate ) {
  629.         return;
  630.     }
  632.     const size_t dim1 = rect.i2 - rect.i1 + 1;
  633.     const size_t dim2 = rect.j2 - rect.j1 + 1;
  634.     const size_t N1   = dim1 + 1;
  635.     const size_t N2   = dim2 + 1;
  637.     vector<TScore> stl_rowV (N2), stl_rowF (N2);
  638.     TScore* rowV    = &stl_rowV [0];
  639.     TScore* rowF    = &stl_rowF [0];
  641.     TScore* pV = rowV + 1;
  643.     const char* seq1 = m_Seq1 + rect.i1;
  644.     const char* seq2 = m_Seq2 + rect.j1;
  646.     const TNCBIScore (*sm) [NCBI_FSM_DIM] = m_ScoreMatrix.s;
  648.     bool bFreeGapRight1  = m_esf_R1 && rect.i2 == m_SeqLen1 - 1;
  649.     bool bFreeGapRight2  = m_esf_R2 && rect.j2 == m_SeqLen2 - 1;
  650.     bool bFreeGapLeft2  =  m_esf_L2 && rect.j1 == 0;
  652.     // progress reporting
  654.     const size_t prg_rep_rate = 100;
  655.     const size_t prg_rep_increment = prg_rep_rate*N2;
  657.     // bottom row
  659.     TScore wg = bFreeGapRight1? 0: m_Wg;
  660.     TScore ws = bFreeGapRight1? 0: m_Ws;
  662.     rowV[N2 - 1] = wg;
  663.     int i, j;
  664.     for (j = N2 - 2; j >= 0; --j) {
  665.         rowV[j] = pV[j] + ws;
  666.         rowF[j] = kInfMinus;
  667.     }
  668.     rowV[N2 - 1] = 0;
  670.     // recurrences
  672.     wg = bFreeGapRight2? 0: m_Wg;
  673.     ws = bFreeGapRight2? 0: m_Ws;
  675.     TScore V  = 0;
  676.     TScore V0 = rb? 0: wg;
  677.     TScore E, G, n0;
  679.     for(i = N1 - 2;  i > 0;  --i) {
  680.         
  681.         V = V0 += ws;
  682.         E = kInfMinus;
  683.         unsigned char ci = seq1[i];
  685.         TScore wg2 = m_Wg, ws2 = m_Ws;
  687.         for (j = N2 - 2; j >= 0; --j) {
  689.             G = pV[j] + sm[ci][(unsigned char)seq2[j]];
  690.             pV[j] = V;
  692.             n0 = V + m_Wg;
  693.             if(E >= n0)
  694.                 E += m_Ws;      // continue the gap
  695.             else
  696.                 E = n0 + m_Ws;  // open a new gap
  698.             if(j == 0 && bFreeGapLeft2) {
  699.                 wg2 = ws2 = 0;
  700.             }
  702.             n0 = rowV[j] + wg2;
  703.             if (rowF[j] > n0)
  704.                 rowF[j] += ws2;
  705.             else
  706.                 rowF[j] = n0 + ws2;
  708.             V = (E >= rowF[j])? (E >= G? E: G): (rowF[j] >= G? rowF[j]: G);
  709.         }
  710.         pV[j] = V;
  712.         if( m_prg_callback && (N1 - i) % prg_rep_rate == 0 ) {
  713. #ifdef NCBI_THREADS
  714.             CFastMutexGuard guard (progress_mutex);
  715. #endif
  716.             m_prg_info.m_iter_done += prg_rep_increment;
  717.             if( m_terminate = m_prg_callback(&m_prg_info) ) {
  718.                 break;
  719.             }
  720.         }
  721.     }
  723.     // the top row (i == 0)
  724.     if(!m_terminate) {
  726.         vF[N2-1] = V = V0 += ws;
  727.         vG[N2-1] = vE[N2-1] = E = kInfMinus;
  728.         trace[N2-1] = kMaskFc;
  729.         unsigned char ci = seq1[i];
  731.         TScore wg2 = m_Wg, ws2 = m_Ws;
  733.         unsigned char tracer;
  734.         for (j = N2 - 2; j >= 0; --j) {
  736.             vG[j] = G = pV[j] + sm[ci][(unsigned char)seq2[j]];
  737.             pV[j] = V;
  739.             n0 = V + m_Wg;
  740.             if(E >= n0) {
  741.                 E += m_Ws;      // continue the gap
  742.                 tracer = kMaskEc;
  743.             }
  744.             else {
  745.                 E = n0 + m_Ws;  // open a new gap
  746.                 tracer = 0;
  747.             }
  748.             vE[j] = E;
  750.             if(j == 0 && bFreeGapLeft2) {
  751.                 wg2 = ws2 = 0;
  752.             }
  754.             n0 = rowV[j] + wg2;
  755.             if(rowF[j] >= n0) {
  756.                 rowF[j] += ws2;
  757.                 tracer |= kMaskFc;
  758.             }
  759.             else {
  760.                 rowF[j] = n0 + ws2;
  761.             }
  762.             vF[j] = rowF[j];
  764.             if (E >= rowF[j]) {
  765.                 if(E >= G) {
  766.                     V = E;
  767.                     tracer |= kMaskE;
  768.                 }
  769.                 else {
  770.                     V = G;
  771.                     tracer |= kMaskD;
  772.                 }
  773.             } else {
  774.                 if(rowF[j] >= G) {
  775.                     V = rowF[j];
  776.                 }
  777.                 else {
  778.                     V = G;
  779.                     tracer |= kMaskD;
  780.                 }
  781.             }
  782.             trace[j] = tracer;
  783.         }
  784.     }
  786.     if( m_prg_callback ) {
  787. #ifdef NCBI_THREADS
  788.         CFastMutexGuard guard (progress_mutex);
  789. #endif
  790.         m_prg_info.m_iter_done += (N1 - i) % prg_rep_rate;
  791.         m_terminate = m_prg_callback(&m_prg_info);
  792.     }
  793. }
  796. CNWAligner::TScore CMMAligner::x_RunTerm(const SCoordRect& rect,
  797.                                          bool left_top, bool right_bottom,
  798.                                          list<ETranscriptSymbol>& subpath)
  799. {
  800.     if( m_terminate ) {
  801.         return 0;
  802.     }
  804.     const size_t N1 = rect.i2 - rect.i1 + 2;
  805.     const size_t N2 = rect.j2 - rect.j1 + 2;
  807.     vector<TScore> stl_rowV (N2), stl_rowF (N2);
  808.     TScore* rowV    = &stl_rowV [0];
  809.     TScore* rowF    = &stl_rowF [0];
  811.     // index calculation: [i,j] = i*n2 + j
  812.     vector<unsigned char> stl_bm (N1*N2);
  813.     unsigned char* backtrace = &stl_bm[0];
  815.     TScore* pV = rowV - 1;
  817.     const char* seq1 = m_Seq1 + rect.i1 - 1;
  818.     const char* seq2 = m_Seq2 + rect.j1 - 1;
  820.     const TNCBIScore (*sm) [NCBI_FSM_DIM] = m_ScoreMatrix.s;
  822.     bool bFreeGapLeft1  = m_esf_L1 && rect.i1 == 0;
  823.     bool bFreeGapRight1 = m_esf_R1 && rect.i2 == m_SeqLen1 - 1;
  824.     bool bFreeGapLeft2  = m_esf_L2 && rect.j1 == 0;
  825.     bool bFreeGapRight2 = m_esf_R2 && rect.j2 == m_SeqLen2 - 1;
  827.     TScore wgleft1   = bFreeGapLeft1? 0: m_Wg;
  828.     TScore wsleft1   = bFreeGapLeft1? 0: m_Ws;
  829.     TScore wg1 = m_Wg, ws1 = m_Ws;
  831.     // first row
  832.     size_t k;
  833.     {
  834.         rowV[0] = wgleft1;
  835.         for (k = 1; k < N2; k++) {
  836.             rowV[k] = pV[k] + wsleft1;
  837.             rowF[k] = kInfMinus;
  838.             backtrace[k] = kMaskE | kMaskEc;
  839.         }
  840.         rowV[0] = 0;
  841.     }
  843.     // recurrences
  844.     TScore wgleft2   = bFreeGapLeft2? 0: m_Wg;
  845.     TScore wsleft2   = bFreeGapLeft2? 0: m_Ws;
  846.     TScore V  = 0;
  847.     TScore V0 = left_top? 0: wgleft2;
  848.     TScore E, G, n0;
  849.     unsigned char tracer;
  851.     size_t i, j;
  852.     for(i = 1;  i < N1;  ++i) {
  853.         
  854.         V = V0 += wsleft2;
  855.         E = kInfMinus;
  856.         backtrace[k++] = kMaskFc;
  857.         unsigned char ci = seq1[i];
  859.         if(i == N1 - 1 && bFreeGapRight1) {
  860.                 wg1 = ws1 = 0;
  861.         }
  863.         TScore wg2 = m_Wg, ws2 = m_Ws;
  865.         for (j = 1; j < N2; ++j, ++k) {
  867.             G = pV[j] + sm[ci][(unsigned char)seq2[j]];
  868.             pV[j] = V;
  870.             n0 = V + wg1;
  871.             if(E >= n0) {
  872.                 E += ws1;      // continue the gap
  873.                 tracer = kMaskEc;
  874.             }
  875.             else {
  876.                 E = n0 + ws1;  // open a new gap
  877.                 tracer = 0;
  878.             }
  880.             if(j == N2 - 1 && bFreeGapRight2) {
  881.                 wg2 = ws2 = 0;
  882.             }
  883.             n0 = rowV[j] + ((right_bottom && j == N2 - 1)? 0: wg2);
  884.             if(rowF[j] >= n0) {
  885.                 rowF[j] += ws2;
  886.                 tracer |= kMaskFc;
  887.             }
  888.             else {
  889.                 rowF[j] = n0 + ws2;
  890.             }
  892.             if (E >= rowF[j]) {
  893.                 if(E >= G) {
  894.                     V = E;
  895.                     tracer |= kMaskE;
  896.                 }
  897.                 else {
  898.                     V = G;
  899.                     tracer |= kMaskD;
  900.                 }
  901.             } else {
  902.                 if(rowF[j] >= G) {
  903.                     V = rowF[j];
  904.                 }
  905.                 else {
  906.                     V = G;
  907.                     tracer |= kMaskD;
  908.                 }
  909.             }
  910.             backtrace[k] = tracer;
  911.         }
  913.         pV[j] = V;
  914.     }
  916.     // fill the subpath
  917.     subpath.clear();
  918.     
  919.     // run backtrace
  920.     k = N1*N2 - 1;
  921.     while (k != 0) {
  922.         unsigned char Key = backtrace[k];
  923.         if (Key & kMaskD) {
  924.             subpath.push_front(eTS_Match);
  925.             k -= N2 + 1;
  926.         }
  927.         else if (Key & kMaskE) {
  928.             subpath.push_front(eTS_Insert); --k;
  929.             while(k > 0 && (Key & kMaskEc)) {
  930.                 subpath.push_front(eTS_Insert);
  931.                 Key = backtrace[k--];
  932.             }
  933.         }
  934.         else {
  935.             subpath.push_front(eTS_Delete);
  936.             k -= N2;
  937.             while(k > 0 && (Key & kMaskFc)) {
  938.                 subpath.push_front(eTS_Delete);
  939.                 Key = backtrace[k];
  940.                 k -= N2;
  941.             }
  942.         }
  943.     }
  945.     return V;
  946. }
  949. bool CMMAligner::x_CheckMemoryLimit()
  950. {
  951.     return true;
  952. }
  956. END_NCBI_SCOPE
  958. /*
  959.  * ===========================================================================
  960.  * $Log: mm_aligner.cpp,v $
  961.  * Revision 1000.1  2004/06/01 18:04:46  gouriano
  962.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.20
  963.  *
  964.  * Revision 1.20  2004/05/21 21:41:01  gorelenk
  965.  * Added PCH ncbi_pch.hpp
  966.  *
  967.  * Revision 1.19  2004/05/18 21:43:40  kapustin
  968.  * Code cleanup
  969.  *
  970.  * Revision 1.18  2004/05/17 14:50:56  kapustin
  971.  * Add/remove/rearrange some includes and object declarations
  972.  *
  973.  * Revision 1.17  2003/09/30 19:50:04  kapustin
  974.  * Make use of standard score matrix interface
  975.  *
  976.  * Revision 1.16  2003/09/26 14:43:18  kapustin
  977.  * Remove exception specifications
  978.  *
  979.  * Revision 1.15  2003/09/04 16:07:38  kapustin
  980.  * Use STL vectors for exception-safe dynamic arrays and matrices
  981.  *
  982.  * Revision 1.14  2003/09/02 22:36:57  kapustin
  983.  * Adjust for transcript symbol name changes
  984.  *
  985.  * Revision 1.13  2003/06/17 17:20:44  kapustin
  986.  * CNWAlignerException -> CAlgoAlignException
  987.  *
  988.  * Revision 1.12  2003/06/17 14:51:04  dicuccio
  989.  * Fixed after algo/ rearragnement
  990.  *
  991.  * Revision 1.11  2003/06/03 18:19:06  rsmith
  992.  * Move static mutex initialization to file from function scope,
  993.  * because of MW compiler choking (wrongly) over complex initialization.
  994.  *
  995.  * Revision 1.10  2003/06/02 14:04:49  kapustin
  996.  * Progress indication-related updates
  997.  *
  998.  * Revision 1.9  2003/05/23 18:26:38  kapustin
  999.  * Use weak comparisons in core recurrences.
  1000.  *
  1001.  * Revision 1.8  2003/04/14 19:01:49  kapustin
  1002.  * Run() -> x_Run()
  1003.  *
  1004.  * Revision 1.7  2003/03/18 15:15:51  kapustin
  1005.  * Implement virtual memory checking function. Allow separate free
  1006.  * end gap specification
  1007.  *
  1008.  * Revision 1.6  2003/03/17 15:30:57  kapustin
  1009.  * Support end-space free alignments
  1010.  *
  1011.  * Revision 1.5  2003/01/30 20:32:06  kapustin
  1012.  * Call x_LoadScoringMatrix() from the base class constructor.
  1013.  *
  1014.  * Revision 1.4  2003/01/28 12:37:40  kapustin
  1015.  * Move m_score to the base class
  1016.  *
  1017.  * Revision 1.3  2003/01/22 13:40:09  kapustin
  1018.  * Implement multithreaded algorithm
  1019.  *
  1020.  * ===========================================================================
  1021.  */