« back to all changes in this revision
Viewing changes to kalign2_dp.c
- browse files at revision 3
- compare with another revision
- download diff
- download tarball
- view history from revision 3
- Committer: Bazaar Package Importer
- Author(s): Charles Plessy
- Date: 2007-02-21 09:07:33 UTC
- mfrom: (1.1.1 upstream)
- Revision ID: james.westby@ubuntu.com-20070221090733-ttdj0bl7yes4to8f
Tags: 2.03-1
* New upstream release.
* Mended the watch file.
* Updated the manpage.
* Mended the watch file.
* Updated the manpage.
- files added:
- Makefile.in
- files removed:
- Makefile
- files modified:
- README
added
removed
kalign2_dp.c
1
/*
2
kalign2_dp.c
3
4
Released under GPL - see the 'COPYING' file
5
6
Copyright (C) 2006 Timo Lassmann <timolassmann@gmail.com>
7
8
This program is free software; you can redistribute it and/or modify
9
it under the terms of the GNU General Public License as published by
10
the Free Software Foundation; either version 2 of the License, or
11
any later version.
12
13
This program is distributed in the hope that it will be useful,
14
but WITHOUT ANY WARRANTY; without even the implied warranty of
15
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16
GNU General Public License for more details.
17
18
You should have received a copy of the GNU General Public License
19
along with this program; if not, write to the Free Software
20
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21
22
Please send bug reports, comments etc. to:
23
timolassmann@gmail.com
24
*/
25
26
#include "kalign2.h"
27
28
int* f_only_pp_dyn(int* path, struct dp_matrix *dp,const float* fprof1,const float* fprof2,const int len_a,const int len_b,int fdim,int stride)
29
{
30
// unsigned int freq[26];
31
32
struct states* s = 0;
33
char** trace = 0;
34
char* tracep = 0;
35
register float pa = 0;
36
register float pga = 0;
37
register float pgb = 0;
38
register float ca = 0;
39
register int i = 0;
40
register int j = 0;
41
register int c = 0;
42
register int f = 0;
43
44
s = dp->s;
45
46
trace = dp->tb;
47
48
trace[len_a][len_b] = 32;
49
50
fprof1 += len_a * stride;
51
52
53
s[len_b].a = 0.0;
54
s[len_b].ga = -FLOATINFTY;
55
s[len_b].gb = -FLOATINFTY;
56
//init of first row;
57
tracep = trace[len_a];
58
59
j = len_b;
60
while(--j){
61
s[j].a = -FLOATINFTY;
62
//s[j].ga = 0;
63
64
s[j].ga = s[j+1].a;//+prof2[29];
65
if (s[j+1].ga > s[j].ga){
66
s[j].ga = s[j+1].ga ;
67
}
68
s[j].gb = -FLOATINFTY;
69
tracep[j] = 8;
70
}
71
72
s[0].a = -FLOATINFTY;
73
s[0].ga = -FLOATINFTY;
74
s[0].gb = -FLOATINFTY;
75
i = len_a;
76
while(--i){
77
78
fprof1 -= stride;
79
80
tracep = trace[i];
81
pa = s[len_b].a;
82
pga = s[len_b].ga;
83
pgb = s[len_b].gb;
84
s[len_b].a = -FLOATINFTY;
85
s[len_b].ga = -FLOATINFTY;
86
//s[len_b].gb = 0;
87
88
s[len_b].gb = pa;//+prof1[29];
89
if(pgb > s[len_b].gb){
90
s[len_b].gb = pgb;
91
}
92
93
tracep[len_b] = 16;
94
95
j = len_b;
96
97
fprof2 += len_b * stride;
98
99
while(--j){
100
101
fprof2 -= stride;
102
ca = s[j].a;
103
104
c = 1;
105
if((pga) > pa){
106
pa = pga;
107
c = 2;
108
}
109
if((pgb) > pa){
110
pa = pgb;
111
c = 4;
112
}
113
114
for (f = 0; f < fdim;f++){
115
// fprintf(stderr,"%d %d: %d\n",i,j,fprof1[pga] * fprof2[pga+fdim]);
116
pa += fprof1[f] * fprof2[f+fdim];
117
}
118
119
s[j].a = pa;
120
121
pga = s[j].ga;
122
123
s[j].ga = s[j+1].a;
124
if (s[j+1].ga > s[j].ga){
125
s[j].ga = s[j+1].ga;
126
c |= 8;
127
}
128
129
pgb = s[j].gb;
130
131
s[j].gb = ca;
132
if(pgb > s[j].gb){
133
s[j].gb = pgb;
134
c |= 16;
135
}
136
tracep[j] = c;
137
pa = ca;
138
139
}
140
141
fprof2 -= stride;
142
//LAST CELL (0)
143
ca = s[0].a;
144
145
c = 1;
146
if((pga) > pa){
147
pa = pga;
148
c = 2;
149
}
150
if((pgb) > pa){
151
pa = pgb;
152
c = 4;
153
}
154
155
for (f = 0; f < fdim;f++){
156
// fprintf(stderr,"%d %d: %d\n",i,j,fprof1[pga] * fprof2[pga+fdim]);
157
pa += fprof1[f] * fprof2[f+fdim];
158
}
159
160
s[0].a = pa;
161
162
s[0].ga = -FLOATINFTY;
163
164
pgb = s[0].gb;
165
s[0].gb = ca;
166
if(pgb> s[0].gb){
167
s[0].gb = pgb;
168
c |= 16;
169
}
170
tracep[0] = c;
171
172
}
173
174
fprof1 -= stride;
175
176
tracep = trace[0];
177
j = len_b;
178
179
fprof2 += len_b * stride;
180
181
pa = s[j].a;
182
pga = s[j].ga;
183
pgb = s[j].gb;
184
s[j].a = -FLOATINFTY;
185
s[j].ga = -FLOATINFTY;
186
//s[j].gb = -INFTY;
187
s[len_b].gb = pa;//+prof1[29];
188
if(pgb > s[len_b].gb){
189
s[len_b].gb = pgb;
190
}
191
192
193
194
while(--j){
195
196
fprof2 -= stride;
197
198
ca = s[j].a;
199
200
c = 1;
201
202
if((pga) > pa){
203
pa = pga;
204
c = 2;
205
}
206
207
if((pgb) > pa){
208
pa = pgb;
209
c = 4;
210
}
211
212
for (f = 0; f < fdim;f++){
213
pa += fprof1[f] * fprof2[f+fdim];
214
// fprintf(stderr,"%d %d: %d\n",i,j,fprof1[pga] * fprof2[pga+fdim]);
215
}
216
217
218
s[j].a = pa;
219
pga = s[j].ga;
220
s[j].ga = s[j+1].a;
221
if (s[j+1].ga > s[j].ga){
222
s[j].ga = s[j+1].ga;
223
c |= 8;
224
}
225
pgb = s[j].gb;
226
s[j].gb = -FLOATINFTY;
227
228
tracep[j] = c;
229
pa = ca;
230
}
231
232
fprof2 -= stride;
233
234
ca = s[0].a;
235
236
c = 1;
237
238
if((pga) > pa){
239
pa = pga;
240
c = 2;
241
}
242
if((pgb) > pa){
243
pa = pgb;
244
c = 4;
245
}
246
247
for (f = 0; f < fdim;f++){
248
pa += fprof1[f] * fprof2[f+fdim];
249
// fprintf(stderr,"%d %d: %d\n",i,j,fprof1[pga] * fprof2[pga+fdim]);
250
}
251
252
s[0].a = pa;
253
254
s[0].ga = s[1].a;
255
if (s[1].ga > s[0].ga){
256
s[0].ga = s[1].ga;
257
c |= 8;
258
}
259
260
pgb = s[0].gb;
261
s[0].gb = ca;
262
if(pgb> s[0].gb){
263
s[0].gb = pgb;
264
c |= 16;
265
}
266
tracep[0] = c;
267
268
pgb = s[0].gb;
269
c = 2;
270
if(s[0].ga > pgb){
271
pgb = s[0].ga;
272
c = 1;
273
}
274
if(s[0].a >= pgb){
275
pgb = s[0].a;
276
c = 0;
277
}
278
279
//fprintf(stderr,"SCORE:%d\n",ca);
280
ca = c;
281
282
i = 0;
283
j = 0;
284
f = 1;
285
while(trace[i][j] < 32){
286
// fprintf(stderr,"%d->%d %d:%d %d:%d\n",c,trace[i][j],i,j,len_a,len_b);
287
switch(f){
288
case 0:
289
if (trace[i][j] & 2){
290
f = 1;
291
if(i+1!= len_a){
292
path[c+1] |= 16;
293
// fprintf(stderr,"GAP_CLOSE\n");
294
}else{
295
path[c+1] |= 32+16;
296
}
297
}else if (trace[i][j] & 4){
298
f = 2;
299
if(j+1!= len_b){
300
path[c+1] |= 16;
301
// fprintf(stderr,"GAP_CLOSE\n");
302
}else{
303
path[c+1] |= 32+16;
304
}
305
}
306
307
//path[c] = 0;
308
i++;
309
j++;
310
break;
311
case 1:
312
if(trace[i][j] & 8){
313
f = 1;
314
if(i!=0 && i!= len_a){
315
// / fprintf(stderr,"GAP_EXT\n");
316
if(!(path[c]&16)){
317
path[c] |= 8;
318
}
319
}else{
320
if(!(path[c]&16)){
321
path[c] |= 32+8;
322
}
323
}
324
}else{
325
f = 0;
326
if(i!=0 && i!= len_a){
327
// fprintf(stderr,"GAP_OPEN\n");
328
path[c] |= 4;
329
}else{
330
path[c] |= 32+4;
331
}
332
}
333
path[c] |= 1;
334
j++;
335
break;
336
case 2:
337
if(trace[i][j] & 16){
338
f = 2;
339
if(j !=0 && j != len_b){
340
// fprintf(stderr,"GAP_EXT\n");
341
if(!(path[c]&16)){
342
path[c] |= 8;
343
}
344
}else{
345
if(!(path[c]&16)){
346
path[c] |= 32+8;
347
}
348
}
349
}else{
350
f = 0;
351
if(j!=0 && j != len_b){
352
// fprintf(stderr,"GAP_OPEN\n");
353
path[c] |= 4;
354
}else{
355
path[c] |= 32+4;
356
}
357
358
}
359
path[c] |= 2;
360
i++;
361
break;
362
}
363
c++;
364
}
365
path[0] = c-1;
366
path[c] = 3;
367
path[c+1] = pgb;
368
return path;
369
}
370
371
372
373
int* fpp_dyn(int* path, struct dp_matrix *dp,const float* prof1,const float* prof2,const float* fprof1,const float* fprof2,const int len_a,const int len_b,int fdim,int stride)
374
{
375
unsigned int freq[26];
376
377
struct states* s = 0;
378
char** trace = 0;
379
char* tracep = 0;
380
register float pa = 0;
381
register float pga = 0;
382
register float pgb = 0;
383
register float ca = 0;
384
register int i = 0;
385
register int j = 0;
386
register int c = 0;
387
register int f = 0;
388
389
390
s = dp->s;
391
392
trace = dp->tb;
393
394
trace[len_a][len_b] = 32;
395
396
prof1 += len_a << 6;
397
398
fprof1 += len_a * stride;
399
400
401
s[len_b].a = 0;
402
s[len_b].ga = -FLOATINFTY;
403
s[len_b].gb = -FLOATINFTY;
404
//init of first row;
405
tracep = trace[len_a];
406
407
j = len_b;
408
while(--j){
409
s[j].a = -FLOATINFTY;
410
//s[j].ga = 0;
411
412
s[j].ga = s[j+1].a+prof2[29];//+prof2[29];
413
if (s[j+1].ga+prof2[29] > s[j].ga){
414
s[j].ga = s[j+1].ga+prof2[29];
415
}
416
s[j].gb = -FLOATINFTY;
417
tracep[j] = 8;
418
}
419
420
s[0].a = -FLOATINFTY;
421
s[0].ga = -FLOATINFTY;
422
s[0].gb = -FLOATINFTY;
423
i = len_a;
424
while(--i){
425
prof1 -= 64;
426
427
fprof1 -= stride;
428
429
c = 1;
430
for (j = 26; j--;){
431
if(prof1[j]){
432
freq[c] = j;
433
c++;
434
}
435
}
436
freq[0] = c;
437
438
tracep = trace[i];
439
pa = s[len_b].a;
440
pga = s[len_b].ga;
441
pgb = s[len_b].gb;
442
s[len_b].a = -FLOATINFTY;
443
s[len_b].ga = -FLOATINFTY;
444
//s[len_b].gb = 0;
445
446
s[len_b].gb = pa+prof1[29];//+prof1[29];
447
if(pgb+prof1[29] > s[len_b].gb){
448
s[len_b].gb = pgb+prof1[29];
449
}
450
451
tracep[len_b] = 16;
452
453
j = len_b;
454
prof2 += len_b << 6;
455
456
fprof2 += len_b * stride;
457
458
while(--j){
459
prof2 -= 64;
460
461
fprof2 -= stride;
462
ca = s[j].a;
463
464
c = 1;
465
if((pga += prof2[91]) > pa){
466
pa = pga;
467
c = 2;
468
}
469
if((pgb += prof1[91]) > pa){
470
pa = pgb;
471
c = 4;
472
}
473
474
prof2 += 32;
475
for (f = freq[0];--f;){
476
pa += prof1[freq[f]]*prof2[freq[f]];
477
}
478
prof2 -= 32;
479
480
for (f = 0; f < fdim;f++){
481
// fprintf(stderr,"%d %d: %d\n",i,j,fprof1[pga] * fprof2[pga+fdim]);
482
pa += fprof1[f] * fprof2[f+fdim];
483
}
484
485
s[j].a = pa;
486
487
pga = s[j].ga;
488
489
s[j].ga = s[j+1].a+prof2[27];
490
if (s[j+1].ga+prof2[28] > s[j].ga){
491
s[j].ga = s[j+1].ga+prof2[28];
492
c |= 8;
493
}
494
495
pgb = s[j].gb;
496
497
s[j].gb = ca+prof1[27];
498
if(pgb+prof1[28] > s[j].gb){
499
s[j].gb = pgb+prof1[28];
500
c |= 16;
501
}
502
tracep[j] = c;
503
pa = ca;
504
505
}
506
507
prof2 -= 64;
508
509
fprof2 -= stride;
510
//LAST CELL (0)
511
ca = s[0].a;
512
513
c = 1;
514
if((pga+=prof2[91]) > pa){
515
pa = pga;
516
c = 2;
517
}
518
if((pgb+=prof1[91]) > pa){
519
pa = pgb;
520
c = 4;
521
}
522
523
prof2 += 32;
524
for (f = freq[0];--f;){
525
pa += prof1[freq[f]]*prof2[freq[f]];
526
}
527
prof2 -= 32;
528
529
for (f = 0; f < fdim;f++){
530
// fprintf(stderr,"%d %d: %d\n",i,j,fprof1[pga] * fprof2[pga+fdim]);
531
pa += fprof1[f] * fprof2[f+fdim];
532
}
533
534
s[0].a = pa;
535
536
s[0].ga = -FLOATINFTY;
537
538
pgb = s[0].gb;
539
s[0].gb = ca+prof1[27]+prof1[29];
540
if(pgb+prof1[29] > s[0].gb){
541
s[0].gb = pgb+prof1[29];
542
c |= 16;
543
}
544
tracep[0] = c;
545
546
}
547
prof1 -= 64;
548
549
fprof1 -= stride;
550
551
c = 1;
552
for (j = 26; j--;){
553
if(prof1[j]){
554
freq[c] = j;
555
c++;
556
}
557
}
558
freq[0] = c;
559
560
tracep = trace[0];
561
j = len_b;
562
prof2 += len_b << 6;
563
564
fprof2 += len_b * stride;
565
566
pa = s[j].a;
567
pga = s[j].ga;
568
pgb = s[j].gb;
569
s[j].a = -FLOATINFTY;
570
s[j].ga = -FLOATINFTY;
571
//s[j].gb = -INFTY;
572
s[len_b].gb = pa+prof1[29];//+prof1[29];
573
if(pgb+prof1[29] > s[len_b].gb){
574
s[len_b].gb = pgb+prof1[29];
575
}
576
577
578
579
while(--j){
580
prof2 -= 64;
581
582
fprof2 -= stride;
583
584
ca = s[j].a;
585
586
c = 1;
587
588
if((pga+=prof2[91]) > pa){
589
pa = pga;
590
c = 2;
591
}
592
593
if((pgb+=prof1[91]) > pa){
594
pa = pgb;
595
c = 4;
596
}
597
598
prof2+=32;
599
for (f = freq[0];--f;){
600
pa += prof1[freq[f]]*prof2[freq[f]];
601
}
602
prof2-=32;
603
604
605
for (f = 0; f < fdim;f++){
606
pa += fprof1[f] * fprof2[f+fdim];
607
// fprintf(stderr,"%d %d: %d\n",i,j,fprof1[pga] * fprof2[pga+fdim]);
608
}
609
610
611
s[j].a = pa;
612
pga = s[j].ga;
613
s[j].ga = s[j+1].a+prof2[27]+prof2[29];
614
if (s[j+1].ga+prof2[29] > s[j].ga){
615
s[j].ga = s[j+1].ga+prof2[29];
616
c |= 8;
617
}
618
pgb = s[j].gb;
619
s[j].gb = -FLOATINFTY;
620
621
tracep[j] = c;
622
pa = ca;
623
}
624
prof2 -= 64;
625
626
fprof2 -= stride;
627
628
ca = s[0].a;
629
630
c = 1;
631
632
if((pga+=prof2[91]) > pa){
633
pa = pga;
634
c = 2;
635
}
636
if((pgb+=prof1[91]) > pa){
637
pa = pgb;
638
c = 4;
639
}
640
prof2+=32;
641
for (f = freq[0];--f;){
642
pa += prof1[freq[f]]*prof2[freq[f]];
643
}
644
prof2-=32;
645
646
for (f = 0; f < fdim;f++){
647
pa += fprof1[f] * fprof2[f+fdim];
648
// fprintf(stderr,"%d %d: %d\n",i,j,fprof1[pga] * fprof2[pga+fdim]);
649
}
650
651
s[0].a = pa;
652
653
s[0].ga = s[1].a+prof2[27]+prof2[29];
654
if (s[1].ga+prof2[29] > s[0].ga){
655
s[0].ga = s[1].ga+prof2[29];
656
c |= 8;
657
}
658
659
pgb = s[0].gb;
660
s[0].gb = ca+prof1[27]+prof1[29];
661
if(pgb +prof1[29]> s[0].gb){
662
s[0].gb = pgb+prof1[29];
663
c |= 16;
664
}
665
tracep[0] = c;
666
667
pgb = s[0].gb;
668
c = 2;
669
if(s[0].ga > pgb){
670
pgb = s[0].ga;
671
c = 1;
672
}
673
if(s[0].a >= pgb){
674
pgb = s[0].a;
675
c = 0;
676
}
677
678
//fprintf(stderr,"SCORE:%d\n",ca);
679
f = c;
680
681
i = 0;
682
j = 0;
683
c = 1;
684
while(trace[i][j] < 32){
685
// fprintf(stderr,"%d->%d %d:%d %d:%d\n",c,trace[i][j],i,j,len_a,len_b);
686
switch(f){
687
case 0:
688
if (trace[i][j] & 2){
689
f = 1;
690
if(i+1!= len_a){
691
path[c+1] |= 16;
692
// fprintf(stderr,"GAP_CLOSE\n");
693
}else{
694
path[c+1] |= 32+16;
695
}
696
}else if (trace[i][j] & 4){
697
f = 2;
698
if(j+1!= len_b){
699
path[c+1] |= 16;
700
// fprintf(stderr,"GAP_CLOSE\n");
701
}else{
702
path[c+1] |= 32+16;
703
}
704
}
705
706
//path[c] = 0;
707
i++;
708
j++;
709
break;
710
case 1:
711
if(trace[i][j] & 8){
712
f = 1;
713
if(i!=0 && i!= len_a){
714
// / fprintf(stderr,"GAP_EXT\n");
715
if(!(path[c]&16)){
716
path[c] |= 8;
717
}
718
}else{
719
if(!(path[c]&16)){
720
path[c] |= 32+8;
721
}
722
}
723
}else{
724
f = 0;
725
if(i!=0 && i!= len_a){
726
// fprintf(stderr,"GAP_OPEN\n");
727
path[c] |= 4;
728
}else{
729
path[c] |= 32+4;
730
}
731
}
732
path[c] |= 1;
733
j++;
734
break;
735
case 2:
736
if(trace[i][j] & 16){
737
f = 2;
738
if(j !=0 && j != len_b){
739
// fprintf(stderr,"GAP_EXT\n");
740
if(!(path[c]&16)){
741
path[c] |= 8;
742
}
743
}else{
744
if(!(path[c]&16)){
745
path[c] |= 32+8;
746
}
747
}
748
}else{
749
f = 0;
750
if(j!=0 && j != len_b){
751
// fprintf(stderr,"GAP_OPEN\n");
752
path[c] |= 4;
753
}else{
754
path[c] |= 32+4;
755
}
756
757
}
758
path[c] |= 2;
759
i++;
760
break;
761
}
762
c++;
763
}
764
path[0] = c-1;
765
path[c] = 3;
766
path[c+1] = pgb;
767
return path;
768
}
769
/*
770
int* dna_pp_dyn(int* path, struct dp_matrix *dp,const int* prof1,const int* prof2,const int len_a,const int len_b)
771
{
772
773
struct states* s = 0;
774
char** trace = 0;
775
char* tracep = 0;
776
register int pa = 0;
777
register int pga = 0;
778
register int pgb = 0;
779
register int ca = 0;
780
register int i = 0;
781
register int j = 0;
782
register int c = 0;
783
784
s = dp->s;
785
786
trace = dp->tb;
787
788
trace[len_a][len_b] = 32;
789
790
prof1 += len_a * 22;
791
792
s[len_b].a = 0;
793
s[len_b].ga = -INFTY;
794
s[len_b].gb = -INFTY;
795
//init of first row;
796
tracep = trace[len_a];
797
798
j = len_b;
799
while(--j){
800
s[j].a = -INFTY;
801
//s[j].ga = 0;
802
803
s[j].ga = s[j+1].a+prof2[10];//+prof2[29];
804
if (s[j+1].ga+prof2[10] > s[j].ga){
805
s[j].ga = s[j+1].ga+prof2[10];
806
}
807
s[j].gb = -INFTY;
808
tracep[j] = 8;
809
}
810
811
s[0].a = -INFTY;
812
s[0].ga = -INFTY;
813
s[0].gb = -INFTY;
814
i = len_a;
815
while(--i){
816
prof1 -= 22;
817
818
tracep = trace[i];
819
pa = s[len_b].a;
820
pga = s[len_b].ga;
821
pgb = s[len_b].gb;
822
s[len_b].a = -INFTY;
823
s[len_b].ga = -INFTY;
824
//s[len_b].gb = 0;
825
826
s[len_b].gb = pa+prof1[10];//+prof1[29];
827
if(pgb+prof1[10] > s[len_b].gb){
828
s[len_b].gb = pgb+prof1[10];
829
}
830
831
tracep[len_b] = 16;
832
833
j = len_b;
834
prof2 += len_b *22;
835
while(--j){
836
prof2 -= 22;
837
ca = s[j].a;
838
839
c = 1;
840
if((pga += prof2[30]) > pa){
841
pa = pga;
842
c = 2;
843
}
844
if((pgb += prof1[30]) > pa){
845
pa = pgb;
846
c = 4;
847
}
848
849
prof2 += 11;
850
for (pga = 8;pga--;){
851
pa += prof1[pga]*prof2[pga];
852
}
853
prof2 -= 11;
854
855
s[j].a = pa;
856
857
pga = s[j].ga;
858
859
s[j].ga = s[j+1].a+prof2[8];
860
if (s[j+1].ga+prof2[9] > s[j].ga){
861
s[j].ga = s[j+1].ga+prof2[9];
862
c |= 8;
863
}
864
865
pgb = s[j].gb;
866
867
s[j].gb = ca+prof1[8];
868
if(pgb+prof1[9] > s[j].gb){
869
s[j].gb = pgb+prof1[9];
870
c |= 16;
871
}
872
tracep[j] = c;
873
pa = ca;
874
875
}
876
877
prof2 -= 22;
878
//LAST CELL (0)
879
ca = s[0].a;
880
881
c = 1;
882
if((pga+=prof2[30]) > pa){
883
pa = pga;
884
c = 2;
885
}
886
if((pgb+=prof1[30]) > pa){
887
pa = pgb;
888
c = 4;
889
}
890
891
prof2 += 11;
892
for (pga = 8;pga--;){
893
pa += prof1[pga]*prof2[pga];
894
}
895
prof2 -= 11;
896
897
s[0].a = pa;
898
899
s[0].ga = -INFTY;
900
901
pgb = s[0].gb;
902
s[0].gb = ca+prof1[8]+prof1[10];
903
if(pgb+prof1[10] > s[0].gb){
904
s[0].gb = pgb+prof1[10];
905
c |= 16;
906
}
907
tracep[0] = c;
908
909
}
910
prof1 -= 22;
911
912
tracep = trace[0];
913
j = len_b;
914
prof2 += len_b *22;
915
pa = s[j].a;
916
pga = s[j].ga;
917
pgb = s[j].gb;
918
s[j].a = -INFTY;
919
s[j].ga = -INFTY;
920
//s[j].gb = -INFTY;
921
s[len_b].gb = pa+prof1[10];//+prof1[29];
922
if(pgb+prof1[10] > s[len_b].gb){
923
s[len_b].gb = pgb+prof1[10];
924
}
925
926
927
928
while(--j){
929
prof2 -= 22;
930
ca = s[j].a;
931
932
c = 1;
933
934
if((pga+=prof2[30]) > pa){
935
pa = pga;
936
c = 2;
937
}
938
939
if((pgb+=prof1[30]) > pa){
940
pa = pgb;
941
c = 4;
942
}
943
944
prof2+=11;
945
946
for (pga = 8;pga--;){
947
pa += prof1[pga]*prof2[pga];
948
}
949
prof2-=11;
950
951
s[j].a = pa;
952
pga = s[j].ga;
953
s[j].ga = s[j+1].a+prof2[2]+prof2[10];
954
if (s[j+1].ga+prof2[10] > s[j].ga){
955
s[j].ga = s[j+1].ga+prof2[10];
956
c |= 8;
957
}
958
pgb = s[j].gb;
959
s[j].gb = -INFTY;
960
961
tracep[j] = c;
962
pa = ca;
963
}
964
prof2 -= 22;
965
966
ca = s[0].a;
967
968
c = 1;
969
970
if((pga+=prof2[30]) > pa){
971
pa = pga;
972
c = 2;
973
}
974
if((pgb+=prof1[30]) > pa){
975
pa = pgb;
976
c = 4;
977
}
978
prof2+=11;
979
for (pga = 8;pga--;){
980
pa += prof1[pga]*prof2[pga];
981
}
982
prof2-=11;
983
984
s[0].a = pa;
985
986
s[0].ga = s[1].a+prof2[8]+prof2[10];
987
if (s[1].ga+prof2[10] > s[0].ga){
988
s[0].ga = s[1].ga+prof2[10];
989
c |= 8;
990
}
991
992
pgb = s[0].gb;
993
s[0].gb = ca+prof1[8]+prof1[10];
994
if(pgb +prof1[10]> s[0].gb){
995
s[0].gb = pgb+prof1[10];
996
c |= 16;
997
}
998
tracep[0] = c;
999
1000
pgb = s[0].gb;
1001
c = 2;
1002
if(s[0].ga > pgb){
1003
pgb = s[0].ga;
1004
c = 1;
1005
}
1006
if(s[0].a >= pgb){
1007
pgb = s[0].a;
1008
c = 0;
1009
}
1010
1011
//fprintf(stderr,"SCORE:%d\n",ca);
1012
ca = c;
1013
1014
i = 0;
1015
j = 0;
1016
c = 1;
1017
while(trace[i][j] < 32){
1018
// fprintf(stderr,"%d->%d %d:%d %d:%d\n",c,trace[i][j],i,j,len_a,len_b);
1019
switch(ca){
1020
case 0:
1021
if (trace[i][j] & 2){
1022
ca = 1;
1023
if(i+1!= len_a){
1024
path[c+1] |= 16;
1025
// fprintf(stderr,"GAP_CLOSE\n");
1026
}else{
1027
path[c+1] |= 32+16;
1028
}
1029
}else if (trace[i][j] & 4){
1030
ca = 2;
1031
if(j+1!= len_b){
1032
path[c+1] |= 16;
1033
// fprintf(stderr,"GAP_CLOSE\n");
1034
}else{
1035
path[c+1] |= 32+16;
1036
}
1037
}
1038
1039
//path[c] = 0;
1040
i++;
1041
j++;
1042
break;
1043
case 1:
1044
if(trace[i][j] & 8){
1045
ca = 1;
1046
if(i!=0 && i!= len_a){
1047
// / fprintf(stderr,"GAP_EXT\n");
1048
if(!(path[c]&16)){
1049
path[c] |= 8;
1050
}
1051
}else{
1052
if(!(path[c]&16)){
1053
path[c] |= 32+8;
1054
}
1055
}
1056
}else{
1057
ca = 0;
1058
if(i!=0 && i!= len_a){
1059
// fprintf(stderr,"GAP_OPEN\n");
1060
path[c] |= 4;
1061
}else{
1062
path[c] |= 32+4;
1063
}
1064
}
1065
path[c] |= 1;
1066
j++;
1067
break;
1068
case 2:
1069
if(trace[i][j] & 16){
1070
ca = 2;
1071
if(j !=0 && j != len_b){
1072
// fprintf(stderr,"GAP_EXT\n");
1073
if(!(path[c]&16)){
1074
path[c] |= 8;
1075
}
1076
}else{
1077
if(!(path[c]&16)){
1078
path[c] |= 32+8;
1079
}
1080
}
1081
}else{
1082
ca = 0;
1083
if(j!=0 && j != len_b){
1084
// fprintf(stderr,"GAP_OPEN\n");
1085
path[c] |= 4;
1086
}else{
1087
path[c] |= 32+4;
1088
}
1089
1090
}
1091
path[c] |= 2;
1092
i++;
1093
break;
1094
}
1095
c++;
1096
}
1097
path[0] = c-1;
1098
path[c] = 3;
1099
path[c+1] = pgb;
1100
return path;
1101
}
1102
1103
1104
int* pp_dyn2(int* path, struct dp_matrix *dp,const int* prof1,const int* prof2,const int len_a,const int len_b)
1105
{
1106
unsigned int freq[26];
1107
1108
struct states* s = 0;
1109
char** trace = 0;
1110
char* tracep = 0;
1111
register int pa = 0;
1112
register int pga = 0;
1113
register int pgb = 0;
1114
register int ca = 0;
1115
register int i = 0;
1116
register int j = 0;
1117
register int c = 0;
1118
1119
s = dp->s;
1120
1121
trace = dp->tb;
1122
1123
trace[len_a][len_b] = 32;
1124
1125
prof1 += len_a << 6;
1126
1127
s[len_b].a = 0;
1128
s[len_b].ga = -INFTY;
1129
s[len_b].gb = -INFTY;
1130
//init of first row;
1131
tracep = trace[len_a];
1132
1133
j = len_b;
1134
while(--j){
1135
s[j].a = -INFTY;
1136
1137
s[j].ga = s[j+1].a+prof2[29];//+prof2[29];
1138
if (s[j+1].ga+prof2[29] > s[j].ga){
1139
s[j].ga = s[j+1].ga+prof2[29];
1140
}
1141
s[j].gb = -INFTY;
1142
tracep[j] = 8;
1143
}
1144
1145
s[0].a = -INFTY;
1146
s[0].ga = -INFTY;
1147
s[0].gb = -INFTY;
1148
i = len_a;
1149
while(--i){
1150
prof1 -= 64;
1151
1152
c = 1;
1153
for (j = 23; j--;){
1154
if(prof1[j]){
1155
freq[c] = j;
1156
c++;
1157
}
1158
}
1159
freq[0] = c;
1160
1161
tracep = trace[i];
1162
pa = s[len_b].a;
1163
pga = s[len_b].ga;
1164
pgb = s[len_b].gb;
1165
s[len_b].a = -INFTY;
1166
s[len_b].ga = -INFTY;
1167
1168
s[len_b].gb = pa+prof1[29];
1169
if(pgb+prof1[29] > s[len_b].gb){
1170
s[len_b].gb = pgb+prof1[29];
1171
}
1172
1173
tracep[len_b] = 16;
1174
1175
j = len_b;
1176
prof2 += len_b << 6;
1177
while(--j){
1178
prof2 -= 64;
1179
ca = s[j].a;
1180
1181
c = 1;
1182
if((pga += prof2[91]) > pa){
1183
pa = pga;
1184
c = 2;
1185
}
1186
if((pgb += prof1[91]) > pa){
1187
pa = pgb;
1188
c = 4;
1189
}
1190
1191
prof2 += 32;
1192
for (pga = freq[0];--pga;){
1193
pgb = freq[pga];
1194
pa += prof1[pgb]*prof2[pgb];
1195
}
1196
prof2 -= 32;
1197
1198
s[j].a = pa;
1199
1200
pga = s[j].ga;
1201
1202
s[j].ga = s[j+1].a+prof2[27];
1203
if (s[j+1].ga+prof2[28] > s[j].ga){
1204
s[j].ga = s[j+1].ga+prof2[28];
1205
c |= 8;
1206
}
1207
1208
pgb = s[j].gb;
1209
1210
s[j].gb = ca+prof1[27];
1211
if(pgb+prof1[28] > s[j].gb){
1212
s[j].gb = pgb+prof1[28];
1213
c |= 16;
1214
}
1215
tracep[j] = c;
1216
pa = ca;
1217
1218
}
1219
1220
prof2 -= 64;
1221
//LAST CELL (0)
1222
ca = s[0].a;
1223
1224
c = 1;
1225
if((pga+=prof2[91]) > pa){
1226
pa = pga;
1227
c = 2;
1228
}
1229
if((pgb+=prof1[91]) > pa){
1230
pa = pgb;
1231
c = 4;
1232
}
1233
1234
prof2 += 32;
1235
for (pga = freq[0];--pga;){
1236
pgb = freq[pga];
1237
pa += prof1[pgb]*prof2[pgb];
1238
}
1239
prof2 -= 32;
1240
1241
s[0].a = pa;
1242
1243
s[0].ga = -INFTY;
1244
1245
pgb = s[0].gb;
1246
s[0].gb = ca+prof1[27]+prof1[29];
1247
if(pgb+prof1[29] > s[0].gb){
1248
s[0].gb = pgb+prof1[29];
1249
c |= 16;
1250
}
1251
tracep[0] = c;
1252
1253
}
1254
prof1 -= 64;
1255
1256
c = 1;
1257
for (j = 23; j--;){
1258
if(prof1[j]){
1259
freq[c] = j;
1260
c++;
1261
}
1262
}
1263
freq[0] = c;
1264
1265
tracep = trace[0];
1266
j = len_b;
1267
prof2 += len_b << 6;
1268
pa = s[j].a;
1269
pga = s[j].ga;
1270
pgb = s[j].gb;
1271
s[j].a = -INFTY;
1272
s[j].ga = -INFTY;
1273
1274
s[len_b].gb = pa+prof1[29];
1275
if(pgb+prof1[29] > s[len_b].gb){
1276
s[len_b].gb = pgb+prof1[29];
1277
}
1278
1279
1280
1281
while(--j){
1282
prof2 -= 64;
1283
ca = s[j].a;
1284
1285
c = 1;
1286
1287
if((pga+=prof2[91]) > pa){
1288
pa = pga;
1289
c = 2;
1290
}
1291
1292
if((pgb+=prof1[91]) > pa){
1293
pa = pgb;
1294
c = 4;
1295
}
1296
1297
prof2+=32;
1298
1299
for (pga = freq[0];--pga;){
1300
pgb = freq[pga];
1301
pa += prof1[pgb]*prof2[pgb];
1302
}
1303
prof2-=32;
1304
1305
s[j].a = pa;
1306
pga = s[j].ga;
1307
s[j].ga = s[j+1].a+prof2[27]+prof2[29];
1308
if (s[j+1].ga+prof2[29] > s[j].ga){
1309
s[j].ga = s[j+1].ga+prof2[29];
1310
c |= 8;
1311
}
1312
pgb = s[j].gb;
1313
s[j].gb = -INFTY;
1314
1315
tracep[j] = c;
1316
pa = ca;
1317
}
1318
prof2 -= 64;
1319
1320
ca = s[0].a;
1321
1322
c = 1;
1323
1324
if((pga+=prof2[91]) > pa){
1325
pa = pga;
1326
c = 2;
1327
}
1328
if((pgb+=prof1[91]) > pa){
1329
pa = pgb;
1330
c = 4;
1331
}
1332
prof2+=32;
1333
for (pga = freq[0];--pga;){
1334
pgb = freq[pga];
1335
pa += prof1[pgb]*prof2[pgb];
1336
}
1337
prof2-=32;
1338
1339
s[0].a = pa;
1340
1341
s[0].ga = s[1].a+prof2[27]+prof2[29];
1342
if (s[1].ga+prof2[29] > s[0].ga){
1343
s[0].ga = s[1].ga+prof2[29];
1344
c |= 8;
1345
}
1346
1347
pgb = s[0].gb;
1348
s[0].gb = ca+prof1[27]+prof1[29];
1349
if(pgb +prof1[29]> s[0].gb){
1350
s[0].gb = pgb+prof1[29];
1351
c |= 16;
1352
}
1353
tracep[0] = c;
1354
1355
pgb = s[0].gb;
1356
c = 2;
1357
if(s[0].ga > pgb){
1358
pgb = s[0].ga;
1359
c = 1;
1360
}
1361
if(s[0].a >= pgb){
1362
pgb = s[0].a;
1363
c = 0;
1364
}
1365
1366
ca = c;
1367
1368
int ga = 1;
1369
int gb = 1;
1370
i = 0;
1371
j = 0;
1372
c = 1;
1373
while(trace[i][j] < 32){
1374
if(i ==0 || j == 0){
1375
path[c] |= 128;
1376
}
1377
if(i ==len_a || j == len_b){
1378
path[c] |= 64;
1379
}
1380
1381
switch(ca){
1382
case 0:
1383
if (trace[i][j] & 2){
1384
ca = 1;
1385
}else if (trace[i][j] & 4){
1386
ca = 2;
1387
}
1388
path[c] = 0;
1389
i++;
1390
j++;
1391
break;
1392
case 1:
1393
if(trace[i][j] & 8){
1394
ca = 1;
1395
}else{
1396
path[c-(gb-1)] |= gb << 16;
1397
gb = 0;
1398
ca = 0;
1399
}
1400
path[c] |= 1;
1401
j++;
1402
gb++;
1403
break;
1404
case 2:
1405
if(trace[i][j] & 16){
1406
ca = 2;
1407
}else{
1408
path[c-(ga-1)] |= ga << 16;
1409
ga = 0;
1410
ca = 0;
1411
}
1412
path[c] |= 2;
1413
i++;
1414
ga++;
1415
break;
1416
}
1417
c++;
1418
}
1419
if (ca == 1){
1420
path[c-(gb-1)] |= (gb-1) << 16;
1421
}
1422
if(ca == 2){
1423
path[c-(ga-1)] |= (ga-1) << 16;
1424
}
1425
path[0] = c-1;
1426
path[c] = 3;
1427
path[c+1] = pgb;
1428
return path;
1429
}
1430
1431
int* ps_dyn2(int* path, struct dp_matrix *dp,const int* prof1,const int* seq2,const int len_a,const int len_b,int sip)
1432
{
1433
1434
struct states* s = 0;
1435
char** trace = 0;
1436
char* tracep = 0;
1437
register int pa = 0;
1438
register int pga = 0;
1439
register int pgb = 0;
1440
register int ca = 0;
1441
register int i = 0;
1442
register int j = 0;
1443
register int c = 0;
1444
1445
const int open = gpo * sip;
1446
const int ext = gpe *sip;
1447
1448
1449
1450
s = dp->s;
1451
1452
trace = dp->tb;
1453
1454
trace[len_a][len_b] = 32;
1455
1456
prof1 += len_a << 6;
1457
1458
s[len_b].a = 0;
1459
s[len_b].ga = -INFTY;
1460
s[len_b].gb = -INFTY;
1461
tracep = trace[len_a];
1462
j = len_b;
1463
1464
1465
while(--j){
1466
s[j].a = -INFTY;
1467
s[j].ga = s[j+1].a-tgpe;
1468
if (s[j+1].ga-tgpe > s[j].ga){
1469
s[j].ga = s[j+1].ga-tgpe;
1470
}
1471
1472
1473
1474
s[j].gb = -INFTY;
1475
tracep[j] = 8;
1476
}
1477
1478
s[0].a = -INFTY;
1479
s[0].ga = -INFTY;
1480
s[0].gb = -INFTY;
1481
i = len_a;
1482
while(--i){
1483
prof1 -= 64;
1484
1485
tracep = trace[i];
1486
pa = s[len_b].a;
1487
pga = s[len_b].ga;
1488
pgb = s[len_b].gb;
1489
s[len_b].a = -INFTY;
1490
s[len_b].ga = -INFTY;
1491
1492
s[len_b].gb = pa+prof1[29];
1493
if(pgb+prof1[29] > s[len_b].gb){
1494
s[len_b].gb = pgb+prof1[29];
1495
}
1496
1497
1498
1499
1500
tracep[len_b] = 16;
1501
1502
j = len_b;
1503
1504
while(--j){
1505
1506
ca = s[j].a;
1507
1508
c = 1;
1509
if((pga -= open) > pa){
1510
pa = pga;
1511
c = 2;
1512
}
1513
if((pgb += prof1[91]) > pa){
1514
pa = pgb;
1515
c = 4;
1516
}
1517
1518
pa += prof1[32 + seq2[j]];
1519
1520
s[j].a = pa;
1521
1522
pga = s[j].ga;
1523
1524
s[j].ga = s[j+1].a-open;
1525
if (s[j+1].ga-ext > s[j].ga){
1526
s[j].ga = s[j+1].ga-ext;
1527
c |= 8;
1528
}
1529
1530
pgb = s[j].gb;
1531
1532
s[j].gb = ca+prof1[27];
1533
if(pgb+prof1[28] > s[j].gb){
1534
s[j].gb = pgb+prof1[28];
1535
c |= 16;
1536
}
1537
tracep[j] = c;
1538
pa = ca;
1539
1540
}
1541
1542
ca = s[0].a;
1543
1544
c = 1;
1545
if((pga-=open) > pa){
1546
pa = pga;
1547
c = 2;
1548
}
1549
if((pgb+=prof1[91]) > pa){
1550
pa = pgb;
1551
c = 4;
1552
}
1553
pa += prof1[32+seq2[0]];
1554
s[0].a = pa;
1555
1556
s[0].ga = -INFTY;
1557
1558
pgb = s[0].gb;
1559
s[0].gb = ca+prof1[27]+prof1[29];
1560
if(pgb+prof1[29] > s[0].gb){
1561
s[0].gb = pgb+prof1[29];
1562
c |= 16;
1563
}
1564
tracep[0] = c;
1565
1566
}
1567
prof1 -= 64;
1568
1569
1570
1571
tracep = trace[0];
1572
j = len_b;
1573
pa = s[j].a;
1574
pga = s[j].ga;
1575
pgb = s[j].gb;
1576
s[j].a = -INFTY;
1577
s[j].ga = -INFTY;
1578
1579
s[len_b].gb = pa+prof1[29];
1580
if(pgb+prof1[29] > s[len_b].gb){
1581
s[len_b].gb = pgb+prof1[29];
1582
}
1583
1584
1585
while(--j){
1586
1587
ca = s[j].a;
1588
1589
c = 1;
1590
1591
if((pga-=open) > pa){
1592
pa = pga;
1593
c = 2;
1594
}
1595
1596
if((pgb+=prof1[91]) > pa){
1597
pa = pgb;
1598
c = 4;
1599
}
1600
pa += prof1[32+seq2[j]];
1601
s[j].a = pa;
1602
pga = s[j].ga;
1603
s[j].ga = s[j+1].a-(open+tgpe);
1604
if (s[j+1].ga-tgpe > s[j].ga){
1605
s[j].ga = s[j+1].ga-tgpe;
1606
c |= 8;
1607
}
1608
pgb = s[j].gb;
1609
s[j].gb = -INFTY;
1610
1611
tracep[j] = c;
1612
pa = ca;
1613
}
1614
1615
1616
ca = s[0].a;
1617
1618
c = 1;
1619
1620
if((pga-=open) > pa){
1621
pa = pga;
1622
c = 2;
1623
}
1624
if((pgb+=prof1[91]) > pa){
1625
pa = pgb;
1626
c = 4;
1627
}
1628
pa += prof1[32+seq2[0]];
1629
s[0].a = pa;
1630
1631
s[0].ga = s[1].a-(open+tgpe);
1632
if (s[1].ga-tgpe > s[0].ga){
1633
s[0].ga = s[1].ga-tgpe;
1634
c |= 8;
1635
}
1636
1637
pgb = s[0].gb;
1638
s[0].gb = ca+prof1[27]+prof1[29];
1639
if(pgb+prof1[29] > s[0].gb){
1640
s[0].gb = pgb+prof1[29];
1641
c |= 16;
1642
}
1643
tracep[0] = c;
1644
1645
1646
pgb = s[0].gb;
1647
c = 2;
1648
if(s[0].ga > pgb){
1649
pgb = s[0].ga;
1650
c = 1;
1651
}
1652
if(s[0].a >= pgb){
1653
pgb = s[0].a;
1654
c = 0;
1655
}
1656
1657
ca = c;
1658
int ga = 1;
1659
int gb = 1;
1660
i = 0;
1661
j = 0;
1662
c = 1;
1663
while(trace[i][j] < 32){
1664
if(i ==0 || j == 0){
1665
path[c] |= 128;
1666
}
1667
if(i ==len_a || j == len_b){
1668
path[c] |= 64;
1669
}
1670
1671
switch(ca){
1672
case 0:
1673
if (trace[i][j] & 2){
1674
ca = 1;
1675
}else if (trace[i][j] & 4){
1676
ca = 2;
1677
}
1678
path[c] = 0;
1679
i++;
1680
j++;
1681
break;
1682
case 1:
1683
if(trace[i][j] & 8){
1684
ca = 1;
1685
}else{
1686
path[c-(gb-1)] |= gb << 16;
1687
gb = 0;
1688
ca = 0;
1689
}
1690
path[c] |= 1;
1691
j++;
1692
gb++;
1693
break;
1694
case 2:
1695
if(trace[i][j] & 16){
1696
ca = 2;
1697
}else{
1698
path[c-(ga-1)] |= ga << 16;
1699
ga = 0;
1700
ca = 0;
1701
}
1702
path[c] |= 2;
1703
i++;
1704
ga++;
1705
break;
1706
}
1707
c++;
1708
}
1709
if (ca == 1){
1710
path[c-(gb-1)] |= (gb-1) << 16;
1711
}
1712
if(ca == 2){
1713
path[c-(ga-1)] |= (ga-1) << 16;
1714
}
1715
path[0] = c-1;
1716
path[c] = 3;
1717
path[c+1] = pgb;
1718
return path;
1719
}
1720
1721
int* ss_dyn2(int**subm,int* path, struct dp_matrix *dp,const int* seq1,const int* seq2,const int len_a,const int len_b)
1722
{
1723
struct states* s = 0;
1724
int *subp = 0;
1725
char** trace = 0;
1726
char* tracep = 0;
1727
register int pa = 0;
1728
register int pga = 0;
1729
register int pgb = 0;
1730
register int ca = 0;
1731
register int i = 0;
1732
register int j = 0;
1733
register int c = 0;
1734
1735
s = dp->s;
1736
1737
trace = dp->tb;
1738
1739
trace[len_a][len_b] = 32;
1740
1741
s[len_b].a = 0;
1742
s[len_b].ga = -INFTY;
1743
s[len_b].gb = -INFTY;
1744
1745
1746
tracep = trace[len_a];
1747
j = len_b;
1748
1749
1750
while(--j){
1751
s[j].a = -INFTY;
1752
1753
s[j].ga = s[j+1].a-tgpe;
1754
if (s[j+1].ga-tgpe > s[j].ga){
1755
s[j].ga = s[j+1].ga-tgpe;
1756
}
1757
1758
s[j].gb = -INFTY;
1759
tracep[j] = 8;
1760
}
1761
1762
s[0].a = -INFTY;
1763
s[0].ga = -INFTY;
1764
s[0].gb = -INFTY;
1765
1766
i = len_a;
1767
while(--i){
1768
1769
tracep = trace[i];
1770
pa = s[len_b].a;
1771
pga = s[len_b].ga;
1772
pgb = s[len_b].gb;
1773
1774
s[len_b].a = -INFTY;
1775
s[len_b].ga = -INFTY;
1776
1777
s[len_b].gb = pa-tgpe;
1778
if(pgb-tgpe > s[len_b].gb){
1779
s[len_b].gb = pgb-tgpe;
1780
}
1781
1782
1783
tracep[len_b] = 16;
1784
j = len_b;
1785
subp = subm[seq1[i]];
1786
while(--j){
1787
ca = s[j].a;
1788
1789
c = 1;
1790
if((pga -= gpo) > pa){
1791
pa = pga;
1792
c = 2;
1793
}
1794
if((pgb -= gpo) > pa){
1795
pa = pgb;
1796
c = 4;
1797
}
1798
1799
pa += subp[seq2[j]];
1800
1801
s[j].a = pa;
1802
1803
pga = s[j].ga;
1804
1805
s[j].ga = s[j+1].a-gpo;
1806
if (s[j+1].ga-gpe > s[j].ga){
1807
s[j].ga = s[j+1].ga-gpe;
1808
c |= 8;
1809
}
1810
1811
pgb = s[j].gb;
1812
1813
s[j].gb = ca-gpo;
1814
if(pgb-gpe > s[j].gb){
1815
s[j].gb = pgb-gpe;
1816
c |= 16;
1817
}
1818
tracep[j] = c;
1819
pa = ca;
1820
1821
}
1822
1823
ca = s[0].a;
1824
1825
c = 1;
1826
if((pga-=gpo) > pa){
1827
pa = pga;
1828
c = 2;
1829
}
1830
if((pgb-=gpo) > pa){
1831
pa = pgb;
1832
c = 4;
1833
}
1834
1835
pa += subp[seq2[0]];
1836
1837
s[0].a = pa;
1838
1839
s[0].ga = -INFTY;
1840
1841
pgb = s[0].gb;
1842
s[0].gb = ca-(gpo+tgpe);
1843
if(pgb-tgpe > s[0].gb){
1844
s[0].gb = pgb-tgpe;
1845
c |= 16;
1846
}
1847
tracep[0] = c;
1848
}
1849
1850
subp = subm[seq1[0]];
1851
tracep = trace[0];
1852
j = len_b;
1853
pa = s[j].a;
1854
pga = s[j].ga;
1855
pgb = s[j].gb;
1856
s[j].a = -INFTY;
1857
s[j].ga = -INFTY;
1858
1859
s[j].gb = pa-tgpe;
1860
if(pgb-tgpe > s[j].gb){
1861
s[j].gb = pgb-tgpe;
1862
}
1863
1864
while(--j){
1865
1866
ca = s[j].a;
1867
1868
c = 1;
1869
1870
if((pga-=gpo) > pa){
1871
pa = pga;
1872
c = 2;
1873
}
1874
1875
if((pgb-=gpo) > pa){
1876
pa = pgb;
1877
c = 4;
1878
}
1879
1880
pa += subp[seq2[j]];
1881
1882
s[j].a = pa;
1883
1884
pga = s[j].ga;
1885
s[j].ga = s[j+1].a-(gpo+tgpe);
1886
if (s[j+1].ga-tgpe > s[j].ga){
1887
s[j].ga = s[j+1].ga-tgpe;
1888
c |= 8;
1889
}
1890
pgb = s[j].gb;
1891
s[j].gb = -INFTY;
1892
tracep[j] = c;
1893
pa = ca;
1894
}
1895
1896
ca = s[0].a;
1897
1898
c = 1;
1899
1900
if((pga-=gpo) > pa){
1901
pa = pga;
1902
c = 2;
1903
}
1904
if((pgb-=gpo) > pa){
1905
pa = pgb;
1906
c = 4;
1907
}
1908
1909
pa += subp[seq2[0]];
1910
1911
s[0].a = pa;
1912
1913
1914
s[0].ga = s[1].a-(gpo+tgpe);
1915
if (s[1].ga-tgpe > s[0].ga){
1916
s[0].ga = s[1].ga-tgpe;
1917
c |= 8;
1918
}
1919
1920
pgb = s[0].gb;
1921
s[0].gb = ca-(gpo+tgpe);
1922
if(pgb-tgpe > s[0].gb){
1923
s[0].gb = pgb-tgpe;
1924
c |= 16;
1925
}
1926
tracep[0] = c;
1927
1928
1929
pgb = s[0].gb;
1930
c = 2;
1931
if(s[0].ga > pgb){
1932
pgb = s[0].ga;
1933
c = 1;
1934
}
1935
if(s[0].a >= pgb){
1936
pgb = s[0].a;
1937
c = 0;
1938
}
1939
1940
ca = c;
1941
1942
int ga = 1;
1943
int gb = 1;
1944
i = 0;
1945
j = 0;
1946
c = 1;
1947
while(trace[i][j] < 32){
1948
if(i ==0 || j == 0){
1949
path[c] |= 128;
1950
}
1951
if(i ==len_a || j == len_b){
1952
path[c] |= 64;
1953
}
1954
1955
switch(ca){
1956
case 0:
1957
if (trace[i][j] & 2){
1958
ca = 1;
1959
}else if (trace[i][j] & 4){
1960
ca = 2;
1961
}
1962
path[c] = 0;
1963
i++;
1964
j++;
1965
break;
1966
case 1:
1967
if(trace[i][j] & 8){
1968
ca = 1;
1969
}else{
1970
path[c-(gb-1)] |= gb << 16;
1971
gb = 0;
1972
ca = 0;
1973
}
1974
path[c] |= 1;
1975
j++;
1976
gb++;
1977
break;
1978
case 2:
1979
if(trace[i][j] & 16){
1980
ca = 2;
1981
}else{
1982
path[c-(ga-1)] |= ga << 16;
1983
ga = 0;
1984
ca = 0;
1985
}
1986
path[c] |= 2;
1987
i++;
1988
ga++;
1989
break;
1990
}
1991
c++;
1992
}
1993
if (ca == 1){
1994
path[c-(gb-1)] |= (gb-1) << 16;
1995
}
1996
if(ca == 2){
1997
path[c-(ga-1)] |= (ga-1) << 16;
1998
}
1999
path[0] = c-1;
2000
path[c] = 3;
2001
path[c+1] = pgb;
2002
return path;
2003
}
2004
2005
2006
2007
int* aapp_dyn(int* path, struct dp_matrix *dp,const int* prof1,const int* prof2,const int len_a,const int len_b,const int mmbonus)
2008
{
2009
unsigned int freq[26];
2010
2011
struct states* s = 0;
2012
char** trace = 0;
2013
char* tracep = 0;
2014
register int pa = 0;
2015
register int pga = 0;
2016
register int pgb = 0;
2017
register int ca = 0;
2018
register int i = 0;
2019
register int j = 0;
2020
register int c = 0;
2021
2022
s = dp->s;
2023
2024
trace = dp->tb;
2025
2026
trace[len_a][len_b] = 32;
2027
2028
prof1 += len_a << 6;
2029
2030
s[len_b].a = 0;
2031
s[len_b].ga = -INFTY;
2032
s[len_b].gb = -INFTY;
2033
//init of first row;
2034
tracep = trace[len_a];
2035
2036
j = len_b;
2037
while(--j){
2038
s[j].a = -INFTY;
2039
2040
s[j].ga = s[j+1].a+prof2[29];
2041
if (s[j+1].ga+prof2[29] > s[j].ga){
2042
s[j].ga = s[j+1].ga+prof2[29];
2043
}
2044
s[j].gb = -INFTY;
2045
tracep[j] = 8;
2046
}
2047
2048
s[0].a = -INFTY;
2049
s[0].ga = -INFTY;
2050
s[0].gb = -INFTY;
2051
i = len_a;
2052
while(--i){
2053
prof1 -= 64;
2054
2055
c = 1;
2056
for (j = 26; j--;){
2057
if(prof1[j]){
2058
freq[c] = j;
2059
c++;
2060
}
2061
}
2062
freq[0] = c;
2063
2064
tracep = trace[i];
2065
pa = s[len_b].a + mmbonus;
2066
pga = s[len_b].ga;
2067
pgb = s[len_b].gb;
2068
s[len_b].a = -INFTY;
2069
s[len_b].ga = -INFTY;
2070
2071
s[len_b].gb = pa+prof1[29];
2072
if(pgb+prof1[29] > s[len_b].gb){
2073
s[len_b].gb = pgb+prof1[29];
2074
}
2075
2076
tracep[len_b] = 16;
2077
2078
j = len_b;
2079
prof2 += len_b << 6;
2080
while(--j){
2081
prof2 -= 64;
2082
ca = s[j].a;
2083
2084
c = 1;
2085
if((pga += prof2[91]) > pa){
2086
pa = pga;
2087
c = 2;
2088
}
2089
if((pgb += prof1[91]) > pa){
2090
pa = pgb;
2091
c = 4;
2092
}
2093
2094
prof2 += 32;
2095
for (pga = freq[0];--pga;){
2096
pgb = freq[pga];
2097
pa += prof1[pgb]*prof2[pgb];
2098
}
2099
prof2 -= 32;
2100
2101
s[j].a = pa;
2102
2103
pga = s[j].ga;
2104
2105
s[j].ga = s[j+1].a+prof2[27];
2106
if (s[j+1].ga+prof2[28] > s[j].ga){
2107
s[j].ga = s[j+1].ga+prof2[28];
2108
c |= 8;
2109
}
2110
2111
pgb = s[j].gb;
2112
2113
s[j].gb = ca+prof1[27];
2114
if(pgb+prof1[28] > s[j].gb){
2115
s[j].gb = pgb+prof1[28];
2116
c |= 16;
2117
}
2118
tracep[j] = c;
2119
pa = ca+ mmbonus;
2120
2121
}
2122
2123
prof2 -= 64;
2124
//LAST CELL (0)
2125
ca = s[0].a;
2126
2127
c = 1;
2128
if((pga+=prof2[91]) > pa){
2129
pa = pga;
2130
c = 2;
2131
}
2132
if((pgb+=prof1[91]) > pa){
2133
pa = pgb;
2134
c = 4;
2135
}
2136
2137
prof2 += 32;
2138
for (pga = freq[0];--pga;){
2139
pgb = freq[pga];
2140
pa += prof1[pgb]*prof2[pgb];
2141
}
2142
prof2 -= 32;
2143
2144
s[0].a = pa;
2145
2146
s[0].ga = -INFTY;
2147
2148
pgb = s[0].gb;
2149
s[0].gb = ca+prof1[27]+prof1[29];
2150
if(pgb+prof1[29] > s[0].gb){
2151
s[0].gb = pgb+prof1[29];
2152
c |= 16;
2153
}
2154
tracep[0] = c;
2155
2156
}
2157
prof1 -= 64;
2158
2159
c = 1;
2160
for (j = 26; j--;){
2161
if(prof1[j]){
2162
freq[c] = j;
2163
c++;
2164
}
2165
}
2166
freq[0] = c;
2167
2168
tracep = trace[0];
2169
j = len_b;
2170
prof2 += len_b << 6;
2171
pa = s[j].a+ mmbonus;
2172
pga = s[j].ga;
2173
pgb = s[j].gb;
2174
s[j].a = -INFTY;
2175
s[j].ga = -INFTY;
2176
2177
s[len_b].gb = pa+prof1[29];
2178
if(pgb+prof1[29] > s[len_b].gb){
2179
s[len_b].gb = pgb+prof1[29];
2180
}
2181
2182
2183
2184
while(--j){
2185
prof2 -= 64;
2186
ca = s[j].a;
2187
2188
c = 1;
2189
2190
if((pga+=prof2[91]) > pa){
2191
pa = pga;
2192
c = 2;
2193
}
2194
2195
if((pgb+=prof1[91]) > pa){
2196
pa = pgb;
2197
c = 4;
2198
}
2199
prof2+=32;
2200
2201
for (pga = freq[0];--pga;){
2202
pgb = freq[pga];
2203
pa += prof1[pgb]*prof2[pgb];
2204
}
2205
prof2-=32;
2206
2207
s[j].a = pa;
2208
pga = s[j].ga;
2209
s[j].ga = s[j+1].a+prof2[27]+prof2[29];
2210
if (s[j+1].ga+prof2[29] > s[j].ga){
2211
s[j].ga = s[j+1].ga+prof2[29];
2212
c |= 8;
2213
}
2214
pgb = s[j].gb;
2215
s[j].gb = -INFTY;
2216
2217
tracep[j] = c;
2218
pa = ca+ mmbonus;
2219
}
2220
prof2 -= 64;
2221
2222
ca = s[0].a;
2223
2224
c = 1;
2225
2226
if((pga+=prof2[91]) > pa){
2227
pa = pga;
2228
c = 2;
2229
}
2230
if((pgb+=prof1[91]) > pa){
2231
pa = pgb;
2232
c = 4;
2233
}
2234
2235
prof2+=32;
2236
for (pga = freq[0];--pga;){
2237
pgb = freq[pga];
2238
pa += prof1[pgb]*prof2[pgb];
2239
}
2240
prof2-=32;
2241
2242
s[0].a = pa;
2243
2244
s[0].ga = s[1].a+prof2[27]+prof2[29];
2245
if (s[1].ga+prof2[29] > s[0].ga){
2246
s[0].ga = s[1].ga+prof2[29];
2247
c |= 8;
2248
}
2249
2250
pgb = s[0].gb;
2251
s[0].gb = ca+prof1[27]+prof1[29];
2252
if(pgb +prof1[29]> s[0].gb){
2253
s[0].gb = pgb+prof1[29];
2254
c |= 16;
2255
}
2256
tracep[0] = c;
2257
2258
pgb = s[0].gb;
2259
c = 2;
2260
if(s[0].ga > pgb){
2261
pgb = s[0].ga;
2262
c = 1;
2263
}
2264
if(s[0].a >= pgb){
2265
pgb = s[0].a;
2266
c = 0;
2267
}
2268
2269
//fprintf(stderr,"SCORE:%d\n",ca);
2270
ca = c;
2271
2272
i = 0;
2273
j = 0;
2274
c = 1;
2275
while(trace[i][j] < 32){
2276
// fprintf(stderr,"%d->%d %d:%d %d:%d\n",c,trace[i][j],i,j,len_a,len_b);
2277
switch(ca){
2278
case 0:
2279
if (trace[i][j] & 2){
2280
ca = 1;
2281
if(i+1!= len_a){
2282
path[c+1] |= 16;
2283
// fprintf(stderr,"GAP_CLOSE\n");
2284
}else{
2285
path[c+1] |= 32+16;
2286
}
2287
}else if (trace[i][j] & 4){
2288
ca = 2;
2289
if(j+1!= len_b){
2290
path[c+1] |= 16;
2291
// fprintf(stderr,"GAP_CLOSE\n");
2292
}else{
2293
path[c+1] |= 32+16;
2294
}
2295
}
2296
2297
//path[c] = 0;
2298
i++;
2299
j++;
2300
break;
2301
case 1:
2302
if(trace[i][j] & 8){
2303
ca = 1;
2304
if(i!=0 && i!= len_a){
2305
// / fprintf(stderr,"GAP_EXT\n");
2306
if(!(path[c]&16)){
2307
path[c] |= 8;
2308
}
2309
}else{
2310
if(!(path[c]&16)){
2311
path[c] |= 32+8;
2312
}
2313
}
2314
}else{
2315
ca = 0;
2316
if(i!=0 && i!= len_a){
2317
// fprintf(stderr,"GAP_OPEN\n");
2318
path[c] |= 4;
2319
}else{
2320
path[c] |= 32+4;
2321
}
2322
}
2323
path[c] |= 1;
2324
j++;
2325
break;
2326
case 2:
2327
if(trace[i][j] & 16){
2328
ca = 2;
2329
if(j !=0 && j != len_b){
2330
// fprintf(stderr,"GAP_EXT\n");
2331
if(!(path[c]&16)){
2332
path[c] |= 8;
2333
}
2334
}else{
2335
if(!(path[c]&16)){
2336
path[c] |= 32+8;
2337
}
2338
}
2339
}else{
2340
ca = 0;
2341
if(j!=0 && j != len_b){
2342
// fprintf(stderr,"GAP_OPEN\n");
2343
path[c] |= 4;
2344
}else{
2345
path[c] |= 32+4;
2346
}
2347
2348
}
2349
path[c] |= 2;
2350
i++;
2351
break;
2352
}
2353
c++;
2354
}
2355
path[0] = c-1;
2356
path[c] = 3;
2357
path[c+1] = pgb;
2358
2359
return path;
2360
}*/
2361
2362
2363
2364
2365
2366
int* pp_dyn(int* path, struct dp_matrix *dp,const float* prof1,const float* prof2,const int len_a,const int len_b)
2367
{
2368
unsigned int freq[26];
2369
2370
struct states* s = 0;
2371
char** trace = 0;
2372
char* tracep = 0;
2373
register float pa = 0;
2374
register float pga = 0;
2375
register float pgb = 0;
2376
register float ca = 0;
2377
register int i = 0;
2378
register int j = 0;
2379
register int c = 0;
2380
register int f = 0;
2381
2382
s = dp->s;
2383
2384
trace = dp->tb;
2385
2386
trace[len_a][len_b] = 32;
2387
2388
prof1 += len_a << 6;
2389
2390
s[len_b].a = 0.0;
2391
s[len_b].ga = -FLOATINFTY;
2392
s[len_b].gb = -FLOATINFTY;
2393
//init of first row;
2394
tracep = trace[len_a];
2395
2396
j = len_b;
2397
while(--j){
2398
s[j].a = -FLOATINFTY;
2399
2400
s[j].ga = s[j+1].a+prof2[29];
2401
if (s[j+1].ga+prof2[29] > s[j].ga){
2402
s[j].ga = s[j+1].ga+prof2[29];
2403
}
2404
s[j].gb = -INFTY;
2405
tracep[j] = 8;
2406
}
2407
2408
s[0].a = -FLOATINFTY;
2409
s[0].ga = -FLOATINFTY;
2410
s[0].gb = -FLOATINFTY;
2411
i = len_a;
2412
while(--i){
2413
prof1 -= 64;
2414
2415
c = 1;
2416
for (j = 26; j--;){
2417
if(prof1[j]){
2418
freq[c] = j;
2419
c++;
2420
}
2421
}
2422
freq[0] = c;
2423
2424
tracep = trace[i];
2425
pa = s[len_b].a;
2426
pga = s[len_b].ga;
2427
pgb = s[len_b].gb;
2428
s[len_b].a = -FLOATINFTY;
2429
s[len_b].ga = -FLOATINFTY;
2430
2431
s[len_b].gb = pa+prof1[29];
2432
if(pgb+prof1[29] > s[len_b].gb){
2433
s[len_b].gb = pgb+prof1[29];
2434
}
2435
2436
tracep[len_b] = 16;
2437
2438
j = len_b;
2439
prof2 += len_b << 6;
2440
while(--j){
2441
prof2 -= 64;
2442
ca = s[j].a;
2443
2444
c = 1;
2445
if((pga += prof2[91]) > pa){
2446
pa = pga;
2447
c = 2;
2448
}
2449
if((pgb += prof1[91]) > pa){
2450
pa = pgb;
2451
c = 4;
2452
}
2453
2454
prof2 += 32;
2455
for (f = freq[0];--f;){
2456
pa += prof1[freq[f]]*prof2[freq[f]];
2457
}
2458
prof2 -= 32;
2459
2460
s[j].a = pa;
2461
2462
pga = s[j].ga;
2463
2464
s[j].ga = s[j+1].a+prof2[27];
2465
if (s[j+1].ga+prof2[28] > s[j].ga){
2466
s[j].ga = s[j+1].ga+prof2[28];
2467
c |= 8;
2468
}
2469
2470
pgb = s[j].gb;
2471
2472
s[j].gb = ca+prof1[27];
2473
if(pgb+prof1[28] > s[j].gb){
2474
s[j].gb = pgb+prof1[28];
2475
c |= 16;
2476
}
2477
tracep[j] = c;
2478
pa = ca;
2479
2480
}
2481
2482
prof2 -= 64;
2483
//LAST CELL (0)
2484
ca = s[0].a;
2485
2486
c = 1;
2487
if((pga+=prof2[91]) > pa){
2488
pa = pga;
2489
c = 2;
2490
}
2491
if((pgb+=prof1[91]) > pa){
2492
pa = pgb;
2493
c = 4;
2494
}
2495
2496
prof2 += 32;
2497
for (f = freq[0];--f;){
2498
pa += prof1[freq[f]]*prof2[freq[f]];
2499
}
2500
prof2 -= 32;
2501
2502
s[0].a = pa;
2503
2504
s[0].ga = -FLOATINFTY;
2505
2506
pgb = s[0].gb;
2507
s[0].gb = ca+prof1[27]+prof1[29];
2508
if(pgb+prof1[29] > s[0].gb){
2509
s[0].gb = pgb+prof1[29];
2510
c |= 16;
2511
}
2512
tracep[0] = c;
2513
2514
}
2515
prof1 -= 64;
2516
2517
c = 1;
2518
for (j = 26; j--;){
2519
if(prof1[j]){
2520
freq[c] = j;
2521
c++;
2522
}
2523
}
2524
freq[0] = c;
2525
2526
tracep = trace[0];
2527
j = len_b;
2528
prof2 += len_b << 6;
2529
pa = s[j].a;
2530
pga = s[j].ga;
2531
pgb = s[j].gb;
2532
s[j].a = -FLOATINFTY;
2533
s[j].ga = -FLOATINFTY;
2534
2535
s[len_b].gb = pa+prof1[29];
2536
if(pgb+prof1[29] > s[len_b].gb){
2537
s[len_b].gb = pgb+prof1[29];
2538
}
2539
2540
2541
2542
while(--j){
2543
prof2 -= 64;
2544
ca = s[j].a;
2545
2546
c = 1;
2547
2548
if((pga+=prof2[91]) > pa){
2549
pa = pga;
2550
c = 2;
2551
}
2552
2553
if((pgb+=prof1[91]) > pa){
2554
pa = pgb;
2555
c = 4;
2556
}
2557
2558
prof2+=32;
2559
2560
for (f = freq[0];--f;){
2561
pa += prof1[freq[f]]*prof2[freq[f]];
2562
}
2563
prof2-=32;
2564
2565
s[j].a = pa;
2566
pga = s[j].ga;
2567
s[j].ga = s[j+1].a+prof2[27]+prof2[29];
2568
if (s[j+1].ga+prof2[29] > s[j].ga){
2569
s[j].ga = s[j+1].ga+prof2[29];
2570
c |= 8;
2571
}
2572
pgb = s[j].gb;
2573
s[j].gb = -FLOATINFTY;
2574
2575
tracep[j] = c;
2576
pa = ca;
2577
}
2578
prof2 -= 64;
2579
2580
ca = s[0].a;
2581
2582
c = 1;
2583
2584
if((pga+=prof2[91]) > pa){
2585
pa = pga;
2586
c = 2;
2587
}
2588
if((pgb+=prof1[91]) > pa){
2589
pa = pgb;
2590
c = 4;
2591
}
2592
prof2+=32;
2593
for (f = freq[0];--f;){
2594
pa += prof1[freq[f]]*prof2[freq[f]];
2595
}
2596
prof2-=32;
2597
2598
s[0].a = pa;
2599
2600
s[0].ga = s[1].a+prof2[27]+prof2[29];
2601
if (s[1].ga+prof2[29] > s[0].ga){
2602
s[0].ga = s[1].ga+prof2[29];
2603
c |= 8;
2604
}
2605
2606
pgb = s[0].gb;
2607
s[0].gb = ca+prof1[27]+prof1[29];
2608
if(pgb +prof1[29]> s[0].gb){
2609
s[0].gb = pgb+prof1[29];
2610
c |= 16;
2611
}
2612
tracep[0] = c;
2613
2614
pgb = s[0].gb;
2615
c = 2;
2616
if(s[0].ga > pgb){
2617
pgb = s[0].ga;
2618
c = 1;
2619
}
2620
if(s[0].a >= pgb){
2621
pgb = s[0].a;
2622
c = 0;
2623
}
2624
2625
//fprintf(stderr,"SCORE:%d\n",ca);
2626
f = c;
2627
2628
i = 0;
2629
j = 0;
2630
c = 1;
2631
while(trace[i][j] < 32){
2632
// fprintf(stderr,"%d->%d %d:%d %d:%d\n",c,trace[i][j],i,j,len_a,len_b);
2633
switch(f){
2634
case 0:
2635
if (trace[i][j] & 2){
2636
f = 1;
2637
if(i+1!= len_a){
2638
path[c+1] |= 16;
2639
// fprintf(stderr,"GAP_CLOSE\n");
2640
}else{
2641
path[c+1] |= 32+16;
2642
}
2643
}else if (trace[i][j] & 4){
2644
f = 2;
2645
if(j+1!= len_b){
2646
path[c+1] |= 16;
2647
// fprintf(stderr,"GAP_CLOSE\n");
2648
}else{
2649
path[c+1] |= 32+16;
2650
}
2651
}
2652
2653
//path[c] = 0;
2654
i++;
2655
j++;
2656
break;
2657
case 1:
2658
if(trace[i][j] & 8){
2659
f = 1;
2660
if(i!=0 && i!= len_a){
2661
// / fprintf(stderr,"GAP_EXT\n");
2662
if(!(path[c]&16)){
2663
path[c] |= 8;
2664
}
2665
}else{
2666
if(!(path[c]&16)){
2667
path[c] |= 32+8;
2668
}
2669
}
2670
}else{
2671
f = 0;
2672
if(i!=0 && i!= len_a){
2673
// fprintf(stderr,"GAP_OPEN\n");
2674
path[c] |= 4;
2675
}else{
2676
path[c] |= 32+4;
2677
}
2678
}
2679
path[c] |= 1;
2680
j++;
2681
break;
2682
case 2:
2683
if(trace[i][j] & 16){
2684
f = 2;
2685
if(j !=0 && j != len_b){
2686
// fprintf(stderr,"GAP_EXT\n");
2687
if(!(path[c]&16)){
2688
path[c] |= 8;
2689
}
2690
}else{
2691
if(!(path[c]&16)){
2692
path[c] |= 32+8;
2693
}
2694
}
2695
}else{
2696
f = 0;
2697
if(j!=0 && j != len_b){
2698
// fprintf(stderr,"GAP_OPEN\n");
2699
path[c] |= 4;
2700
}else{
2701
path[c] |= 32+4;
2702
}
2703
2704
}
2705
path[c] |= 2;
2706
i++;
2707
break;
2708
}
2709
c++;
2710
}
2711
path[0] = c-1;
2712
path[c] = 3;
2713
path[c+1] = pgb;
2714
return path;
2715
}
2716
2717
2718
int* ps_dyn(int* path, struct dp_matrix *dp,const float* prof1,const int* seq2,const int len_a,const int len_b,int sip)
2719
{
2720
struct states* s = 0;
2721
char** trace = 0;
2722
char* tracep = 0;
2723
register float pa = 0;
2724
register float pga = 0;
2725
register float pgb = 0;
2726
register float ca = 0;
2727
register int i = 0;
2728
register int j = 0;
2729
register int c = 0;
2730
register int f = 0;
2731
2732
const float open = gpo * sip;
2733
const float ext = gpe *sip;
2734
2735
s = dp->s;
2736
2737
trace = dp->tb;
2738
2739
trace[len_a][len_b] = 32;
2740
2741
prof1 += len_a << 6;
2742
2743
s[len_b].a = 0.0;
2744
s[len_b].ga = -FLOATINFTY;
2745
s[len_b].gb = -FLOATINFTY;
2746
//init of first row;
2747
tracep = trace[len_a];
2748
j = len_b;
2749
2750
2751
while(--j){
2752
s[j].a = -FLOATINFTY;
2753
//s[j].ga = 0;
2754
2755
s[j].ga = s[j+1].a-tgpe;//-topen;
2756
if (s[j+1].ga-tgpe > s[j].ga){
2757
s[j].ga = s[j+1].ga-tgpe;
2758
}
2759
2760
s[j].gb = -FLOATINFTY;
2761
tracep[j] = 8;
2762
}
2763
2764
s[0].a = -FLOATINFTY;
2765
s[0].ga = -FLOATINFTY;
2766
s[0].gb = -FLOATINFTY;
2767
i = len_a;
2768
while(--i){
2769
prof1 -= 64;
2770
2771
tracep = trace[i];
2772
pa = s[len_b].a;
2773
pga = s[len_b].ga;
2774
pgb = s[len_b].gb;
2775
s[len_b].a = -FLOATINFTY;
2776
s[len_b].ga = -FLOATINFTY;
2777
//s[len_b].gb = 0;
2778
s[len_b].gb = pa+prof1[29];//+prof1[29];
2779
if(pgb+prof1[29] > s[len_b].gb){
2780
s[len_b].gb = pgb+prof1[29];
2781
}
2782
2783
tracep[len_b] = 16;
2784
2785
j = len_b;
2786
2787
while(--j){
2788
2789
ca = s[j].a;
2790
2791
c = 1;
2792
if((pga -= open) > pa){
2793
pa = pga;
2794
c = 2;
2795
}
2796
if((pgb += prof1[91]) > pa){
2797
pa = pgb;
2798
c = 4;
2799
}
2800
2801
pa += prof1[32 + seq2[j]];
2802
2803
s[j].a = pa;
2804
2805
pga = s[j].ga;
2806
2807
s[j].ga = s[j+1].a-open;
2808
if (s[j+1].ga-ext > s[j].ga){
2809
s[j].ga = s[j+1].ga-ext;
2810
c |= 8;
2811
}
2812
2813
pgb = s[j].gb;
2814
2815
s[j].gb = ca+prof1[27];
2816
if(pgb+prof1[28] > s[j].gb){
2817
s[j].gb = pgb+prof1[28];
2818
c |= 16;
2819
}
2820
tracep[j] = c;
2821
pa = ca;
2822
2823
}
2824
2825
//LAST CELL (0)
2826
ca = s[0].a;
2827
2828
c = 1;
2829
if((pga-=open) > pa){
2830
pa = pga;
2831
c = 2;
2832
}
2833
if((pgb+=prof1[91]) > pa){
2834
pa = pgb;
2835
c = 4;
2836
}
2837
pa += prof1[32+seq2[0]];
2838
s[0].a = pa;
2839
2840
s[0].ga = -FLOATINFTY;
2841
2842
pgb = s[0].gb;
2843
s[0].gb = ca+prof1[27]+prof1[29];
2844
if(pgb+prof1[29] > s[0].gb){
2845
s[0].gb = pgb+prof1[29];
2846
c |= 16;
2847
}
2848
tracep[0] = c;
2849
2850
}
2851
prof1 -= 64;
2852
2853
2854
2855
tracep = trace[0];
2856
j = len_b;
2857
pa = s[j].a;
2858
pga = s[j].ga;
2859
pgb = s[j].gb;
2860
s[j].a = -FLOATINFTY;
2861
s[j].ga = -FLOATINFTY;
2862
//s[j].gb = -INFTY;
2863
s[len_b].gb = pa+prof1[29];//+prof1[29];
2864
if(pgb+prof1[29] > s[len_b].gb){
2865
s[len_b].gb = pgb+prof1[29];
2866
}
2867
2868
2869
while(--j){
2870
2871
ca = s[j].a;
2872
2873
c = 1;
2874
2875
if((pga-=open) > pa){
2876
pa = pga;
2877
c = 2;
2878
}
2879
2880
if((pgb+=prof1[91]) > pa){
2881
pa = pgb;
2882
c = 4;
2883
}
2884
pa += prof1[32+seq2[j]];
2885
s[j].a = pa;
2886
pga = s[j].ga;
2887
s[j].ga = s[j+1].a-(open+tgpe);
2888
if (s[j+1].ga-tgpe > s[j].ga){
2889
s[j].ga = s[j+1].ga-tgpe;
2890
c |= 8;
2891
}
2892
pgb = s[j].gb;
2893
s[j].gb = -INFTY;
2894
2895
tracep[j] = c;
2896
pa = ca;
2897
}
2898
2899
2900
ca = s[0].a;
2901
2902
c = 1;
2903
2904
if((pga-=open) > pa){
2905
pa = pga;
2906
c = 2;
2907
}
2908
if((pgb+=prof1[91]) > pa){
2909
pa = pgb;
2910
c = 4;
2911
}
2912
pa += prof1[32+seq2[0]];
2913
s[0].a = pa;
2914
2915
s[0].ga = s[1].a-(open+tgpe);
2916
if (s[1].ga-tgpe > s[0].ga){
2917
s[0].ga = s[1].ga-tgpe;
2918
c |= 8;
2919
}
2920
2921
pgb = s[0].gb;
2922
s[0].gb = ca+prof1[27]+prof1[29];
2923
if(pgb+prof1[29] > s[0].gb){
2924
s[0].gb = pgb+prof1[29];
2925
c |= 16;
2926
}
2927
tracep[0] = c;
2928
2929
2930
pgb = s[0].gb;
2931
c = 2;
2932
if(s[0].ga > pgb){
2933
pgb = s[0].ga;
2934
c = 1;
2935
}
2936
if(s[0].a >= pgb){
2937
pgb = s[0].a;
2938
c = 0;
2939
}
2940
2941
//fprintf(stderr,"SCORE:%d\n",ca);
2942
f = c;
2943
2944
i = 0;
2945
j = 0;
2946
c = 1;
2947
while(trace[i][j] < 32){
2948
// fprintf(stderr,"%d->%d %d:%d %d:%d\n",c,trace[i][j],i,j,len_a,len_b);
2949
switch(f){
2950
case 0:
2951
if (trace[i][j] & 2){
2952
f = 1;
2953
if(i+1!= len_a){
2954
path[c+1] |= 16;
2955
// fprintf(stderr,"GAP_CLOSE\n");
2956
}else{
2957
path[c+1] |= 32+16;
2958
}
2959
}else if (trace[i][j] & 4){
2960
f = 2;
2961
if(j+1!= len_b){
2962
path[c+1] |= 16;
2963
// fprintf(stderr,"GAP_CLOSE\n");
2964
}else{
2965
path[c+1] |= 32+16;
2966
}
2967
}
2968
2969
//path[c] = 0;
2970
i++;
2971
j++;
2972
break;
2973
case 1:
2974
if(trace[i][j] & 8){
2975
f = 1;
2976
if(i!=0 && i!= len_a){
2977
// / fprintf(stderr,"GAP_EXT\n");
2978
if(!(path[c]&16)){
2979
path[c] |= 8;
2980
}
2981
}else{
2982
if(!(path[c]&16)){
2983
path[c] |= 32+8;
2984
}
2985
}
2986
}else{
2987
f = 0;
2988
if(i!=0 && i!= len_a){
2989
// fprintf(stderr,"GAP_OPEN\n");
2990
path[c] |= 4;
2991
}else{
2992
path[c] |= 32+4;
2993
}
2994
}
2995
path[c] |= 1;
2996
j++;
2997
break;
2998
case 2:
2999
if(trace[i][j] & 16){
3000
f = 2;
3001
if(j !=0 && j != len_b){
3002
// fprintf(stderr,"GAP_EXT\n");
3003
if(!(path[c]&16)){
3004
path[c] |= 8;
3005
}
3006
}else{
3007
if(!(path[c]&16)){
3008
path[c] |= 32+8;
3009
}
3010
}
3011
}else{
3012
f = 0;
3013
if(j!=0 && j != len_b){
3014
// fprintf(stderr,"GAP_OPEN\n");
3015
path[c] |= 4;
3016
}else{
3017
path[c] |= 32+4;
3018
}
3019
3020
}
3021
path[c] |= 2;
3022
i++;
3023
break;
3024
}
3025
c++;
3026
}
3027
path[0] = c-1;
3028
path[c] = 3;
3029
path[c+1] = pgb;
3030
return path;
3031
}
3032
3033
int* ss_dyn(float**subm,int* path, struct dp_matrix *dp,const int* seq1,const int* seq2,const int len_a,const int len_b)
3034
{
3035
struct states* s = 0;
3036
const float *subp = 0;
3037
char** trace = 0;
3038
char* tracep = 0;
3039
register float pa = 0;
3040
register float pga = 0;
3041
register float pgb = 0;
3042
register float ca = 0;
3043
register int i = 0;
3044
register int j = 0;
3045
register int c = 0;
3046
register int f = 0;
3047
3048
s = dp->s;
3049
3050
trace = dp->tb;
3051
3052
trace[len_a][len_b] = 32;
3053
3054
s[len_b].a = 0.0;
3055
s[len_b].ga = -FLOATINFTY;
3056
s[len_b].gb = -FLOATINFTY;
3057
3058
//init of first row;
3059
tracep = trace[len_a];
3060
j = len_b;
3061
3062
3063
while(--j){
3064
s[j].a = -FLOATINFTY;
3065
//s[j].ga = 0;
3066
s[j].ga = s[j+1].a-tgpe;//-gpo;
3067
if (s[j+1].ga-tgpe > s[j].ga){
3068
s[j].ga = s[j+1].ga-tgpe;
3069
}
3070
3071
3072
3073
s[j].gb = -FLOATINFTY;
3074
tracep[j] = 8;
3075
}
3076
3077
s[0].a = -FLOATINFTY;
3078
s[0].ga = -FLOATINFTY;
3079
s[0].gb = -FLOATINFTY;
3080
3081
i = len_a;
3082
while(--i){
3083
3084
tracep = trace[i];
3085
pa = s[len_b].a;
3086
pga = s[len_b].ga;
3087
pgb = s[len_b].gb;
3088
3089
s[len_b].a = -FLOATINFTY;
3090
s[len_b].ga = -FLOATINFTY;
3091
//s[len_b].gb = 0;
3092
3093
s[len_b].gb = pa-tgpe;//-gpo;
3094
if(pgb-tgpe > s[len_b].gb){
3095
s[len_b].gb = pgb-tgpe;
3096
}
3097
3098
3099
tracep[len_b] = 16;
3100
j = len_b;
3101
subp = subm[seq1[i]];
3102
while(--j){
3103
ca = s[j].a;
3104
3105
c = 1;
3106
if((pga -= gpo) > pa){
3107
pa = pga;
3108
c = 2;
3109
}
3110
if((pgb -= gpo) > pa){
3111
pa = pgb;
3112
c = 4;
3113
}
3114
3115
pa += subp[seq2[j]];
3116
3117
s[j].a = pa;
3118
3119
pga = s[j].ga;
3120
3121
s[j].ga = s[j+1].a-gpo;
3122
if (s[j+1].ga-gpe > s[j].ga){
3123
s[j].ga = s[j+1].ga-gpe;
3124
c |= 8;
3125
}
3126
3127
pgb = s[j].gb;
3128
3129
s[j].gb = ca-gpo;
3130
if(pgb-gpe > s[j].gb){
3131
s[j].gb = pgb-gpe;
3132
c |= 16;
3133
}
3134
tracep[j] = c;
3135
pa = ca;
3136
3137
}
3138
3139
//LAST CELL (0)
3140
ca = s[0].a;
3141
3142
c = 1;
3143
if((pga-=gpo) > pa){
3144
pa = pga;
3145
c = 2;
3146
}
3147
if((pgb-=gpo) > pa){
3148
pa = pgb;
3149
c = 4;
3150
}
3151
3152
pa += subp[seq2[0]];
3153
3154
s[0].a = pa;
3155
3156
s[0].ga = -FLOATINFTY;
3157
3158
pgb = s[0].gb;
3159
s[0].gb = ca-(gpo+tgpe);
3160
if(pgb-tgpe > s[0].gb){
3161
s[0].gb = pgb-tgpe;
3162
c |= 16;
3163
}
3164
tracep[0] = c;
3165
}
3166
3167
subp = subm[seq1[0]];
3168
tracep = trace[0];
3169
j = len_b;
3170
pa = s[j].a;
3171
pga = s[j].ga;
3172
pgb = s[j].gb;
3173
s[j].a = -FLOATINFTY;
3174
s[j].ga = -FLOATINFTY;
3175
3176
s[j].gb = pa-tgpe;//-gpo;
3177
if(pgb-tgpe > s[j].gb){
3178
s[j].gb = pgb-tgpe;
3179
}
3180
3181
//s[j].gb = -INFTY;
3182
while(--j){
3183
3184
ca = s[j].a;
3185
3186
c = 1;
3187
3188
if((pga-=gpo) > pa){
3189
pa = pga;
3190
c = 2;
3191
}
3192
3193
if((pgb-=gpo) > pa){
3194
pa = pgb;
3195
c = 4;
3196
}
3197
3198
pa += subp[seq2[j]];
3199
3200
s[j].a = pa;
3201
3202
pga = s[j].ga;
3203
s[j].ga = s[j+1].a-(gpo+tgpe);
3204
if (s[j+1].ga-tgpe > s[j].ga){
3205
s[j].ga = s[j+1].ga-tgpe;
3206
c |= 8;
3207
}
3208
pgb = s[j].gb;
3209
s[j].gb = -FLOATINFTY;
3210
tracep[j] = c;
3211
pa = ca;
3212
}
3213
3214
ca = s[0].a;
3215
3216
c = 1;
3217
3218
if((pga-=gpo) > pa){
3219
pa = pga;
3220
c = 2;
3221
}
3222
if((pgb-=gpo) > pa){
3223
pa = pgb;
3224
c = 4;
3225
}
3226
3227
pa += subp[seq2[0]];
3228
3229
s[0].a = pa;
3230
3231
3232
s[0].ga = s[1].a-(gpo+tgpe);
3233
if (s[1].ga-tgpe > s[0].ga){
3234
s[0].ga = s[1].ga-tgpe;
3235
c |= 8;
3236
}
3237
3238
pgb = s[0].gb;
3239
s[0].gb = ca-(gpo+tgpe);
3240
if(pgb-tgpe > s[0].gb){
3241
s[0].gb = pgb-tgpe;
3242
c |= 16;
3243
}
3244
tracep[0] = c;
3245
3246
3247
pgb = s[0].gb;
3248
c = 2;
3249
if(s[0].ga > pgb){
3250
pgb = s[0].ga;
3251
c = 1;
3252
}
3253
if(s[0].a >= pgb){
3254
pgb = s[0].a;
3255
c = 0;
3256
}
3257
3258
f = c;
3259
3260
i = 0;
3261
j = 0;
3262
c = 1;
3263
while(trace[i][j] < 32){
3264
// fprintf(stderr,"%d->%d %d:%d %d:%d\n",c,trace[i][j],i,j,len_a,len_b);
3265
switch(f){
3266
case 0:
3267
if (trace[i][j] & 2){
3268
f = 1;
3269
if(i+1!= len_a){
3270
path[c+1] |= 16;
3271
// fprintf(stderr,"GAP_CLOSE\n");
3272
}else{
3273
path[c+1] |= 32+16;
3274
}
3275
}else if (trace[i][j] & 4){
3276
f = 2;
3277
if(j+1!= len_b){
3278
path[c+1] |= 16;
3279
// fprintf(stderr,"GAP_CLOSE\n");
3280
}else{
3281
path[c+1] |= 32+16;
3282
}
3283
}
3284
3285
//path[c] = 0;
3286
i++;
3287
j++;
3288
break;
3289
case 1:
3290
if(trace[i][j] & 8){
3291
f = 1;
3292
if(i!=0 && i!= len_a){
3293
// / fprintf(stderr,"GAP_EXT\n");
3294
if(!(path[c]&16)){
3295
path[c] |= 8;
3296
}
3297
}else{
3298
if(!(path[c]&16)){
3299
path[c] |= 32+8;
3300
}
3301
}
3302
}else{
3303
f = 0;
3304
if(i!=0 && i!= len_a){
3305
// fprintf(stderr,"GAP_OPEN\n");
3306
path[c] |= 4;
3307
}else{
3308
path[c] |= 32+4;
3309
}
3310
}
3311
path[c] |= 1;
3312
j++;
3313
break;
3314
case 2:
3315
if(trace[i][j] & 16){
3316
f = 2;
3317
if(j !=0 && j != len_b){
3318
// fprintf(stderr,"GAP_EXT\n");
3319
if(!(path[c]&16)){
3320
path[c] |= 8;
3321
}
3322
}else{
3323
if(!(path[c]&16)){
3324
path[c] |= 32+8;
3325
}
3326
}
3327
}else{
3328
f = 0;
3329
if(j!=0 && j != len_b){
3330
// fprintf(stderr,"GAP_OPEN\n");
3331
path[c] |= 4;
3332
}else{
3333
path[c] |= 32+4;
3334
}
3335
3336
}
3337
path[c] |= 2;
3338
i++;
3339
break;
3340
}
3341
c++;
3342
}
3343
path[0] = c-1;
3344
path[c] = 3;
3345
path[c+1] = pgb;
3346
return path;
3347
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1