1 #include "mltaln.h"
2 #include "dp.h"
3
4 #define DEBUG 0
5 #define DEBUG2 0
6 #define XXXXXXX 0
7 #define USE_PENALTY_EX 1
8
9
10 static TLS int localstop; // 060910
11
12 #if 1
match_calc_mtx(double ** mtx,double * match,char ** s1,char ** s2,int i1,int lgth2)13 static void match_calc_mtx( double **mtx, double *match, char **s1, char **s2, int i1, int lgth2 )
14 {
15 char *seq2 = s2[0];
16 double *doubleptr = mtx[(int)s1[0][i1]];
17
18 while( lgth2-- )
19 *match++ = doubleptr[(int)*seq2++];
20 }
21 #else
match_calc(double * match,char ** s1,char ** s2,int i1,int lgth2)22 static void match_calc( double *match, char **s1, char **s2, int i1, int lgth2 )
23 {
24 int j;
25
26 for( j=0; j<lgth2; j++ )
27 match[j] = amino_dis[(*s1)[i1]][(*s2)[j]];
28 }
29 #endif
30
31 #if 0
32 static void match_calc_bk( double *match, double **cpmx1, double **cpmx2, int i1, int lgth2, double **doublework, int **intwork, int initialize )
33 {
34 int j, k, l;
35 double scarr[nalphabets];
36 double **cpmxpd = doublework;
37 int **cpmxpdn = intwork;
38 int count = 0;
39
40 if( initialize )
41 {
42 for( j=0; j<lgth2; j++ )
43 {
44 count = 0;
45 for( l=0; l<nalphabets; l++ )
46 {
47 if( cpmx2[l][j] )
48 {
49 cpmxpd[count][j] = cpmx2[l][j];
50 cpmxpdn[count][j] = l;
51 count++;
52 }
53 }
54 cpmxpdn[count][j] = -1;
55 }
56 }
57
58 for( l=0; l<nalphabets; l++ )
59 {
60 scarr[l] = 0.0;
61 for( k=0; k<nalphabets; k++ )
62 scarr[l] += n_dis[k][l] * cpmx1[k][i1];
63 }
64 #if 0 /* �����Ȥ��Ȥ���doublework�Υ������Ȥ�դˤ��� */
65 {
66 double *fpt, **fptpt, *fpt2;
67 int *ipt, **iptpt;
68 fpt2 = match;
69 iptpt = cpmxpdn;
70 fptpt = cpmxpd;
71 while( lgth2-- )
72 {
73 *fpt2 = 0.0;
74 ipt=*iptpt,fpt=*fptpt;
75 while( *ipt > -1 )
76 *fpt2 += scarr[*ipt++] * *fpt++;
77 fpt2++,iptpt++,fptpt++;
78 }
79 }
80 #else
81 for( j=0; j<lgth2; j++ )
82 {
83 match[j] = 0.0;
84 for( k=0; cpmxpdn[k][j]>-1; k++ )
85 match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
86 }
87 #endif
88 }
89 #endif
90
Ltracking(double * lasthorizontalw,double * lastverticalw,char ** seq1,char ** seq2,char ** mseq1,char ** mseq2,int ** ijp,int * off1pt,int * off2pt,int endi,int endj,int * warpis,int * warpjs,int warpbase)91 static double Ltracking( double *lasthorizontalw, double *lastverticalw,
92 char **seq1, char **seq2,
93 char **mseq1, char **mseq2,
94 int **ijp, int *off1pt, int *off2pt, int endi, int endj,
95 int *warpis, int *warpjs, int warpbase )
96 {
97 int i, j, l, iin, jin, lgth1, lgth2, k, limk;
98 int ifi=0, jfi=0; // by D.Mathog, a guess
99 // char gap[] = "-";
100 char *gap;
101 gap = newgapstr;
102 lgth1 = strlen( seq1[0] );
103 lgth2 = strlen( seq2[0] );
104
105 #if 0
106 for( i=0; i<lgth1; i++ )
107 {
108 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
109 }
110 #endif
111
112 for( i=0; i<lgth1+1; i++ )
113 {
114 ijp[i][0] = localstop;
115 }
116 for( j=0; j<lgth2+1; j++ )
117 {
118 ijp[0][j] = localstop;
119 }
120
121 mseq1[0] += lgth1+lgth2;
122 *mseq1[0] = 0;
123 mseq2[0] += lgth1+lgth2;
124 *mseq2[0] = 0;
125 iin = endi; jin = endj;
126 limk = lgth1+lgth2;
127 for( k=0; k<=limk; k++ )
128 {
129 if( ijp[iin][jin] >= warpbase )
130 {
131 // fprintf( stderr, "WARP!\n" );
132 ifi = warpis[ijp[iin][jin]-warpbase];
133 jfi = warpjs[ijp[iin][jin]-warpbase];
134 }
135 else if( ijp[iin][jin] < 0 )
136 {
137 ifi = iin-1; jfi = jin+ijp[iin][jin];
138 }
139 else if( ijp[iin][jin] > 0 )
140 {
141 ifi = iin-ijp[iin][jin]; jfi = jin-1;
142 }
143 else
144 {
145 ifi = iin-1; jfi = jin-1;
146 }
147
148
149 #if 1 // sentou de warp?
150 if( ifi == -warpbase && jfi == -warpbase )
151 {
152 l = iin;
153 while( --l >= 0 )
154 {
155 *--mseq1[0] = seq1[0][l];
156 *--mseq2[0] = *gap;
157 k++;
158 }
159 l= jin;
160 while( --l >= 0 )
161 {
162 *--mseq1[0] = *gap;
163 *--mseq2[0] = seq2[0][l];
164 k++;
165 }
166 break;
167 }
168 else
169 #endif
170 {
171 l = iin - ifi;
172 while( --l > 0 )
173 {
174 *--mseq1[0] = seq1[0][ifi+l];
175 *--mseq2[0] = *gap;
176 k++;
177 }
178 l= jin - jfi;
179 while( --l > 0 )
180 {
181 *--mseq1[0] = *gap;
182 *--mseq2[0] = seq2[0][jfi+l];
183 k++;
184 }
185 }
186
187
188 if( iin <= 0 || jin <= 0 ) break;
189 *--mseq1[0] = seq1[0][ifi];
190 *--mseq2[0] = seq2[0][jfi];
191 if( ijp[ifi][jfi] == localstop ) break;
192 k++;
193 iin = ifi; jin = jfi;
194 }
195 if( ifi == -1 ) *off1pt = 0; else *off1pt = ifi;
196 if( jfi == -1 ) *off2pt = 0; else *off2pt = jfi;
197
198 // fprintf( stderr, "ifn = %d, jfn = %d\n", ifi, jfi );
199 // fprintf( stderr, "\n" );
200 // fprintf( stderr, "%s\n", mseq1[0] );
201 // fprintf( stderr, "%s\n", mseq2[0] );
202
203
204 return( 0.0 );
205 }
206
207
L__align11(double ** n_dynamicmtx,double scoreoffset,char ** seq1,char ** seq2,int alloclen,int * off1pt,int * off2pt)208 double L__align11( double **n_dynamicmtx, double scoreoffset, char **seq1, char **seq2, int alloclen, int *off1pt, int *off2pt )
209 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
210 {
211 // int k;
212 int i, j;
213 int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
214 int lgth1, lgth2;
215 int resultlen;
216 double wm = 0.0; /* int ?????? */
217 double g;
218 double *currentw, *previousw;
219 #if 1
220 double *wtmp;
221 int *ijppt;
222 double *mjpt, *prept, *curpt;
223 int *mpjpt;
224 #endif
225 static TLS double mi, *m;
226 static TLS int **ijp;
227 static TLS int mpi, *mp;
228 static TLS double *w1, *w2;
229 static TLS double *match;
230 static TLS double *initverticalw; /* kufuu sureba iranai */
231 static TLS double *lastverticalw; /* kufuu sureba iranai */
232 static TLS char **mseq1;
233 static TLS char **mseq2;
234 static TLS char **mseq;
235 // static TLS int **intwork;
236 // static TLS double **doublework;
237 static TLS int orlgth1 = 0, orlgth2 = 0;
238 static TLS double **amino_dynamicmtx = NULL; // ??
239 double maxwm;
240 int endali = 0, endalj = 0; // by D.Mathog, a guess
241 // int endali, endalj;
242 double localthr = -offset + scoreoffset * 600; // 2013/12/13
243 double localthr2 = -offset + scoreoffset * 600; // 2013/12/13
244 // double localthr = -offset;
245 // double localthr2 = -offset;
246 double fpenalty = (double)penalty;
247 double fpenalty_ex = (double)penalty_ex;
248 double fpenalty_shift = (double)penalty_shift;
249 double fpenalty_tmp; // atode kesu
250
251 int *warpis = NULL;
252 int *warpjs = NULL;
253 int *warpi = NULL;
254 int *warpj = NULL;
255 int *prevwarpi = NULL;
256 int *prevwarpj = NULL;
257 double *wmrecords = NULL;
258 double *prevwmrecords = NULL;
259 int warpn = 0;
260 int warpbase;
261 double curm = 0.0;
262 double *wmrecordspt, *wmrecords1pt, *prevwmrecordspt;
263 int *warpipt, *warpjpt;
264
265
266
267 if( seq1 == NULL )
268 {
269 if( orlgth1 > 0 && orlgth2 > 0 )
270 {
271 orlgth1 = 0;
272 orlgth2 = 0;
273 free( mseq1 );
274 free( mseq2 );
275 FreeFloatVec( w1 );
276 FreeFloatVec( w2 );
277 FreeFloatVec( match );
278 FreeFloatVec( initverticalw );
279 FreeFloatVec( lastverticalw );
280
281 FreeFloatVec( m );
282 FreeIntVec( mp );
283
284 FreeCharMtx( mseq );
285 if( amino_dynamicmtx ) FreeDoubleMtx( amino_dynamicmtx ); amino_dynamicmtx = NULL;
286
287 }
288 return( 0.0 );
289 }
290
291
292 if( orlgth1 == 0 )
293 {
294 mseq1 = AllocateCharMtx( njob, 0 );
295 mseq2 = AllocateCharMtx( njob, 0 );
296 }
297
298
299 lgth1 = strlen( seq1[0] );
300 lgth2 = strlen( seq2[0] );
301
302
303 warpbase = lgth1 + lgth2;
304 warpis = NULL;
305 warpjs = NULL;
306 warpn = 0;
307 if( trywarp )
308 {
309 wmrecords = AllocateFloatVec( lgth2+1 );
310 warpi = AllocateIntVec( lgth2+1 );
311 warpj = AllocateIntVec( lgth2+1 );
312 prevwmrecords = AllocateFloatVec( lgth2+1 );
313 prevwarpi = AllocateIntVec( lgth2+1 );
314 prevwarpj = AllocateIntVec( lgth2+1 );
315 for( i=0; i<lgth2+1; i++ ) prevwmrecords[i] = 0.0;
316 for( i=0; i<lgth2+1; i++ ) wmrecords[i] = 0.0;
317 for( i=0; i<lgth2+1; i++ ) prevwarpi[i] = -warpbase;
318 for( i=0; i<lgth2+1; i++ ) prevwarpj[i] = -warpbase;
319 for( i=0; i<lgth2+1; i++ ) warpi[i] = -warpbase;
320 for( i=0; i<lgth2+1; i++ ) warpj[i] = -warpbase;
321 }
322
323
324 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
325 {
326 int ll1, ll2;
327
328 if( orlgth1 > 0 && orlgth2 > 0 )
329 {
330 FreeFloatVec( w1 );
331 FreeFloatVec( w2 );
332 FreeFloatVec( match );
333 FreeFloatVec( initverticalw );
334 FreeFloatVec( lastverticalw );
335
336 FreeFloatVec( m );
337 FreeIntVec( mp );
338
339 FreeCharMtx( mseq );
340 if( amino_dynamicmtx ) FreeDoubleMtx( amino_dynamicmtx ); amino_dynamicmtx = NULL;
341
342
343 // FreeFloatMtx( doublework );
344 // FreeIntMtx( intwork );
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
370 #if DEBUG
371 fprintf( stderr, "succeeded\n" );
372 #endif
373 amino_dynamicmtx = AllocateDoubleMtx( 0x80, 0x80 );
374 orlgth1 = ll1 - 100;
375 orlgth2 = ll2 - 100;
376 }
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 match_calc_mtx( amino_dynamicmtx, initverticalw, seq2, seq1, 0, lgth1 );
423
424 match_calc_mtx( amino_dynamicmtx, currentw, seq1, seq2, 0, lgth2 );
425
426
427 lasti = lgth2+1;
428 for( j=1; j<lasti; ++j )
429 {
430 m[j] = currentw[j-1]; mp[j] = 0;
431 #if 0
432 if( m[j] < localthr ) m[j] = localthr2;
433 #endif
434 }
435
436 lastverticalw[0] = currentw[lgth2-1];
437
438 lasti = lgth1+1;
439
440 #if 0
441 fprintf( stderr, "currentw = \n" );
442 for( i=0; i<lgth1+1; i++ )
443 {
444 fprintf( stderr, "%5.2f ", currentw[i] );
445 }
446 fprintf( stderr, "\n" );
447 fprintf( stderr, "initverticalw = \n" );
448 for( i=0; i<lgth2+1; i++ )
449 {
450 fprintf( stderr, "%5.2f ", initverticalw[i] );
451 }
452 fprintf( stderr, "\n" );
453 #endif
454 #if DEBUG2
455 fprintf( stderr, "\n" );
456 fprintf( stderr, " " );
457 for( j=0; j<lgth2; j++ )
458 fprintf( stderr, "%c ", seq2[0][j] );
459 fprintf( stderr, "\n" );
460 #endif
461
462 localstop = lgth1+lgth2+1;
463 maxwm = -999999999.9;
464 #if DEBUG2
465 fprintf( stderr, "\n" );
466 fprintf( stderr, "%c ", seq1[0][0] );
467
468 for( j=0; j<lgth2+1; j++ )
469 fprintf( stderr, "%5.0f ", currentw[j] );
470 fprintf( stderr, "\n" );
471 #endif
472
473 for( i=1; i<lasti; i++ )
474 {
475 wtmp = previousw;
476 previousw = currentw;
477 currentw = wtmp;
478
479 previousw[0] = initverticalw[i-1];
480
481 match_calc_mtx( amino_dynamicmtx, currentw, seq1, seq2, i, lgth2 );
482 #if DEBUG2
483 fprintf( stderr, "%c ", seq1[0][i] );
484 fprintf( stderr, "%5.0f ", currentw[0] );
485 #endif
486
487 #if XXXXXXX
488 fprintf( stderr, "\n" );
489 fprintf( stderr, "i=%d\n", i );
490 fprintf( stderr, "currentw = \n" );
491 for( j=0; j<lgth2; j++ )
492 {
493 fprintf( stderr, "%5.2f ", currentw[j] );
494 }
495 fprintf( stderr, "\n" );
496 #endif
497 #if XXXXXXX
498 fprintf( stderr, "\n" );
499 fprintf( stderr, "i=%d\n", i );
500 fprintf( stderr, "currentw = \n" );
501 for( j=0; j<lgth2; j++ )
502 {
503 fprintf( stderr, "%5.2f ", currentw[j] );
504 }
505 fprintf( stderr, "\n" );
506 #endif
507 currentw[0] = initverticalw[i];
508
509 mi = previousw[0]; mpi = 0;
510
511 #if 0
512 if( mi < localthr ) mi = localthr2;
513 #endif
514
515 ijppt = ijp[i] + 1;
516 mjpt = m + 1;
517 prept = previousw;
518 curpt = currentw + 1;
519 mpjpt = mp + 1;
520 lastj = lgth2+1;
521
522 if( trywarp )
523 {
524 prevwmrecordspt = prevwmrecords;
525 wmrecordspt = wmrecords+1;
526 wmrecords1pt = wmrecords;
527 warpipt = warpi + 1;
528 warpjpt = warpj + 1;
529 }
530 for( j=1; j<lastj; j++ )
531 {
532 wm = *prept;
533 *ijppt = 0;
534
535 #if 0
536 fprintf( stderr, "%5.0f->", wm );
537 #endif
538 #if 0
539 fprintf( stderr, "%5.0f?", g );
540 #endif
541 if( (g=mi+fpenalty) > wm )
542 {
543 wm = g;
544 *ijppt = -( j - mpi );
545 }
546 if( *prept > mi )
547 {
548 mi = *prept;
549 mpi = j-1;
550 }
551
552 #if USE_PENALTY_EX
553 mi += fpenalty_ex;
554 #endif
555
556 #if 0
557 fprintf( stderr, "%5.0f?", g );
558 #endif
559 if( (g=*mjpt+fpenalty) > wm )
560 {
561 wm = g;
562 *ijppt = +( i - *mpjpt );
563 }
564 if( *prept > *mjpt )
565 {
566 *mjpt = *prept;
567 *mpjpt = i-1;
568 }
569 #if USE_PENALTY_EX
570 *mjpt += fpenalty_ex;
571 #endif
572
573 if( maxwm < wm )
574 {
575 maxwm = wm;
576 endali = i;
577 endalj = j;
578 }
579 #if 1
580 if( wm < localthr )
581 {
582 // fprintf( stderr, "stop i=%d, j=%d, curpt=%f, localthr = %f\n", i, j, *curpt, localthr );
583 *ijppt = localstop;
584 wm = localthr2;
585 }
586 #endif
587 #if 0
588 fprintf( stderr, "%5.0f ", *curpt );
589 #endif
590 #if 0
591 fprintf( stderr, "wm (%d,%d) = %5.0f\n", i, j, wm );
592 // fprintf( stderr, "%c-%c *ijppt = %d, localstop = %d\n", seq1[0][i], seq2[0][j], *ijppt, localstop );
593 #endif
594 if( trywarp )
595 {
596 fpenalty_tmp = fpenalty_shift + fpenalty_ex * ( i - prevwarpi[j-1] + j - prevwarpj[j-1] );
597 // fprintf( stderr, "fpenalty_shift = %f\n", fpenalty_tmp );
598
599 // 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 );
600 // if( (g=prevwmrecords[j-1] + fpenalty_shift )> wm )
601 if( ( g=*prevwmrecordspt++ + fpenalty_tmp )> wm ) // naka ha osokute kamawanai
602 {
603 // 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 );
604 if( warpn && prevwarpi[j-1] == warpis[warpn-1] && prevwarpj[j-1] == warpjs[warpn-1] )
605 {
606 *ijppt = warpbase + warpn - 1;
607 }
608 else
609 {
610 *ijppt = warpbase + warpn;
611 warpis = realloc( warpis, sizeof(int) * ( warpn+1 ) );
612 warpjs = realloc( warpjs, sizeof(int) * ( warpn+1 ) );
613 warpis[warpn] = prevwarpi[j-1];
614 warpjs[warpn] = prevwarpj[j-1];
615 warpn++;
616 }
617 wm = g;
618 }
619 else
620 {
621 }
622
623 curm = *curpt + wm;
624
625 // fprintf( stderr, "###### curm = %f at %c-%c, i=%d, j=%d\n", curm, seq1[0][i], seq2[0][j], i, j );
626
627 // fprintf( stderr, "copy from i, j-1? %f > %f?\n", wmrecords[j-1], curm );
628 // if( wmrecords[j-1] > wmrecords[j] )
629 if( *wmrecords1pt > *wmrecordspt )
630 {
631 // fprintf( stderr, "yes\n" );
632 // wmrecords[j] = wmrecords[j-1];
633 *wmrecordspt = *wmrecords1pt;
634 // warpi[j] = warpi[j-1];
635 // warpj[j] = warpj[j-1];
636 *warpipt = *(warpipt-1);
637 *warpjpt = *(warpjpt-1);
638 // fprintf( stderr, "warpi[j]=%d, warpj[j]=%d wmrecords[j] = %f\n", warpi[j], warpj[j], wmrecords[j] );
639 }
640 // else
641 // {
642 // fprintf( stderr, "no\n" );
643 // }
644
645 // fprintf( stderr, " curm = %f at %c-%c\n", curm, seq1[0][i], seq2[0][j] );
646 // fprintf( stderr, " wmrecords[%d] = %f\n", j, wmrecords[j] );
647 // fprintf( stderr, "replace?\n" );
648
649 // if( curm > wmrecords[j] )
650 if( curm > *wmrecordspt )
651 {
652 // 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 );
653 // wmrecords[j] = curm;
654 *wmrecordspt = curm;
655 // warpi[j] = i;
656 // warpj[j] = j;
657 *warpipt = i;
658 *warpjpt = j;
659 }
660 // else
661 // {
662 // fprintf( stderr, "No! warpi[j]=%d, warpj[j]=%d wmrecords[j] = %f\n", warpi[j], warpj[j], wmrecords[j] );
663 // }
664 // fprintf( stderr, "%d-%d (%c-%c) curm = %5.0f, wmrecords[j]=%f\n", i, j, seq1[0][i], seq2[0][j], curm, wmrecords[j] );
665 wmrecordspt++;
666 wmrecords1pt++;
667 warpipt++;
668 warpjpt++;
669 }
670
671 *curpt++ += wm;
672 ijppt++;
673 mjpt++;
674 prept++;
675 mpjpt++;
676 }
677 #if DEBUG2
678 fprintf( stderr, "\n" );
679 #endif
680
681 lastverticalw[i] = currentw[lgth2-1];
682 if( trywarp )
683 {
684 fltncpy( prevwmrecords, wmrecords, lastj );
685 intncpy( prevwarpi, warpi, lastj );
686 intncpy( prevwarpj, warpj, lastj );
687 }
688
689 }
690 // fprintf( stderr, "\nwm = %f\n", wm );
691 if( trywarp )
692 {
693 // if( warpn ) fprintf( stderr, "warpn = %d\n", warpn );
694 free( wmrecords );
695 free( prevwmrecords );
696 free( warpi );
697 free( warpj );
698 free( prevwarpi );
699 free( prevwarpj );
700 }
701
702 #if 0
703 fprintf( stderr, "maxwm = %f\n", maxwm );
704 fprintf( stderr, "endali = %d\n", endali );
705 fprintf( stderr, "endalj = %d\n", endalj );
706 #endif
707
708 if( ijp[endali][endalj] == localstop )
709 {
710 strcpy( seq1[0], "" );
711 strcpy( seq2[0], "" );
712 *off1pt = *off2pt = 0;
713 fprintf( stderr, "maxwm <- 0.0 \n" );
714 return( 0.0 );
715 }
716
717 Ltracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, off1pt, off2pt, endali, endalj, warpis, warpjs, warpbase );
718 if( warpis ) free( warpis );
719 if( warpjs ) free( warpjs );
720
721
722 resultlen = strlen( mseq1[0] );
723 if( alloclen < resultlen || resultlen > N )
724 {
725 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
726 ErrorExit( "LENGTH OVER!\n" );
727 }
728
729
730 strcpy( seq1[0], mseq1[0] );
731 strcpy( seq2[0], mseq2[0] );
732
733 #if 0
734 fprintf( stderr, "wm=%f\n", wm );
735 fprintf( stderr, ">\n%s\n", mseq1[0] );
736 fprintf( stderr, ">\n%s\n", mseq2[0] );
737
738 fprintf( stderr, "*off1pt = %d, *off2pt = %d\n", *off1pt, *off2pt );
739
740 fprintf( stderr, "maxwm = %f\n", maxwm );
741 fprintf( stderr, " wm = %f\n", wm );
742 #endif
743
744 return( maxwm );
745 }
746
747
L__align11_noalign(double ** n_dynamicmtx,char ** seq1,char ** seq2)748 double L__align11_noalign( double **n_dynamicmtx, char **seq1, char **seq2 )
749 // warp mitaiou
750 {
751 // int k;
752 int i, j;
753 int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
754 int lgth1, lgth2;
755 // int resultlen;
756 double wm = 0.0; /* int ?????? */
757 double g;
758 double *currentw, *previousw;
759 #if 1
760 double *wtmp;
761 // int *ijppt;
762 double *mjpt, *prept, *curpt;
763 // int *mpjpt;
764 #endif
765 static TLS double mi, *m;
766 // static TLS int **ijp;
767 // static TLS int mpi, *mp;
768 static TLS double *w1, *w2;
769 static TLS double *match;
770 static TLS double *initverticalw; /* kufuu sureba iranai */
771 static TLS double *lastverticalw; /* kufuu sureba iranai */
772 // static TLS char **mseq1;
773 // static TLS char **mseq2;
774 // static TLS char **mseq;
775 // static TLS int **intwork;
776 // static TLS double **doublework;
777 static TLS int orlgth1 = 0, orlgth2 = 0;
778 static TLS double **amino_dynamicmtx = NULL; // ??
779 double maxwm;
780 // int endali = 0, endalj = 0; // by D.Mathog, a guess
781 // int endali, endalj;
782 double localthr = -offset;
783 double localthr2 = -offset;
784 // double localthr = 100;
785 // double localthr2 = 100;
786 double fpenalty = (double)penalty;
787 double fpenalty_ex = (double)penalty_ex;
788
789 if( seq1 == NULL )
790 {
791 if( orlgth1 > 0 && orlgth2 > 0 )
792 {
793 orlgth1 = 0;
794 orlgth2 = 0;
795 // free( mseq1 );
796 // free( mseq2 );
797 FreeFloatVec( w1 );
798 FreeFloatVec( w2 );
799 FreeFloatVec( match );
800 FreeFloatVec( initverticalw );
801 FreeFloatVec( lastverticalw );
802
803 FreeFloatVec( m );
804 // FreeIntVec( mp );
805
806 // FreeCharMtx( mseq );
807 if( amino_dynamicmtx ) FreeDoubleMtx( amino_dynamicmtx ); amino_dynamicmtx = NULL;
808
809 }
810 return( 0.0 );
811 }
812
813
814 // if( orlgth1 == 0 )
815 // {
816 // mseq1 = AllocateCharMtx( njob, 0 );
817 // mseq2 = AllocateCharMtx( njob, 0 );
818 // }
819
820
821 lgth1 = strlen( seq1[0] );
822 lgth2 = strlen( seq2[0] );
823
824 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
825 {
826 int ll1, ll2;
827
828 if( orlgth1 > 0 && orlgth2 > 0 )
829 {
830 FreeFloatVec( w1 );
831 FreeFloatVec( w2 );
832 FreeFloatVec( match );
833 FreeFloatVec( initverticalw );
834 FreeFloatVec( lastverticalw );
835
836 FreeFloatVec( m );
837 // FreeIntVec( mp );
838
839 // FreeCharMtx( mseq );
840
841
842
843 // FreeFloatMtx( doublework );
844 // FreeIntMtx( intwork );
845 if( amino_dynamicmtx ) FreeDoubleMtx( amino_dynamicmtx ); amino_dynamicmtx = NULL;
846 }
847
848 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
849 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
850
851 #if DEBUG
852 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
853 #endif
854
855 w1 = AllocateFloatVec( ll2+2 );
856 w2 = AllocateFloatVec( ll2+2 );
857 match = AllocateFloatVec( ll2+2 );
858
859 initverticalw = AllocateFloatVec( ll1+2 );
860 lastverticalw = AllocateFloatVec( ll1+2 );
861
862 m = AllocateFloatVec( ll2+2 );
863 // mp = AllocateIntVec( ll2+2 );
864
865 // mseq = AllocateCharMtx( njob, ll1+ll2 );
866
867
868 // doublework = AllocateFloatMtx( nalphabets, MAX( ll1, ll2 )+2 );
869 // intwork = AllocateIntMtx( nalphabets, MAX( ll1, ll2 )+2 );
870
871 #if DEBUG
872 fprintf( stderr, "succeeded\n" );
873 #endif
874 amino_dynamicmtx = AllocateDoubleMtx( 0x80, 0x80 );
875 orlgth1 = ll1 - 100;
876 orlgth2 = ll2 - 100;
877 }
878
879 for( i=0; i<nalphabets; i++) for( j=0; j<nalphabets; j++ )
880 amino_dynamicmtx[(int)amino[i]][(int)amino[j]] = (double)n_dynamicmtx[i][j];
881
882
883
884 // mseq1[0] = mseq[0];
885 // mseq2[0] = mseq[1];
886
887
888 // if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
889 // {
890 // int ll1, ll2;
891 //
892 // if( commonAlloc1 && commonAlloc2 )
893 // {
894 // FreeIntMtx( commonIP );
895 // }
896 //
897 // ll1 = MAX( orlgth1, commonAlloc1 );
898 // ll2 = MAX( orlgth2, commonAlloc2 );
899
900 #if DEBUG
901 // fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
902 #endif
903
904 // commonIP = AllocateIntMtx( ll1+10, ll2+10 );
905
906 #if DEBUG
907 // fprintf( stderr, "succeeded\n\n" );
908 #endif
909
910 // commonAlloc1 = ll1;
911 // commonAlloc2 = ll2;
912 // }
913 // ijp = commonIP;
914
915
916 #if 0
917 for( i=0; i<lgth1; i++ )
918 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
919 #endif
920
921 currentw = w1;
922 previousw = w2;
923
924 match_calc_mtx( amino_dynamicmtx, initverticalw, seq2, seq1, 0, lgth1 );
925
926 match_calc_mtx( amino_dynamicmtx, currentw, seq1, seq2, 0, lgth2 );
927
928
929 lasti = lgth2+1;
930 for( j=1; j<lasti; ++j )
931 {
932 m[j] = currentw[j-1];
933 // mp[j] = 0;
934 #if 0
935 if( m[j] < localthr ) m[j] = localthr2;
936 #endif
937 }
938
939 lastverticalw[0] = currentw[lgth2-1];
940
941 lasti = lgth1+1;
942
943 #if 0
944 fprintf( stderr, "currentw = \n" );
945 for( i=0; i<lgth1+1; i++ )
946 {
947 fprintf( stderr, "%5.2f ", currentw[i] );
948 }
949 fprintf( stderr, "\n" );
950 fprintf( stderr, "initverticalw = \n" );
951 for( i=0; i<lgth2+1; i++ )
952 {
953 fprintf( stderr, "%5.2f ", initverticalw[i] );
954 }
955 fprintf( stderr, "\n" );
956 #endif
957 #if DEBUG2
958 fprintf( stderr, "\n" );
959 fprintf( stderr, " " );
960 for( j=0; j<lgth2; j++ )
961 fprintf( stderr, "%c ", seq2[0][j] );
962 fprintf( stderr, "\n" );
963 #endif
964
965 localstop = lgth1+lgth2+1;
966 maxwm = -999999999.9;
967 #if DEBUG2
968 fprintf( stderr, "\n" );
969 fprintf( stderr, "%c ", seq1[0][0] );
970
971 for( j=0; j<lgth2+1; j++ )
972 fprintf( stderr, "%5.0f ", currentw[j] );
973 fprintf( stderr, "\n" );
974 #endif
975
976 for( i=1; i<lasti; i++ )
977 {
978 wtmp = previousw;
979 previousw = currentw;
980 currentw = wtmp;
981
982 previousw[0] = initverticalw[i-1];
983
984 match_calc_mtx( amino_dynamicmtx, currentw, seq1, seq2, i, lgth2 );
985 #if DEBUG2
986 fprintf( stderr, "%c ", seq1[0][i] );
987 fprintf( stderr, "%5.0f ", currentw[0] );
988 #endif
989
990 #if XXXXXXX
991 fprintf( stderr, "\n" );
992 fprintf( stderr, "i=%d\n", i );
993 fprintf( stderr, "currentw = \n" );
994 for( j=0; j<lgth2; j++ )
995 {
996 fprintf( stderr, "%5.2f ", currentw[j] );
997 }
998 fprintf( stderr, "\n" );
999 #endif
1000 #if XXXXXXX
1001 fprintf( stderr, "\n" );
1002 fprintf( stderr, "i=%d\n", i );
1003 fprintf( stderr, "currentw = \n" );
1004 for( j=0; j<lgth2; j++ )
1005 {
1006 fprintf( stderr, "%5.2f ", currentw[j] );
1007 }
1008 fprintf( stderr, "\n" );
1009 #endif
1010 currentw[0] = initverticalw[i];
1011
1012 mi = previousw[0];
1013 // mpi = 0;
1014
1015 #if 0
1016 if( mi < localthr ) mi = localthr2;
1017 #endif
1018
1019 // ijppt = ijp[i] + 1;
1020 mjpt = m + 1;
1021 prept = previousw;
1022 curpt = currentw + 1;
1023 // mpjpt = mp + 1;
1024 lastj = lgth2+1;
1025 for( j=1; j<lastj; j++ )
1026 {
1027 wm = *prept;
1028 // *ijppt = 0;
1029
1030 #if 0
1031 fprintf( stderr, "%5.0f->", wm );
1032 #endif
1033 #if 0
1034 fprintf( stderr, "%5.0f?", g );
1035 #endif
1036 if( (g=mi+fpenalty) > wm )
1037 {
1038 wm = g;
1039 // *ijppt = -( j - mpi );
1040 }
1041 if( *prept > mi )
1042 {
1043 mi = *prept;
1044 // mpi = j-1;
1045 }
1046
1047 #if USE_PENALTY_EX
1048 mi += fpenalty_ex;
1049 #endif
1050
1051 #if 0
1052 fprintf( stderr, "%5.0f?", g );
1053 #endif
1054 if( (g=*mjpt+fpenalty) > wm )
1055 {
1056 wm = g;
1057 // *ijppt = +( i - *mpjpt );
1058 }
1059 if( *prept > *mjpt )
1060 {
1061 *mjpt = *prept;
1062 // *mpjpt = i-1;
1063 }
1064 #if USE_PENALTY_EX
1065 *mjpt += fpenalty_ex;
1066 #endif
1067
1068 if( maxwm < wm )
1069 {
1070 maxwm = wm;
1071 // endali = i;
1072 // endalj = j;
1073 }
1074 #if 1
1075 if( wm < localthr )
1076 {
1077 // fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
1078 // *ijppt = localstop;
1079 wm = localthr2;
1080 }
1081 #endif
1082 #if 0
1083 fprintf( stderr, "%5.0f ", *curpt );
1084 #endif
1085 #if DEBUG2
1086 fprintf( stderr, "%5.0f ", wm );
1087 // fprintf( stderr, "%c-%c *ijppt = %d, localstop = %d\n", seq1[0][i], seq2[0][j], *ijppt, localstop );
1088 #endif
1089
1090 *curpt++ += wm;
1091 // ijppt++;
1092 mjpt++;
1093 prept++;
1094 // mpjpt++;
1095 }
1096 #if DEBUG2
1097 fprintf( stderr, "\n" );
1098 #endif
1099
1100 lastverticalw[i] = currentw[lgth2-1];
1101 }
1102
1103
1104 #if 0
1105 fprintf( stderr, "maxwm = %f\n", maxwm );
1106 fprintf( stderr, "endali = %d\n", endali );
1107 fprintf( stderr, "endalj = %d\n", endalj );
1108 #endif
1109
1110
1111 #if 0 // IRUKAMO!!!!
1112 if( ijp[endali][endalj] == localstop )
1113 {
1114 strcpy( seq1[0], "" );
1115 strcpy( seq2[0], "" );
1116 *off1pt = *off2pt = 0;
1117 fprintf( stderr, "maxwm <- 0.0 \n" );
1118 return( 0.0 );
1119 }
1120 #else
1121 if( maxwm < localthr )
1122 {
1123 fprintf( stderr, "maxwm <- 0.0 \n" );
1124 return( 0.0 );
1125 }
1126 #endif
1127
1128 // Ltracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, off1pt, off2pt, endali, endalj );
1129
1130
1131 // resultlen = strlen( mseq1[0] );
1132 // if( alloclen < resultlen || resultlen > N )
1133 // {
1134 // fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1135 // ErrorExit( "LENGTH OVER!\n" );
1136 // }
1137
1138
1139 // strcpy( seq1[0], mseq1[0] );
1140 // strcpy( seq2[0], mseq2[0] );
1141
1142 #if 0
1143 fprintf( stderr, "wm=%f\n", wm );
1144 fprintf( stderr, ">\n%s\n", mseq1[0] );
1145 fprintf( stderr, ">\n%s\n", mseq2[0] );
1146
1147 fprintf( stderr, "maxwm = %f\n", maxwm );
1148 fprintf( stderr, " wm = %f\n", wm );
1149 #endif
1150
1151 return( maxwm );
1152 }
1153