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

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

生物技术

开发平台：

C/C++

splign_hitparser.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: splign_hitparser.cpp,v $
  4.  * PRODUCTION Revision 1000.0  2004/06/01 18:12:53  gouriano
  5.  * PRODUCTION PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.7
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: splign_hitparser.cpp,v 1000.0 2004/06/01 18:12:53 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. * Author:  Yuri Kapustin
  36. *
  37. * File Description:
  38. *   Temporary code - will be part of the hit filtering library
  39. *
  40. */
  43. #include <ncbi_pch.hpp>
  44. #include "splign_hitparser.hpp"
  45. #include <algo/align/align_exception.hpp>
  47. #include <corelib/ncbi_limits.hpp>
  48. #include <corelib/ncbistre.hpp>
  50. #include <algorithm>
  51. #include <numeric>
  54. BEGIN_NCBI_SCOPE
  56. // auxiliary structure utilized in calculation of collisions and other routines
  57. struct SNode {
  58.     int x; // coordinate
  59.     bool type; // true = open, false = close
  60.     vector<CHit>::const_iterator pHit; // owner
  61.     SNode(int n, bool b, vector<CHit>::const_iterator ph):
  62.         x(n), type(b), pHit(ph) {};
  63.     bool operator>(const SNode& rhs) const { return x > rhs.x; }
  64.     bool operator==(const SNode& rhs) const { return x == rhs.x; }
  65.     bool operator<(const SNode& rhs) const { return x < rhs.x; }
  66. };
  68. //
  69. // general single-linkage clustering (SLC) algorithm 
  70. //
  72. struct SNodePair
  73. {
  74.     SNodePair(int i, int j): n1(i), n2(j) {}
  75.     int n1, n2;
  76. };
  79. template<class T, class PClustCrit, class CClustID_Set, class CClustID_Get>
  80. void DoSingleLinkageClustering (
  81.    vector<T>& vt,          // vector of objects to be clustered
  82.    const PClustCrit& pr2,        // clustering criterion (two-argument predicate)
  83.    const CClustID_Set& CID_Set,  // function that assigns ClustID to objects
  84.    const CClustID_Get& CID_Get   // function that retrieves ClustID from objects
  85. )
  87. {
  88.     vector<SNodePair> vnodes;    // auxiliary vector
  89.     size_t nSize = vt.size(), j = 0, i;
  91.     for( j = 0; j < nSize; ++j ) {
  93.         CID_Set( vt[j], j );
  94.         for( i = 0; i < j; ++i ) {
  96.             if( pr2(vt[i], vt[j] ) )
  97.                 vnodes.push_back( SNodePair( i, j ) );
  98.          }
  99.     }
  100.  
  101.     vector<SNodePair>::const_iterator kk, kend;
  102.     for( kk = vnodes.begin(), kend = vnodes.end(); kk != kend; ++kk ) {
  104.         int n1 = CID_Get( vt[kk->n1] ), n2 = CID_Get( vt[kk->n2] );
  105.         if( n1 != n2 ) {
  106.             NON_CONST_ITERATE (typename vector<T>, jj, vt) {
  107.                 if( CID_Get( *jj ) == n1 )
  108.                     CID_Set( *jj, n2 );
  109.             }
  110.         }
  111.     }
  112. }
  115. CHitParser::CHitParser()
  116. {
  117.     x_Set2Defaults();
  118. }
  120. CHitParser::CHitParser(const vector<CHit>& vh, int& nGroupID)
  121. {
  122.     Init(vh, nGroupID);
  123. }
  125. void CHitParser::Init(const vector<CHit>& vh, int& nGroupID)
  126. {
  127.     m_Out = vh;
  128.     x_RemoveEqual();
  129.     x_Set2Defaults();
  130.     nGroupID++;
  131.     m_GroupID = &nGroupID;
  132. }
  134. void CHitParser::x_Set2Defaults()
  135. {
  136.     m_SameOrder = true;
  137.     m_Strand = eStrand_Auto;
  138.     m_MaxHitDistQ = kMax_Int;
  139.     m_MaxHitDistS = kMax_Int;
  140.     m_Method = eMethod_MaxScore;
  141.     m_SplitQuery   = eSplitMode_MaxScore;
  142.     m_SplitSubject = eSplitMode_MaxScore;
  143.     m_Prot2Nucl = false;
  145.     m_CombinationProximity_pre = m_CombinationProximity_post = -1.;
  147.     m_group_identification_method = eNone;
  148.     m_CovStep = 90.;
  150.     m_MinQueryCoverage = 0.;
  151.     m_MinSubjCoverage = 0.;
  153.     m_OutputAllGroups = false;
  155.     m_Query = m_Subj = "";
  157.     m_GroupID = 0;
  158.     
  159.     x_CalcGlobalEnvelop();
  160. }
  162. void CHitParser::x_CalcGlobalEnvelop()
  163. {
  164.     const size_t dim = m_Out.size();
  165.     if(dim == 0) {
  166.         m_ange[0] = m_ange[1] = m_ange[2] = m_ange[3] = 0;
  167.         return;
  168.     }
  169.     m_ange[0] = m_ange[2] = kMax_Int;
  170.     m_ange[1] = m_ange[3] = kMin_Int;
  171.     vector<CHit>::iterator ii, iend;
  172.     for(ii = m_Out.begin(), iend = m_Out.end(); ii != iend; ii++)
  173.     {
  174.         if(m_ange[0] > ii->m_ai[0]) m_ange[0] = ii->m_ai[0];
  175.         if(m_ange[2] > ii->m_ai[2]) m_ange[2] = ii->m_ai[2];
  176.         if(m_ange[1] < ii->m_ai[1]) m_ange[1] = ii->m_ai[1];
  177.         if(m_ange[3] < ii->m_ai[3]) m_ange[3] = ii->m_ai[3];
  178.     }
  179. }
  181. CHitParser::~CHitParser()
  182. {
  183. }
  186. bool  CHitParser::x_DetectInclusion(const CHit& h1, const CHit& h2) const
  187. {
  188.   if( h2.m_ai[0] <= h1.m_ai[0] && h1.m_ai[1] <= h2.m_ai[1]) {
  189.     return true;
  190.   }
  191.   if( h2.m_ai[2] <= h1.m_ai[2] && h1.m_ai[3] <= h2.m_ai[3]) {
  192.     return true;
  193.   }
  194.   if( h1.m_ai[0] <= h2.m_ai[0] && h2.m_ai[1] <= h1.m_ai[1]) {
  195.     return true;
  196.   }
  197.   if( h1.m_ai[2] <= h2.m_ai[2] && h2.m_ai[3] <= h1.m_ai[3]) {
  198.     return true;
  199.   }
  200.   return false;
  201. }
  204. // implements half splitmode
  205. bool CHitParser::x_RunAltSplitMode(char cSide, CHit& hit1, CHit& hit2)
  206. {
  207.     int nb = cSide*2;
  209.     // detect intersection
  210.     if(hit1.m_ai[nb] <= hit2.m_ai[nb] && hit2.m_ai[nb] <= hit1.m_ai[nb+1])
  211.     {
  212.         int n1 = hit1.m_ai[nb + 1] + 1;
  213.         int n2 = hit2.m_ai[nb] - 1;
  214.         if(m_Prot2Nucl)
  215.         {  // frame must be preserved
  216.             int nShift = n1 - hit2.m_ai[nb];
  217.             while(nShift % 3) { nShift++; n1++; }
  218.             nShift = hit1.m_ai[nb + 1] - n2;
  219.             while(nShift % 3) { nShift++; n2--; }
  220.         }
  221.         hit1.Move(nb + 1, n2, m_Prot2Nucl);
  222.         hit2.Move(nb, n1, m_Prot2Nucl);
  223.         hit1.UpdateAttributes();
  224.         hit2.UpdateAttributes();
  225.     }
  226.     else
  227.     if(hit2.m_ai[nb] <= hit1.m_ai[nb] && hit1.m_ai[nb] <= hit2.m_ai[nb+1])
  228.     {
  229.         int n1 = hit1.m_ai[nb] - 1;
  230.         int n2 = hit2.m_ai[nb + 1] + 1;
  231.         if(m_Prot2Nucl)
  232.         {  // frame must be preserved
  233.             int nShift = hit2.m_ai[nb + 1] - n1;
  234.             while(nShift % 3) { nShift++; n1--; }
  235.             nShift = n2 - hit1.m_ai[nb];
  236.             while(nShift % 3) { nShift++; n2++; }
  237.         }
  238.         hit1.Move(nb, n2, m_Prot2Nucl);
  239.         hit2.Move(nb + 1, n1, m_Prot2Nucl);
  240.         hit1.UpdateAttributes();
  241.         hit2.UpdateAttributes();
  242.     }
  244.     return true;
  245. }
  247. int GetIntersectionSize (int a0, int b0, int c0, int d0) {
  248.     int a = min(a0, b0), b = max(a0, b0), c = min(c0, d0), d = max(c0, d0);
  249.     int an1[] = {a, b};
  250.     int an2[] = {c, d};
  251.     int* apn[] = {an1, an2};
  252.     if(c <= a) {
  253.         apn[0] = an2; apn[1] = an1;
  254.     }
  255.     int n = (apn[0][1] - apn[1][0]) + 1;
  256.     return n > 0? n: 0;
  257. }
  260. // accumulate coverage on query or subject
  261. class AcmCvg: public binary_function<int, const CHit*, int>
  262. {
  263. public:
  264.     AcmCvg(char where): m_where(where), m_nFinish(-1) {
  265.         m_i1 = m_where == 'q'? 0: 2;
  266.         m_i2 = m_where == 'q'? 1: 3;
  267.     }
  268.     int operator() (int iVal, const CHit* ph)
  269.     {
  270.         const CHit& h = *ph;
  271.         if(h.m_ai[m_i1] > m_nFinish)
  272.             return iVal + (m_nFinish = h.m_ai[m_i2]) - h.m_ai[m_i1] + 1;
  273.         else
  274.         {
  275.             int nFinish0 = m_nFinish;
  276.             return (h.m_ai[m_i2] > nFinish0)?
  277.                 (iVal + (m_nFinish = h.m_ai[m_i2]) - nFinish0): iVal;
  278.         }
  279.     }
  280. private:
  281.     char m_where;
  282.     int  m_nFinish;
  283.     unsigned char m_i1, m_i2;
  284. };
  287. int CHitParser::CalcCoverage(vector<CHit>::const_iterator ib,
  288.                              vector<CHit>::const_iterator ie,
  289.                              char where)
  290. {
  291.     vector<const CHit*> vh (ie - ib, (const CHit*)0);
  292.     for(int i = 0; i < ie - ib; ++i)  vh[i] = &(*ib) + i;
  293.     if(where == 'q') {
  294.         stable_sort(vh.begin(), vh.end(), CHit::PPrecedeQ_ptr);
  295.     }
  296.     else {
  297.         stable_sort(vh.begin(), vh.end(), CHit::PPrecedeS_ptr);
  298.     }
  299.     int s = accumulate(vh.begin(), vh.end(), 0, AcmCvg(where));
  300.     return s;
  301. }
  303. //////// Maximum Score method:
  304. // 1. Sort the vector in desc order starting from the current hit
  305. // 2. Consider max score hit (the current)
  306. // 3. Process the other hits, so the current hit remains intact
  307. //    and overlapping free
  308. // 4. Go to step 1.
  309. int CHitParser::x_RunMaxScore()
  310. {
  311.     // special pre-processing for prot2nucl mode:
  312.     // remove those hits that intersect on the genomic
  313.     // side having different reading frame
  314.     if(m_Prot2Nucl && false)  // -- currently disabled
  315.     {
  316.         list<vector<CHit>::iterator > l0;
  317.         vector<CHit>::iterator ivh;
  318.         for(ivh = m_Out.begin(); ivh != m_Out.end(); ivh++)
  319.             l0.push_back(ivh);
  320.         list<vector<CHit>::iterator >::iterator ii;
  321.         bool bDel = false;
  322.         for(ii = l0.begin(); ii != l0.end(); )
  323.         {
  324.             list<vector<CHit>::iterator >::iterator jj = ii;
  325.             for(++jj; jj != l0.end(); jj++)
  326.             {
  327.                 int& a = (*ii)->m_ai[2];
  328.                 int& b = (*ii)->m_ai[3];                
  329.                 int& c = (*jj)->m_ai[2];
  330.                 int& d = (*jj)->m_ai[3];
  331.                 bool bSameFrame = a%3 == c%3;
  332.                 if( !bSameFrame && GetIntersectionSize(a, b, c, d) > 0)
  333.                 {
  334.                     l0.erase(jj);
  335.                     bDel = true;
  336.                     break;
  337.                 }
  338.             }
  339.             if(bDel)
  340.             {
  341.                 list<vector<CHit>::iterator >::iterator jj0 = ii; jj0++;
  342.                 l0.erase(ii);
  343.                 ii = jj0;
  344.                 bDel = false;
  345.             }
  346.             else
  347.                 ii++;
  348.         }
  349.         vector<CHit> vh1;
  350.         for(ii = l0.begin(); ii != l0.end(); ii++)
  351.             vh1.push_back(**ii);
  352.         m_Out = vh1;
  353.     }
  355.     size_t nQuerySize = 0;
  356.     if(m_MinQueryCoverage > 0.) {
  357.         NCBI_THROW( CAlgoAlignException,
  358.                     eInternal,
  359.                     "Query coverage filtering not supported" );
  360.     }
  362.     size_t nSubjSize = 0;
  363.     if(m_MinSubjCoverage > 0.) {
  364.         NCBI_THROW( CAlgoAlignException,
  365.                     eInternal,
  366.                     "Subj coverage filtering not supported" );
  367.     }
  369.     // main processing
  370.     int i0 = 0;
  371.     bool bRestart = true;
  372.     while(bRestart) {
  373.         bRestart = false;
  374.         for(i0 = 0; i0 < int(m_Out.size())-1; ++i0) {
  376.             vector<CHit>::iterator ii = m_Out.begin() + i0;
  377.             // Sort by score (& priority) starting from the current hit
  378.             sort(ii, m_Out.end(), hit_greater());
  379.             // Consider max score hit (the current)
  380.             CHit& hit = *ii;
  381.             // Process the other hits, so the current hit
  382.             // remains intact and overlapping free
  383.             int j0 = 0;
  384.             for(j0 = int(m_Out.size())-1; j0 > i0; --j0) {
  386.                 vector<CHit>::iterator jj = m_Out.begin() + j0;
  388. if(m_SplitQuery == eSplitMode_Clear
  389.    || m_SplitSubject == eSplitMode_Clear) {
  390.   if(x_DetectInclusion(hit, *jj)) {
  391.     m_Out.erase(jj);
  392.     continue;
  393.   }
  394. }
  396.                 if(m_SplitQuery == eSplitMode_Clear) {
  397.                     if(x_RunAltSplitMode(0, hit, *jj) == false) {
  398.                         m_Out.erase(jj);
  399.                         continue;
  400.                     }
  401.                 }
  402.                 if(m_SplitSubject == eSplitMode_Clear) {
  403.                     if(x_RunAltSplitMode(1, hit, *jj) == false) {
  404.                         m_Out.erase(jj);
  405.                         continue;
  406.                     }
  407.                 }
  409.                 // We have two possibilities:
  410.                 // (1) there is at least one end of *jj inside the current hit;
  411.                 // (2) when *jj "embraces" the current hit;
  413.                 // (1) - proceed it in cyclic manner
  414.                 bool b [4];
  415.                 bool bFirstLoop = true;
  416.                 bool bContinue = false;
  417.                 do {
  418.                     if(!bFirstLoop) {
  419.                         size_t i = 0;
  420.                         for(i = 0; i < 4; ++i) {
  421.                             if(b[i]) {
  422.                                 switch(i) {
  424.                                 case 0: {
  425.                                     jj->Move(0,
  426.                                              hit.m_ai[1] + 1,
  427.                                              m_Prot2Nucl);
  428.                                     for(size_t d = 2; d < 4; ++d) {
  429.                                         b[d] = hit.m_ai[2] <= jj->m_ai[d] &&
  430.                                                jj->m_ai[d] <= hit.m_ai[3];
  431.                                     }
  432.                                 }
  433.                                 break;
  435.                                 case 1: {
  436.                                     jj->Move(1,
  437.                                              hit.m_ai[0] - 1,
  438.                                              m_Prot2Nucl);
  439.                                     for(size_t d = 2; d < 4; ++d) {
  440.                                         b[d] = hit.m_ai[2] <= jj->m_ai[d] &&
  441.                                                jj->m_ai[d] <= hit.m_ai[3];
  442.                                     }
  443.                                 }
  444.                                 break;
  446.                                 case 2: {
  447.                                     jj->Move(2,
  448.                                              hit.m_ai[3] + 1,
  449.                                              m_Prot2Nucl);
  450.                                     for(size_t d = 0; d < 2; ++d) {
  451.                                         b[d] = hit.m_ai[0] <= jj->m_ai[d] &&
  452.                                                jj->m_ai[d] <= hit.m_ai[1];
  453.                                     }
  454.                                 }
  455.                                 break;
  457.                                 case 3: {
  458.                                     jj->Move(3,
  459.                                              hit.m_ai[2] - 1,
  460.                                              m_Prot2Nucl);
  461.                                     for(size_t d = 0; d < 2; ++d) {
  462.                                         b[d] = hit.m_ai[0] <= jj->m_ai[d] &&
  463.                                                jj->m_ai[d] <= hit.m_ai[1];
  464.                                     }
  465.                                 }
  466.                                 break;
  467.                                 }
  468.                             }
  469.                             
  470.                             // check if *jj still exists
  471.                             if(!jj->IsConsistent())
  472.                             {
  473.                                 m_Out.erase(jj);
  474.                                 bContinue = true;
  475.                                 break;
  476.                             }
  477.                         }
  478.                         if(bContinue) break;
  479.                     }
  480.                     else
  481.                         bFirstLoop = false;
  483.                     size_t i = 0;
  484.                     for(i = 0; i < 2; ++i) {
  485.                         b[i] = hit.m_ai[0] <= jj->m_ai[i] &&
  486.                             jj->m_ai[i] <= hit.m_ai[1];
  487.                     }
  488.                     for(i = 2; i < 4; ++i) {
  489.                         b[i] = hit.m_ai[2] <= jj->m_ai[i] &&
  490.                             jj->m_ai[i] <= hit.m_ai[3];
  491.                     }
  492.                     if(b[0] && b[1] || b[2] && b[3]) {
  493.                         m_Out.erase(jj);
  494.                         bContinue = true;
  495.                         break;
  496.                     }
  497.                 }
  498.                 while(b[0] || b[1] || b[2] || b[3]);
  499.             
  500.                 if(bContinue)
  501.                     continue;    // the hit has dissapeared
  503.                 // (2)
  504.                 double ad[] = {
  505.                     hit.m_ai[0] - jj->m_ai[0], jj->m_ai[1] - hit.m_ai[1],
  506.                     hit.m_ai[2] - jj->m_ai[2], jj->m_ai[3] - hit.m_ai[3] };
  507.                 bool ab[] = {ad[0] > 0, ad[1] > 0, ad[2] > 0, ad[3] > 0 };
  508.                 bool bNewHit = false;
  509.                 CHit hit2;
  510.                 if(ab[0] && ab[1] || ab[2] && ab[3]) { // split *jj
  511.                     int n1, n2, n1x, n2x;
  512.                     if(jj->IsStraight()) {
  513.                         if(!(ab[2] && ab[3])) {
  514.                             n1 = 1; n2 = 0;
  515.                             n1x = hit.m_ai[0] - 1;
  516.                             n2x = hit.m_ai[1] + 1;
  517.                         }
  518.                         else if(!(ab[0] && ab[1])) {
  519.                             n1 = 3; n2 = 2;
  520.                             n1x = hit.m_ai[2] - 1;
  521.                             n2x = hit.m_ai[3] + 1;
  522.                         }
  523.                         else {
  524.                             n1 = (ad[0] < ad[2])? 1: 3;
  525.                             n2 = (ad[1] < ad[3])? 0: 2;
  526.                             n1x = (ad[0] < ad[2])?
  527.                                 hit.m_ai[0] - 1: hit.m_ai[2] - 1;
  528.                             n2x = (ad[1] < ad[3])?
  529.                                 hit.m_ai[1] + 1: hit.m_ai[3] + 1;
  530.                         }
  531.                     }
  532.                     else {
  533.                         if(!(ab[2] && ab[3])) {
  534.                             n1 = 1; n2 = 0;
  535.                             n1x = hit.m_ai[0] - 1; n2x = hit.m_ai[1] + 1;
  536.                         }
  537.                         else if(!(ab[0] && ab[1])) {
  538.                             n1 = 2; n2 = 3;
  539.                             n1x = hit.m_ai[3] + 1; n2x = hit.m_ai[2] - 1;
  540.                         }
  541.                         else {
  542.                             n1 = (ad[0] < ad[3])? 1: 2;
  543.                             n2 = (ad[1] < ad[2])? 0: 3;
  544.                             n1x = (ad[0] < ad[3])?
  545.                                 hit.m_ai[0] - 1: hit.m_ai[3] + 1;
  546.                             n2x = (ad[1] < ad[2])?
  547.                                 hit.m_ai[1] + 1: hit.m_ai[2] - 1;
  548.                         }
  549.                     }
  550.                     hit2 = *jj;
  551.                     bNewHit = true;
  552.                     jj->Move(n1, n1x, m_Prot2Nucl);
  553.                     if(!jj->IsConsistent()) {
  554.                         m_Out.erase(jj);
  555.                     }
  556.                     else { // some trapezoidal hits may still overlap
  557.                         bRestart = true;
  558.                     }
  559.                     hit2.Move(n2, n2x, m_Prot2Nucl);
  560.                     if(hit2.IsConsistent()) {
  561.                         // some trapezoidal hits may still overlap
  562.                         hit2.UpdateAttributes();
  563.                         bRestart = true;
  564.                     }
  565.                     else {
  566.                         bNewHit = false;
  567.                     }
  568.                 }
  570.                 jj->UpdateAttributes();
  571.                 if(bNewHit) {
  572.                     m_Out.push_back(hit2);
  573.                 }
  574.                 if(bRestart) break;
  575.             }
  577.             if(m_SameOrder) {
  578.                 x_FilterByOrder(i0, false);
  579.             }
  581.             // step 4: re-assess groups and continue the cycle
  582.             if(bRestart) break;
  583.         }
  584.     }
  586.     if(m_MinQueryCoverage > 0.) {
  587.         // calculate query coverage
  588.         double dQCov = CalcCoverage(m_Out.begin(), m_Out.end(), 'q');
  589.         if(dQCov < m_MinQueryCoverage* nQuerySize)
  590.             m_Out.clear();
  591.     }
  593.     if(m_MinSubjCoverage > 0.) {
  594.         // calculate subj coverage
  595.         double dSCov = CalcCoverage(m_Out.begin(), m_Out.end(), 's');
  596.         if(dSCov < m_MinSubjCoverage* nSubjSize)
  597.             m_Out.clear();
  598.     }
  599.     return 0;
  600. }
  603. void CHitParser::x_FilterByOrder(size_t offset, bool use_chainfilter)
  604. {
  605.     size_t dim = m_Out.size();
  606.     if(offset >= dim) return;
  607.     vector<CHit>::iterator ii_beg = m_Out.begin() + offset;
  609.     if(use_chainfilter) {  // groupwise chain filtering
  610.       /*
  611.         sort(ii_beg, m_Out.end(), hit_groupid_less());
  612.         int nGroupID = m_Out[0].m_GroupID;
  613.         size_t start = 0, fin = 0, target = 0;
  614.         while(fin < dim) {
  615.             for(; fin < dim && nGroupID == m_Out[fin].m_GroupID; ++fin);
  616.             start = fin;
  617.             nGroupID = m_Out[fin].m_GroupID;
  618.             vector<CHit> vh (ii_beg + start, ii_beg + fin);
  619.             CChainFilter cf (vh);
  620.             cf.Run();
  621.             copy(vh.begin(), vh.end(), ii_beg + target);
  622.             target += vh.size();
  623.         }
  624.         m_Out.resize(target);
  625.       */
  626.     }
  627.     else { // greedy filtering
  628.         vector<CHit>::iterator ii2 = ii_beg, iend2;
  629.         //for(iend2 = m_Out.end(); ii2 < iend2; ++ii2) {
  630.         //    iend2 = remove_if(ii2 + 1, iend2,
  631.         //                     bind1st(hit_not_same_order(), *ii2));
  632.         //}
  633.         iend2 = remove_if(ii2 + 1, m_Out.end(),
  634.                              bind1st(hit_not_same_order(), *ii2));
  635.         m_Out.erase(iend2, m_Out.end());
  636.     }
  637. }
  640. // Calculates proximity btw any two hits.
  641. // Return value can range from 0 (similar) to 1 (different)
  642. double CHitParser::x_CalculateProximity(const CHit& h1, const CHit& h2) const
  643. {
  644.     //    if(m_SameOrder)
  645.     //    {
  646.     //        hit_not_same_order hnso;
  647.     //        if(hnso(h1,h2)) return 1.;
  648.     //    }
  649.     double adm [] = { m_ange[1] - m_ange[0] + 1, m_ange[3] - m_ange[2] + 1};
  650.     double dm = max(adm[0], adm[1]);
  651.     double ad [2];
  652.     for(int i = 0; i < 4; i += 2) {
  653.         int i0 = i, i1 = i + 1;
  654.         if(h1.m_ai[i1] <= h2.m_ai[i0])  // no intersection
  655.             ad[i/2] = double(h2.m_ai[i0] - h1.m_ai[i1] - 1) / dm;// adm[i/2];
  656.         else
  657.         if(h2.m_ai[i1] <= h1.m_ai[i0])  // no intersection
  658.             ad[i/2] = double(h1.m_ai[i0] - h2.m_ai[i1] - 1) / dm;//adm[i/2];
  659.         else // inclusion
  660.         if(h1.m_ai[i0] <= h2.m_ai[i0] && h2.m_ai[i1] <= h1.m_ai[i1] ||
  661.              h2.m_ai[i0] <= h1.m_ai[i0] && h1.m_ai[i1] <= h2.m_ai[i1] )
  662.             ad[i/2] = 1.0;
  663.         else // intersection
  664.         if(h2.m_ai[i0] <= h1.m_ai[i0] && h1.m_ai[i0] <= h2.m_ai[i1])
  665.             ad[i/2] = double(h2.m_ai[i1] - h1.m_ai[i0] + 1) /
  666.                 (h1.m_ai[i1] - h2.m_ai[i0] + 1);
  667.         else
  668.             ad[i/2] = double(h1.m_ai[i1] - h2.m_ai[i0] + 1) /
  669.                 (h2.m_ai[i1] - h1.m_ai[i0] + 1);
  670.     }
  671.     return max(ad[0],ad[1]);
  672. }
  675. int CHitParser::x_RunMSGS(bool bSelectGroupsOnly, int& nGroupId)
  676. {
  677.     if( m_group_identification_method != eNone || bSelectGroupsOnly )
  678.     {
  679.         if(m_group_identification_method == eQueryCoverage) {
  680.             x_GroupsIdentifyByCoverage(0, m_Out.size(), nGroupId, 0, 'q');
  681.         }
  682.         else {
  683.             x_GroupsIdentifyByCoverage(0, m_Out.size(), nGroupId, 0, 's');
  684.         }
  685.         x_SyncGroupsByMaxDist(nGroupId);
  686.     }
  688.     if(m_OutputAllGroups && !bSelectGroupsOnly && m_Out.size())
  689.     {
  690.         // make copies to m_Out here
  691.         // new subj string must reflect both number and strand
  692.         // like [s4], [i8]
  694.         // first we count the number of groups
  695.         size_t nGroupCount = 1;
  696.         {{
  697.         sort(m_Out.begin(), m_Out.end(), hit_groupid_less());
  698.         vector<CHit>::iterator ibegin = m_Out.begin(), iend = m_Out.end(),
  699.             ii = ibegin;
  700.         int ngid = ii->m_GroupID;
  701.         for(++ii; ii != iend; ++ii)
  702.         {
  703.             if(ii->m_GroupID != ngid)
  704.             {
  705.                 ngid = ii->m_GroupID;
  706.                 ++nGroupCount;
  707.             }
  708.         }
  709.         }}
  711.         {{
  712.             // determine strands of the groups
  713.             // 1 = straight, 0 = inverse, -1 = both / undefined
  714.             vector<char> vStrands (nGroupCount, -1);
  715.             vector<CHit>::iterator ibegin = m_Out.begin(),
  716.                 ii = ibegin, iend = m_Out.end();
  717.             for(size_t ig = 0; ig < nGroupCount; ++ig)
  718.             {
  719.                 vStrands[ig] = ii->IsStraight()? 1: 0;
  720.                 int ngid = ii->m_GroupID;
  721.                 for(; ii != iend; ii++)
  722.                 {
  723.                     if(ii->m_GroupID != ngid) break;
  724.                     char cs = ii->IsStraight()? 1: 0;
  725.                     if(cs != vStrands[ig]) vStrands[ig] = -1;
  726.                 }
  727.             }
  728.             
  729.             // now separate the groups
  730.             ii = ibegin;
  731.             string strSubject0 = ii->m_Subj;
  732.             int ngc_s = 0, ngc_i = 0;
  733.             for(size_t ig = 0; ig < nGroupCount; ig++)
  734.             {
  735.                 int ngid = ii->m_GroupID;
  736. CNcbiOstrstream oss; // new subject
  737.                 char c0 = (vStrands[ig] == -1)? 'x':
  738.                     (vStrands[ig] == 0)? 'm': 'p';
  739.                 int ngc = (c0 == 'i')? ++ngc_i: ++ngc_s;
  740. oss << strSubject0 << "_[" << c0 << ngc << ']';
  741. const string new_subj = CNcbiOstrstreamToString(oss);
  742.                 for(; ii != iend; ++ii)
  743.                 {
  744.                     if(ii->m_GroupID != ngid) break;
  745.                     ii->m_Subj = new_subj;
  746.                 }
  747.             }
  748.         }}
  750.         return 0; // we don't resolve ambiguities in this mode
  751.     }
  753.     if(bSelectGroupsOnly) return 0;
  755.     vector<CHit> vh (m_Out);
  756.     vector<CHit> vh_out;
  758.     double max_group_score = 0;
  759.     size_t i = 0, dim = vh.size();
  760.     while(i < dim) {
  761.         int groupid0 = vh[i].m_GroupID;
  762.         int istart = i;
  763.         for(; i < dim && groupid0 == vh[i].m_GroupID; ++i);
  764.         int ifin = i;
  765.         m_Out.resize(ifin - istart);
  766.         copy(vh.begin() + istart, vh.begin() + ifin, m_Out.begin());
  767.         x_RunMaxScore();
  768.         double group_score = 0;
  769.         size_t dim_group = m_Out.size();
  770.         for(size_t k = 0; k < dim_group; ++k) group_score += m_Out[k].m_Score;
  771.         
  772.         if(m_OutputAllGroups) {
  773.             copy(m_Out.begin(), m_Out.end(), back_inserter(vh_out));
  774.         }
  775.         else if(group_score > max_group_score) {
  776.             max_group_score = group_score;
  777.             vh_out.clear();
  778.             copy(m_Out.begin(), m_Out.end(), back_inserter(vh_out));
  779.         }
  780.     }
  782.     m_Out = vh_out;
  784.     return 0;
  785. }
  788. void CHitParser::x_GroupsIdentifyByCoverage(int iStart, int iStop, int& iGroupId,
  789.                                           double dTotalCoverage, char where)
  790. {
  791.     bool bTopLevelCall = (iStart == 0) && (size_t(iStop) == m_Out.size());
  793.     if(bTopLevelCall) { // top-level call
  794.         x_CalcGlobalEnvelop();
  795.         if(where == 'q') {
  796.             stable_sort(m_Out.begin() + iStart, m_Out.begin() + iStop,
  797.                         CHit::PPrecedeS);
  798.         }
  799.         else {
  800.             stable_sort(m_Out.begin() + iStart, m_Out.begin() + iStop,
  801.                         CHit::PPrecedeQ);
  802.         }
  803.         dTotalCoverage = CalcCoverage(m_Out.begin() + iStart,
  804.                                       m_Out.begin() + iStop, where);
  805.     }
  807.     // assign new groupid
  808.     {{
  809.         ++iGroupId;
  810.         for(int i = iStart; i < iStop; ++i)
  811.             m_Out[i].m_GroupID = iGroupId;
  812.     }}
  815.     if(iStop - iStart <= 1) return;
  817.     // choose the best location
  818.     double dMax = dTotalCoverage;
  819.     int iMax = iStop;
  820.     unsigned char idx1, idx2;
  821.     if(where == 'q') {
  822.         idx1 = 2;
  823.         idx2 = 3;
  824.     }
  825.     else {
  826.         idx1 = 0;
  827.         idx2 = 1;
  828.     }
  829.     // find all gaps, sort them by size, then 
  830.     // proceed from larger to lower
  831.     typedef multimap<int, int> mm_t;
  832.     mm_t mm_gapsize2index;
  833.     for(int i = iStart; i < iStop - 1; i++)
  834.     {
  835.         int gap = m_Out[i+1].m_ai[idx1] - m_Out[i].m_ai[idx2];
  836.         if(gap > 0) {
  837.             mm_gapsize2index.insert(mm_t::value_type(gap, i));
  838.         }
  839.     }
  841.     mm_t::reverse_iterator mmirb = mm_gapsize2index.rbegin();
  842.     mm_t::reverse_iterator mmire = mm_gapsize2index.rend();
  843.     mm_t::reverse_iterator mmir;
  844.     bool match = false;
  845.     for(mmir = mmirb; mmir != mmire; ++mmir) {
  846.         int i = mmir->second;
  847.         int n1 = CalcCoverage(m_Out.begin() + iStart,
  848.                               m_Out.begin() + i + 1, where);
  849.         int n2 = CalcCoverage(m_Out.begin() + i + 1,
  850.                               m_Out.begin() + iStop, where);
  851.         double d = (n1 + n2);
  852.         match = (d - dTotalCoverage) / dTotalCoverage  >= m_CovStep;
  853.         if(match) {
  854.              dMax = d;
  855.              iMax = i;
  856.              break;
  857.          }
  858.     }
  860.     if(match)
  861.     {
  862.         x_GroupsIdentifyByCoverage(iStart, iMax + 1, iGroupId,
  863.    dTotalCoverage, where);
  864.         x_GroupsIdentifyByCoverage(iMax + 1, iStop, iGroupId,
  865.    dTotalCoverage, where);
  866.     }
  867. }
  870. size_t CHitParser::x_RemoveEqual()
  871. {
  872.     typedef map<size_t, vector<size_t> > map_t;
  873.     map_t m;  // hit checksum to hits
  875.     vector<size_t> vhd;  // hits to delete
  876.     size_t dim = m_Out.size();
  877.     for(size_t i = 0; i < dim ; ++i) {
  878.         const CHit& h = m_Out[i];
  879.         size_t checksum = ( h.m_an[0] + h.m_an[1] + h.m_an[2] + h.m_an[3] )
  880.                           % 1000;
  881.         map_t::iterator mi = m.find(checksum);
  882.         if(mi == m.end()) {
  883.             m[checksum].push_back(i);
  884.         }
  885.         else {
  886.             vector<size_t>& vph = mi->second;
  887.             size_t dim_hits = vph.size(), k = 0;
  888.             for(; k < dim_hits; ++k) {
  889.                 const CHit& h2 = m_Out[vph[k]];
  890.                 if( h.m_an[0] == h2.m_an[0] && h.m_an[1] == h2.m_an[1] &&
  891.                     h.m_an[2] == h2.m_an[2] && h.m_an[3] == h2.m_an[3]    ) {
  893.                     vhd.push_back(i);
  894.                     break;
  895.                 }
  896.             }
  897.             if(k == dim_hits) {
  898.                     m[checksum].push_back(i);
  899.             }
  900.         }
  901.     }
  903.     size_t del_count = 0;
  904.     if(vhd.size()) {
  905.         sort(vhd.begin(), vhd.end());
  906.         size_t dim_vhd = vhd.size();
  907.         size_t vhd_i = 0;
  908.         for(size_t k = vhd[vhd_i]; k < dim - del_count; ++k) {
  909.             bool b = false;
  910.             if(vhd_i < dim_vhd && k == vhd[vhd_i] - del_count) {
  911.                 ++del_count;
  912.                 ++vhd_i;
  913.                 b = true;
  914.             }
  915.             size_t k1 = k + del_count;
  916.             if(k1 < dim) {
  917.                 m_Out[k] = m_Out[k1];
  918.             }
  919.             else  break;
  920.             if(b) --k;
  921.         }
  922.         m_Out.resize(dim - del_count);
  923.     }
  925.     return del_count;
  926. }
  929. int CHitParser::Run(EMode erm)
  930. {
  931.     int nCode = x_CheckParameters();
  932.     if(nCode) return nCode;
  934.     vector<CHit> vhInvStrand;      // used only if strand mode = auto
  935.     int& nGroupId = *m_GroupID;  // auxiliary variable
  936.     size_t dim = m_Out.size();
  937.     if(dim)
  938.     {
  939.         x_TransformCoordinates(true);
  940.         x_Combine(m_CombinationProximity_pre);
  942.         {{
  944.         // filter by strand or separate strands, if applicable
  945.         if(m_Strand == eStrand_Plus || m_Strand == eStrand_Minus ||
  946.            m_Strand == eStrand_Auto)
  947.         {
  948.             vector<CHit>::iterator ii;
  949.             for(ii = m_Out.end()-1; ii >= m_Out.begin(); ii--)
  950.             {
  951.                 bool bStraight = ii->IsStraight();
  952.                 switch(m_Strand)
  953.                 {
  954.                     case eStrand_Plus:
  955.                         if (bStraight) continue; else m_Out.erase(ii);
  956.                         break;
  957.                     case eStrand_Minus:
  958.                         if (!bStraight) continue; else m_Out.erase(ii);
  959.                         break;
  960.                     case eStrand_Auto: //  separate strands before processing
  961.                         if (bStraight) continue;
  962.                         else
  963.                         {
  964.                             vhInvStrand.push_back(*ii);
  965.                             m_Out.erase(ii);
  966.                         }
  967.                         break;
  968.                 }
  969.             }
  970.         }
  972.         x_CalcGlobalEnvelop();
  974.         switch(erm)
  975.         {
  976.             case eMode_Combine: // already combined
  977.                 break;
  978.             case eMode_GroupSelect:
  979.                 x_RunMSGS(true, nGroupId);
  980.                 break;
  981.             case eMode_Normal:
  982.                 x_IdentifyMaxDistGroups();
  983.                 switch(m_Method)
  984.                 {
  985.                     case eMethod_MaxScore:
  986.                         x_RunMaxScore(); break;
  987.                     case eMethod_MaxScore_GroupSelect:
  988.                         x_RunMSGS(false, nGroupId); break;
  989.                     default: return 1; //  not supported
  990.                 }
  991.                 if(!m_OutputAllGroups) {
  992.                     // we need to identify again because "connecting"
  993.                     // hits might might have perished
  994.                     x_IdentifyMaxDistGroups();  
  995.                     x_FilterByMaxDist();
  996.                 }
  997.                 break;
  998.         }
  1001.         // if strand mode is auto then we run it again for inverse strand
  1002.         // and then compare the strands
  1003.         if(m_Strand == 3)
  1004.         {
  1005.             vector<CHit> vh0 = m_Out;
  1006.             m_Out = vhInvStrand;
  1007.             x_CalcGlobalEnvelop();
  1008.             switch(erm)
  1009.             {
  1010.                 case eMode_Combine:
  1011.                     break;
  1012.                 case eMode_GroupSelect:
  1013.                     x_RunMSGS(true, nGroupId);
  1014.                     break;
  1015.                 case eMode_Normal:
  1016.                     x_IdentifyMaxDistGroups();
  1017.                     switch(m_Method)
  1018.                     {
  1019.                         case eMethod_MaxScore:
  1020.                             x_RunMaxScore(); break;
  1021.                         case eMethod_MaxScore_GroupSelect:
  1022.                             x_RunMSGS(false, nGroupId); break;
  1023.                         default: return 1; //  not supported
  1024.                     }
  1025.                     x_IdentifyMaxDistGroups();
  1026.                     x_FilterByMaxDist();
  1027.                     break;
  1028.             }
  1031.             if(erm == eMode_Normal && !m_OutputAllGroups)
  1032.             { // leave only the best strand
  1033.                 double dScoreStraight = 0., dScoreInverse = 0.;
  1034.                 vector<CHit>::iterator ii = vh0.begin(), iend = vh0.end();
  1035.                 for(; ii != iend; ii++) dScoreStraight += ii->m_Score;
  1036.                 ii = m_Out.begin(); iend = m_Out.end();
  1037.                 for(; ii != iend; ii++) dScoreInverse += ii->m_Score;
  1038.                 if(dScoreStraight > dScoreInverse) m_Out = vh0;
  1039.             }
  1040.             else
  1041.             { // merge the strands since we don't actually want to
  1042.               // eliminate ambiguities here
  1043.               copy(vh0.begin(), vh0.end(), back_inserter(m_Out));
  1044.             }
  1045.         }
  1047.         x_Combine(m_CombinationProximity_post);
  1049.         }}
  1051.         x_TransformCoordinates(false);
  1052.     }
  1054.     return nCode;
  1055. }
  1058. int CHitParser::x_CheckParameters() const
  1059. {
  1060.     switch (m_Method)
  1061.     {
  1062.         case eMethod_MaxScore:
  1063.         case eMethod_MaxScore_GroupSelect:
  1064.         break;
  1065.         default: return 1; // unsupported
  1066.     };
  1067.     switch (m_Strand)
  1068.     {
  1069.         case 0: case 1: case 3: break;
  1070.             // same order and different strands are incompatible
  1071.         case 2: if(m_SameOrder) return 3;
  1072.                         else break;
  1073.         default: return 2; // unsupported
  1074.     };
  1075.     return 0;
  1076. }
  1079. // counts the number of collisions
  1080. // pvh0== 0 => count at m_Out
  1081. int CHitParser::GetCollisionCount (int& nq, int& ns, const vector<CHit>& vh)
  1082. {
  1083.     nq = ns = 0;
  1084.     vector<SNode> vi_q, vi_s;
  1085.     vector<CHit>::const_iterator phi0 = vh.begin(),
  1086.       phi1 = vh.end(), phi = phi0;
  1087.     for(phi = phi0; phi != phi1; ++phi)
  1088.     {
  1089.         SNode Nodes [] = {
  1090.             SNode(phi->m_ai[0],true,phi), SNode(phi->m_ai[1],false,phi),
  1091.             SNode(phi->m_ai[2],true,phi), SNode(phi->m_ai[3],false,phi)
  1092.         };
  1093.         vi_q.push_back(Nodes[0]); vi_q.push_back(Nodes[1]);
  1094.         vi_s.push_back(Nodes[2]); vi_s.push_back(Nodes[3]);
  1095.     }
  1096.     sort(vi_q.begin(),vi_q.end());
  1097.     sort(vi_s.begin(),vi_s.end());
  1098.     int nLevel = 0;
  1099.     vector<SNode>::iterator ii;
  1100.     for(ii = vi_q.begin(); ii != vi_q.end(); ii++)
  1101.         if(ii->type) { if(++nLevel == 2) nq++; } else nLevel--;
  1102.     nLevel = 0;
  1103.     vector<SNode>::iterator jj;
  1104.     for(jj = vi_s.begin(); jj != vi_s.end(); jj++)
  1105.         if(jj->type) { if(++nLevel == 2) ns++; } else nLevel--;
  1106.     return nq + ns;
  1107. }
  1110. void CHitParser::x_TransformCoordinates(bool bDir)
  1111. {
  1112.     if(bDir)
  1113.     {
  1114.         // transform hits
  1115.         m_Origin[0] = m_ange[0];
  1116.         m_Origin[1] = m_ange[2];
  1117.         vector<CHit>::iterator ii, iend;
  1118.         for(ii = m_Out.begin(), iend = m_Out.end(); ii != iend; ii++)
  1119.         {
  1120.             for(int k=0; k<4; k++)
  1121.             {
  1122.                 ii->m_an[k] -= m_Origin[k/2];
  1123.                 ii->m_ai[k] -= m_Origin[k/2];
  1124.             }
  1125.         }
  1126.     }
  1127.     else
  1128.     {
  1129.         vector<CHit>::iterator ii, iend;
  1130.         vector<vector<CHit>::iterator> vhi_del;
  1131.         for(ii = m_Out.begin(), iend = m_Out.end(); ii != iend; ii++)
  1132.         {
  1133.             for(int k=0; k<4; k++)
  1134.             {
  1135.                 ii->m_an[k] += m_Origin[k/2];
  1136.                 ii->m_ai[k] += m_Origin[k/2];
  1137.             }
  1138.             if(ii->m_ai[1] == ii->m_ai[0] || ii->m_ai[3] == ii->m_ai[2])
  1139.                 vhi_del.push_back(ii);
  1140.         }
  1141.         // delete the slack
  1142.         vector<vector<CHit>::iterator>::reverse_iterator i1, i1e;
  1143.         for(i1 = vhi_del.rbegin(), i1e = vhi_del.rend(); i1 != i1e; i1++)
  1144.             m_Out.erase(*i1);
  1145.     }
  1146.     x_CalcGlobalEnvelop();
  1147. }
  1150. void CHitParser::x_FilterByMaxDist()
  1151. {
  1152.     //  sum and compare the scores
  1153.     map<int,double> mgr2sc;
  1154.     int nSize = m_Out.size();
  1155.     for(int j = 0; j < nSize; j++)
  1156.     {
  1157.         const int nClustID = CMaxDistClIdGet()(m_Out[j]);
  1158.         mgr2sc[nClustID] += m_Out[j].m_Score;
  1159.     }
  1160.     double dmax = -1e38;
  1161.     map<int,double>::iterator im, im1;
  1162.     for(im1 = im = mgr2sc.begin(); im != mgr2sc.end(); im++)
  1163.     {
  1164.         if(im->second > dmax)
  1165.         {
  1166.             dmax = im->second;
  1167.             im1 = im;
  1168.         }
  1169.     }
  1171.     // leave only the top-score max dist group
  1172.     vector<CHit>::iterator ie = remove_if(
  1173.         m_Out.begin(), m_Out.end(), bind2nd(hit_out_of_group(),im1->first));
  1174.     m_Out.erase(ie,m_Out.end());
  1175. }
  1178. void CHitParser::x_IdentifyMaxDistGroups()
  1179. {
  1180.     DoSingleLinkageClustering ( m_Out,
  1181.         CMaxDistClustPred(m_MaxHitDistQ, m_MaxHitDistS),
  1182.         CMaxDistClIdSet(),
  1183.         CMaxDistClIdGet() );
  1184. }
  1187. void CHitParser::x_SyncGroupsByMaxDist(int& nGroupID)
  1188. {
  1189.     if(m_Out.size() == 0) return;
  1191.     // order by GroupID
  1192.     stable_sort(m_Out.begin(), m_Out.end(), hit_groupid_less());
  1194.     // renumber the groups to make sure that
  1195.     // we have same nClustID throughout the group
  1196.     vector<CHit>::iterator ibegin = m_Out.begin(), iend = m_Out.end(),
  1197.         ii = ibegin;
  1198.     int ngid = ii->m_GroupID, nmdcid = ii->m_MaxDistClustID;
  1199.     ii->m_GroupID = ++nGroupID;
  1200.     for(++ii; ii != iend; ++ii)
  1201.     {
  1202.         if(ngid == ii->m_GroupID && nmdcid != ii->m_MaxDistClustID)
  1203.         {
  1204.             ++nGroupID;
  1205.             nmdcid = ii->m_MaxDistClustID;
  1206.         }
  1207.         if(ngid != ii->m_GroupID)
  1208.         {
  1209.             ++nGroupID;
  1210.             ngid = ii->m_GroupID;
  1211.         }
  1212.         ii->m_GroupID = nGroupID;
  1213.     }
  1214. }
  1217. size_t CHitParser::GetSeqLength(const string& accession)
  1218. {
  1219.     return 0;
  1221.     /*
  1222.   using namespace objects;
  1224.   map<string, size_t>::const_iterator i1 = m_seqsizes.find(accession),
  1225.     i2 = m_seqsizes.end();
  1226.   if(i1 != i2) {
  1227.     size_t s = i1->second;
  1228.     return s;
  1229.   }
  1230.   
  1231.   CObjectManager om;
  1232.   om.RegisterDataLoader ( *new CGBDataLoader("ID", new CId1Reader, 2),
  1233.   CObjectManager::eDefault );
  1234.   CScope scope(om);
  1235.   scope.AddDefaults();
  1236.     
  1237.   CSeq_id seq_id (accession);
  1238.   
  1239.   CBioseq_Handle bioseq_handle = scope.GetBioseqHandle(seq_id);
  1240.   
  1241.   if ( bioseq_handle ) {
  1242.     const CBioseq& bioseq = bioseq_handle.GetBioseq();
  1243.     const CSeq_inst& inst = bioseq.GetInst();
  1244.     if(inst.IsSetLength()) {
  1245.       size_t s = inst.GetLength();
  1246.       m_seqsizes[accession] = s;
  1247.       return s;
  1248.     }
  1249.   }
  1250.   return 0;
  1252.     */
  1253. }
  1256. bool CHitGroup::CheckSameOrder()
  1257. {
  1258.     if(m_vHits.size() == 0) return true;
  1259.     vector<int>::const_iterator ivh0 = m_vHits.begin(),
  1260.         ivh1 = m_vHits.end(), ii, jj;
  1262.     hit_not_same_order hnso;
  1263.     for(ii = ivh0; ii < ivh1; ++ii) {
  1264.         for(jj = ivh0; jj < ii; ++jj) {
  1265.             if(hnso( (*m_pHits)[*ii], (*m_pHits)[*jj] ))
  1266.                 return false;
  1267.         }
  1268.     }
  1270.     return true;
  1271. }
  1274. void CHitParser::x_Combine(double dProximity)
  1275. {
  1276.     if(dProximity < 0 || m_Out.size() == 0) {
  1277.         return;
  1278.     }
  1280.     x_CalcGlobalEnvelop();
  1282.     size_t total_hits = m_Out.size();
  1283.     vector<CHit*> vh (total_hits);
  1284.     CHit* phit = &(m_Out[0]);
  1285.     for(size_t i = 0; i < total_hits; ++i) {
  1286.         vh[i] = phit++;
  1287.     }
  1288.     typedef vector<CHit*>::iterator TVHitPtrI;
  1290.     // separate strands
  1291.     sort(vh.begin(), vh.end(),  CHit::PPrecedeStrand_ptr);
  1292.     TVHitPtrI strand_minus_start = 
  1293.         find_if(vh.begin(), vh.end(), bind2nd(hit_ptr_strand_is(), true));
  1295.     typedef pair<TVHitPtrI, TVHitPtrI> THitRange;
  1296.     THitRange strands [2];
  1297.     strands[0].first = vh.begin();
  1298.     strands[0].second = strands[1].first = strand_minus_start;
  1299.     strands[1].second = strands[0].first + total_hits;
  1300.  
  1301.     for(size_t strand = 0; strand < 2; ++strand) {
  1302.         TVHitPtrI start  = strands[strand].first;
  1303.         TVHitPtrI finish = strands[strand].second;
  1304.         
  1305.         size_t dim = finish - start;
  1306.         if(dim < 2) continue;
  1307.         size_t count;
  1308.         do {
  1309.             count = 0;
  1310.             stable_sort(start, finish, CHit::PGreaterScore_ptr);
  1312.             for(TVHitPtrI i1 = start; i1 < finish - 1; ++i1) {
  1313.                 if(*i1 == 0) continue;
  1314.                 for(TVHitPtrI j1 = i1 + 1; j1 < finish; ++j1) {
  1315.                     if(*j1 == 0) continue;
  1316.                     double d = x_CalculateProximity(**i1, **j1);
  1317.                     if(d <= dProximity) {
  1318.                         CHit hc (**i1, **j1);
  1319.                         **i1 = hc;
  1320.                         *j1 = 0;
  1321.                         ++count;
  1322.                     }
  1323.                 }
  1324.             }
  1325.         }
  1326.         while(count);
  1327.     }
  1329.     // flush deletions, compress pointers vector, fill original hit vector
  1330.     TVHitPtrI iv = remove(vh.begin(), vh.end(), (CHit*) 0);
  1331.     vh.erase(iv, vh.end());
  1332.     TVHitPtrI ib = vh.begin(), ie = vh.end();
  1333.     sort(ib, ie, CHit::PGreaterScore_ptr);
  1334.     vector<CHit> vout (vh.size());
  1335.     for(iv = ib; iv != ie; ++iv) {
  1336.         vout[iv - ib] = **iv;
  1337.     }
  1338.     m_Out = vout;
  1340.     return;
  1341. }
  1343. END_NCBI_SCOPE
  1345. /*
  1346. * ===========================================================================
  1347. *
  1348. * $Log: splign_hitparser.cpp,v $
  1349. * Revision 1000.0  2004/06/01 18:12:53  gouriano
  1350. * PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.7
  1351. *
  1352. * Revision 1.7  2004/05/24 16:13:57  gorelenk
  1353. * Added PCH ncbi_pch.hpp
  1354. *
  1355. * Revision 1.6  2004/05/18 21:43:40  kapustin
  1356. * Code cleanup
  1357. *
  1358. * Revision 1.5  2004/05/03 21:53:57  kapustin
  1359. * Eliminate OM-dependant code
  1360. *
  1361. * Revision 1.4  2004/04/26 16:52:44  ucko
  1362. * Add an explicit "typename" annotation required by GCC 3.4, and adjust
  1363. * the code to take advantage of the ITERATE macro.
  1364. *
  1365. * Revision 1.3  2004/04/23 14:37:44  kapustin
  1366. * *** empty log message ***
  1367. *
  1368. * 
  1369. * ===========================================================================
  1370. */