1 #include "mltaln.h"
2 #include "dp.h"
3
4 #define DEBUG 0
5 #define XXXXXXX 0
6 #define USE_PENALTY_EX 1
7
8 #if 1
match_calc_mtx(double ** mtx,double * match,char ** s1,char ** s2,int i1,int lgth2)9 static void match_calc_mtx( double **mtx, double *match, char **s1, char **s2, int i1, int lgth2 )
10 {
11 char *seq2 = s2[0];
12 double *doubleptr = mtx[(int)s1[0][i1]];
13
14 while( lgth2-- )
15 *match++ = doubleptr[(int)*seq2++];
16 }
17 #else
match_calc(double * match,char ** s1,char ** s2,int i1,int lgth2)18 static void match_calc( double *match, char **s1, char **s2, int i1, int lgth2 )
19 {
20 int j;
21
22 for( j=0; j<lgth2; j++ )
23 match[j] = amino_dis[(*s1)[i1]][(*s2)[j]];
24 }
25 #endif
26
Atracking(double * lasthorizontalw,double * lastverticalw,char ** seq1,char ** seq2,char ** mseq1,char ** mseq2,int ** ijp,int tailgp,int * warpis,int * warpjs,int warpbase)27 static double Atracking( double *lasthorizontalw, double *lastverticalw,
28 char **seq1, char **seq2,
29 char **mseq1, char **mseq2,
30 int **ijp,
31 int tailgp,
32 int *warpis, int *warpjs, int warpbase )
33 {
34 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, limk;
35 // char gap[] = "-";
36 char *gap;
37 gap = newgapstr;
38 lgth1 = strlen( seq1[0] );
39 lgth2 = strlen( seq2[0] );
40 double wm;
41
42
43 #if 0
44 for( i=0; i<lgth1; i++ )
45 {
46 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
47 }
48 #endif
49
50 for( i=0; i<lgth1+1; i++ )
51 {
52 ijp[i][0] = i + 1;
53 }
54 for( j=0; j<lgth2+1; j++ )
55 {
56 ijp[0][j] = -( j + 1 );
57 }
58
59 // if( tailgp == 1 || ijp[lgth1][lgth2] >= warpbase )
60 if( tailgp == 1 )
61 ;
62 else
63 {
64 wm = lastverticalw[0];
65 for( i=0; i<lgth1; i++ )
66 {
67 if( lastverticalw[i] >= wm )
68 {
69 wm = lastverticalw[i];
70 iin = i; jin = lgth2-1;
71 ijp[lgth1][lgth2] = +( lgth1 - i );
72 }
73 }
74 for( j=0; j<lgth2; j++ )
75 {
76 if( lasthorizontalw[j] >= wm )
77 {
78 wm = lasthorizontalw[j];
79 iin = lgth1-1; jin = j;
80 ijp[lgth1][lgth2] = -( lgth2 - j );
81 }
82 }
83 }
84
85
86
87 mseq1[0] += lgth1+lgth2;
88 *mseq1[0] = 0;
89 mseq2[0] += lgth1+lgth2;
90 *mseq2[0] = 0;
91
92
93
94 iin = lgth1; jin = lgth2;
95 limk = lgth1+lgth2 + 1;
96 for( k=0; k<limk; k++ )
97 {
98 if( ijp[iin][jin] >= warpbase )
99 {
100 // fprintf( stderr, "WARP!\n" );
101 ifi = warpis[ijp[iin][jin]-warpbase];
102 jfi = warpjs[ijp[iin][jin]-warpbase];
103 }
104 else if( ijp[iin][jin] < 0 )
105 {
106 ifi = iin-1; jfi = jin+ijp[iin][jin];
107 }
108 else if( ijp[iin][jin] > 0 )
109 {
110 ifi = iin-ijp[iin][jin]; jfi = jin-1;
111 }
112 else
113 {
114 ifi = iin-1; jfi = jin-1;
115 }
116
117 if( ifi == -warpbase && jfi == -warpbase )
118 {
119 l = iin;
120 while( --l >= 0 )
121 {
122 *--mseq1[0] = seq1[0][l];
123 *--mseq2[0] = *gap;
124 k++;
125 }
126 l= jin;
127 while( --l >= 0 )
128 {
129 *--mseq1[0] = *gap;
130 *--mseq2[0] = seq2[0][l];
131 k++;
132 }
133 break;
134 }
135 else
136 {
137 l = iin - ifi;
138 while( --l > 0 )
139 {
140 *--mseq1[0] = seq1[0][ifi+l];
141 *--mseq2[0] = *gap;
142 k++;
143 }
144 l= jin - jfi;
145 while( --l > 0 )
146 {
147 *--mseq1[0] = *gap;
148 *--mseq2[0] = seq2[0][jfi+l];
149 k++;
150 }
151 }
152 if( iin <= 0 || jin <= 0 ) break;
153 *--mseq1[0] = seq1[0][ifi];
154 *--mseq2[0] = seq2[0][jfi];
155 k++;
156 iin = ifi; jin = jfi;
157 }
158
159 // fprintf( stderr, "%s\n", mseq1[0] );
160 // fprintf( stderr, "%s\n", mseq2[0] );
161 return( 0.0 );
162 }
163
164
G__align11(double ** n_dynamicmtx,char ** seq1,char ** seq2,int alloclen,int headgp,int tailgp)165 double G__align11( double **n_dynamicmtx, char **seq1, char **seq2, int alloclen, int headgp, int tailgp )
166 {
167 // int k;
168 register int i, j;
169 int lasti; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
170 int lastj;
171 int lgth1, lgth2;
172 int resultlen;
173 double wm; /* int ?????? */
174 double g;
175 double *currentw, *previousw;
176 double fpenalty = (double)penalty;
177 double fpenalty_shift = (double)penalty_shift;
178 double fpenalty_tmp;
179 #if USE_PENALTY_EX
180 double fpenalty_ex = (double)penalty_ex;
181 #endif
182 #if 1
183 double *wtmp;
184 int *ijppt;
185 double *mjpt, *prept, *curpt;
186 int *mpjpt;
187 #endif
188 static TLS double mi = 0.0;
189 static TLS double *m = NULL;
190 static TLS int **ijp = NULL;
191 static TLS int mpi = 0;
192 static TLS int *mp = NULL;
193 static TLS double *w1 = NULL;
194 static TLS double *w2 = NULL;
195 static TLS double *match = NULL;
196 static TLS double *initverticalw = NULL; /* kufuu sureba iranai */
197 static TLS double *lastverticalw = NULL; /* kufuu sureba iranai */
198 static TLS char **mseq1 = NULL;
199 static TLS char **mseq2 = NULL;
200 static TLS char **mseq = NULL;
201 static TLS int **intwork = NULL;
202 static TLS double **doublework = NULL;
203 static TLS int orlgth1 = 0, orlgth2 = 0;
204 static TLS double **amino_dynamicmtx = NULL; // ??
205
206 int *warpis = NULL;
207 int *warpjs = NULL;
208 int *warpi = NULL;
209 int *warpj = NULL;
210 int *prevwarpi = NULL;
211 int *prevwarpj = NULL;
212 double *wmrecords = NULL;
213 double *prevwmrecords = NULL;
214 int warpn = 0;
215 int warpbase;
216 double curm = 0.0;
217 double *wmrecordspt, *wmrecords1pt, *prevwmrecordspt;
218 int *warpipt, *warpjpt;
219
220
221 if( seq1 == NULL )
222 {
223 if( orlgth1 > 0 && orlgth2 > 0 )
224 {
225 orlgth1 = 0;
226 orlgth2 = 0;
227 if( mseq1 ) free( mseq1 ); mseq1 = NULL;
228 if( mseq2 ) free( mseq2 ); mseq2 = NULL;
229 if( w1 ) FreeFloatVec( w1 ); w1 = NULL;
230 if( w2 ) FreeFloatVec( w2 ); w2 = NULL;
231 if( match ) FreeFloatVec( match ); match = NULL;
232 if( initverticalw ) FreeFloatVec( initverticalw ); initverticalw = NULL;
233 if( lastverticalw ) FreeFloatVec( lastverticalw ); lastverticalw = NULL;
234
235 if( m ) FreeFloatVec( m ); m = NULL;
236 if( mp ) FreeIntVec( mp ); mp = NULL;
237
238 if( mseq ) FreeCharMtx( mseq ); mseq = NULL;
239
240
241
242 if( doublework ) FreeFloatMtx( doublework ); doublework = NULL;
243 if( intwork ) FreeIntMtx( intwork ); intwork = NULL;
244
245 if( amino_dynamicmtx ) FreeDoubleMtx( amino_dynamicmtx ); amino_dynamicmtx = NULL;
246 }
247 orlgth1 = 0;
248 orlgth2 = 0;
249 return( 0.0 );
250 }
251
252
253
254 lgth1 = strlen( seq1[0] );
255 lgth2 = strlen( seq2[0] );
256
257 warpbase = lgth1 + lgth2;
258 warpis = NULL;
259 warpjs = NULL;
260 warpn = 0;
261 if( trywarp )
262 {
263 // fprintf( stderr, "IN G__align11\n" );
264 if( headgp == 0 || tailgp == 0 )
265 {
266 fprintf( stderr, "At present, headgp and tailgp must be 1.\n" );
267 exit( 1 );
268 }
269
270 wmrecords = AllocateFloatVec( lgth2+1 );
271 warpi = AllocateIntVec( lgth2+1 );
272 warpj = AllocateIntVec( lgth2+1 );
273 prevwmrecords = AllocateFloatVec( lgth2+1 );
274 prevwarpi = AllocateIntVec( lgth2+1 );
275 prevwarpj = AllocateIntVec( lgth2+1 );
276 for( i=0; i<lgth2+1; i++ ) prevwmrecords[i] = 0.0;
277 for( i=0; i<lgth2+1; i++ ) wmrecords[i] = 0.0;
278 for( i=0; i<lgth2+1; i++ ) prevwarpi[i] = -warpbase;
279 for( i=0; i<lgth2+1; i++ ) prevwarpj[i] = -warpbase;
280 for( i=0; i<lgth2+1; i++ ) warpi[i] = -warpbase;
281 for( i=0; i<lgth2+1; i++ ) warpj[i] = -warpbase;
282 }
283
284 #if 0
285 if( lgth1 <= 0 || lgth2 <= 0 )
286 {
287 fprintf( stderr, "WARNING (g11): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
288 }
289 #endif
290
291 #if 1
292 if( lgth1 == 0 && lgth2 == 0 )
293 return( 0.0 );
294
295 if( lgth1 == 0 )
296 {
297 seq1[0][lgth2] = 0;
298 while( lgth2 ) seq1[0][--lgth2] = *newgapstr;
299 // fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
300 return( 0.0 );
301 }
302
303 if( lgth2 == 0 )
304 {
305 seq2[0][lgth1] = 0;
306 while( lgth1 ) seq2[0][--lgth1] = *newgapstr;
307 // fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
308 return( 0.0 );
309 }
310 #endif
311
312
313 wm = 0.0;
314
315 if( orlgth1 == 0 )
316 {
317 mseq1 = AllocateCharMtx( njob, 0 );
318 mseq2 = AllocateCharMtx( njob, 0 );
319 }
320
321
322
323 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
324 {
325 int ll1, ll2;
326
327 if( orlgth1 > 0 && orlgth2 > 0 )
328 {
329 FreeFloatVec( w1 );
330 FreeFloatVec( w2 );
331 FreeFloatVec( match );
332 FreeFloatVec( initverticalw );
333 FreeFloatVec( lastverticalw );
334
335 FreeFloatVec( m );
336 FreeIntVec( mp );
337
338 FreeCharMtx( mseq );
339
340
341
342 FreeFloatMtx( doublework );
343 FreeIntMtx( intwork );
344 FreeDoubleMtx( amino_dynamicmtx );
345 }
346
347 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
348 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
349
350 #if DEBUG
351 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
352 #endif
353
354 w1 = AllocateFloatVec( ll2+2 );
355 w2 = AllocateFloatVec( ll2+2 );
356 match = AllocateFloatVec( ll2+2 );
357
358 initverticalw = AllocateFloatVec( ll1+2 );
359 lastverticalw = AllocateFloatVec( ll1+2 );
360
361 m = AllocateFloatVec( ll2+2 );
362 mp = AllocateIntVec( ll2+2 );
363
364 mseq = AllocateCharMtx( njob, ll1+ll2 );
365
366
367 doublework = AllocateFloatMtx( nalphabets, MAX( ll1, ll2 )+2 );
368 intwork = AllocateIntMtx( nalphabets, MAX( ll1, ll2 )+2 );
369 amino_dynamicmtx = AllocateDoubleMtx( 0x80, 0x80 );
370
371 #if DEBUG
372 fprintf( stderr, "succeeded\n" );
373 #endif
374
375 orlgth1 = ll1 - 100;
376 orlgth2 = ll2 - 100;
377 }
378 for( i=0; i<nalphabets; i++) for( j=0; j<nalphabets; j++ )
379 amino_dynamicmtx[(int)amino[i]][(int)amino[j]] = (double)n_dynamicmtx[i][j];
380
381
382 mseq1[0] = mseq[0];
383 mseq2[0] = mseq[1];
384
385
386 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
387 {
388 int ll1, ll2;
389
390 if( commonAlloc1 && commonAlloc2 )
391 {
392 FreeIntMtx( commonIP );
393 }
394
395 ll1 = MAX( orlgth1, commonAlloc1 );
396 ll2 = MAX( orlgth2, commonAlloc2 );
397
398 #if DEBUG
399 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
400 #endif
401
402 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
403
404 #if DEBUG
405 fprintf( stderr, "succeeded\n\n" );
406 #endif
407
408 commonAlloc1 = ll1;
409 commonAlloc2 = ll2;
410 }
411 ijp = commonIP;
412
413
414 #if 0
415 for( i=0; i<lgth1; i++ )
416 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
417 #endif
418
419 currentw = w1;
420 previousw = w2;
421
422
423 match_calc_mtx( amino_dynamicmtx, initverticalw, seq2, seq1, 0, lgth1 );
424 match_calc_mtx( amino_dynamicmtx, currentw, seq1, seq2, 0, lgth2 );
425
426 if( headgp == 1 )
427 {
428 for( i=1; i<lgth1+1; i++ )
429 {
430 initverticalw[i] += fpenalty;
431 }
432 for( j=1; j<lgth2+1; j++ )
433 {
434 currentw[j] += fpenalty;
435 }
436 }
437
438 for( j=1; j<lgth2+1; ++j )
439 {
440 m[j] = currentw[j-1]; mp[j] = 0;
441 }
442
443 if( lgth2 == 0 )
444 lastverticalw[0] = 0.0; // lgth2==0 no toki error
445 else
446 lastverticalw[0] = currentw[lgth2-1]; // lgth2==0 no toki error
447
448 if( tailgp ) lasti = lgth1+1; else lasti = lgth1;
449 lastj = lgth2+1;
450
451 #if XXXXXXX
452 fprintf( stderr, "currentw = \n" );
453 for( i=0; i<lgth1+1; i++ )
454 {
455 fprintf( stderr, "%5.2f ", currentw[i] );
456 }
457 fprintf( stderr, "\n" );
458 fprintf( stderr, "initverticalw = \n" );
459 for( i=0; i<lgth2+1; i++ )
460 {
461 fprintf( stderr, "%5.2f ", initverticalw[i] );
462 }
463 fprintf( stderr, "\n" );
464 #endif
465
466 for( i=1; i<lasti; i++ )
467 {
468 wtmp = previousw;
469 previousw = currentw;
470 currentw = wtmp;
471
472 previousw[0] = initverticalw[i-1];
473
474 match_calc_mtx( amino_dynamicmtx, currentw, seq1, seq2, i, lgth2 );
475 #if XXXXXXX
476 fprintf( stderr, "\n" );
477 fprintf( stderr, "i=%d\n", i );
478 fprintf( stderr, "currentw = \n" );
479 for( j=0; j<lgth2; j++ )
480 {
481 fprintf( stderr, "%5.2f ", currentw[j] );
482 }
483 fprintf( stderr, "\n" );
484 #endif
485 #if XXXXXXX
486 fprintf( stderr, "\n" );
487 fprintf( stderr, "i=%d\n", i );
488 fprintf( stderr, "currentw = \n" );
489 for( j=0; j<lgth2; j++ )
490 {
491 fprintf( stderr, "%5.2f ", currentw[j] );
492 }
493 fprintf( stderr, "\n" );
494 #endif
495 currentw[0] = initverticalw[i];
496
497 mi = previousw[0]; mpi = 0;
498
499 ijppt = ijp[i] + 1;
500 mjpt = m + 1;
501 prept = previousw;
502 curpt = currentw + 1;
503 mpjpt = mp + 1;
504 if( trywarp )
505 {
506 prevwmrecordspt = prevwmrecords;
507 wmrecordspt = wmrecords+1;
508 wmrecords1pt = wmrecords;
509 warpipt = warpi + 1;
510 warpjpt = warpj + 1;
511 }
512
513 for( j=1; j<lastj; j++ )
514 {
515
516 wm = *prept;
517 *ijppt = 0;
518
519 #if 0
520 fprintf( stderr, "%5.0f->", wm );
521 #endif
522 #if 0
523 fprintf( stderr, "%5.0f?", g );
524 #endif
525 if( (g=mi+fpenalty) > wm )
526 {
527 wm = g;
528 *ijppt = -( j - mpi );
529 }
530 if( (g=*prept) >= mi )
531 {
532 mi = g;
533 mpi = j-1;
534 }
535 #if USE_PENALTY_EX
536 mi += fpenalty_ex;
537 #endif
538
539 #if 0
540 fprintf( stderr, "%5.0f?", g );
541 #endif
542 if( (g=*mjpt + fpenalty) > wm )
543 {
544 wm = g;
545 *ijppt = +( i - *mpjpt );
546 }
547 if( (g=*prept) >= *mjpt )
548 {
549 *mjpt = g;
550 *mpjpt = i-1;
551 }
552 #if USE_PENALTY_EX
553 m[j] += fpenalty_ex;
554 #endif
555 #if 1
556 if( trywarp )
557 {
558 fpenalty_tmp = fpenalty_shift + fpenalty_ex * ( i - prevwarpi[j-1] + j - prevwarpj[j-1] );
559 // fprintf( stderr, "fpenalty_shift = %f\n", fpenalty_tmp );
560
561 // fprintf( stderr, "\n\n\nwarp to %c-%c (%d-%d) from %c-%c (%d-%d) ? prevwmrecords[%d] = %f + %f <- wm = %f\n", seq1[0][prevwarpi[j-1]], seq2[0][prevwarpj[j-1]], prevwarpi[j-1], prevwarpj[j-1], seq1[0][i], seq2[0][j], i, j, j, prevwmrecords[j-1], fpenalty_tmp, wm );
562 // if( (g=prevwmrecords[j-1] + fpenalty_shift )> wm )
563 if( ( g=*prevwmrecordspt++ + fpenalty_tmp )> wm ) // naka ha osokute kamawanai
564 {
565 // fprintf( stderr, "Yes! Warp!! from %d-%d (%c-%c) to %d-%d (%c-%c) fpenalty_tmp = %f! warpn = %d\n", i, j, seq1[0][i], seq2[0][j-1], prevwarpi[j-1], prevwarpj[j-1],seq1[0][prevwarpi[j-1]], seq2[0][prevwarpj[j-1]], fpenalty_tmp, warpn );
566 if( warpn && prevwarpi[j-1] == warpis[warpn-1] && prevwarpj[j-1] == warpjs[warpn-1] )
567 {
568 *ijppt = warpbase + warpn - 1;
569 }
570 else
571 {
572 *ijppt = warpbase + warpn;
573 warpis = realloc( warpis, sizeof(int) * ( warpn+1 ) );
574 warpjs = realloc( warpjs, sizeof(int) * ( warpn+1 ) );
575 warpis[warpn] = prevwarpi[j-1];
576 warpjs[warpn] = prevwarpj[j-1];
577 warpn++;
578 }
579 wm = g;
580 }
581 else
582 {
583 }
584
585 curm = *curpt + wm;
586
587 // fprintf( stderr, "###### curm = %f at %c-%c, i=%d, j=%d\n", curm, seq1[0][i], seq2[0][j], i, j );
588
589 // fprintf( stderr, "copy from i, j-1? %f > %f?\n", wmrecords[j-1], curm );
590 // if( wmrecords[j-1] > wmrecords[j] )
591 if( *wmrecords1pt > *wmrecordspt )
592 {
593 // fprintf( stderr, "yes\n" );
594 // wmrecords[j] = wmrecords[j-1];
595 *wmrecordspt = *wmrecords1pt;
596 // warpi[j] = warpi[j-1];
597 // warpj[j] = warpj[j-1];
598 *warpipt = *(warpipt-1);
599 *warpjpt = *(warpjpt-1);
600 // fprintf( stderr, "warpi[j]=%d, warpj[j]=%d wmrecords[j] = %f\n", warpi[j], warpj[j], wmrecords[j] );
601 }
602 // else
603 // {
604 // fprintf( stderr, "no\n" );
605 // }
606
607 // fprintf( stderr, " curm = %f at %c-%c\n", curm, seq1[0][i], seq2[0][j] );
608 // fprintf( stderr, " wmrecords[%d] = %f\n", j, wmrecords[j] );
609 // fprintf( stderr, "replace?\n" );
610
611 // if( curm > wmrecords[j] )
612 if( curm > *wmrecordspt )
613 {
614 // fprintf( stderr, "yes at %d-%d (%c-%c), replaced warp: warpi[j]=%d, warpj[j]=%d warpn=%d, wmrecords[j] = %f -> %f\n", i, j, seq1[0][i], seq2[0][j], i, j, warpn, wmrecords[j], curm );
615 // wmrecords[j] = curm;
616 *wmrecordspt = curm;
617 // warpi[j] = i;
618 // warpj[j] = j;
619 *warpipt = i;
620 *warpjpt = j;
621 }
622 // else
623 // {
624 // fprintf( stderr, "No! warpi[j]=%d, warpj[j]=%d wmrecords[j] = %f\n", warpi[j], warpj[j], wmrecords[j] );
625 // }
626 // fprintf( stderr, "%d-%d (%c-%c) curm = %5.0f, wmrecords[j]=%f\n", i, j, seq1[0][i], seq2[0][j], curm, wmrecords[j] );
627 wmrecordspt++;
628 wmrecords1pt++;
629 warpipt++;
630 warpjpt++;
631 }
632 #endif
633
634 *curpt++ += wm;
635 ijppt++;
636 mjpt++;
637 prept++;
638 mpjpt++;
639 }
640 lastverticalw[i] = currentw[lgth2-1]; // lgth2==0 no toki error
641
642 if( trywarp )
643 {
644 fltncpy( prevwmrecords, wmrecords, lastj );
645 intncpy( prevwarpi, warpi, lastj );
646 intncpy( prevwarpj, warpj, lastj );
647 }
648 }
649
650 if( trywarp )
651 {
652 // fprintf( stderr, "\nwm = %f\n", wm );
653 // fprintf( stderr, "warpn = %d\n", warpn );
654 free( wmrecords );
655 free( prevwmrecords );
656 free( warpi );
657 free( warpj );
658 free( prevwarpi );
659 free( prevwarpj );
660 }
661
662 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, tailgp, warpis, warpjs, warpbase );
663 if( warpis ) free( warpis );
664 if( warpjs ) free( warpjs );
665
666
667 resultlen = strlen( mseq1[0] );
668 if( alloclen < resultlen || resultlen > N )
669 {
670 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
671 ErrorExit( "LENGTH OVER!\n" );
672 }
673
674
675 strcpy( seq1[0], mseq1[0] );
676 strcpy( seq2[0], mseq2[0] );
677 #if 0
678 fprintf( stderr, "\n" );
679 fprintf( stderr, ">\n%s\n", mseq1[0] );
680 fprintf( stderr, ">\n%s\n", mseq2[0] );
681 fprintf( stderr, "wm = %f\n", wm );
682 #endif
683
684 return( wm );
685 }
686
G__align11_noalign(double ** n_dynamicmtx,int penal,int penal_ex,char ** seq1,char ** seq2,int alloclen)687 double G__align11_noalign( double **n_dynamicmtx, int penal, int penal_ex, char **seq1, char **seq2, int alloclen )
688 /* warp mitaiou */
689 {
690 // int k;
691 register int i, j;
692 int lasti; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
693 int lgth1, lgth2;
694 // int resultlen;
695 double wm; /* int ?????? */
696 double g;
697 double *currentw, *previousw;
698 double fpenalty = (double)penal;
699 #if USE_PENALTY_EX
700 double fpenalty_ex = (double)penal_ex;
701 #endif
702 #if 1
703 double *wtmp;
704 double *mjpt, *prept, *curpt;
705 // int *mpjpt;
706 #endif
707 static TLS double mi, *m;
708 static TLS double *w1, *w2;
709 static TLS double *match;
710 static TLS double *initverticalw; /* kufuu sureba iranai */
711 static TLS double *lastverticalw; /* kufuu sureba iranai */
712 static TLS int **intwork;
713 static TLS double **doublework;
714 static TLS int orlgth1 = 0, orlgth2 = 0;
715 static TLS double **amino_dynamicmtx;
716
717 if( seq1 == NULL )
718 {
719 if( orlgth1 > 0 && orlgth2 > 0 )
720 {
721 orlgth1 = 0;
722 orlgth2 = 0;
723 FreeFloatVec( w1 );
724 FreeFloatVec( w2 );
725 FreeFloatVec( match );
726 FreeFloatVec( initverticalw );
727 FreeFloatVec( lastverticalw );
728 free( m );
729
730
731 FreeFloatMtx( doublework );
732 FreeIntMtx( intwork );
733 FreeDoubleMtx( amino_dynamicmtx );
734 }
735 return( 0.0 );
736 }
737
738
739 wm = 0.0;
740
741
742
743 lgth1 = strlen( seq1[0] );
744 lgth2 = strlen( seq2[0] );
745
746
747
748 #if 0
749 if( lgth1 <= 0 || lgth2 <= 0 )
750 {
751 fprintf( stderr, "WARNING (g11): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
752 }
753 #endif
754
755 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
756 {
757 int ll1, ll2;
758
759 if( orlgth1 > 0 && orlgth2 > 0 )
760 {
761 FreeFloatVec( w1 );
762 FreeFloatVec( w2 );
763 FreeFloatVec( match );
764 FreeFloatVec( initverticalw );
765 FreeFloatVec( lastverticalw );
766
767 FreeFloatVec( m );
768
769
770
771
772 FreeFloatMtx( doublework );
773 FreeIntMtx( intwork );
774
775 FreeDoubleMtx( amino_dynamicmtx );
776 }
777
778 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
779 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
780
781 #if DEBUG
782 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
783 #endif
784
785 w1 = AllocateFloatVec( ll2+2 );
786 w2 = AllocateFloatVec( ll2+2 );
787 match = AllocateFloatVec( ll2+2 );
788
789 initverticalw = AllocateFloatVec( ll1+2 );
790 lastverticalw = AllocateFloatVec( ll1+2 );
791
792 m = AllocateFloatVec( ll2+2 );
793
794
795
796 doublework = AllocateFloatMtx( nalphabets, MAX( ll1, ll2 )+2 );
797 intwork = AllocateIntMtx( nalphabets, MAX( ll1, ll2 )+2 );
798
799
800 amino_dynamicmtx = AllocateDoubleMtx( 0x80, 0x80 );
801 #if DEBUG
802 fprintf( stderr, "succeeded\n" );
803 #endif
804
805 orlgth1 = ll1 - 100;
806 orlgth2 = ll2 - 100;
807 }
808
809
810 for( i=0; i<nalphabets; i++) for( j=0; j<nalphabets; j++ )
811 amino_dynamicmtx[(int)amino[i]][(int)amino[j]] = (double)n_dynamicmtx[i][j];
812
813
814
815
816 #if 0
817 for( i=0; i<lgth1; i++ )
818 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
819 #endif
820
821 currentw = w1;
822 previousw = w2;
823
824
825 match_calc_mtx( amino_dynamicmtx, initverticalw, seq2, seq1, 0, lgth1 );
826
827
828 match_calc_mtx( amino_dynamicmtx, currentw, seq1, seq2, 0, lgth2 );
829
830 if( 1 ) // tsuneni outgap-1
831 {
832 for( i=1; i<lgth1+1; i++ )
833 {
834 initverticalw[i] += fpenalty;
835 }
836 for( j=1; j<lgth2+1; j++ )
837 {
838 currentw[j] += fpenalty;
839 }
840 }
841
842 for( j=1; j<lgth2+1; ++j )
843 {
844 m[j] = currentw[j-1];
845 }
846
847 if( lgth2 == 0 )
848 lastverticalw[0] = 0.0; // lgth2==0 no toki error
849 else
850 lastverticalw[0] = currentw[lgth2-1]; // lgth2==0 no toki error
851
852 if( 1 ) lasti = lgth1+1; else lasti = lgth1; // tsuneni outgap-1
853
854 #if XXXXXXX
855 fprintf( stderr, "currentw = \n" );
856 for( i=0; i<lgth1+1; i++ )
857 {
858 fprintf( stderr, "%5.2f ", currentw[i] );
859 }
860 fprintf( stderr, "\n" );
861 fprintf( stderr, "initverticalw = \n" );
862 for( i=0; i<lgth2+1; i++ )
863 {
864 fprintf( stderr, "%5.2f ", initverticalw[i] );
865 }
866 fprintf( stderr, "\n" );
867 #endif
868
869 for( i=1; i<lasti; i++ )
870 {
871 wtmp = previousw;
872 previousw = currentw;
873 currentw = wtmp;
874
875 previousw[0] = initverticalw[i-1];
876
877 match_calc_mtx( amino_dynamicmtx, currentw, seq1, seq2, i, lgth2 );
878 #if XXXXXXX
879 fprintf( stderr, "\n" );
880 fprintf( stderr, "i=%d\n", i );
881 fprintf( stderr, "currentw = \n" );
882 for( j=0; j<lgth2; j++ )
883 {
884 fprintf( stderr, "%5.2f ", currentw[j] );
885 }
886 fprintf( stderr, "\n" );
887 #endif
888 #if XXXXXXX
889 fprintf( stderr, "\n" );
890 fprintf( stderr, "i=%d\n", i );
891 fprintf( stderr, "currentw = \n" );
892 for( j=0; j<lgth2; j++ )
893 {
894 fprintf( stderr, "%5.2f ", currentw[j] );
895 }
896 fprintf( stderr, "\n" );
897 #endif
898 currentw[0] = initverticalw[i];
899
900 mi = previousw[0];
901
902 mjpt = m + 1;
903 prept = previousw;
904 curpt = currentw + 1;
905 for( j=1; j<lgth2+1; j++ )
906 {
907 wm = *prept;
908
909 #if 0
910 fprintf( stderr, "%5.0f->", wm );
911 #endif
912 #if 0
913 fprintf( stderr, "%5.0f?", g );
914 #endif
915 if( (g=mi+fpenalty) > wm )
916 {
917 wm = g;
918 }
919 if( (g=*prept) >= mi )
920 {
921 mi = g;
922 }
923 #if USE_PENALTY_EX
924 mi += fpenalty_ex;
925 #endif
926
927 #if 0
928 fprintf( stderr, "%5.0f?", g );
929 #endif
930 if( (g=*mjpt + fpenalty) > wm )
931 {
932 wm = g;
933 }
934 if( (g=*prept) >= *mjpt )
935 {
936 *mjpt = g;
937 }
938 #if USE_PENALTY_EX
939 m[j] += fpenalty_ex;
940 #endif
941
942 #if 0
943 fprintf( stderr, "%5.0f ", wm );
944 #endif
945 *curpt++ += wm;
946 mjpt++;
947 prept++;
948 }
949 lastverticalw[i] = currentw[lgth2-1]; // lgth2==0 no toki error
950 }
951
952 #if 0
953 fprintf( stderr, "\n" );
954 fprintf( stderr, ">\n%s\n", mseq1[0] );
955 fprintf( stderr, ">\n%s\n", mseq2[0] );
956 fprintf( stderr, "wm = %f\n", wm );
957 #endif
958
959 return( wm );
960 }
961