?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
modified: makefile
[GalaxyCodeBases.git] / BGI / soap_src / soap_builder / DNACount.c
blobd2178e25b9f09a438fe7d0de4c95449e8ea1c362
1 /*
2
3 DNACount.c DNA Count
4
5 This module contains DNA occurrence counting functions. The DNA must be
6 in word-packed format.
7
8 Copyright (C) 2004, Wong Chi Kwong.
9
10 This program is free software; you can redistribute it and/or
11 modify it under the terms of the GNU General Public License
12 as published by the Free Software Foundation; either version 2
13 of the License, or (at your option) any later version.
14
15 This program is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU General Public License for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with this program; if not, write to the Free Software
22 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
23
24 */
25
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include "DNACount.h"
29 #include "MiscUtilities.h"
30
31
32 void GenerateDNAOccCountTable(unsigned int *dnaDecodeTable) {
33
34 unsigned int i, j, c, t;
35
36 for (i=0; i<DNA_OCC_CNT_TABLE_SIZE_IN_WORD; i++) {
37 dnaDecodeTable[i] = 0;
38 c = i;
39 for (j=0; j<8; j++) {
40 t = c & 0x00000003;
41 dnaDecodeTable[i] += 1 << (t * 8);
42 c >>= 2;
43 }
44 }
45
46 }
47
48 unsigned int ForwardDNAOccCount(const unsigned int* dna, const unsigned int index, const unsigned int character, const unsigned int* dnaDecodeTable) {
49
50 static const unsigned int truncateRightMask[16] = { 0x00000000, 0xC0000000, 0xF0000000, 0xFC000000,
51 0xFF000000, 0xFFC00000, 0xFFF00000, 0xFFFC0000,
52 0xFFFF0000, 0xFFFFC000, 0xFFFFF000, 0xFFFFFC00,
53 0xFFFFFF00, 0xFFFFFFC0, 0xFFFFFFF0, 0xFFFFFFFC };
54
55 unsigned int wordToCount, charToCount;
56 unsigned int i, c;
57 unsigned int sum = 0;
58
59 #ifdef DEBUG
60 if (index >= 256) {
61 fprintf(stderr, "ForwardDNAOccCount() : index >= 256!\n");
62 exit(1);
63 }
64 #endif
65
66 wordToCount = index / 16;
67 charToCount = index - wordToCount * 16;
68
69 for (i=0; i<wordToCount; i++) {
70 sum += dnaDecodeTable[dna[i] >> 16];
71 sum += dnaDecodeTable[dna[i] & 0x0000FFFF];
72 }
73
74 if (charToCount > 0) {
75 c = dna[i] & truncateRightMask[charToCount]; // increase count of 'a' by 16 - c;
76 sum += dnaDecodeTable[c >> 16];
77 sum += dnaDecodeTable[c & 0xFFFF];
78 sum += charToCount - 16; // decrease count of 'a' by 16 - positionToProcess
79 }
80
81 return (sum >> (character * 8)) & 0x000000FF;
82
83 }
84
85 unsigned int BackwardDNAOccCount(const unsigned int* dna, const unsigned int index, const unsigned int character, const unsigned int* dnaDecodeTable) {
86
87 static const unsigned int truncateLeftMask[16] = { 0x00000000, 0x00000003, 0x0000000F, 0x0000003F,
88 0x000000FF, 0x000003FF, 0x00000FFF, 0x00003FFF,
89 0x0000FFFF, 0x0003FFFF, 0x000FFFFF, 0x003FFFFF,
90 0x00FFFFFF, 0x03FFFFFF, 0x0FFFFFFF, 0x3FFFFFFF };
91
92 unsigned int wordToCount, charToCount;
93 unsigned int i, c;
94 unsigned int sum = 0;
95
96 #ifdef DEBUG
97 if (index >= 256) {
98 fprintf(stderr, "ForwardDNAOccCount() : index >= 256!\n");
99 exit(1);
100 }
101 #endif
102
103 wordToCount = index / 16;
104 charToCount = index - wordToCount * 16;
105
106 dna -= wordToCount + 1;
107
108 if (charToCount > 0) {
109 c = *dna & truncateLeftMask[charToCount]; // increase count of 'a' by 16 - c;
110 sum += dnaDecodeTable[c >> 16];
111 sum += dnaDecodeTable[c & 0xFFFF];
112 sum += charToCount - 16; // decrease count of 'a' by 16 - positionToProcess
113 }
114
115 for (i=0; i<wordToCount; i++) {
116 dna++;
117 sum += dnaDecodeTable[*dna >> 16];
118 sum += dnaDecodeTable[*dna & 0x0000FFFF];
119 }
120
121 return (sum >> (character * 8)) & 0x000000FF;
122
123 }
124
125 void ForwardDNAAllOccCount(const unsigned int* dna, const unsigned int index, unsigned int* __restrict occCount, const unsigned int* dnaDecodeTable) {
126
127 static const unsigned int truncateRightMask[16] = { 0x00000000, 0xC0000000, 0xF0000000, 0xFC000000,
128 0xFF000000, 0xFFC00000, 0xFFF00000, 0xFFFC0000,
129 0xFFFF0000, 0xFFFFC000, 0xFFFFF000, 0xFFFFFC00,
130 0xFFFFFF00, 0xFFFFFFC0, 0xFFFFFFF0, 0xFFFFFFFC };
131
132 unsigned int wordToCount, charToCount;
133 unsigned int i, c;
134 unsigned int sum = 0;
135
136 #ifdef DEBUG
137 if (index >= 256) {
138 fprintf(stderr, "ForwardDNAOccCount() : index >= 256!\n");
139 exit(1);
140 }
141 #endif
142
143 wordToCount = index / 16;
144 charToCount = index - wordToCount * 16;
145
146 for (i=0; i<wordToCount; i++) {
147 sum += dnaDecodeTable[dna[i] >> 16];
148 sum += dnaDecodeTable[dna[i] & 0x0000FFFF];
149 }
150
151 if (charToCount > 0) {
152 c = dna[i] & truncateRightMask[charToCount]; // increase count of 'a' by 16 - c;
153 sum += dnaDecodeTable[c >> 16];
154 sum += dnaDecodeTable[c & 0xFFFF];
155 sum += charToCount - 16; // decrease count of 'a' by 16 - positionToProcess
156 }
157
158 occCount[0] = sum & 0x000000FF; sum >>= 8;
159 occCount[1] = sum & 0x000000FF; sum >>= 8;
160 occCount[2] = sum & 0x000000FF; sum >>= 8;
161 occCount[3] = sum;
162
163 }
164
165 void BackwardDNAAllOccCount(const unsigned int* dna, const unsigned int index, unsigned int* __restrict occCount, const unsigned int* dnaDecodeTable) {
166
167 static const unsigned int truncateLeftMask[16] = { 0x00000000, 0x00000003, 0x0000000F, 0x0000003F,
168 0x000000FF, 0x000003FF, 0x00000FFF, 0x00003FFF,
169 0x0000FFFF, 0x0003FFFF, 0x000FFFFF, 0x003FFFFF,
170 0x00FFFFFF, 0x03FFFFFF, 0x0FFFFFFF, 0x3FFFFFFF };
171
172 unsigned int wordToCount, charToCount;
173 unsigned int i, c;
174 unsigned int sum = 0;
175
176 #ifdef DEBUG
177 if (index >= 256) {
178 fprintf(stderr, "ForwardDNAOccCount() : index >= 256!\n");
179 exit(1);
180 }
181 #endif
182
183 wordToCount = index / 16;
184 charToCount = index - wordToCount * 16;
185
186 dna -= wordToCount + 1;
187
188 if (charToCount > 0) {
189 c = *dna & truncateLeftMask[charToCount]; // increase count of 'a' by 16 - c;
190 sum += dnaDecodeTable[c >> 16];
191 sum += dnaDecodeTable[c & 0xFFFF];
192 sum += charToCount - 16; // decrease count of 'a' by 16 - positionToProcess
193 }
194
195 for (i=0; i<wordToCount; i++) {
196 dna++;
197 sum += dnaDecodeTable[*dna >> 16];
198 sum += dnaDecodeTable[*dna & 0x0000FFFF];
199 }
200
201 occCount[0] = sum & 0x000000FF; sum >>= 8;
202 occCount[1] = sum & 0x000000FF; sum >>= 8;
203 occCount[2] = sum & 0x000000FF; sum >>= 8;
204 occCount[3] = sum;
205
206 }
207
208 unsigned int Forward1OccCount(const unsigned int* bitVector, const unsigned int index, const unsigned int* dnaDecodeTable) {
209
210 static const unsigned int truncateRightMask[32] = { 0x00000000, 0x80000000, 0xC0000000, 0xE0000000, 0xF0000000, 0xF8000000, 0xFC000000, 0xFE000000,
211 0xFF000000, 0xFF800000, 0xFFC00000, 0xFFE00000, 0xFFF00000, 0xFFF80000, 0xFFFC0000, 0xFFFE0000,
212 0xFFFF0000, 0xFFFF8000, 0xFFFFC000, 0xFFFFE000, 0xFFFFF000, 0xFFFFF800, 0xFFFFFC00, 0xFFFFFE00,
213 0xFFFFFF00, 0xFFFFFF80, 0xFFFFFFC0, 0xFFFFFFE0, 0xFFFFFFF0, 0xFFFFFFF8, 0xFFFFFFFC, 0xFFFFFFFE};
214
215 unsigned int wordToCount, bitToCount;
216 unsigned int i, c;
217 unsigned int sum = 0;
218 unsigned int numberOf1;
219
220 #ifdef DEBUG
221 if (index >= 256) {
222 fprintf(stderr, "Forward1OccCount() : index >= 256!\n");
223 exit(1);
224 }
225 #endif
226
227 wordToCount = index / 32;
228 bitToCount = index - wordToCount * 32;
229
230 for (i=0; i<wordToCount; i++) {
231 sum += dnaDecodeTable[bitVector[i] >> 16];
232 sum += dnaDecodeTable[bitVector[i] & 0x0000FFFF];
233 }
234
235 if (bitToCount > 0) {
236 c = bitVector[i] & truncateRightMask[bitToCount];
237 sum += dnaDecodeTable[c >> 16];
238 sum += dnaDecodeTable[c & 0x0000FFFF];
239 }
240
241 sum >>= 8;
242 numberOf1 = sum & 0x000000FF;
243 sum >>= 8;
244 numberOf1 = sum & 0x000000FF;
245 sum >>= 8;
246 numberOf1 = sum * 2;
247
248 return numberOf1;
249
250 }
251
252 unsigned int Backward1OccCount(const unsigned int* bitVector, const unsigned int index, const unsigned int* dnaDecodeTable) {
253
254 static const unsigned int truncateLeftMask[32] = { 0x00000000, 0x00000001, 0x00000003, 0x00000007, 0x0000000F, 0x0000001F, 0x0000003F, 0x0000007F,
255 0x000000FF, 0x000001FF, 0x000003FF, 0x000007FF, 0x00000FFF, 0x00001FFF, 0x00003FFF, 0x00007FFF,
256 0x0000FFFF, 0x0001FFFF, 0x0003FFFF, 0x0007FFFF, 0x000FFFFF, 0x001FFFFF, 0x003FFFFF, 0x007FFFFF,
257 0x00FFFFFF, 0x01FFFFFF, 0x03FFFFFF, 0x07FFFFFF, 0x0FFFFFFF, 0x1FFFFFFF, 0x3FFFFFFF, 0x7FFFFFFF};
258
259 unsigned int wordToCount, bitToCount;
260 unsigned int i, c;
261 unsigned int sum = 0;
262 unsigned int numberOf1;
263
264 #ifdef DEBUG
265 if (index >= 256) {
266 fprintf(stderr, "ForwardDNAOccCount() : index >= 256!\n");
267 exit(1);
268 }
269 #endif
270
271 wordToCount = index / 32;
272 bitToCount = index - wordToCount * 32;
273
274 bitVector -= wordToCount + 1;
275
276 if (bitToCount > 0) {
277 c = *bitVector & truncateLeftMask[bitToCount];
278 sum += dnaDecodeTable[c >> 16];
279 sum += dnaDecodeTable[c & 0xFFFF];
280 }
281
282 for (i=0; i<wordToCount; i++) {
283 bitVector++;
284 sum += dnaDecodeTable[*bitVector >> 16];
285 sum += dnaDecodeTable[*bitVector & 0x0000FFFF];
286 }
287
288 sum >>= 8;
289 numberOf1 = sum & 0x000000FF;
290 sum >>= 8;
291 numberOf1 = sum & 0x000000FF;
292 sum >>= 8;
293 numberOf1 = sum * 2;
294
295 return numberOf1;
296
297 }
298
299 unsigned int ForwardDNAOccCountNoLimit(const unsigned int* dna, const unsigned int index, const unsigned int character, const unsigned int* dnaDecodeTable) {
300
301 static const unsigned int truncateRightMask[16] = { 0x00000000, 0xC0000000, 0xF0000000, 0xFC000000,
302 0xFF000000, 0xFFC00000, 0xFFF00000, 0xFFFC0000,
303 0xFFFF0000, 0xFFFFC000, 0xFFFFF000, 0xFFFFFC00,
304 0xFFFFFF00, 0xFFFFFFC0, 0xFFFFFFF0, 0xFFFFFFFC };
305
306 unsigned int iteration, wordToCount, charToCount;
307 unsigned int i, j, c;
308 unsigned int sum;
309 unsigned int occCount = 0;
310
311 iteration = index / 256;
312 wordToCount = (index - iteration * 256) / 16;
313 charToCount = index - iteration * 256 - wordToCount * 16;
314
315 for (i=0; i<iteration; i++) {
316
317 sum = 0;
318 for (j=0; j<16; j++) {
319 sum += dnaDecodeTable[*dna >> 16];
320 sum += dnaDecodeTable[*dna & 0x0000FFFF];
321 dna++;
322 }
323 if (!DNA_OCC_SUM_EXCEPTION(sum)) {
324 occCount += (sum >> (character * 8)) & 0x000000FF;
325 } else {
326 // only some or all of the 3 bits are on
327 // in reality, only one of the four cases are possible
328 if (sum == 0x00000100) {
329 if (character == 0) {
330 occCount += 256;
331 }
332 } else if (sum == 0x00010000) {
333 if (character == 1) {
334 occCount += 256;
335 }
336 } else if (sum == 0x01000000) {
337 if (character == 2) {
338 occCount += 256;
339 }
340 } else if (sum == 0x00000000) {
341 if (character == 3) {
342 occCount += 256;
343 }
344 } else {
345 fprintf(stderr, "ForwardDNAOccCountNoLimit(): DNA occ sum exception!\n");
346 exit(1);
347 }
348 }
349
350 }
351
352 sum = 0;
353 for (j=0; j<wordToCount; j++) {
354 sum += dnaDecodeTable[*dna >> 16];
355 sum += dnaDecodeTable[*dna & 0x0000FFFF];
356 dna++;
357 }
358
359 if (charToCount > 0) {
360 c = *dna & truncateRightMask[charToCount]; // increase count of 'a' by 16 - c;
361 sum += dnaDecodeTable[c >> 16];
362 sum += dnaDecodeTable[c & 0xFFFF];
363 sum += charToCount - 16; // decrease count of 'a' by 16 - positionToProcess
364 }
365
366 occCount += (sum >> (character * 8)) & 0x000000FF;
367
368 return occCount;
369
370 }
371
372 unsigned int BackwardDNAOccCountNoLimit(const unsigned int* dna, const unsigned int index, const unsigned int character, const unsigned int* dnaDecodeTable) {
373
374 static const unsigned int truncateLeftMask[16] = { 0x00000000, 0x00000003, 0x0000000F, 0x0000003F,
375 0x000000FF, 0x000003FF, 0x00000FFF, 0x00003FFF,
376 0x0000FFFF, 0x0003FFFF, 0x000FFFFF, 0x003FFFFF,
377 0x00FFFFFF, 0x03FFFFFF, 0x0FFFFFFF, 0x3FFFFFFF };
378
379 unsigned int iteration, wordToCount, charToCount;
380 unsigned int i, j, c;
381 unsigned int sum = 0;
382 unsigned int occCount;
383
384 dna -= index / 16 + 1;
385
386 iteration = index / 256;
387 wordToCount = (index - iteration * 256) / 16;
388 charToCount = index - iteration * 256 - wordToCount * 16;
389
390 if (charToCount > 0) {
391 c = *dna & truncateLeftMask[charToCount]; // increase count of 'a' by 16 - c;
392 sum += dnaDecodeTable[c >> 16];
393 sum += dnaDecodeTable[c & 0xFFFF];
394 sum += charToCount - 16; // decrease count of 'a' by 16 - positionToProcess
395 }
396
397 for (j=0; j<wordToCount; j++) {
398 dna++;
399 sum += dnaDecodeTable[*dna >> 16];
400 sum += dnaDecodeTable[*dna & 0x0000FFFF];
401 }
402
403 occCount = (sum >> (character * 8)) & 0x000000FF;
404
405
406 for (i=0; i<iteration; i++) {
407
408 sum = 0;
409 for (j=0; j<16; j++) {
410 dna++;
411 sum += dnaDecodeTable[*dna >> 16];
412 sum += dnaDecodeTable[*dna & 0x0000FFFF];
413 }
414 if (!DNA_OCC_SUM_EXCEPTION(sum)) {
415 occCount += (sum >> (character * 8)) & 0x000000FF;
416 } else {
417 // only some or all of the 3 bits are on
418 // in reality, only one of the four cases are possible
419 if (sum == 0x00000100) {
420 if (character == 0) {
421 occCount += 256;
422 }
423 } else if (sum == 0x00010000) {
424 if (character == 1) {
425 occCount += 256;
426 }
427 } else if (sum == 0x01000000) {
428 if (character == 2) {
429 occCount += 256;
430 }
431 } else if (sum == 0x00000000) {
432 if (character == 3) {
433 occCount += 256;
434 }
435 } else {
436 fprintf(stderr, "BackwardDNAOccCountNoLimit(): DNA occ sum exception!\n");
437 exit(1);
438 }
439 }
440
441 }
442
443 return occCount;
444
445 }
446
447 void ForwardDNAAllOccCountNoLimit(const unsigned int* dna, const unsigned int index, unsigned int* __restrict occCount, const unsigned int* dnaDecodeTable) {
448
449 static const unsigned int truncateRightMask[16] = { 0x00000000, 0xC0000000, 0xF0000000, 0xFC000000,
450 0xFF000000, 0xFFC00000, 0xFFF00000, 0xFFFC0000,
451 0xFFFF0000, 0xFFFFC000, 0xFFFFF000, 0xFFFFFC00,
452 0xFFFFFF00, 0xFFFFFFC0, 0xFFFFFFF0, 0xFFFFFFFC };
453
454 unsigned int iteration, wordToCount, charToCount;
455 unsigned int i, j, c;
456 unsigned int sum;
457
458 occCount[0] = 0;
459 occCount[1] = 0;
460 occCount[2] = 0;
461 occCount[3] = 0;
462
463 iteration = index / 256;
464 wordToCount = (index - iteration * 256) / 16;
465 charToCount = index - iteration * 256 - wordToCount * 16;
466
467 for (i=0; i<iteration; i++) {
468
469 sum = 0;
470 for (j=0; j<16; j++) {
471 sum += dnaDecodeTable[*dna >> 16];
472 sum += dnaDecodeTable[*dna & 0x0000FFFF];
473 dna++;
474 }
475 if (!DNA_OCC_SUM_EXCEPTION(sum)) {
476 occCount[0] += sum & 0x000000FF; sum >>= 8;
477 occCount[1] += sum & 0x000000FF; sum >>= 8;
478 occCount[2] += sum & 0x000000FF; sum >>= 8;
479 occCount[3] += sum;
480 } else {
481 // only some or all of the 3 bits are on
482 // in reality, only one of the four cases are possible
483 if (sum == 0x00000100) {
484 occCount[0] += 256;
485 } else if (sum == 0x00010000) {
486 occCount[1] += 256;
487 } else if (sum == 0x01000000) {
488 occCount[2] += 256;
489 } else if (sum == 0x00000000) {
490 occCount[3] += 256;
491 } else {
492 fprintf(stderr, "ForwardDNAAllOccCountNoLimit(): DNA occ sum exception!\n");
493 exit(1);
494 }
495 }
496
497 }
498
499 sum = 0;
500 for (j=0; j<wordToCount; j++) {
501 sum += dnaDecodeTable[*dna >> 16];
502 sum += dnaDecodeTable[*dna & 0x0000FFFF];
503 dna++;
504 }
505
506 if (charToCount > 0) {
507 c = *dna & truncateRightMask[charToCount]; // increase count of 'a' by 16 - c;
508 sum += dnaDecodeTable[c >> 16];
509 sum += dnaDecodeTable[c & 0xFFFF];
510 sum += charToCount - 16; // decrease count of 'a' by 16 - positionToProcess
511 }
512
513 occCount[0] += sum & 0x000000FF; sum >>= 8;
514 occCount[1] += sum & 0x000000FF; sum >>= 8;
515 occCount[2] += sum & 0x000000FF; sum >>= 8;
516 occCount[3] += sum;
517
518 }
519
520 void BackwardDNAAllOccCountNoLimit(const unsigned int* dna, const unsigned int index, unsigned int* __restrict occCount, const unsigned int* dnaDecodeTable) {
521
522 static const unsigned int truncateLeftMask[16] = { 0x00000000, 0x00000003, 0x0000000F, 0x0000003F,
523 0x000000FF, 0x000003FF, 0x00000FFF, 0x00003FFF,
524 0x0000FFFF, 0x0003FFFF, 0x000FFFFF, 0x003FFFFF,
525 0x00FFFFFF, 0x03FFFFFF, 0x0FFFFFFF, 0x3FFFFFFF };
526
527 unsigned int iteration, wordToCount, charToCount;
528 unsigned int i, j, c;
529 unsigned int sum;
530
531 dna -= index / 16 + 1;
532
533 iteration = index / 256;
534 wordToCount = (index - iteration * 256) / 16;
535 charToCount = index - iteration * 256 - wordToCount * 16;
536
537 sum = 0;
538 if (charToCount > 0) {
539 c = *dna & truncateLeftMask[charToCount]; // increase count of 'a' by 16 - c;
540 sum += dnaDecodeTable[c >> 16];
541 sum += dnaDecodeTable[c & 0xFFFF];
542 sum += charToCount - 16; // decrease count of 'a' by 16 - positionToProcess
543 }
544
545 for (j=0; j<wordToCount; j++) {
546 dna++;
547 sum += dnaDecodeTable[*dna >> 16];
548 sum += dnaDecodeTable[*dna & 0x0000FFFF];
549 }
550
551 occCount[0] = sum & 0x000000FF; sum >>= 8;
552 occCount[1] = sum & 0x000000FF; sum >>= 8;
553 occCount[2] = sum & 0x000000FF; sum >>= 8;
554 occCount[3] = sum;
555
556 for (i=0; i<iteration; i++) {
557
558 sum = 0;
559 for (j=0; j<16; j++) {
560 dna++;
561 sum += dnaDecodeTable[*dna >> 16];
562 sum += dnaDecodeTable[*dna & 0x0000FFFF];
563 }
564 if (!DNA_OCC_SUM_EXCEPTION(sum)) {
565 occCount[0] += sum & 0x000000FF; sum >>= 8;
566 occCount[1] += sum & 0x000000FF; sum >>= 8;
567 occCount[2] += sum & 0x000000FF; sum >>= 8;
568 occCount[3] += sum;
569 } else {
570 // only some or all of the 3 bits are on
571 // in reality, only one of the four cases are possible
572 if (sum == 0x00000100) {
573 occCount[0] += 256;
574 } else if (sum == 0x00010000) {
575 occCount[1] += 256;
576 } else if (sum == 0x01000000) {
577 occCount[2] += 256;
578 } else if (sum == 0x00000000) {
579 occCount[3] += 256;
580 } else {
581 fprintf(stderr, "BackwardDNAAllOccCountNoLimit(): DNA occ sum exception!\n");
582 exit(1);
583 }
584 }
585
586 }
587
588 }
589
590
591 void GenerateDNA_NOccCountTable(unsigned int *dnaDecodeTable) {
592
593 unsigned int i, j, c, t;
594
595 for (i=0; i<DNA_N_OCC_CNT_TABLE_SIZE_IN_WORD; i++) {
596 dnaDecodeTable[i] = 0;
597 c = i;
598 for (j=0; j<5; j++) {
599 t = c & 0x00000007;
600 if (t != 4) {
601 // Count of 'n' is to be derived from acgt
602 dnaDecodeTable[i] += 1 << (t * 8);
603 }
604 c >>= 3;
605 }
606 }
607
608 }
609
610 unsigned int ForwardDNA_NOccCount(const unsigned int* dna, const unsigned int index, const unsigned int character, const unsigned int* dnaDecodeTable) {
611
612 static const unsigned int truncateRightMask[10] = { 0x00000000, 0xE0000000, 0xFC000000, 0xFF800000,
613 0xFFF00000, 0xFFFE0000, 0xFFFFC000, 0xFFFFF800,
614 0xFFFFFF00, 0xFFFFFFE0};
615
616 unsigned int wordToCount, charToCount;
617 unsigned int i, c;
618 unsigned int sum = 0;
619 unsigned int occCount;
620
621 #ifdef DEBUG
622 if (index > 250) {
623 fprintf(stderr, "ForwardDNA_NOccCount() : index > 250!\n");
624 exit(1);
625 }
626 #endif
627
628 wordToCount = index / 10;
629 charToCount = index - wordToCount * 10;
630
631 for (i=0; i<wordToCount; i++) {
632 sum += dnaDecodeTable[dna[i] >> 17];
633 sum += dnaDecodeTable[(dna[i] >> 2) & 0x00007FFF];
634 }
635
636 if (charToCount > 0) {
637 c = dna[i] & truncateRightMask[charToCount]; // increase count of 'a' by 10 - charToCount;
638 sum += dnaDecodeTable[c >> 17];
639 sum += dnaDecodeTable[(c >> 2) & 0x00007FFF];
640 sum += charToCount - 10; // decrease count of 'a' by 10 - charToCount
641 }
642
643 if (character != 4) {
644 occCount = (sum >> (character * 8)) & 0x000000FF;
645 } else {
646 occCount = index;
647 occCount -= sum & 0x000000FF; sum >>= 8;
648 occCount -= sum & 0x000000FF; sum >>= 8;
649 occCount -= sum & 0x000000FF; sum >>= 8;
650 occCount -= sum;
651 }
652
653 return occCount;
654
655 }
656
657 unsigned int BackwardDNA_NOccCount(const unsigned int* dna, const unsigned int index, const unsigned int character, const unsigned int* dnaDecodeTable) {
658
659 static const unsigned int truncateLeftMask[10] = { 0x00000000, 0x0000001C, 0x000000FC, 0x000007FC,
660 0x00003FFC, 0x0001FFFC, 0x000FFFFC, 0x007FFFFC,
661 0x03FFFFFC, 0x1FFFFFFC};
662
663 unsigned int wordToCount, charToCount;
664 unsigned int j, c;
665 unsigned int sum = 0;
666 unsigned int occCount;
667
668 #ifdef DEBUG
669 if (index > 250) {
670 fprintf(stderr, "BackwardDNA_NOccCount() : index >= 250!\n");
671 exit(1);
672 }
673 #endif
674
675 wordToCount = index / 10;
676 charToCount = index - wordToCount * 10;
677
678 dna -= wordToCount + 1;
679
680 if (charToCount > 0) {
681 c = *dna & truncateLeftMask[charToCount]; // increase count of 'a' by 10 - charToCount;
682 sum += dnaDecodeTable[c >> 17];
683 sum += dnaDecodeTable[(c >> 2) & 0x00007FFF];
684 sum += charToCount - 10; // decrease count of 'a' by 10 - charToCount
685 }
686
687 for (j=0; j<wordToCount; j++) {
688 dna++;
689 sum += dnaDecodeTable[*dna >> 17];
690 sum += dnaDecodeTable[(*dna >> 2) & 0x00007FFF];
691 }
692
693 if (character != 4) {
694 occCount = (sum >> (character * 8)) & 0x000000FF;
695 } else {
696 occCount = index;
697 occCount -= sum & 0x000000FF; sum >>= 8;
698 occCount -= sum & 0x000000FF; sum >>= 8;
699 occCount -= sum & 0x000000FF; sum >>= 8;
700 occCount -= sum;
701 }
702
703 #ifdef DEBUG
704 if (occCount > index + 1) {
705 fprintf(stderr, "BackwardDNA_NOccCount() : occCount > index + 1!\n");
706 exit(1);
707 }
708 #endif
709
710 return occCount;
711
712 }
713
714 void ForwardDNA_NAllOccCount(const unsigned int* dna, const unsigned int index, unsigned int* __restrict occCount, const unsigned int* dnaDecodeTable) {
715
716 static const unsigned int truncateRightMask[10] = { 0x00000000, 0xE0000000, 0xFC000000, 0xFF800000,
717 0xFFF00000, 0xFFFE0000, 0xFFFFC000, 0xFFFFF800,
718 0xFFFFFF00, 0xFFFFFFE0};
719
720 unsigned int wordToCount, charToCount;
721 unsigned int i, c;
722 unsigned int sum = 0;
723
724 #ifdef DEBUG
725 if (index > 250) {
726 fprintf(stderr, "ForwardDNA_NAllOccCount() : index >= 250!\n");
727 exit(1);
728 }
729 #endif
730
731 wordToCount = index / 10;
732 charToCount = index - wordToCount * 10;
733
734 for (i=0; i<wordToCount; i++) {
735 sum += dnaDecodeTable[dna[i] >> 17];
736 sum += dnaDecodeTable[(dna[i] >> 2) & 0x00007FFF];
737 }
738
739 if (charToCount > 0) {
740 c = dna[i] & truncateRightMask[charToCount]; // increase count of 'a' by 10 - charToCount;
741 sum += dnaDecodeTable[c >> 17];
742 sum += dnaDecodeTable[(c >> 2) & 0x00007FFF];
743 sum += charToCount - 10; // decrease count of 'a' by 10 - charToCount
744 }
745
746 occCount[0] = sum & 0x000000FF; sum >>= 8;
747 occCount[1] = sum & 0x000000FF; sum >>= 8;
748 occCount[2] = sum & 0x000000FF; sum >>= 8;
749 occCount[3] = sum;
750
751 }
752
753 void BackwardDNA_NAllOccCount(const unsigned int* dna, const unsigned int index, unsigned int* __restrict occCount, const unsigned int* dnaDecodeTable) {
754
755 static const unsigned int truncateLeftMask[10] = { 0x00000000, 0x0000001C, 0x000000FC, 0x000007FC,
756 0x00003FFC, 0x0001FFFC, 0x000FFFFC, 0x007FFFFC,
757 0x03FFFFFC, 0x1FFFFFFC};
758
759 unsigned int wordToCount, charToCount;
760 unsigned int j, c;
761 unsigned int sum = 0;
762
763 #ifdef DEBUG
764 if (index > 250) {
765 fprintf(stderr, "BackwardDNA_NAllOccCount() : index >= 250!\n");
766 exit(1);
767 }
768 #endif
769
770 wordToCount = index / 10;
771 charToCount = index - wordToCount * 10;
772
773 dna -= wordToCount + 1;
774
775 if (charToCount > 0) {
776 c = *dna & truncateLeftMask[charToCount]; // increase count of 'a' by 10 - charToCount;
777 sum += dnaDecodeTable[c >> 17];
778 sum += dnaDecodeTable[(c >> 2) & 0x00007FFF];
779 sum += charToCount - 10; // decrease count of 'a' by 16 - charToCount
780 }
781
782 for (j=0; j<wordToCount; j++) {
783 dna++;
784 sum += dnaDecodeTable[*dna >> 17];
785 sum += dnaDecodeTable[(*dna >> 2) & 0x00007FFF];
786 }
787
788 occCount[0] = sum & 0x000000FF; sum >>= 8;
789 occCount[1] = sum & 0x000000FF; sum >>= 8;
790 occCount[2] = sum & 0x000000FF; sum >>= 8;
791 occCount[3] = sum;
792
793 }
794
795 unsigned int ForwardDNA_NOccCountNoLimit(const unsigned int* dna, const unsigned int index, const unsigned int character, const unsigned int* dnaDecodeTable) {
796
797 static const unsigned int truncateRightMask[10] = { 0x00000000, 0xE0000000, 0xFC000000, 0xFF800000,
798 0xFFF00000, 0xFFFE0000, 0xFFFFC000, 0xFFFFF800,
799 0xFFFFFF00, 0xFFFFFFE0};
800
801 unsigned int iteration, wordToCount, charToCount;
802 unsigned int i, j, c;
803 unsigned int sum;
804 unsigned int occCount = 0;
805
806 iteration = index / 250;
807 wordToCount = (index - iteration * 250) / 10;
808 charToCount = index - iteration * 250 - wordToCount * 10;
809
810 for (i=0; i<iteration; i++) {
811
812 sum = 0;
813 for (j=0; j<25; j++) {
814 sum += dnaDecodeTable[*dna >> 17];
815 sum += dnaDecodeTable[(*dna >> 2) & 0x00007FFF];
816 dna++;
817 }
818 if (character != 4) {
819 occCount += (sum >> (character * 8)) & 0x000000FF;
820 } else {
821 occCount -= sum & 0x000000FF; sum >>= 8;
822 occCount -= sum & 0x000000FF; sum >>= 8;
823 occCount -= sum & 0x000000FF; sum >>= 8;
824 occCount -= sum;
825 }
826 }
827
828 sum = 0;
829 for (j=0; j<wordToCount; j++) {
830 sum += dnaDecodeTable[*dna >> 17];
831 sum += dnaDecodeTable[(*dna >> 2) & 0x00007FFF];
832 dna++;
833 }
834
835 if (charToCount > 0) {
836 c = *dna & truncateRightMask[charToCount]; // increase count of 'a' by 10 - charToCount;
837 sum += dnaDecodeTable[c >> 17];
838 sum += dnaDecodeTable[(c >> 2) & 0x00007FFF];
839 sum += charToCount - 10; // decrease count of 'a' by 10 - charToCount
840 }
841
842 if (character != 4) {
843 occCount += (sum >> (character * 8)) & 0x000000FF;
844 } else {
845 occCount += index;
846 occCount -= sum & 0x000000FF; sum >>= 8;
847 occCount -= sum & 0x000000FF; sum >>= 8;
848 occCount -= sum & 0x000000FF; sum >>= 8;
849 occCount -= sum;
850 }
851
852 return occCount;
853
854 }
855
856 unsigned int BackwardDNA_NOccCountNoLimit(const unsigned int* dna, const unsigned int index, const unsigned int character, const unsigned int* dnaDecodeTable) {
857
858 static const unsigned int truncateLeftMask[10] = { 0x00000000, 0x0000001C, 0x000000FC, 0x000007FC,
859 0x00003FFC, 0x0001FFFC, 0x000FFFFC, 0x007FFFFC,
860 0x03FFFFFC, 0x1FFFFFFC};
861
862 unsigned int iteration, wordToCount, charToCount;
863 unsigned int i, j, c;
864 unsigned int sum = 0;
865 unsigned int occCount = 0;
866
867 dna -= index / 10 + 1;
868
869 iteration = index / 250;
870 wordToCount = (index - iteration * 250) / 10;
871 charToCount = index - iteration * 250 - wordToCount * 10;
872
873 if (charToCount > 0) {
874 c = *dna & truncateLeftMask[charToCount]; // increase count of 'a' by 10 - charToCount;
875 sum += dnaDecodeTable[c >> 17];
876 sum += dnaDecodeTable[(c >> 2) & 0x00007FFF];
877 sum += charToCount - 10; // decrease count of 'a' by 10 - charToCount
878 }
879
880 for (j=0; j<wordToCount; j++) {
881 dna++;
882 sum += dnaDecodeTable[*dna >> 17];
883 sum += dnaDecodeTable[(*dna >> 2) & 0x00007FFF];
884 }
885
886 if (character != 4) {
887 occCount = (sum >> (character * 8)) & 0x000000FF;
888 } else {
889 occCount = index;
890 occCount -= sum & 0x000000FF; sum >>= 8;
891 occCount -= sum & 0x000000FF; sum >>= 8;
892 occCount -= sum & 0x000000FF; sum >>= 8;
893 occCount -= sum;
894 }
895
896 for (i=0; i<iteration; i++) {
897
898 sum = 0;
899 for (j=0; j<25; j++) {
900 dna++;
901 sum += dnaDecodeTable[*dna >> 17];
902 sum += dnaDecodeTable[(*dna >> 2) & 0x00007FFF];
903 }
904 if (character != 4) {
905 occCount += (sum >> (character * 8)) & 0x000000FF;
906 } else {
907 occCount -= sum & 0x000000FF; sum >>= 8;
908 occCount -= sum & 0x000000FF; sum >>= 8;
909 occCount -= sum & 0x000000FF; sum >>= 8;
910 occCount -= sum;
911 }
912 }
913
914 return occCount;
915
916 }
917
918 void ForwardDNA_NAllOccCountNoLimit(const unsigned int* dna, const unsigned int index, unsigned int* __restrict occCount, const unsigned int* dnaDecodeTable) {
919
920 static const unsigned int truncateRightMask[10] = { 0x00000000, 0xE0000000, 0xFC000000, 0xFF800000,
921 0xFFF00000, 0xFFFE0000, 0xFFFFC000, 0xFFFFF800,
922 0xFFFFFF00, 0xFFFFFFE0};
923
924 unsigned int iteration, wordToCount, charToCount;
925 unsigned int i, j, c;
926 unsigned int sum;
927
928 occCount[0] = 0;
929 occCount[1] = 0;
930 occCount[2] = 0;
931 occCount[3] = 0;
932
933 iteration = index / 250;
934 wordToCount = (index - iteration * 250) / 10;
935 charToCount = index - iteration * 250 - wordToCount * 10;
936
937 for (i=0; i<iteration; i++) {
938
939 sum = 0;
940 for (j=0; j<25; j++) {
941 sum += dnaDecodeTable[*dna >> 17];
942 sum += dnaDecodeTable[(*dna >> 2) & 0x00007FFF];
943 dna++;
944 }
945 occCount[0] += sum & 0x000000FF; sum >>= 8;
946 occCount[1] += sum & 0x000000FF; sum >>= 8;
947 occCount[2] += sum & 0x000000FF; sum >>= 8;
948 occCount[3] += sum;
949 }
950
951 sum = 0;
952 for (j=0; j<wordToCount; j++) {
953 sum += dnaDecodeTable[*dna >> 17];
954 sum += dnaDecodeTable[(*dna >> 2) & 0x00007FFF];
955 dna++;
956 }
957
958 if (charToCount > 0) {
959 c = *dna & truncateRightMask[charToCount]; // increase count of 'a' by 10 - charToCount;
960 sum += dnaDecodeTable[c >> 17];
961 sum += dnaDecodeTable[(c >> 2) & 0x00007FFF];
962 sum += charToCount - 10; // decrease count of 'a' by 10 - charToCount
963 }
964
965 occCount[0] += sum & 0x000000FF; sum >>= 8;
966 occCount[1] += sum & 0x000000FF; sum >>= 8;
967 occCount[2] += sum & 0x000000FF; sum >>= 8;
968 occCount[3] += sum;
969
970 }
971
972 void BackwardDNA_NAllOccCountNoLimit(const unsigned int* dna, const unsigned int index, unsigned int* __restrict occCount, const unsigned int* dnaDecodeTable) {
973
974 static const unsigned int truncateLeftMask[10] = { 0x00000000, 0x0000001C, 0x000000FC, 0x000007FC,
975 0x00003FFC, 0x0001FFFC, 0x000FFFFC, 0x007FFFFC,
976 0x03FFFFFC, 0x1FFFFFFC};
977
978 unsigned int iteration, wordToCount, charToCount;
979 unsigned int i, j, c;
980 unsigned int sum = 0;
981
982 dna -= index / 10 + 1;
983
984 iteration = index / 250;
985 wordToCount = (index - iteration * 250) / 10;
986 charToCount = index - iteration * 250 - wordToCount * 10;
987
988 if (charToCount > 0) {
989 c = *dna & truncateLeftMask[charToCount]; // increase count of 'a' by 10 - charToCount;
990 sum += dnaDecodeTable[c >> 17];
991 sum += dnaDecodeTable[(c >> 2) & 0x00007FFF];
992 sum += charToCount - 10; // decrease count of 'a' by 16 - charToCount
993 }
994
995 for (j=0; j<wordToCount; j++) {
996 dna++;
997 sum += dnaDecodeTable[*dna >> 17];
998 sum += dnaDecodeTable[(*dna >> 2) & 0x00007FFF];
999 }
1000
1001 occCount[0] = sum & 0x000000FF; sum >>= 8;
1002 occCount[1] = sum & 0x000000FF; sum >>= 8;
1003 occCount[2] = sum & 0x000000FF; sum >>= 8;
1004 occCount[3] = sum;
1005
1006 for (i=0; i<iteration; i++) {
1007
1008 sum = 0;
1009 for (j=0; j<25; j++) {
1010 dna++;
1011 sum += dnaDecodeTable[*dna >> 17];
1012 sum += dnaDecodeTable[(*dna >> 2) & 0x00007FFF];
1013 }
1014 occCount[0] += sum & 0x000000FF; sum >>= 8;
1015 occCount[1] += sum & 0x000000FF; sum >>= 8;
1016 occCount[2] += sum & 0x000000FF; sum >>= 8;
1017 occCount[3] += sum;
1018 }
1019
1020 }
1021
1022 unsigned int ForwardOccCount(const unsigned int* packed, const unsigned int index, const unsigned int character, const unsigned int alphabetSize) {
1023
1024 unsigned int wordToCount, charToCount;
1025 unsigned int bitPerChar, charPerWord;
1026 unsigned int i, j, c;
1027 unsigned int occCount = 0;
1028
1029 bitPerChar = ceilLog2(alphabetSize);
1030 charPerWord = BITS_IN_WORD / bitPerChar;
1031
1032 wordToCount = index / charPerWord;
1033 charToCount = index - wordToCount * charPerWord;
1034
1035 for (i=0; i<wordToCount; i++) {
1036 c = packed[i];
1037 for (j=0; j<charPerWord; j++) {
1038 if (c >> (BITS_IN_WORD - bitPerChar) == character) {
1039 occCount++;
1040 }
1041 c <<= bitPerChar;
1042 }
1043 }
1044 if (charToCount > 0) {
1045 c = packed[i];
1046 for (j=0; j<charToCount; j++) {
1047 if (c >> (BITS_IN_WORD - bitPerChar) == character) {
1048 occCount++;
1049 }
1050 c <<= bitPerChar;
1051 }
1052 }
1053
1054 return occCount;
1055
1056 }
1057
1058
1059 unsigned int BackwardOccCount(const unsigned int* packed, const unsigned int index, const unsigned int character, const unsigned int alphabetSize) {
1060
1061 unsigned int wordToCount, charToCount;
1062 unsigned int bitPerChar, charPerWord;
1063 unsigned int i, j, c;
1064 unsigned int occCount = 0;
1065
1066 bitPerChar = ceilLog2(alphabetSize);
1067 charPerWord = BITS_IN_WORD / bitPerChar;
1068
1069 wordToCount = index / charPerWord;
1070 charToCount = index - wordToCount * charPerWord;
1071
1072 packed -= wordToCount + 1;
1073
1074 if (charToCount > 0) {
1075 c = *packed << (bitPerChar * (charPerWord - charToCount));
1076 for (j=0; j<charToCount; j++) {
1077 if (c >> (BITS_IN_WORD - bitPerChar) == character) {
1078 occCount++;
1079 }
1080 c <<= bitPerChar;
1081 }
1082 }
1083
1084 for (i=1; i<=wordToCount; i++) {
1085 packed++;
1086 c = *packed;
1087 for (j=0; j<charPerWord; j++) {
1088 if (c >> (BITS_IN_WORD - bitPerChar) == character) {
1089 occCount++;
1090 }
1091 c <<= bitPerChar;
1092 }
1093 }
1094
1095 return occCount;
1096 }
1097
1098 void ForwardAllOccCount(const unsigned int* packed, const unsigned int index, const unsigned int alphabetSize, unsigned int* occCount) {
1099
1100 unsigned int wordToCount, charToCount;
1101 unsigned int bitPerChar, charPerWord;
1102 unsigned int i, j, c;
1103
1104 bitPerChar = ceilLog2(alphabetSize);
1105 charPerWord = BITS_IN_WORD / bitPerChar;
1106
1107 wordToCount = index / charPerWord;
1108 charToCount = index - wordToCount * charPerWord;
1109
1110 for (i=0; i<wordToCount; i++) {
1111 c = packed[i];
1112 for (j=0; j<charPerWord; j++) {
1113 occCount[c >> (BITS_IN_WORD - bitPerChar)]++;
1114 c <<= bitPerChar;
1115 }
1116 }
1117 if (charToCount > 0) {
1118 c = packed[i];
1119 for (j=0; j<charToCount; j++) {
1120 occCount[c >> (BITS_IN_WORD - bitPerChar)]++;
1121 c <<= bitPerChar;
1122 }
1123 }
1124
1125 }
1126
1127 void BackwardAllOccCount(const unsigned int* packed, const unsigned int index, const unsigned int alphabetSize, unsigned int* occCount) {
1128
1129 unsigned int wordToCount, charToCount;
1130 unsigned int bitPerChar, charPerWord;
1131 unsigned int i, j, c;
1132
1133 bitPerChar = ceilLog2(alphabetSize);
1134 charPerWord = BITS_IN_WORD / bitPerChar;
1135
1136 wordToCount = index / charPerWord;
1137 charToCount = index - wordToCount * charPerWord;
1138
1139 packed -= wordToCount + 1;
1140
1141 if (charToCount > 0) {
1142 c = *packed << (bitPerChar * (charPerWord - charToCount));
1143 for (j=0; j<charToCount; j++) {
1144 occCount[c >> (BITS_IN_WORD - bitPerChar)]++;
1145 c <<= bitPerChar;
1146 }
1147 }
1148
1149 for (i=1; i<=wordToCount; i++) {
1150 packed++;
1151 c = *packed;
1152 for (j=0; j<charPerWord; j++) {
1153 occCount[c >> (BITS_IN_WORD - bitPerChar)]++;
1154 c <<= bitPerChar;
1155 }
1156 }
1157
1158 }
1159
A place to share codes with mates