~ubuntu-branches/ubuntu/precise/ncbi-tools6/precise

« back to all changes in this revision

Viewing changes to algo/blast/core/blast_inline.h

  • Committer: Bazaar Package Importer
  • Author(s): Aaron M. Ucko
  • Date: 2005-03-27 12:00:15 UTC
  • mfrom: (2.1.2 hoary)
  • Revision ID: james.westby@ubuntu.com-20050327120015-embhesp32nj73p9r
Tags: 6.1.20041020-3
* Fix FTBFS under GCC 4.0 caused by inconsistent use of "static" on
  functions.  (Closes: #295110.)
* Add a watch file, now that we can.  (Upstream's layout needs version=3.)

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/* $Id: blast_inline.h,v 1.4 2004/09/01 13:19:43 madden Exp $
 
2
 * ===========================================================================
 
3
 *
 
4
 *                            PUBLIC DOMAIN NOTICE
 
5
 *               National Center for Biotechnology Information
 
6
 *
 
7
 *  This software/database is a "United States Government Work" under the
 
8
 *  terms of the United States Copyright Act.  It was written as part of
 
9
 *  the author's offical duties as a United States Government employee and
 
10
 *  thus cannot be copyrighted.  This software/database is freely available
 
11
 *  to the public for use. The National Library of Medicine and the U.S.
 
12
 *  Government have not placed any restriction on its use or reproduction.
 
13
 *
 
14
 *  Although all reasonable efforts have been taken to ensure the accuracy
 
15
 *  and reliability of the software and data, the NLM and the U.S.
 
16
 *  Government do not and cannot warrant the performance or results that
 
17
 *  may be obtained by using this software or data. The NLM and the U.S.
 
18
 *  Government disclaim all warranties, express or implied, including
 
19
 *  warranties of performance, merchantability or fitness for any particular
 
20
 *  purpose.
 
21
 *
 
22
 *  Please cite the author in any work or product based on this material.
 
23
 *
 
24
 * ===========================================================================
 
25
 *
 
26
 */
 
27
 
 
28
/** @file blast_inline.h
 
29
 * @todo FIXME needs file description
 
30
 */
 
31
 
 
32
#include <algo/blast/core/blast_util.h>
 
33
#include <algo/blast/core/blast_lookup.h>
 
34
#include <algo/blast/core/mb_lookup.h>
 
35
 
 
36
/** Given a word packed into an integer, compute a discontiguous word lookup 
 
37
 *  index.
 
38
 * @param subject Pointer to the next byte of the sequence after the end of 
 
39
 *        the word (needed when word template is longer than 16 bases) [in]
 
40
 * @param word A piece of the sequence packed into an integer [in]
 
41
 * @param template_type What type of discontiguous word template to use [in]
 
42
 * @return The lookup table index of the discontiguous word [out]
 
43
 */
 
44
static NCBI_INLINE Int4 ComputeDiscontiguousIndex(Uint1* subject, Int4 word,
 
45
                  Uint1 template_type)
 
46
{
 
47
   Int4 index;
 
48
   Int4 extra_code;   
 
49
 
 
50
   switch (template_type) {
 
51
   case TEMPL_11_16:
 
52
      index = GET_WORD_INDEX_11_16(word);
 
53
      break;
 
54
   case TEMPL_12_16:
 
55
      index = GET_WORD_INDEX_12_16(word);
 
56
      break;
 
57
   case TEMPL_11_16_OPT:
 
58
      index = GET_WORD_INDEX_11_16_OPT(word);
 
59
      break;
 
60
   case TEMPL_12_16_OPT:
 
61
      index = GET_WORD_INDEX_12_16_OPT(word);
 
62
      break;
 
63
   case TEMPL_11_18: 
 
64
     extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_18(subject);
 
65
     index = (GET_WORD_INDEX_11_18(word) | extra_code);
 
66
     break;
 
67
   case TEMPL_12_18: 
 
68
      extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_18(subject);
 
69
      index = (GET_WORD_INDEX_12_18(word) | extra_code);
 
70
      break;
 
71
   case TEMPL_11_18_OPT: 
 
72
      extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_18_OPT(subject);
 
73
      index = (GET_WORD_INDEX_11_18_OPT(word) | extra_code);
 
74
      break;
 
75
   case TEMPL_12_18_OPT:
 
76
      extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_18_OPT(subject);
 
77
      index = (GET_WORD_INDEX_12_18_OPT(word) | extra_code);
 
78
      break;
 
79
   case TEMPL_11_21: 
 
80
      extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_21(subject);
 
81
      index = (GET_WORD_INDEX_11_21(word) | extra_code);
 
82
      break;
 
83
   case TEMPL_12_21:
 
84
      extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_21(subject);
 
85
      index = (GET_WORD_INDEX_12_21(word) | extra_code);
 
86
      break;
 
87
   case TEMPL_11_21_OPT: 
 
88
      extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_21_OPT(subject);
 
89
      index = (GET_WORD_INDEX_11_21_OPT(word) | extra_code);
 
90
      break;
 
91
   case TEMPL_12_21_OPT:
 
92
      extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_21_OPT(subject);
 
93
      index = (GET_WORD_INDEX_12_21_OPT(word) | extra_code);
 
94
         break;
 
95
   default: 
 
96
      extra_code = 0; 
 
97
      index = 0;
 
98
      break;
 
99
   }
 
100
#ifdef USE_HASH_TABLE
 
101
   hash_buf = (Uint1*)&index;
 
102
   CRC32(crc, hash_buf);
 
103
   index = (crc>>hash_shift) & hash_mask;
 
104
#endif
 
105
 
 
106
   return index;
 
107
}
 
108
 
 
109
/** Compute the lookup table index for the first word template, given a word 
 
110
 * position, template type and previous value of the word, in case of 
 
111
 * one-base (2 bit) database scanning.
 
112
 * @param word_start Pointer to the start of a word in the sequence [in]
 
113
 * @param word The word packed into an integer value [in]
 
114
 * @param sequence_bit By how many bits the real word start is shifted within 
 
115
 *        a compressed sequence byte [in]
 
116
 * @param template_type What discontiguous word template to use for index 
 
117
 *        computation [in]
 
118
 * @return The lookup index for the discontiguous word.
 
119
*/
 
120
static NCBI_INLINE Int4 ComputeDiscontiguousIndex_1b(const Uint1* word_start, 
 
121
                      Int4 word, Uint1 sequence_bit, Uint1 template_type)
 
122
{
 
123
   Int4 index;
 
124
   Uint1* subject = (Uint1 *) word_start;
 
125
   Uint1 bit;
 
126
   Int4 extra_code, tmpval;   
 
127
 
 
128
   /* Prepare auxiliary variables for extra code calculation */
 
129
   tmpval = 0;
 
130
   extra_code = 0;
 
131
   /* The bits in an integer byte are counted in a reverse order than in a
 
132
      sequence byte */
 
133
   bit = 6 - sequence_bit;
 
134
 
 
135
   switch (template_type) {
 
136
   case TEMPL_11_16:
 
137
      index = GET_WORD_INDEX_11_16(word);
 
138
      break;
 
139
   case TEMPL_12_16:
 
140
      index = GET_WORD_INDEX_12_16(word);
 
141
      break;
 
142
   case TEMPL_11_16_OPT:
 
143
      index = GET_WORD_INDEX_11_16_OPT(word);
 
144
      break;
 
145
   case TEMPL_12_16_OPT:
 
146
      index = GET_WORD_INDEX_12_16_OPT(word);
 
147
      break;
 
148
   case TEMPL_11_18: 
 
149
      GET_EXTRA_CODE_PACKED_18(subject, bit, tmpval, extra_code);
 
150
      index = (GET_WORD_INDEX_11_18(word) | extra_code);
 
151
      break;
 
152
   case TEMPL_12_18: 
 
153
      GET_EXTRA_CODE_PACKED_18(subject, bit, tmpval, extra_code);
 
154
      index = (GET_WORD_INDEX_12_18(word) | extra_code);
 
155
      break;
 
156
   case TEMPL_11_18_OPT: 
 
157
      GET_EXTRA_CODE_PACKED_18_OPT(subject, bit, tmpval, extra_code);
 
158
      index = (GET_WORD_INDEX_11_18_OPT(word) | extra_code);
 
159
      break;
 
160
   case TEMPL_12_18_OPT:
 
161
      GET_EXTRA_CODE_PACKED_18_OPT(subject, bit, tmpval, extra_code);
 
162
      index = (GET_WORD_INDEX_12_18_OPT(word) | extra_code);
 
163
      break;
 
164
   case TEMPL_11_21: 
 
165
      GET_EXTRA_CODE_PACKED_21(subject, bit, tmpval, extra_code);
 
166
      index = (GET_WORD_INDEX_11_21(word) | extra_code);
 
167
      break;
 
168
   case TEMPL_12_21:
 
169
      GET_EXTRA_CODE_PACKED_21(subject, bit, tmpval, extra_code);
 
170
      index = (GET_WORD_INDEX_12_21(word) | extra_code);
 
171
      break;
 
172
   case TEMPL_11_21_OPT: 
 
173
      GET_EXTRA_CODE_PACKED_21_OPT(subject, bit, tmpval, extra_code);
 
174
      index = (GET_WORD_INDEX_11_21_OPT(word) | extra_code);
 
175
      break;
 
176
   case TEMPL_12_21_OPT:
 
177
      GET_EXTRA_CODE_PACKED_21_OPT(subject, bit, tmpval, extra_code);
 
178
      index = (GET_WORD_INDEX_12_21_OPT(word) | extra_code);
 
179
         break;
 
180
   default: 
 
181
      extra_code = 0; 
 
182
      index = 0;
 
183
      break;
 
184
   }
 
185
#ifdef USE_HASH_TABLE
 
186
   hash_buf = (Uint1*)&index;
 
187
   CRC32(crc, hash_buf);
 
188
   index = (crc>>hash_shift) & hash_mask;
 
189
#endif
 
190
  
 
191
   return index;
 
192
}
 
193
 
 
194
static NCBI_INLINE void  _ComputeIndex(Int4 wordsize,
 
195
                                  Int4 charsize,
 
196
                                  Int4 mask,
 
197
                                  const Uint1* word,
 
198
                                  Int4* index);
 
199
 
 
200
static NCBI_INLINE void  _ComputeIndexIncremental(Int4 wordsize,
 
201
                                             Int4 charsize,
 
202
                                             Int4 mask,
 
203
                                             const Uint1* word,
 
204
                                             Int4* index);
 
205
 
 
206
/** Given a word, compute its index value from scratch.
 
207
 *
 
208
 * @param wordsize length of the word, in residues [in]
 
209
 * @param charsize length of one residue, in bits [in]
 
210
 * @param mask value used to mask the index so that only the bottom wordsize * charsize bits remain [in]
 
211
 * @param word pointer to the beginning of the word [in]
 
212
 * @param index the computed index value [out] 
 
213
 */
 
214
 
 
215
static NCBI_INLINE void  _ComputeIndex(Int4 wordsize,
 
216
                                  Int4 charsize,
 
217
                                  Int4 mask,
 
218
                                  const Uint1* word,
 
219
                                  Int4* index)
 
220
{
 
221
  Int4 i;
 
222
 
 
223
  *index = 0;
 
224
 
 
225
  for(i=0;i<wordsize;i++)
 
226
{
 
227
    *index = ((*index << charsize) | word[i]) & mask;
 
228
}
 
229
 
 
230
  return;
 
231
}
 
232
 
 
233
/** Given a word, compute its index value, reusing a previously 
 
234
 *  computed index value.
 
235
 *
 
236
 * @param wordsize length of the word - 1, in residues [in]
 
237
 * @param charsize length of one residue, in bits [in]
 
238
 * @param mask value used to mask the index so that only the bottom wordsize * charsize bits remain [in]
 
239
 * @param word pointer to the beginning of the word [in]
 
240
 * @param index the computed index value [in/out]
 
241
 */
 
242
 
 
243
static NCBI_INLINE void  _ComputeIndexIncremental(Int4 wordsize,
 
244
                                             Int4 charsize,
 
245
                                             Int4 mask,
 
246
                                             const Uint1* word,
 
247
                                             Int4* index)
 
248
{
 
249
  *index = ((*index << charsize) | word[wordsize - 1]) & mask;
 
250
  return;
 
251
}
 
252
 
 
253
 
 
254
/* Given a starting position of a word in a compressed nucleotide sequence, 
 
255
 * compute this word's lookup table index
 
256
 */
 
257
static NCBI_INLINE Uint1* BlastNaLookupInitIndex(Int4 length,
 
258
                          const Uint1* word, Int4* index)
 
259
{
 
260
   Int4 i;
 
261
   
 
262
   *index = 0;
 
263
   for (i = 0; i < length; ++i)
 
264
      *index = ((*index)<<FULL_BYTE_SHIFT) | word[i];
 
265
   return (Uint1 *) (word + length);
 
266
}
 
267
 
 
268
/* Recompute the word index given its previous value and the new location 
 
269
 *  of the last byte of the word
 
270
 */
 
271
static NCBI_INLINE Int4 BlastNaLookupComputeIndex(Int4 scan_shift, Int4 mask, 
 
272
                      const Uint1* word, Int4 index)
 
273
{
 
274
   return (((index)<<scan_shift) & mask) | *(word);  
 
275
   
 
276
}
 
277
 
 
278
/* Given a word computed from full bytes of a compressed sequence, 
 
279
 *  shift it by 0-3 bases 
 
280
 */
 
281
static NCBI_INLINE Int4 BlastNaLookupAdjustIndex(Uint1* s, Int4 index, 
 
282
                      Int4 mask, Uint1 bit)
 
283
{
 
284
   return (((index)<<bit) & mask) | ((*s)>>(FULL_BYTE_SHIFT-bit));
 
285
 
 
286
}
 
287
 
 
288
 
 
289
#define BLAST2NA_MASK 0xfc
 
290
 
 
291
/** Compute the lookup index for a word in an uncompressed sequence, without 
 
292
 *  using any previous index information.
 
293
 * @param lookup Pointer to the traditional BLASTn lookup table structure [in]
 
294
 * @param word Pointer to the start of the word [in]
 
295
 * @param index The lookup index [out]
 
296
 */
 
297
static NCBI_INLINE Int2
 
298
Na_LookupComputeIndex(BlastLookupTable* lookup, Uint1* word, Int4* index)
 
299
{
 
300
   Int4 i;
 
301
   Int4 wordsize = lookup->reduced_wordsize*COMPRESSION_RATIO; /* i.e. 8 or 4 */
 
302
 
 
303
   *index = 0;
 
304
   for (i = 0; i < wordsize; ++i) {
 
305
      if ((word[i] & BLAST2NA_MASK) != 0) {
 
306
         *index = 0;
 
307
         return -1;
 
308
      } else {
 
309
         *index = (((*index)<<lookup->charsize) & lookup->mask) | word[i];
 
310
      }
 
311
   }
 
312
   return 0;
 
313
}
 
314
 
 
315
/** Pack 4 sequence bytes into a one byte integer, assuming sequence contains
 
316
 *  no ambiguities.
 
317
 */
 
318
#define PACK_WORD(q) ((q[0]<<6) + (q[1]<<4) + (q[2]<<2) + q[3]) 
 
319
 
 
320
/** Compare a given number of bytes of an compressed subject sequence with
 
321
 *  the non-compressed query sequence.
 
322
 * @param q Pointer to the first byte to be compared in the query sequence [in]
 
323
 * @param s Pointer to the first byte to be compared in the subject 
 
324
 *        sequence [in]
 
325
 * @param extra_bytes Number of compressed bytes to compare [in]
 
326
 * @return TRUE if sequences are identical, FALSE if mismatch is found.
 
327
 */
 
328
static NCBI_INLINE Boolean BlastNaCompareExtraBytes(Uint1* q, Uint1* s, 
 
329
                                               Int4 extra_bytes)
 
330
{
 
331
   Int4 index;
 
332
   
 
333
   for (index = 0; index < extra_bytes; ++index) {
 
334
      if (*s++ != PACK_WORD(q))
 
335
         return FALSE;
 
336
      q += COMPRESSION_RATIO;
 
337
   }
 
338
   return TRUE;
 
339
}
 
340
 
 
341
/** Perform mini extension (up to max_left <= 4 bases) to the left; 
 
342
 * @param q Pointer to the query base right after the ones to be extended [in]
 
343
 * @param s Pointer to a byte in the compressed subject sequence that is to be 
 
344
 *        tested for extension [in]
 
345
 * @param max_left Maximal number of bits to compare [in]
 
346
 * @return Number of matched bases
 
347
*/
 
348
static NCBI_INLINE Uint1 
 
349
BlastNaMiniExtendLeft(Uint1* q, const Uint1* s, Uint1 max_left)
 
350
{
 
351
   Uint1 left = 0;
 
352
 
 
353
   for (left = 0; left < max_left; ++left) {
 
354
      if (NCBI2NA_UNPACK_BASE(*s, left) != *--q) {
 
355
         break;
 
356
      }
 
357
   }
 
358
   return left;
 
359
}
 
360
 
 
361
/** Perform mini extension (up to max_right <= 4 bases) to the right; 
 
362
 * @param q Pointer to the start of the extension in the query [in]
 
363
 * @param s Pointer to a byte in the compressed subject sequence that is to be 
 
364
 *        tested for extension [in]
 
365
 * @param max_right Maximal number of bits to compare [in]
 
366
 * @return Number of matched bases
 
367
*/
 
368
static NCBI_INLINE Uint1 
 
369
BlastNaMiniExtendRight(Uint1* q, const Uint1* s, Uint1 max_right)
 
370
{
 
371
   Uint1 right;
 
372
   Uint1 index = 3;
 
373
   
 
374
   for (right = 0; right < max_right; ++right, --index) {
 
375
      if (NCBI2NA_UNPACK_BASE(*s, index) != *q++) {
 
376
         break;
 
377
      }
 
378
   }
 
379
   return right;
 
380
}
 
381
 
 
382
/** Returns the index in the MaskLoc given a context number for the query.
 
383
 * If the query is nucleotide
 
384
 *
 
385
 * @param is_na the query is nucleotide
 
386
 * @param context offset in the QueryInfo array
 
387
 * @return index in the maskloc
 
388
 */
 
389
static NCBI_INLINE Int2 BlastGetMaskLocIndexFromContext(Boolean is_na, Int2 context)
 
390
{
 
391
     return (is_na ? context / 2 : context);
 
392
}
 
393
 
 
394
/** Determines whether this is a nucleotide query and whether this a minus strand or not
 
395
 *
 
396
 * @param is_na the query is nucleotide
 
397
 * @param context offset in the QueryInfo array
 
398
 * @return TRUE if this is minus strand
 
399
 */
 
400
static NCBI_INLINE Boolean BlastIsReverseStrand(Boolean is_na, Int2 context)
 
401
{
 
402
     return (is_na && ((context & 1) != 0));
 
403
 
 
404
}