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

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

生物技术

开发平台：

C/C++

mb_lookup.c：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: mb_lookup.c,v $
  4.  * PRODUCTION Revision 1000.4  2004/06/01 18:08:17  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.33
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: mb_lookup.c,v 1000.4 2004/06/01 18:08:17 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 offical 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: Ilya Dondoshansky
  36.  *
  37.  */
  39. /** @file mb_lookup.c
  40.  * Functions responsible for the creation of a lookup table
  41.  * @todo FIXME @sa mb_lookup.h
  42.  */
  44. static char const rcsid[] = 
  45.     "$Id: mb_lookup.c,v 1000.4 2004/06/01 18:08:17 gouriano Exp $";
  47. #include <algo/blast/core/blast_options.h>
  48. #include <algo/blast/core/blast_def.h>
  49. #include <algo/blast/core/mb_lookup.h>
  50. #include "blast_inline.h"
  53. MBLookupTable* MBLookupTableDestruct(MBLookupTable* mb_lt)
  54. {
  55.    if (!mb_lt)
  56.       return NULL;
  57.    if (mb_lt->hashtable)
  58.       sfree(mb_lt->hashtable);
  59.    if (mb_lt->next_pos)
  60.       sfree(mb_lt->next_pos);
  61.    if (mb_lt->next_pos2)
  62.       sfree(mb_lt->next_pos2);
  63.    sfree(mb_lt->pv_array);
  65.    sfree(mb_lt);
  66.    return mb_lt;
  67. }
  69. /** Convert weight, template length and template type from input options into
  70.     an MBTemplateType enum
  71. */
  72. static DiscTemplateType GetDiscTemplateType(Int2 weight, Uint1 length, 
  73.                                      DiscWordType type)
  74. {
  75.    if (weight == 11) {
  76.       if (length == 16) {
  77.          if (type == MB_WORD_CODING || type == MB_TWO_TEMPLATES)
  78.             return TEMPL_11_16;
  79.          else if (type == MB_WORD_OPTIMAL)
  80.             return TEMPL_11_16_OPT;
  81.       } else if (length == 18) {
  82.          if (type == MB_WORD_CODING || type == MB_TWO_TEMPLATES)
  83.             return TEMPL_11_18;
  84.          else if (type == MB_WORD_OPTIMAL)
  85.             return TEMPL_11_18_OPT;
  86.       } else if (length == 21) {
  87.          if (type == MB_WORD_CODING || type == MB_TWO_TEMPLATES)
  88.             return TEMPL_11_21;
  89.          else if (type == MB_WORD_OPTIMAL)
  90.             return TEMPL_11_21_OPT;
  91.       }
  92.    } else if (weight == 12) {
  93.       if (length == 16) {
  94.          if (type == MB_WORD_CODING || type == MB_TWO_TEMPLATES)
  95.             return TEMPL_12_16;
  96.          else if (type == MB_WORD_OPTIMAL)
  97.             return TEMPL_12_16_OPT;
  98.       } else if (length == 18) {
  99.          if (type == MB_WORD_CODING || type == MB_TWO_TEMPLATES)
  100.             return TEMPL_12_18;
  101.          else if (type == MB_WORD_OPTIMAL)
  102.             return TEMPL_12_18_OPT;
  103.       } else if (length == 21) {
  104.          if (type == MB_WORD_CODING || type == MB_TWO_TEMPLATES)
  105.             return TEMPL_12_21;
  106.          else if (type == MB_WORD_OPTIMAL)
  107.             return TEMPL_12_21_OPT;
  108.       }
  109.    }
  110.    return TEMPL_CONTIGUOUS; /* All unsupported cases default to 0 */
  111. }
  113. #define SMALL_QUERY_CUTOFF 15000
  114. #define LARGE_QUERY_CUTOFF 800000
  116. /* Documentation in mb_lookup.h */
  117. Int2 MB_LookupTableNew(BLAST_SequenceBlk* query, ListNode* location,
  118.         MBLookupTable** mb_lt_ptr,
  119.         const LookupTableOptions* lookup_options)
  120. {
  121.    Int4 query_length;
  122.    Uint1* seq,* pos;
  123.    Int4 index;
  124.    Int4 ecode;
  125.    Int4 mask;
  126.    Int4 ecode1, ecode2;
  127.    Uint1 val, nuc_mask = 0xfc;
  128.    MBLookupTable* mb_lt;
  129.    Int4 masked_word_count = 0;
  130.    PV_ARRAY_TYPE *pv_array=NULL;
  131.    Int4 pv_shift, pv_array_bts, pv_size, pv_index;
  132.    Int2 word_length, extra_length;
  133.    Int4 last_offset;
  134.    Int4 table_entries;
  135. #ifdef USE_HASH_TABLE
  136.    Int4 hash_shift, hash_mask, crc, size, length;
  137.    Uint1* hash_buf;
  138. #endif
  139.    DiscTemplateType template_type;
  140.    Boolean amb_cond;
  141.    Boolean two_templates = 
  142.       (lookup_options->mb_template_type == MB_TWO_TEMPLATES);
  143.    Int2 bytes_in_word;
  144.    Uint1 width;
  145.    Int4 from, to;
  146.    ListNode* loc;
  147.    Int4 longest_chain = 0;
  148.    Int4 pcount, pcount1=0;
  149.    
  150.    if (!location) {
  151.      /* Empty sequence location provided */
  152.      *mb_lt_ptr = NULL;
  153.      return -1;
  154.    }
  156.    template_type = GetDiscTemplateType(lookup_options->word_size,
  157.                       lookup_options->mb_template_length, 
  158.                       (DiscWordType)lookup_options->mb_template_type);
  159.    query_length = query->length;
  160.    mb_lt = (MBLookupTable*) calloc(1, sizeof(MBLookupTable));
  161.     
  162.    bytes_in_word = (lookup_options->word_size + 1)/ 4;
  163.    if (bytes_in_word < 3) 
  164.       width = 2;
  165.    else
  166.       width = 3;
  168.    mb_lt->template_type = template_type;
  169.    mb_lt->two_templates = two_templates;
  170.    if (two_templates)
  171.       /* For now leave only one possibility for the second template */
  172.       mb_lt->second_template_type = template_type + 1;
  174. #ifdef USE_HASH_TABLE
  175.    /* Make hash table size the smallest power of 2 greater than query length */
  176.    for (length=query_length,size=1; length>0; length=length>>1,size++);
  177.    mb_lt->hashsize = 1<<size;
  178.    hash_shift = (32 - size)/2;
  179.    hash_mask = mb_lt->hashsize - 1;
  180.    pv_shift = 0;
  182. #else
  184.    if (width == 2) {
  185.       mb_lt->hashsize = 1 << 16;
  186.       pv_shift = 0;
  187.    }
  188.    else {
  190.       /* determine the approximate number of hashtable entries */
  191.       table_entries = 0;
  192.       for (loc = location; loc; loc = loc->next) {
  193.          from = ((SSeqRange*) loc->ptr)->left;
  194.          to = ((SSeqRange*) loc->ptr)->right;
  195.          table_entries += (to - from);
  196.       }
  198.       /* To fit in the external cache of latter-day micro-
  199.          processors, the PV array must be compressed. pv_shift
  200.  below is the power of two that the array size is
  201.  divided by. The target PV array size is 128 kBytes.
  203.  If the query is too small or too large, the compression 
  204.  should be higher. Small queries don't reuse the PV array,
  205.  and large queries saturate it. In either case, cache
  206.  is better used on other data. */
  208.       if (lookup_options->word_size == 11 && lookup_options->mb_template_length > 0) {
  209.          mb_lt->hashsize = 1 << 22;
  210.          if( table_entries <= SMALL_QUERY_CUTOFF ||
  211.      table_entries >= LARGE_QUERY_CUTOFF )
  212.             pv_shift = 3;
  213.  else
  214.             pv_shift = 2;
  215.       }
  216.       else {
  217.          mb_lt->hashsize = 1 << 24;
  218.          if( table_entries <= SMALL_QUERY_CUTOFF ||
  219.      table_entries >= LARGE_QUERY_CUTOFF )
  220.             pv_shift = 5;
  221.  else
  222.             pv_shift = 4;
  223.       }
  224.    }
  225. #endif
  227.    if (lookup_options->mb_template_length > 0) {
  228.       mb_lt->discontiguous = TRUE;
  229.       mb_lt->compressed_wordsize = width + 1;
  230.       word_length = mb_lt->word_length = COMPRESSION_RATIO*(width+1);
  231.       mb_lt->template_length = lookup_options->mb_template_length;
  232.       extra_length = mb_lt->template_length - word_length;
  233.       mb_lt->mask = (Int4) (-1);
  234.       mask = (1 << (8*(width+1) - 2)) - 1;
  235.    } else {
  236.       mb_lt->compressed_wordsize = width;
  237.       word_length = COMPRESSION_RATIO*width;
  238.       mb_lt->word_length = lookup_options->word_size;
  239.       mb_lt->template_length = COMPRESSION_RATIO*width;
  240.       mb_lt->mask = mb_lt->hashsize - 1;
  241.       mask = (1 << (8*width - 2)) - 1;
  242.       extra_length = 0;
  243.    }
  245.    /* Scan step must be set in lookup_options; even if it is not, still 
  246.       try to set it reasonably in mb_lt - it can't remain 0 */
  247.    if (lookup_options->scan_step)
  248.       mb_lt->scan_step = lookup_options->scan_step;
  249.    else
  250.       mb_lt->scan_step = COMPRESSION_RATIO;
  252.    mb_lt->full_byte_scan = (mb_lt->scan_step % COMPRESSION_RATIO == 0);
  254.    if ((mb_lt->hashtable = (Int4*) 
  255.         calloc(mb_lt->hashsize, sizeof(Int4))) == NULL) {
  256.       MBLookupTableDestruct(mb_lt);
  257.       return -1;
  258.    }
  260.    if ((mb_lt->next_pos = (Int4*) 
  261.         calloc((query_length+1), sizeof(Int4))) == NULL) {
  262.       MBLookupTableDestruct(mb_lt);
  263.       return -1;
  264.    }
  266.    if (two_templates) {
  267.       if ((mb_lt->hashtable2 = (Int4*) 
  268.            calloc(mb_lt->hashsize, sizeof(Int4))) == NULL) {
  269.          MBLookupTableDestruct(mb_lt);
  270.          return -1;
  271.       }
  272.       if ((mb_lt->next_pos2 = (Int4*) 
  273.            calloc((query_length+1), sizeof(Int4))) == NULL) {
  274.          MBLookupTableDestruct(mb_lt);
  275.          return -1;
  276.       }
  277.    }
  279.    for (loc = location; loc; loc = loc->next) {
  280.       /* We want index to be always pointing to the start of the word.
  281.          Since sequence pointer points to the end of the word, subtract
  282.          word length from the loop boundaries. 
  283.       */
  284.       from = ((SSeqRange*) loc->ptr)->left;
  285.       to = ((SSeqRange*) loc->ptr)->right - word_length;
  287.       seq = query->sequence_start + from;
  288.       pos = seq + word_length;
  289.       
  290.       ecode = 0;
  291.       /* Also add 1 to all indices, because lookup table indices count 
  292.          from 1. */
  293.       from -= word_length - 2;
  294.       last_offset = to - extra_length + 2;
  295.       amb_cond = TRUE;
  297.       for (index = from; index <= last_offset; index++) {
  298.          val = *++seq;
  299.          if ((val & nuc_mask) != 0) { /* ambiguity or gap */
  300.             ecode = 0;
  301.             pos = seq + word_length;
  302.          } else {
  303.             /* get next base */
  304.             ecode = ((ecode & mask) << 2) + val;
  305.             if (seq >= pos) {
  306.                if (extra_length) {
  307.                   switch (template_type) {
  308.                   case TEMPL_11_18: case TEMPL_12_18:
  309.                      amb_cond = !GET_AMBIG_CONDITION_18(seq);
  310.                      break;
  311.                   case TEMPL_11_18_OPT: case TEMPL_12_18_OPT:
  312.                      amb_cond = !GET_AMBIG_CONDITION_18_OPT(seq);
  313.                      break;
  314.                   case TEMPL_11_21: case TEMPL_12_21:
  315.                      amb_cond = !GET_AMBIG_CONDITION_21(seq);
  316.                      break;
  317.                   case TEMPL_11_21_OPT: case TEMPL_12_21_OPT:
  318.                      amb_cond = !GET_AMBIG_CONDITION_21_OPT(seq);
  319.                      break;
  320.                   default:
  321.                      break;
  322.                   }
  323.                }
  324.                      
  325.                if (amb_cond) {
  326.                   switch (template_type) {
  327.                   case TEMPL_11_16:
  328.                      ecode1 = GET_WORD_INDEX_11_16(ecode);
  329.                      break;
  330.                   case TEMPL_12_16:
  331.                      ecode1 = GET_WORD_INDEX_12_16(ecode);
  332.                      break;
  333.                   case TEMPL_11_16_OPT:
  334.                      ecode1 = GET_WORD_INDEX_11_16_OPT(ecode);
  335.                      break;
  336.                   case TEMPL_12_16_OPT:
  337.                      ecode1 = GET_WORD_INDEX_12_16_OPT(ecode);
  338.                      break;
  339.                   case TEMPL_11_18:
  340.                      ecode1 = (GET_WORD_INDEX_11_18(ecode)) |
  341.                                (GET_EXTRA_CODE_18(seq));
  342.                      break;
  343.                   case TEMPL_12_18:
  344.                      ecode1 = (GET_WORD_INDEX_12_18(ecode)) |
  345.                                (GET_EXTRA_CODE_18(seq));
  346.                      break;
  347.                   case TEMPL_11_18_OPT:
  348.                      ecode1 = (GET_WORD_INDEX_11_18_OPT(ecode)) |
  349.                                (GET_EXTRA_CODE_18_OPT(seq));
  350.                      break;
  351.                   case TEMPL_12_18_OPT:
  352.                      ecode1 = (GET_WORD_INDEX_12_18_OPT(ecode)) |
  353.                                (GET_EXTRA_CODE_18_OPT(seq));
  354.                      break;
  355.                   case TEMPL_11_21:
  356.                      ecode1 = (GET_WORD_INDEX_11_21(ecode)) |
  357.                                (GET_EXTRA_CODE_21(seq));
  358.                      break;
  359.                   case TEMPL_12_21:
  360.                      ecode1 = (GET_WORD_INDEX_12_21(ecode)) |
  361.                                (GET_EXTRA_CODE_21(seq));
  362.                      break;
  363.                   case TEMPL_11_21_OPT:
  364.                      ecode1 = (GET_WORD_INDEX_11_21_OPT(ecode)) |
  365.                                (GET_EXTRA_CODE_21_OPT(seq));
  366.                      break;
  367.                   case TEMPL_12_21_OPT:
  368.                      ecode1 = (GET_WORD_INDEX_12_21_OPT(ecode)) |
  369.                                (GET_EXTRA_CODE_21_OPT(seq));
  370.                      break;
  371.                   default: /* Contiguous word */
  372.                      ecode1 = ecode;
  373.                      break;
  374.                   }
  375.                      
  376. #ifdef USE_HASH_TABLE                     
  377.                   hash_buf = (Uint1*)&ecode1;
  378.                   CRC32(crc, hash_buf);
  379.                   ecode1 = (crc>>hash_shift) & hash_mask;
  380. #endif
  382.                   if (two_templates) {
  383.                      switch (template_type) {
  384.                      case TEMPL_11_16:
  385.                         ecode2 = GET_WORD_INDEX_11_16_OPT(ecode);
  386.                         break;
  387.                      case TEMPL_12_16:
  388.                         ecode2 = GET_WORD_INDEX_12_16_OPT(ecode);
  389.                         break;
  390.                      case TEMPL_11_18:
  391.                         ecode2 = (GET_WORD_INDEX_11_18_OPT(ecode)) |
  392.                            (GET_EXTRA_CODE_18_OPT(seq));
  393.                         break;
  394.                      case TEMPL_12_18:
  395.                         ecode2 = (GET_WORD_INDEX_12_18_OPT(ecode)) |
  396.                            (GET_EXTRA_CODE_18_OPT(seq));
  397.                         break;
  398.                      case TEMPL_11_21:
  399.                         ecode2 = (GET_WORD_INDEX_11_21_OPT(ecode)) |
  400.                            (GET_EXTRA_CODE_21_OPT(seq));
  401.                         break;
  402.                      case TEMPL_12_21:
  403.                         ecode2 = (GET_WORD_INDEX_12_21_OPT(ecode)) |
  404.                            (GET_EXTRA_CODE_21_OPT(seq));
  405.                         break;
  406.                      default:
  407.                         ecode2 = 0; break;
  408.                      }
  409.                      
  410. #ifdef USE_HASH_TABLE                     
  411.                      hash_buf = (Uint1*)&ecode2;
  412.                      CRC32(crc, hash_buf);
  413.                      ecode2 = (crc>>hash_shift) & hash_mask;
  414. #endif
  415.                      if (mb_lt->hashtable[ecode1] == 0 || 
  416.                          mb_lt->hashtable2[ecode2] == 0)
  417.                         mb_lt->num_unique_pos_added++;
  418.                      mb_lt->next_pos2[index] = mb_lt->hashtable2[ecode2];
  419.                      mb_lt->hashtable2[ecode2] = index;
  420.                   } else {
  421.                      if (mb_lt->hashtable[ecode1] == 0)
  422.                         mb_lt->num_unique_pos_added++;
  423.                   }
  424.                   mb_lt->next_pos[index] = mb_lt->hashtable[ecode1];
  425.                   mb_lt->hashtable[ecode1] = index;
  426.                }
  427.             }
  428.          }
  429.       }
  430.    }
  432.    mb_lt->pv_array_bts = pv_array_bts = PV_ARRAY_BTS + pv_shift; 
  434. #ifdef QUESTION_PV_ARRAY_USE
  435.    if (mb_lt->num_unique_pos_added < 
  436.        (PV_ARRAY_FACTOR*(mb_lt->hashsize>>pv_shift))) {
  437. #endif
  438.       pv_size = mb_lt->hashsize>>pv_array_bts;
  439.       /* Size measured in PV_ARRAY_TYPE's */
  440.       pv_array = calloc(PV_ARRAY_BYTES, pv_size);
  442.       for (index=0; index<mb_lt->hashsize; index++) {
  443.          if (mb_lt->hashtable[index] != 0 ||
  444.              (two_templates && mb_lt->hashtable2[index] != 0))
  445.             pv_array[index>>pv_array_bts] |= 
  446.                (((PV_ARRAY_TYPE) 1)<<(index&PV_ARRAY_MASK));
  447.       }
  448.       mb_lt->pv_array = pv_array;
  449. #ifdef QUESTION_PV_ARRAY_USE
  450.    }
  451. #endif
  453.    /* Now remove the hash entries that have too many positions */
  454. #if 0
  455.    if (lookup_options->max_positions>0) {
  456. #endif
  457.       for (pv_index = 0; pv_index < pv_size; ++pv_index) {
  458.          if (pv_array[pv_index]) {
  459.             for (ecode = pv_index<<pv_array_bts; 
  460.                  ecode < (pv_index+1)<<pv_array_bts; ++ecode) {
  461.                for (index=mb_lt->hashtable[ecode], pcount=0; 
  462.                     index>0; index=mb_lt->next_pos[index], pcount++);
  463.                if (two_templates) {
  464.                   for (index=mb_lt->hashtable2[ecode], pcount1=0; 
  465.                        index>0; index=mb_lt->next_pos2[index], pcount1++);
  466.                }
  467.                if ((pcount=MAX(pcount,pcount1))>lookup_options->max_positions) {
  468.                   mb_lt->hashtable[ecode] = 0;
  469.                   if (two_templates)
  470.                      mb_lt->hashtable2[ecode] = 0;
  471. #ifdef ERR_POST_EX_DEFINED
  472.                   ErrPostEx(SEV_WARNING, 0, 0, "%lx - %d", ecode, pcount);
  473. #endif
  474.                   masked_word_count++;
  475.                }
  476.                longest_chain = MAX(longest_chain, pcount);
  477.             }
  478.          }
  479.       }
  480.       if (masked_word_count) {
  481. #ifdef ERR_POST_EX_DEFINED
  482.  ErrPostEx(SEV_WARNING, 0, 0, "Masked %d %dmers", masked_word_count,
  483.    4*width);
  484. #endif
  485.       }
  486. #if 0
  487.    }
  488. #endif
  490.    mb_lt->longest_chain = longest_chain;
  492.    *mb_lt_ptr = mb_lt;
  494.    return 0;
  495. }
  497. Int4 MB_AG_ScanSubject(const LookupTableWrap* lookup_wrap,
  498.        const BLAST_SequenceBlk* subject, Int4 start_offset,
  499.        Uint4* q_offsets, Uint4* s_offsets, Int4 max_hits,  
  500.        Int4* end_offset)
  501. {
  502.    MBLookupTable* mb_lt = (MBLookupTable*) lookup_wrap->lut;
  503.    Uint1* s;
  504.    Uint1* abs_start;
  505.    Int4  index=0, s_off;
  506.    
  507.    Int4 q_off;
  508.    PV_ARRAY_TYPE *pv_array = mb_lt->pv_array;
  509.    Int4 total_hits = 0;
  510.    Int4 compressed_wordsize, compressed_scan_step, word_size;
  511.    Boolean full_byte_scan = mb_lt->full_byte_scan;
  512.    Uint1 pv_array_bts = mb_lt->pv_array_bts;
  513.    
  514.    /* Since the test for number of hits here is done after adding them, 
  515.       subtract the longest chain length from the allowed offset array size. */
  516.    max_hits -= mb_lt->longest_chain;
  518.    abs_start = subject->sequence;
  519.    s = abs_start + start_offset/COMPRESSION_RATIO;
  520.    compressed_scan_step = mb_lt->scan_step / COMPRESSION_RATIO;
  521.    compressed_wordsize = mb_lt->compressed_wordsize;
  522.    word_size = compressed_wordsize*COMPRESSION_RATIO;
  524.    index = 0;
  525.    
  526.    /* NB: s in this function always points to the start of the word!
  527.     */
  528.    if (full_byte_scan) {
  529.       Uint1* s_end = abs_start + subject->length/COMPRESSION_RATIO - 
  530.          compressed_wordsize;
  531.       for ( ; s <= s_end; s += compressed_scan_step) {
  532.          BlastNaLookupInitIndex(compressed_wordsize, s, &index);
  533.          
  534.          if (NA_PV_TEST(pv_array, index, pv_array_bts)) {
  535.             q_off = mb_lt->hashtable[index];
  536.             s_off = (s - abs_start)*COMPRESSION_RATIO;
  537.             if (q_off && (total_hits >= max_hits))
  538.                break;
  539.             while (q_off) {
  540.                q_offsets[total_hits] = q_off - 1;
  541.                s_offsets[total_hits++] = s_off;
  542.                q_off = mb_lt->next_pos[q_off];
  543.             }
  544.          }
  545.       }
  546.       *end_offset = (s - abs_start)*COMPRESSION_RATIO;
  547.    } else {
  548.       Int4 last_offset = subject->length - word_size;
  549.       Uint1 bit;
  550.       Int4 adjusted_index;
  552.       for (s_off = start_offset; s_off <= last_offset; 
  553.            s_off += mb_lt->scan_step) {
  554.          s = abs_start + (s_off / COMPRESSION_RATIO);
  555.          bit = 2*(s_off % COMPRESSION_RATIO);
  556.          /* Compute index for a word made of full bytes */
  557.          s = BlastNaLookupInitIndex(compressed_wordsize, s, &index);
  558.          /* Adjust the word index by the base within a byte */
  559.          adjusted_index = BlastNaLookupAdjustIndex(s, index, mb_lt->mask,
  560.                                                    bit);
  561.          if (NA_PV_TEST(pv_array, adjusted_index, pv_array_bts)) {
  562.             q_off = mb_lt->hashtable[adjusted_index];
  563.             if (q_off && (total_hits >= max_hits))
  564.                break;
  565.             while (q_off) {
  566.                q_offsets[total_hits] = q_off - 1;
  567.                s_offsets[total_hits++] = s_off; 
  568.                q_off = mb_lt->next_pos[q_off];
  569.             }
  570.          }
  571.       }
  572.       *end_offset = s_off;
  573.    }
  575.    return total_hits;
  576. }
  578. Int4 MB_ScanSubject(const LookupTableWrap* lookup,
  579.        const BLAST_SequenceBlk* subject, Int4 start_offset, 
  580.        Uint4* q_offsets, Uint4* s_offsets, Int4 max_hits,
  581.        Int4* end_offset) 
  582. {
  583.    Uint1* abs_start,* s_end;
  584.    Uint1* s;
  585.    Int4 hitsfound = 0;
  586.    Uint4 query_offset, subject_offset;
  587.    Int4 index;
  588.    MBLookupTable* mb_lt = (MBLookupTable*) lookup->lut;
  589.    Uint4* q_ptr = q_offsets,* s_ptr = s_offsets;
  590.    PV_ARRAY_TYPE *pv_array = mb_lt->pv_array;
  591.    Uint1 pv_array_bts = mb_lt->pv_array_bts;
  592.    Int4 compressed_wordsize = mb_lt->compressed_wordsize;
  594. #ifdef DEBUG_LOG
  595.    FILE *logfp0 = fopen("new0.log", "a");
  596. #endif   
  598.    /* Since the test for number of hits here is done after adding them, 
  599.       subtract the longest chain length from the allowed offset array size.
  600.    */
  601.    max_hits -= mb_lt->longest_chain;
  603.    abs_start = subject->sequence;
  604.    s = abs_start + start_offset/COMPRESSION_RATIO;
  605.    s_end = abs_start + (*end_offset)/COMPRESSION_RATIO;
  607.    s = BlastNaLookupInitIndex(mb_lt->compressed_wordsize, s, &index);
  609.    /* In the following code, s points to the byte right after the end of 
  610.       the word. */
  611.    while (s <= s_end) {
  612.       if (NA_PV_TEST(pv_array, index, pv_array_bts)) {
  613.          query_offset = mb_lt->hashtable[index];
  614.          subject_offset = 
  615.             ((s - abs_start) - compressed_wordsize)*COMPRESSION_RATIO;
  616.          if (query_offset && (hitsfound >= max_hits))
  617.             break;
  618.          while (query_offset) {
  619. #ifdef DEBUG_LOG
  620.             fprintf(logfp0, "%ldt%ldt%ldn", query_offset, 
  621.                     subject_offset, index);
  622. #endif
  623.             *(q_ptr++) = query_offset - 1;
  624.             *(s_ptr++) = subject_offset;
  625.             ++hitsfound;
  626.             query_offset = mb_lt->next_pos[query_offset];
  627.          }
  628.       }
  629.       /* Compute the next value of the lookup index 
  630.          (shifting sequence by a full byte) */
  631.       index = BlastNaLookupComputeIndex(FULL_BYTE_SHIFT, mb_lt->mask, 
  632.                                         s++, index);
  633.    }
  635.    *end_offset = 
  636.      ((s - abs_start) - compressed_wordsize)*COMPRESSION_RATIO;
  637. #ifdef DEBUG_LOG
  638.    fclose(logfp0);
  639. #endif
  641.    return hitsfound;
  642. }
  644. Int4 MB_DiscWordScanSubject(const LookupTableWrap* lookup, 
  645.        const BLAST_SequenceBlk* subject, Int4 start_offset,
  646.        Uint4* q_offsets, Uint4* s_offsets, Int4 max_hits, 
  647.        Int4* end_offset)
  648. {
  649.    Uint1* s;
  650.    Uint1* s_start,* abs_start,* s_end;
  651.    Int4 hitsfound = 0;
  652.    Uint4 query_offset, subject_offset;
  653.    Int4 word, index, index2=0;
  654.    MBLookupTable* mb_lt = (MBLookupTable*) lookup->lut;
  655.    Uint4* q_ptr = q_offsets,* s_ptr = s_offsets;
  656.    Boolean full_byte_scan = mb_lt->full_byte_scan;
  657.    Boolean two_templates = mb_lt->two_templates;
  658.    Uint1 template_type = mb_lt->template_type;
  659.    Uint1 second_template_type = mb_lt->second_template_type;
  660.    PV_ARRAY_TYPE *pv_array = mb_lt->pv_array;
  661.    Uint1 pv_array_bts = mb_lt->pv_array_bts;
  662.    Int4 compressed_wordsize = mb_lt->compressed_wordsize;
  664. #ifdef DEBUG_LOG
  665.    FILE *logfp0 = fopen("new0.log", "a");
  666. #endif   
  667.    
  668.    /* Since the test for number of hits here is done after adding them, 
  669.       subtract the longest chain length from the allowed offset array size. */
  670.    max_hits -= mb_lt->longest_chain;
  672.    abs_start = subject->sequence;
  673.    s_start = abs_start + start_offset/COMPRESSION_RATIO;
  674.    s_end = abs_start + (*end_offset)/COMPRESSION_RATIO;
  676.    s = BlastNaLookupInitIndex(mb_lt->compressed_wordsize, s_start, &word);
  678.    /* s now points to the byte right after the end of the current word */
  679.    if (full_byte_scan) {
  681.      while (s <= s_end) {
  682.        index = ComputeDiscontiguousIndex(s, word, template_type);
  684.        if (two_templates) {
  685.           index2 = ComputeDiscontiguousIndex(s, word, second_template_type);
  686.        }
  687.        if (NA_PV_TEST(pv_array, index, pv_array_bts)) {
  688.           query_offset = mb_lt->hashtable[index];
  689.           subject_offset = 
  690.              ((s - abs_start) - compressed_wordsize)*COMPRESSION_RATIO;
  691.           if (query_offset && (hitsfound >= max_hits))
  692.              break;
  693.           while (query_offset) {
  694. #ifdef DEBUG_LOG
  695.              fprintf(logfp0, "%ldt%ldt%ldn", query_offset, 
  696.                      subject_offset, index);
  697. #endif
  698.              *(q_ptr++) = query_offset - 1;
  699.              *(s_ptr++) = subject_offset;
  700.              ++hitsfound;
  701.              query_offset = mb_lt->next_pos[query_offset];
  702.           }
  703.        }
  704.        if (two_templates && NA_PV_TEST(pv_array, index2, pv_array_bts)) {
  705.           query_offset = mb_lt->hashtable2[index2];
  706.           subject_offset = 
  707.              ((s - abs_start) - compressed_wordsize)*COMPRESSION_RATIO;
  708.           if (query_offset && (hitsfound >= max_hits))
  709.              break;
  710.           while (query_offset) {
  711.              q_offsets[hitsfound] = query_offset - 1;
  712.              s_offsets[hitsfound++] = subject_offset;
  713.              query_offset = mb_lt->next_pos2[query_offset];
  714.           }
  715.        }
  716.        word = BlastNaLookupComputeIndex(FULL_BYTE_SHIFT, 
  717.                            mb_lt->mask, s++, word);
  718.      }
  719.    } else {
  720.       Int4 scan_shift = 2*mb_lt->scan_step;
  721.       Uint1 bit = 2*(start_offset % COMPRESSION_RATIO);
  722.       Int4 adjusted_word;
  724.       while (s <= s_end) {
  725.          /* Adjust the word index by the base within a byte */
  726.          adjusted_word = BlastNaLookupAdjustIndex(s, word, mb_lt->mask,
  727.                                          bit);
  728.    
  729.          index = ComputeDiscontiguousIndex_1b(s, adjusted_word, bit,
  730.                     template_type);
  731.        
  732.          if (two_templates)
  733.             index2 = ComputeDiscontiguousIndex_1b(s, adjusted_word, bit,
  734.                         second_template_type);
  736.          if (NA_PV_TEST(pv_array, index, pv_array_bts)) {
  737.             query_offset = mb_lt->hashtable[index];
  738.             subject_offset = 
  739.                ((s - abs_start) - compressed_wordsize)*COMPRESSION_RATIO
  740.                + bit/2;
  741.             if (query_offset && (hitsfound >= max_hits))
  742.                break;
  743.             while (query_offset) {
  744. #ifdef DEBUG_LOG
  745.                fprintf(logfp0, "%ldt%ldt%ldn", query_offset, 
  746.                        subject_offset, index);
  747. #endif
  748.                *(q_ptr++) = query_offset - 1;
  749.                *(s_ptr++) = subject_offset;
  750.                ++hitsfound;
  751.                query_offset = mb_lt->next_pos[query_offset];
  752.             }
  753.          }
  754.          if (two_templates && NA_PV_TEST(pv_array, index2, pv_array_bts)) {
  755.             query_offset = mb_lt->hashtable2[index2];
  756.             subject_offset = 
  757.                ((s - abs_start) - compressed_wordsize)*COMPRESSION_RATIO
  758.                + bit/2;
  759.             if (query_offset && (hitsfound >= max_hits))
  760.                break;
  761.             while (query_offset) {
  762.                q_offsets[hitsfound] = query_offset - 1;
  763.                s_offsets[hitsfound++] = subject_offset;
  764.                query_offset = mb_lt->next_pos2[query_offset];
  765.             }
  766.          }
  767.          bit += scan_shift;
  768.          if (bit >= FULL_BYTE_SHIFT) {
  769.             /* Advance to the next full byte */
  770.             bit -= FULL_BYTE_SHIFT;
  771.             word = BlastNaLookupComputeIndex(FULL_BYTE_SHIFT, 
  772.                                 mb_lt->mask, s++, word);
  773.          }
  774.       }
  775.    }
  776.    *end_offset = 
  777.      ((s - abs_start) - compressed_wordsize)*COMPRESSION_RATIO;
  778. #ifdef DEBUG_LOG
  779.    fclose(logfp0);
  780. #endif
  782.    return hitsfound;
  783. }