1 #include "mltaln.h"
2 #include "dp.h"
3
4 #define MEMSAVE 1
5
6 #define DEBUG 0
7 #define USE_PENALTY_EX 0
8 #define STOREWM 1
9
10 #define DPTANNI 10
11
12 #define LOCAL 0
13
14 static int reccycle = 0;
15
16 static double localthr;
17
match_ribosum(double * match,double ** cpmx1,double ** cpmx2,int i1,int lgth2,double ** doublework,int ** intwork,int initialize)18 static void match_ribosum( double *match, double **cpmx1, double **cpmx2, int i1, int lgth2, double **doublework, int **intwork, int initialize )
19 {
20 int j, k, l;
21 double scarr[38];
22 double **cpmxpd = doublework;
23 int **cpmxpdn = intwork;
24 int count = 0;
25 double *matchpt;
26 double **cpmxpdpt;
27 int **cpmxpdnpt;
28 int cpkd;
29
30 if( initialize )
31 {
32 for( j=0; j<lgth2; j++ )
33 {
34 count = 0;
35 for( l=0; l<37; l++ )
36 {
37 if( cpmx2[j][l] )
38 {
39 cpmxpd[j][count] = cpmx2[j][l];
40 cpmxpdn[j][count] = l;
41 count++;
42 }
43 }
44 cpmxpdn[j][count] = -1;
45 }
46 }
47
48 for( l=0; l<37; l++ )
49 {
50 scarr[l] = 0.0;
51 for( k=0; k<37; k++ )
52 {
53 scarr[l] += ribosumdis[k][l] * cpmx1[i1][k];
54 }
55 }
56 #if 0 /* �����Ȥ��Ȥ���doublework�Υ������Ȥ�դˤ��� */
57 {
58 double *fpt, **fptpt, *fpt2;
59 int *ipt, **iptpt;
60 fpt2 = match;
61 iptpt = cpmxpdn;
62 fptpt = cpmxpd;
63 while( lgth2-- )
64 {
65 *fpt2 = 0.0;
66 ipt=*iptpt,fpt=*fptpt;
67 while( *ipt > -1 )
68 *fpt2 += scarr[*ipt++] * *fpt++;
69 fpt2++,iptpt++,fptpt++;
70 }
71 }
72 for( j=0; j<lgth2; j++ )
73 {
74 match[j] = 0.0;
75 for( k=0; cpmxpdn[j][k]>-1; k++ )
76 match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
77 }
78 #else
79 matchpt = match;
80 cpmxpdnpt = cpmxpdn;
81 cpmxpdpt = cpmxpd;
82 while( lgth2-- )
83 {
84 *matchpt = 0.0;
85 for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
86 *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
87 matchpt++;
88 cpmxpdnpt++;
89 cpmxpdpt++;
90 }
91 #endif
92 }
93
match_calc(double * match,double ** cpmx1,double ** cpmx2,int i1,int lgth2,double ** doublework,int ** intwork,int initialize)94 static void match_calc( double *match, double **cpmx1, double **cpmx2, int i1, int lgth2, double **doublework, int **intwork, int initialize )
95 {
96 int j, k, l;
97 // double scarr[26];
98 double **cpmxpd = doublework;
99 int **cpmxpdn = intwork;
100 int count = 0;
101 double *matchpt;
102 double **cpmxpdpt;
103 int **cpmxpdnpt;
104 int cpkd;
105 double *scarr;
106 scarr = calloc( nalphabets, sizeof( double ) );
107
108 if( initialize )
109 {
110 for( j=0; j<lgth2; j++ )
111 {
112 count = 0;
113 for( l=0; l<nalphabets; l++ )
114 {
115 if( cpmx2[j][l] )
116 {
117 cpmxpd[j][count] = cpmx2[j][l];
118 cpmxpdn[j][count] = l;
119 count++;
120 }
121 }
122 cpmxpdn[j][count] = -1;
123 }
124 }
125
126 for( l=0; l<nalphabets; l++ )
127 {
128 scarr[l] = 0.0;
129 for( k=0; k<nalphabets; k++ )
130 {
131 scarr[l] += (n_dis[k][l]-RNAthr) * cpmx1[i1][k];
132 }
133 }
134 #if 0 /* �����Ȥ��Ȥ���doublework�Υ������Ȥ�դˤ��� */
135 {
136 double *fpt, **fptpt, *fpt2;
137 int *ipt, **iptpt;
138 fpt2 = match;
139 iptpt = cpmxpdn;
140 fptpt = cpmxpd;
141 while( lgth2-- )
142 {
143 *fpt2 = 0.0;
144 ipt=*iptpt,fpt=*fptpt;
145 while( *ipt > -1 )
146 *fpt2 += scarr[*ipt++] * *fpt++;
147 fpt2++,iptpt++,fptpt++;
148 }
149 }
150 for( j=0; j<lgth2; j++ )
151 {
152 match[j] = 0.0;
153 for( k=0; cpmxpdn[j][k]>-1; k++ )
154 match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
155 }
156 #else
157 matchpt = match;
158 cpmxpdnpt = cpmxpdn;
159 cpmxpdpt = cpmxpd;
160 while( lgth2-- )
161 {
162 *matchpt = 0.0;
163 for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
164 *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
165 matchpt++;
166 cpmxpdnpt++;
167 cpmxpdpt++;
168 }
169 #endif
170 free( scarr );
171 }
172
173 #if 0
174 static void match_add( double *match, double **cpmx1, double **cpmx2, int i1, int lgth2, double **doublework, int **intwork, int initialize )
175 {
176 int j, k, l;
177 double scarr[nalphabets];
178 double **cpmxpd = doublework;
179 int **cpmxpdn = intwork;
180 int count = 0;
181 double *matchpt;
182 double **cpmxpdpt;
183 int **cpmxpdnpt;
184 int cpkd;
185
186
187 if( initialize )
188 {
189 for( j=0; j<lgth2; j++ )
190 {
191 count = 0;
192 for( l=0; l<nalphabets; l++ )
193 {
194 if( cpmx2[j][l] )
195 {
196 cpmxpd[j][count] = cpmx2[j][l];
197 cpmxpdn[j][count] = l;
198 count++;
199 }
200 }
201 cpmxpdn[j][count] = -1;
202 }
203 }
204
205 for( l=0; l<nalphabets; l++ )
206 {
207 scarr[l] = 0.0;
208 for( k=0; k<nalphabets; k++ )
209 {
210 scarr[l] += n_dis[k][l] * cpmx1[i1][k];
211 }
212 }
213 #if 0 /* �����Ȥ��Ȥ���doublework�Υ������Ȥ�դˤ��� */
214 {
215 double *fpt, **fptpt, *fpt2;
216 int *ipt, **iptpt;
217 fpt2 = match;
218 iptpt = cpmxpdn;
219 fptpt = cpmxpd;
220 while( lgth2-- )
221 {
222 *fpt2 = 0.0;
223 ipt=*iptpt,fpt=*fptpt;
224 while( *ipt > -1 )
225 *fpt2 += scarr[*ipt++] * *fpt++;
226 fpt2++,iptpt++,fptpt++;
227 }
228 }
229 for( j=0; j<lgth2; j++ )
230 {
231 match[j] = 0.0;
232 for( k=0; cpmxpdn[j][k]>-1; k++ )
233 match[j] += scarr[cpmxpdn[j][k]] * cpmxpd[j][k];
234 }
235 #else
236 matchpt = match;
237 cpmxpdnpt = cpmxpdn;
238 cpmxpdpt = cpmxpd;
239 while( lgth2-- )
240 {
241 // *matchpt = 0.0; // add dakara
242 for( k=0; (cpkd=(*cpmxpdnpt)[k])>-1; k++ )
243 *matchpt += scarr[cpkd] * (*cpmxpdpt)[k];
244 matchpt++;
245 cpmxpdnpt++;
246 cpmxpdpt++;
247 }
248 #endif
249 }
250 #endif
251
252 #if 0
253 static double Atracking(
254 char **seq1, char **seq2,
255 char **mseq1, char **mseq2,
256 int **ijp, int icyc, int jcyc,
257 int ist, int ien, int jst, int jen )
258 {
259 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, klim;
260 char *gaptable1, *gt1bk;
261 char *gaptable2, *gt2bk;
262 lgth1 = ien-ist+1;
263 lgth2 = jen-jst+1;
264
265 gt1bk = AllocateCharVec( lgth1+lgth2+1 );
266 gt2bk = AllocateCharVec( lgth1+lgth2+1 );
267
268 #if 0
269 for( i=0; i<lgth1; i++ )
270 {
271 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
272 }
273 #endif
274
275
276 // fprintf( stderr, "in Atracking, lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
277
278 for( i=0; i<lgth1+1; i++ )
279 {
280 ijp[i][0] = i + 1;
281 }
282 for( j=0; j<lgth2+1; j++ )
283 {
284 ijp[0][j] = -( j + 1 );
285 }
286
287
288 gaptable1 = gt1bk + lgth1+lgth2;
289 *gaptable1 = 0;
290 gaptable2 = gt2bk + lgth1+lgth2;
291 *gaptable2 = 0;
292
293 // if( lgth2 == 1 ) fprintf( stderr, "in Atracking, mseq1 = %s, mseq2 = %s\n", mseq1[0], mseq2[0] );
294
295 iin = lgth1; jin = lgth2;
296 klim = lgth1+lgth2;
297 for( k=0; k<=klim; k++ )
298 {
299 if( ijp[iin][jin] < 0 )
300 {
301 ifi = iin-1; jfi = jin+ijp[iin][jin];
302 }
303 else if( ijp[iin][jin] > 0 )
304 {
305 ifi = iin-ijp[iin][jin]; jfi = jin-1;
306 }
307 else
308 {
309 ifi = iin-1; jfi = jin-1;
310 }
311 l = iin - ifi;
312 while( --l )
313 {
314 *--gaptable1 = 'o';
315 *--gaptable2 = '-';
316 k++;
317 }
318 l= jin - jfi;
319 while( --l )
320 {
321 *--gaptable1 = '-';
322 *--gaptable2 = 'o';
323 k++;
324 }
325 if( iin <= 0 || jin <= 0 ) break;
326 *--gaptable1 = 'o';
327 *--gaptable2 = 'o';
328 k++;
329 iin = ifi; jin = jfi;
330
331 }
332
333 for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i]+ist, gaptable1 );
334 for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j]+jst, gaptable2 );
335
336 free( gt1bk );
337 free( gt2bk );
338
339 // fprintf( stderr, "in Atracking (owari), mseq1 = %s\n", mseq1[0] );
340 // fprintf( stderr, "in Atracking (owari), mseq2 = %s\n", mseq2[0] );
341 return( 0.0 );
342 }
343 #endif
344
345
MSalign2m2m_rec(int icyc,int jcyc,double * eff1,double * eff2,char ** seq1,char ** seq2,double ** cpmx1,double ** cpmx2,int ist,int ien,int jst,int jen,int alloclen,char ** mseq1,char ** mseq2,int depth,double ** gapinfo,double ** map)346 static double MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, char **seq1, char **seq2, double **cpmx1, double **cpmx2, int ist, int ien, int jst, int jen, int alloclen, char **mseq1, char **mseq2, int depth, double **gapinfo, double **map )
347 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
348 {
349 double value = 0.0;
350 register int i, j;
351 char **aseq1, **aseq2;
352 int ll1, ll2;
353 int lasti, lastj, imid, jmid = 0;
354 double wm = 0.0; /* int ?????? */
355 double g;
356 double *currentw, *previousw;
357 #if USE_PENALTY_EX
358 double fpenalty_ex = (double)penalty_ex;
359 #endif
360 // double fpenalty = (double)penalty;
361 double *wtmp;
362 // short *ijppt;
363 int *mpjpt;
364 // short **ijp;
365 int *mp;
366 int mpi;
367 double *mjpt, *prept, *curpt;
368 double mi;
369 double *m;
370 double *w1, *w2;
371 // double *match;
372 double *initverticalw; /* kufuu sureba iranai */
373 double *lastverticalw; /* kufuu sureba iranai */
374 int **intwork;
375 double **doublework;
376 // short **shortmtx;
377 #if STOREWM
378 double **WMMTX;
379 double **WMMTX2;
380 #endif
381 double *midw;
382 double *midm;
383 double *midn;
384 int lgth1, lgth2;
385 double maxwm = 0.0;
386 int *jumpforwi;
387 int *jumpforwj;
388 int *jumpbacki;
389 int *jumpbackj;
390 int *jumpdummi; //muda
391 int *jumpdummj; //muda
392 int jumpi, jumpj = 0;
393 char *gaps;
394 int ijpi, ijpj;
395 double *ogcp1;
396 double *fgcp1;
397 double *ogcp2;
398 double *fgcp2;
399 double firstm;
400 int firstmp;
401 #if 0
402 static char ttt1[50000];
403 static char ttt2[50000];
404 #endif
405
406 localthr = -offset + 500; // 0?
407
408 ogcp1 = gapinfo[0] + ist;
409 fgcp1 = gapinfo[1] + ist;
410 ogcp2 = gapinfo[2] + jst;
411 fgcp2 = gapinfo[3] + jst;
412
413 depth++;
414 reccycle++;
415
416 lgth1 = ien-ist+1;
417 lgth2 = jen-jst+1;
418
419 // if( lgth1 < 5 )
420 // fprintf( stderr, "\nWARNING: lgth1 = %d\n", lgth1 );
421 // if( lgth2 < 5 )
422 // fprintf( stderr, "\nWARNING: lgth2 = %d\n", lgth2 );
423 //
424
425 #if 0
426 fprintf( stderr, "==== MSalign (depth=%d, reccycle=%d), ist=%d, ien=%d, jst=%d, jen=%d\n", depth, reccycle, ist, ien, jst, jen );
427 strncpy( ttt1, seq1[0]+ist, lgth1 );
428 strncpy( ttt2, seq2[0]+jst, lgth2 );
429 ttt1[lgth1] = 0;
430 ttt2[lgth2] = 0;
431 fprintf( stderr, "seq1 = %s\n", ttt1 );
432 fprintf( stderr, "seq2 = %s\n", ttt2 );
433 #endif
434 if( lgth2 <= 0 ) // lgth1 <= 0 ha?
435 {
436 // fprintf( stderr, "\n\n==== jimei\n\n" );
437 // exit( 1 );
438 for( i=0; i<icyc; i++ )
439 {
440 strncpy( mseq1[i], seq1[i]+ist, lgth1 );
441 mseq1[i][lgth1] = 0;
442 }
443 for( i=0; i<jcyc; i++ )
444 {
445 mseq2[i][0] = 0;
446 for( j=0; j<lgth1; j++ )
447 strcat( mseq2[i], "-" );
448 }
449
450 // fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
451 // fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
452
453 return( 0.0 );
454 }
455
456 #if MEMSAVE
457 aseq1 = AllocateCharMtx( icyc, 0 );
458 aseq2 = AllocateCharMtx( jcyc, 0 );
459 for( i=0; i<icyc; i++ ) aseq1[i] = mseq1[i];
460 for( i=0; i<jcyc; i++ ) aseq2[i] = mseq2[i];
461 #else
462 aseq1 = AllocateCharMtx( icyc, lgth1+lgth2+100 );
463 aseq2 = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
464 #endif
465
466 // if( lgth1 < DPTANNI && lgth2 < DPTANNI ) // & dato lgth ==1 no kanousei ga arunode yokunai
467 // if( lgth1 < DPTANNI ) // kore mo lgth2 ga mijikasugiru kanousei ari
468 if( lgth1 < DPTANNI || lgth2 < DPTANNI ) // zettai ni anzen ka?
469 {
470 // fprintf( stderr, "==== Going to _tanni\n" );
471
472 // value = MSalignmm_tanni( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ien, jst, jen, alloclen, aseq1, aseq2, gapinfo );
473
474
475 #if MEMSAVE
476 free( aseq1 );
477 free( aseq2 );
478 #else
479 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
480 for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
481
482 FreeCharMtx( aseq1 );
483 FreeCharMtx( aseq2 );
484 #endif
485
486 // fprintf( stderr, "value = %f\n", value );
487
488 return( value );
489 }
490 // fprintf( stderr, "Trying to divide the mtx\n" );
491
492 ll1 = ( (int)(lgth1) ) + 100;
493 ll2 = ( (int)(lgth2) ) + 100;
494
495 // fprintf( stderr, "ll1,ll2=%d,%d\n", ll1, ll2 );
496
497 w1 = AllocateFloatVec( ll2+2 );
498 w2 = AllocateFloatVec( ll2+2 );
499 // match = AllocateFloatVec( ll2+2 );
500 midw = AllocateFloatVec( ll2+2 );
501 midn = AllocateFloatVec( ll2+2 );
502 midm = AllocateFloatVec( ll2+2 );
503 jumpbacki = AllocateIntVec( ll2+2 );
504 jumpbackj = AllocateIntVec( ll2+2 );
505 jumpforwi = AllocateIntVec( ll2+2 );
506 jumpforwj = AllocateIntVec( ll2+2 );
507 jumpdummi = AllocateIntVec( ll2+2 );
508 jumpdummj = AllocateIntVec( ll2+2 );
509
510 initverticalw = AllocateFloatVec( ll1+2 );
511 lastverticalw = AllocateFloatVec( ll1+2 );
512
513 m = AllocateFloatVec( ll2+2 );
514 mp = AllocateIntVec( ll2+2 );
515 gaps = AllocateCharVec( MAX( ll1, ll2 ) + 2 );
516
517 doublework = AllocateFloatMtx( MAX( ll1, ll2 )+2, nalphabets );
518 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, nalphabets );
519
520 #if DEBUG
521 fprintf( stderr, "succeeded\n" );
522 #endif
523
524 #if STOREWM
525 WMMTX = AllocateFloatMtx( ll1, ll2 );
526 WMMTX2 = AllocateFloatMtx( ll1, ll2 );
527 #endif
528 #if 0
529 shortmtx = AllocateShortMtx( ll1, ll2 );
530
531 #if DEBUG
532 fprintf( stderr, "succeeded\n\n" );
533 #endif
534
535 ijp = shortmtx;
536 #endif
537
538 currentw = w1;
539 previousw = w2;
540
541 match_ribosum( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, doublework, intwork, 1 );
542
543 match_ribosum( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, doublework, intwork, 1 );
544
545 for( i=1; i<lgth1+1; i++ )
546 {
547 initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] );
548 }
549 for( j=1; j<lgth2+1; j++ )
550 {
551 currentw[j] += ( ogcp2[0] + fgcp2[j-1] );
552 }
553
554 #if STOREWM
555 WMMTX[0][0] = initverticalw[0];
556 for( i=1; i<lgth1+1; i++ )
557 {
558 WMMTX[i][0] = initverticalw[i];
559 }
560 for( j=1; j<lgth2+1; j++ )
561 {
562 WMMTX[0][j] = currentw[j];
563 }
564 #endif
565
566
567 for( j=1; j<lgth2+1; ++j )
568 {
569 m[j] = currentw[j-1] + ogcp1[1];
570 // m[j] = currentw[j-1];
571 mp[j] = 0;
572 }
573
574 lastverticalw[0] = currentw[lgth2-1];
575
576 imid = lgth1 * 0.5;
577
578 jumpi = 0; // atode kawaru.
579 lasti = lgth1+1;
580 #if STOREWM
581 for( i=1; i<lasti; i++ )
582 #else
583 for( i=1; i<=imid; i++ )
584 #endif
585 {
586 wtmp = previousw;
587 previousw = currentw;
588 currentw = wtmp;
589
590 previousw[0] = initverticalw[i-1];
591
592 match_ribosum( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, doublework, intwork, 0 );
593 currentw[0] = initverticalw[i];
594
595 m[0] = ogcp1[i];
596 #if STOREM
597 WMMTX2[i][0] = m[0];
598 #endif
599 if( i == imid ) midm[0] = m[0];
600
601 mi = previousw[0] + ogcp2[1];
602 // mi = previousw[0];
603 mpi = 0;
604
605
606 // ijppt = ijp[i] + 1;
607 mjpt = m + 1;
608 prept = previousw;
609 curpt = currentw + 1;
610 mpjpt = mp + 1;
611
612
613 lastj = lgth2+1;
614 for( j=1; j<lastj; j++ )
615 {
616
617 wm = *prept;
618
619 #if 0
620 fprintf( stderr, "%5.0f->", wm );
621 #endif
622 g = mi + fgcp2[j-1];
623 // g = mi + fpenalty;
624 #if 0
625 fprintf( stderr, "%5.0f?", g );
626 #endif
627 if( g > wm )
628 {
629 wm = g;
630 // *ijppt = -( j - mpi );
631 }
632 g = *prept + ogcp2[j];
633 // g = *prept;
634 if( g >= mi )
635 {
636 mi = g;
637 mpi = j-1;
638 }
639 #if USE_PENALTY_EX
640 mi += fpenalty_ex;
641 #endif
642
643 g = *mjpt + fgcp1[i-1];
644 // g = *mjpt + fpenalty;
645 #if 0
646 fprintf( stderr, "%5.0f?", g );
647 #endif
648 if( g > wm )
649 {
650 wm = g;
651 // *ijppt = +( i - *mpjpt );
652 }
653
654
655 g = *prept + ogcp1[i];
656 // g = *prept;
657 if( g >= *mjpt )
658 {
659 *mjpt = g;
660 *mpjpt = i-1;
661 }
662 #if USE_PENALTY_EX
663 m[j] += fpenalty_ex;
664 #endif
665 #if LOCAL
666 if( wm < localthr )
667 {
668 // fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
669 wm = 0;
670 }
671 #endif
672
673 #if 0
674 fprintf( stderr, "%5.0f ", wm );
675 #endif
676 *curpt += wm;
677
678
679 #if STOREWM
680 WMMTX[i][j] = *curpt;
681 WMMTX2[i][j] = *mjpt;
682 #endif
683
684 if( i == imid ) //muda
685 {
686 jumpbackj[j] = *mpjpt; // muda atode matomeru
687 jumpbacki[j] = mpi; // muda atode matomeru
688 // fprintf( stderr, "jumpbackj[%d] in forward dp is %d\n", j, *mpjpt );
689 // fprintf( stderr, "jumpbacki[%d] in forward dp is %d\n", j, mpi );
690 midw[j] = *curpt;
691 midm[j] = *mjpt;
692 midn[j] = mi;
693 }
694
695 // fprintf( stderr, "m[%d] = %f\n", j, m[j] );
696 mjpt++;
697 prept++;
698 mpjpt++;
699 curpt++;
700
701 }
702 lastverticalw[i] = currentw[lgth2-1];
703
704 #if STOREWM
705 WMMTX2[i][lgth2] = m[lgth2-1];
706 #endif
707
708 #if 0 // ue
709 if( i == imid )
710 {
711 for( j=0; j<lgth2; j++ ) midw[j] = currentw[j];
712 for( j=0; j<lgth2; j++ ) midm[j] = m[j];
713 }
714 #endif
715 }
716 // for( j=0; j<lgth2; j++ ) midw[j] = WMMTX[imid][j];
717 // for( j=0; j<lgth2; j++ ) midm[j] = WMMTX2[imid][j];
718
719 #if 0
720 for( i=0; i<lgth1; i++ )
721 {
722 for( j=0; j<lgth2; j++ )
723 {
724 fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
725 }
726 fprintf( stderr, "\n" );
727 }
728 fprintf( stderr, "\n" );
729 fprintf( stderr, "WMMTX2 = \n" );
730 for( i=0; i<lgth1; i++ )
731 {
732 for( j=0; j<lgth2; j++ )
733 {
734 fprintf( stderr, "% 10.2f ", WMMTX2[i][j] );
735 }
736 fprintf( stderr, "\n" );
737 }
738 fprintf( stderr, "\n" );
739 #endif
740
741 // gyakudp
742
743 match_ribosum( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, doublework, intwork, 1 );
744 match_ribosum( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, doublework, intwork, 1 );
745
746 for( i=0; i<lgth1-1; i++ )
747 {
748 initverticalw[i] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
749 // initverticalw[i] += fpenalty;
750 }
751 for( j=0; j<lgth2-1; j++ )
752 {
753 currentw[j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
754 // currentw[j] += fpenalty;
755 }
756
757 #if STOREWM
758 for( i=0; i<lgth1-1; i++ )
759 {
760 WMMTX[i][lgth2-1] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
761 // fprintf( stderr, "fgcp1[lgth1-1] + ogcp1[i+1] = %f\n", fgcp1[lgth1-1] + ogcp1[i+1] );
762 }
763 for( j=0; j<lgth2-1; j++ )
764 {
765 WMMTX[lgth1-1][j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
766 // fprintf( stderr, "fgcp2[lgth2-1] + ogcp2[j+1] = %f\n", fgcp2[lgth2-1] + ogcp2[j+1] );
767 }
768 #endif
769
770
771
772
773
774
775 for( j=lgth2-1; j>0; --j )
776 {
777 m[j-1] = currentw[j] + fgcp2[lgth2-2];
778 // m[j-1] = currentw[j];
779 mp[j] = lgth1-1;
780 }
781
782 // for( j=0; j<lgth2; j++ ) m[j] = 0.0;
783 // m[lgth2-1] ha irunoka?
784
785
786 // for( i=lgth1-2; i>=imid; i-- )
787 firstm = -9999999.9;
788 firstmp = lgth1-1;
789 for( i=lgth1-2; i>-1; i-- )
790 {
791 wtmp = previousw;
792 previousw = currentw;
793 currentw = wtmp;
794 previousw[lgth2-1] = initverticalw[i+1];
795 // match_calc( currentw, seq1, seq2, i, lgth2 );
796 match_ribosum( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, doublework, intwork, 0 );
797
798 currentw[lgth2-1] = initverticalw[i];
799
800 // m[lgth2] = fgcp1[i];
801 // WMMTX2[i][lgth2] += m[lgth2];
802 // fprintf( stderr, "m[] = %f\n", m[lgth2] );
803
804 mi = previousw[lgth2-1] + fgcp2[lgth2-2];
805 // mi = previousw[lgth2-1];
806 mpi = lgth2 - 1;
807
808 mjpt = m + lgth2 - 2;
809 prept = previousw + lgth2 - 1;
810 curpt = currentw + lgth2 - 2;
811 mpjpt = mp + lgth2 - 2;
812
813
814 for( j=lgth2-2; j>-1; j-- )
815 {
816 wm = *prept;
817 ijpi = i+1;
818 ijpj = j+1;
819
820 g = mi + ogcp2[j+1];
821 // g = mi + fpenalty;
822 if( g > wm )
823 {
824 wm = g;
825 ijpj = mpi;
826 ijpi = i+1;
827 }
828
829 g = *prept + fgcp2[j];
830 // g = *prept;
831 if( g >= mi )
832 {
833 // fprintf( stderr, "i,j=%d,%d - renewed! mpi = %d\n", i, j, j+1 );
834 mi = g;
835 mpi = j + 1;
836 }
837
838 #if USE_PENALTY_EX
839 mi += fpenalty_ex;
840 #endif
841
842 // fprintf( stderr, "i,j=%d,%d *mpjpt = %d\n", i, j, *mpjpt );
843 g = *mjpt + ogcp1[i+1];
844 // g = *mjpt + fpenalty;
845 if( g > wm )
846 {
847 wm = g;
848 ijpi = *mpjpt;
849 ijpj = j+1;
850 }
851
852 // if( i == imid )fprintf( stderr, "i,j=%d,%d \n", i, j );
853 g = *prept + fgcp1[i];
854 // g = *prept;
855 if( g >= *mjpt )
856 {
857 *mjpt = g;
858 *mpjpt = i + 1;
859 }
860
861 #if USE_PENALTY_EX
862 m[j] += fpenalty_ex;
863 #endif
864
865 if( i == jumpi || i == imid - 1 )
866 {
867 jumpforwi[j] = ijpi; //muda
868 jumpforwj[j] = ijpj; //muda
869 // fprintf( stderr, "jumpfori[%d] = %d\n", j, ijpi );
870 // fprintf( stderr, "jumpforj[%d] = %d\n", j, ijpj );
871 }
872 if( i == imid ) // muda
873 {
874 midw[j] += wm;
875 // midm[j+1] += *mjpt + fpenalty; //??????
876 midm[j+1] += *mjpt; //??????
877 }
878 if( i == imid - 1 )
879 {
880 // midn[j] += mi + fpenalty; //????
881 midn[j] += mi; //????
882 }
883 #if LOCAL
884 if( wm < localthr )
885 {
886 // fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
887 wm = 0;
888 }
889 #endif
890
891 #if STOREWM
892 WMMTX[i][j] += wm;
893 // WMMTX2[i][j+1] += *mjpt + fpenalty;
894 WMMTX2[i][j] += *curpt;
895 #endif
896 *curpt += wm;
897
898 mjpt--;
899 prept--;
900 mpjpt--;
901 curpt--;
902 }
903 // fprintf( stderr, "adding *mjpt (=%f) to WMMTX2[%d][%d]\n", *mjpt, i, j+1 );
904 g = *prept + fgcp1[i];
905 if( firstm < g )
906 {
907 firstm = g;
908 firstmp = i + 1;
909 }
910 #if STOREWM
911 // WMMTX2[i][j+1] += firstm;
912 #endif
913 if( i == imid ) midm[j+1] += firstm;
914
915 if( i == imid - 1 )
916 {
917 maxwm = midw[1];
918 jmid = 0;
919 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
920 for( j=2; j<lgth2-1; j++ )
921 {
922 wm = midw[j];
923 if( wm > maxwm )
924 {
925 jmid = j;
926 maxwm = wm;
927 }
928 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
929 }
930 for( j=0; j<lgth2+1; j++ )
931 {
932 wm = midm[j];
933 if( wm > maxwm )
934 {
935 jmid = j;
936 maxwm = wm;
937 }
938 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
939 }
940
941 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
942
943
944 // fprintf( stderr, "### imid=%d, jmid=%d\n", imid, jmid );
945 wm = midw[jmid];
946 jumpi = imid-1;
947 jumpj = jmid-1;
948 if( jmid > 0 && midn[jmid-1] > wm ) //060413
949 {
950 jumpi = imid-1;
951 jumpj = jumpbacki[jmid];
952 wm = midn[jmid-1];
953 // fprintf( stderr, "rejump (n)\n" );
954 }
955 if( midm[jmid] > wm )
956 {
957 jumpi = jumpbackj[jmid];
958 jumpj = jmid-1;
959 wm = midm[jmid];
960 // fprintf( stderr, "rejump (m) jumpi=%d\n", jumpi );
961 }
962
963
964 // fprintf( stderr, "--> imid=%d, jmid=%d\n", imid, jmid );
965 // fprintf( stderr, "--> jumpi=%d, jumpj=%d\n", jumpi, jumpj );
966 #if 0
967 fprintf( stderr, "imid = %d\n", imid );
968 fprintf( stderr, "midn = \n" );
969 for( j=0; j<lgth2; j++ )
970 {
971 fprintf( stderr, "% 7.1f ", midn[j] );
972 }
973 fprintf( stderr, "\n" );
974 fprintf( stderr, "midw = \n" );
975 for( j=0; j<lgth2; j++ )
976 {
977 fprintf( stderr, "% 7.1f ", midw[j] );
978 }
979 fprintf( stderr, "\n" );
980 fprintf( stderr, "midm = \n" );
981 for( j=0; j<lgth2; j++ )
982 {
983 fprintf( stderr, "% 7.1f ", midm[j] );
984 }
985 fprintf( stderr, "\n" );
986 #endif
987 // fprintf( stderr, "maxwm = %f\n", maxwm );
988 }
989 if( i == jumpi ) //saki?
990 {
991 // fprintf( stderr, "imid, jumpi = %d,%d\n", imid, jumpi );
992 // fprintf( stderr, "jmid, jumpj = %d,%d\n", jmid, jumpj );
993 if( jmid == 0 )
994 {
995 // fprintf( stderr, "CHUI2!\n" );
996 jumpj = 0; jmid = 1;
997 jumpi = firstmp - 1;
998 imid = firstmp;
999 }
1000
1001 #if 0
1002 else if( jmid == lgth2 )
1003 {
1004 fprintf( stderr, "CHUI1!\n" );
1005 jumpi=0; jumpj=0;
1006 imid=jumpforwi[0]; jmid=lgth2-1;
1007 }
1008 #else // 060414
1009 else if( jmid >= lgth2 )
1010 {
1011 // fprintf( stderr, "CHUI1!\n" );
1012 jumpi=imid-1; jmid=lgth2;
1013 jumpj = lgth2-1;
1014 }
1015 #endif
1016 else
1017 {
1018 imid = jumpforwi[jumpj];
1019 jmid = jumpforwj[jumpj];
1020 }
1021 #if 0
1022 fprintf( stderr, "jumpi -> %d\n", jumpi );
1023 fprintf( stderr, "jumpj -> %d\n", jumpj );
1024 fprintf( stderr, "imid -> %d\n", imid );
1025 fprintf( stderr, "jmid -> %d\n", jmid );
1026 #endif
1027
1028 #if STOREWM
1029 // break;
1030 #else
1031 break;
1032 #endif
1033 }
1034 }
1035 #if 0
1036 jumpi=0; jumpj=0;
1037 imid=lgth1-1; jmid=lgth2-1;
1038 }
1039 #endif
1040
1041 // fprintf( stderr, "imid = %d, but jumpi = %d\n", imid, jumpi );
1042 // fprintf( stderr, "jmid = %d, but jumpj = %d\n", jmid, jumpj );
1043
1044 // for( j=0; j<lgth2; j++ ) midw[j] += currentw[j];
1045 // for( j=0; j<lgth2; j++ ) midm[j] += m[j+1];
1046 // for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1047 // for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1048
1049
1050
1051 for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
1052 map[i][j] = WMMTX[i][j] / maxwm;
1053 // map[i][j] = WMMTX2[i][j] / maxwm;
1054
1055 #if STOREWM
1056
1057 #if 0
1058 for( i=0; i<lgth1; i++ )
1059 {
1060 double maxpairscore = -9999.9;
1061 double tmpscore;
1062
1063 for( j=0; j<lgth2; j++ )
1064 {
1065 if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1066 {
1067 map12[i].pos = j;
1068 map12[i].score = tmpscore;
1069 maxpairscore = tmpscore;
1070 }
1071 }
1072
1073 for( k=0; k<lgth1; k++ )
1074 {
1075 if( i == k ) continue;
1076 if( map12[i].score <= WMMTX[k][map12[i].pos] )
1077 break;
1078 }
1079 if( k != lgth1 )
1080 {
1081 map12[i].pos = -1;
1082 map12[i].score = -1.0;
1083 }
1084 fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", i, map12[i].pos, map12[i].score, seq1[0][i], seq2[0][map12[i].pos] );
1085 }
1086 for( j=0; j<lgth2; j++ )
1087 {
1088 double maxpairscore = -9999.9;
1089 double tmpscore;
1090
1091 for( i=0; i<lgth1; i++ )
1092 {
1093 if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1094 {
1095 map21[j].pos = i;
1096 map21[j].score = tmpscore;
1097 maxpairscore = tmpscore;
1098 }
1099 }
1100
1101 for( k=0; k<lgth2; k++ )
1102 {
1103 if( j == k ) continue;
1104 if( map21[j].score <= WMMTX[map21[j].pos][k] )
1105 break;
1106 }
1107 if( k != lgth2 )
1108 {
1109 map21[j].pos = -1;
1110 map21[j].score = -1.0;
1111 }
1112 fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", j, map21[j].pos, map21[j].score, seq2[0][j], seq1[0][map21[j].pos] );
1113 }
1114
1115 for( i=0; i<lgth1; i++ )
1116 {
1117 if( map12[i].pos != -1 ) if( map21[map12[i].pos].pos != i )
1118 fprintf( stderr, "ERROR i=%d, but map12[i].pos=%d and map21[map12[i].pos]=%d\n", i, map12[i].pos, map21[map12[i].pos].pos );
1119 }
1120 #endif
1121
1122 #if 0
1123 fprintf( stderr, "WMMTX = \n" );
1124 for( i=0; i<lgth1; i++ )
1125 {
1126 fprintf( stderr, "%d ", i );
1127 for( j=0; j<lgth2; j++ )
1128 {
1129 fprintf( stderr, "% 7.2f ", WMMTX[i][j] );
1130 }
1131 fprintf( stderr, "\n" );
1132 }
1133 // fprintf( stderr, "WMMTX2 = (p = %f)\n", fpenalty );
1134 for( i=0; i<lgth1; i++ )
1135 {
1136 fprintf( stderr, "%d ", i );
1137 for( j=0; j<lgth2+1; j++ )
1138 {
1139 fprintf( stderr, "% 7.2f ", WMMTX2[i][j] );
1140 }
1141 fprintf( stderr, "\n" );
1142 }
1143 #endif
1144
1145 #if 0
1146 fprintf( stdout, "#WMMTX = \n" );
1147 for( i=0; i<lgth1; i++ )
1148 {
1149 // fprintf( stdout, "%d ", i );
1150 for( j=0; j<lgth2; j++ )
1151 {
1152 // if( WMMTX[i][j] > amino_dis['a']['g'] -1 )
1153 fprintf( stdout, "%d %d %8.1f", i, j, WMMTX[i][j] );
1154 if( WMMTX[i][j] == maxwm )
1155 fprintf( stdout, "selected \n" );
1156 else
1157 fprintf( stdout, "\n" );
1158 }
1159 fprintf( stdout, "\n" );
1160 }
1161 #endif
1162
1163 #if 0
1164
1165 fprintf( stderr, "jumpbacki = \n" );
1166 for( j=0; j<lgth2; j++ )
1167 {
1168 fprintf( stderr, "% 7d ", jumpbacki[j] );
1169 }
1170 fprintf( stderr, "\n" );
1171 fprintf( stderr, "jumpbackj = \n" );
1172 for( j=0; j<lgth2; j++ )
1173 {
1174 fprintf( stderr, "% 7d ", jumpbackj[j] );
1175 }
1176 fprintf( stderr, "\n" );
1177 fprintf( stderr, "jumpforwi = \n" );
1178 for( j=0; j<lgth2; j++ )
1179 {
1180 fprintf( stderr, "% 7d ", jumpforwi[j] );
1181 }
1182 fprintf( stderr, "\n" );
1183 fprintf( stderr, "jumpforwj = \n" );
1184 for( j=0; j<lgth2; j++ )
1185 {
1186 fprintf( stderr, "% 7d ", jumpforwj[j] );
1187 }
1188 fprintf( stderr, "\n" );
1189 #endif
1190
1191
1192 #endif
1193
1194
1195 // Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1196
1197 #if 0 // irukamo
1198 resultlen = strlen( mseq1[0] );
1199 if( alloclen < resultlen || resultlen > N )
1200 {
1201 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1202 ErrorExit( "LENGTH OVER!\n" );
1203 }
1204 #endif
1205
1206
1207
1208 #if 0
1209 fprintf( stderr, "jumpi = %d, imid = %d\n", jumpi, imid );
1210 fprintf( stderr, "jumpj = %d, jmid = %d\n", jumpj, jmid );
1211
1212 fprintf( stderr, "imid = %d\n", imid );
1213 fprintf( stderr, "jmid = %d\n", jmid );
1214 #endif
1215
1216
1217 FreeFloatVec( w1 );
1218 FreeFloatVec( w2 );
1219 FreeFloatVec( initverticalw );
1220 FreeFloatVec( lastverticalw );
1221 FreeFloatVec( midw );
1222 FreeFloatVec( midm );
1223 FreeFloatVec( midn );
1224
1225 FreeIntVec( jumpbacki );
1226 FreeIntVec( jumpbackj );
1227 FreeIntVec( jumpforwi );
1228 FreeIntVec( jumpforwj );
1229 FreeIntVec( jumpdummi );
1230 FreeIntVec( jumpdummj );
1231
1232 FreeFloatVec( m );
1233 FreeIntVec( mp );
1234
1235 FreeFloatMtx( doublework );
1236 FreeIntMtx( intwork );
1237
1238 #if STOREWM
1239 FreeFloatMtx( WMMTX );
1240 FreeFloatMtx( WMMTX2 );
1241 #endif
1242
1243 return( value );
1244
1245 }
MSalignmm_rec(int icyc,int jcyc,double * eff1,double * eff2,char ** seq1,char ** seq2,double ** cpmx1,double ** cpmx2,int ist,int ien,int jst,int jen,int alloclen,char ** mseq1,char ** mseq2,int depth,double ** gapinfo,double ** map)1246 static double MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char **seq1, char **seq2, double **cpmx1, double **cpmx2, int ist, int ien, int jst, int jen, int alloclen, char **mseq1, char **mseq2, int depth, double **gapinfo, double **map )
1247 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1248 {
1249 double value = 0.0;
1250 register int i, j;
1251 char **aseq1, **aseq2;
1252 int ll1, ll2;
1253 int lasti, lastj, imid, jmid=0;
1254 double wm = 0.0; /* int ?????? */
1255 double g;
1256 double *currentw, *previousw;
1257 #if USE_PENALTY_EX
1258 double fpenalty_ex = (double)RNApenalty_ex;
1259 #endif
1260 // double fpenalty = (double)penalty;
1261 double *wtmp;
1262 // short *ijppt;
1263 int *mpjpt;
1264 // short **ijp;
1265 int *mp;
1266 int mpi;
1267 double *mjpt, *prept, *curpt;
1268 double mi;
1269 double *m;
1270 double *w1, *w2;
1271 // double *match;
1272 double *initverticalw; /* kufuu sureba iranai */
1273 double *lastverticalw; /* kufuu sureba iranai */
1274 int **intwork;
1275 double **doublework;
1276 // short **shortmtx;
1277 #if STOREWM
1278 double **WMMTX;
1279 double **WMMTX2;
1280 #endif
1281 double *midw;
1282 double *midm;
1283 double *midn;
1284 int lgth1, lgth2;
1285 double maxwm = 0.0;
1286 int *jumpforwi;
1287 int *jumpforwj;
1288 int *jumpbacki;
1289 int *jumpbackj;
1290 int *jumpdummi; //muda
1291 int *jumpdummj; //muda
1292 int jumpi, jumpj = 0;
1293 char *gaps;
1294 int ijpi, ijpj;
1295 double *ogcp1;
1296 double *fgcp1;
1297 double *ogcp2;
1298 double *fgcp2;
1299 double firstm;
1300 int firstmp;
1301 #if 0
1302 static char ttt1[50000];
1303 static char ttt2[50000];
1304 #endif
1305
1306 localthr = -offset + 500; // 0?
1307
1308 ogcp1 = gapinfo[0] + ist;
1309 fgcp1 = gapinfo[1] + ist;
1310 ogcp2 = gapinfo[2] + jst;
1311 fgcp2 = gapinfo[3] + jst;
1312
1313 depth++;
1314 reccycle++;
1315
1316 lgth1 = ien-ist+1;
1317 lgth2 = jen-jst+1;
1318
1319 // if( lgth1 < 5 )
1320 // fprintf( stderr, "\nWARNING: lgth1 = %d\n", lgth1 );
1321 // if( lgth2 < 5 )
1322 // fprintf( stderr, "\nWARNING: lgth2 = %d\n", lgth2 );
1323 //
1324
1325 #if 0
1326 fprintf( stderr, "==== MSalign (depth=%d, reccycle=%d), ist=%d, ien=%d, jst=%d, jen=%d\n", depth, reccycle, ist, ien, jst, jen );
1327 strncpy( ttt1, seq1[0]+ist, lgth1 );
1328 strncpy( ttt2, seq2[0]+jst, lgth2 );
1329 ttt1[lgth1] = 0;
1330 ttt2[lgth2] = 0;
1331 fprintf( stderr, "seq1 = %s\n", ttt1 );
1332 fprintf( stderr, "seq2 = %s\n", ttt2 );
1333 #endif
1334 if( lgth2 <= 0 ) // lgth1 <= 0 ha?
1335 {
1336 // fprintf( stderr, "\n\n==== jimei\n\n" );
1337 // exit( 1 );
1338 for( i=0; i<icyc; i++ )
1339 {
1340 strncpy( mseq1[i], seq1[i]+ist, lgth1 );
1341 mseq1[i][lgth1] = 0;
1342 }
1343 for( i=0; i<jcyc; i++ )
1344 {
1345 mseq2[i][0] = 0;
1346 for( j=0; j<lgth1; j++ )
1347 strcat( mseq2[i], "-" );
1348 }
1349
1350 // fprintf( stderr, "==== mseq1[0] (%d) = %s\n", depth, mseq1[0] );
1351 // fprintf( stderr, "==== mseq2[0] (%d) = %s\n", depth, mseq2[0] );
1352
1353 return( 0.0 );
1354 }
1355
1356 #if MEMSAVE
1357 aseq1 = AllocateCharMtx( icyc, 0 );
1358 aseq2 = AllocateCharMtx( jcyc, 0 );
1359 for( i=0; i<icyc; i++ ) aseq1[i] = mseq1[i];
1360 for( i=0; i<jcyc; i++ ) aseq2[i] = mseq2[i];
1361 #else
1362 aseq1 = AllocateCharMtx( icyc, lgth1+lgth2+100 );
1363 aseq2 = AllocateCharMtx( jcyc, lgth1+lgth2+100 );
1364 #endif
1365
1366 // if( lgth1 < DPTANNI && lgth2 < DPTANNI ) // & dato lgth ==1 no kanousei ga arunode yokunai
1367 // if( lgth1 < DPTANNI ) // kore mo lgth2 ga mijikasugiru kanousei ari
1368 if( lgth1 < DPTANNI || lgth2 < DPTANNI ) // zettai ni anzen ka?
1369 {
1370 // fprintf( stderr, "==== Going to _tanni\n" );
1371
1372 // value = MSalignmm_tanni( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ien, jst, jen, alloclen, aseq1, aseq2, gapinfo );
1373
1374
1375 #if MEMSAVE
1376 free( aseq1 );
1377 free( aseq2 );
1378 #else
1379 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
1380 for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
1381
1382 FreeCharMtx( aseq1 );
1383 FreeCharMtx( aseq2 );
1384 #endif
1385
1386 // fprintf( stderr, "value = %f\n", value );
1387
1388 return( value );
1389 }
1390 // fprintf( stderr, "Trying to divide the mtx\n" );
1391
1392 ll1 = ( (int)(lgth1) ) + 100;
1393 ll2 = ( (int)(lgth2) ) + 100;
1394
1395 // fprintf( stderr, "ll1,ll2=%d,%d\n", ll1, ll2 );
1396
1397 w1 = AllocateFloatVec( ll2+2 );
1398 w2 = AllocateFloatVec( ll2+2 );
1399 // match = AllocateFloatVec( ll2+2 );
1400 midw = AllocateFloatVec( ll2+2 );
1401 midn = AllocateFloatVec( ll2+2 );
1402 midm = AllocateFloatVec( ll2+2 );
1403 jumpbacki = AllocateIntVec( ll2+2 );
1404 jumpbackj = AllocateIntVec( ll2+2 );
1405 jumpforwi = AllocateIntVec( ll2+2 );
1406 jumpforwj = AllocateIntVec( ll2+2 );
1407 jumpdummi = AllocateIntVec( ll2+2 );
1408 jumpdummj = AllocateIntVec( ll2+2 );
1409
1410 initverticalw = AllocateFloatVec( ll1+2 );
1411 lastverticalw = AllocateFloatVec( ll1+2 );
1412
1413 m = AllocateFloatVec( ll2+2 );
1414 mp = AllocateIntVec( ll2+2 );
1415 gaps = AllocateCharVec( MAX( ll1, ll2 ) + 2 );
1416
1417 doublework = AllocateFloatMtx( MAX( ll1, ll2 )+2, nalphabets );
1418 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, nalphabets );
1419
1420 #if DEBUG
1421 fprintf( stderr, "succeeded\n" );
1422 #endif
1423
1424 #if STOREWM
1425 WMMTX = AllocateFloatMtx( ll1, ll2 );
1426 WMMTX2 = AllocateFloatMtx( ll1, ll2 );
1427 #endif
1428 #if 0
1429 shortmtx = AllocateShortMtx( ll1, ll2 );
1430
1431 #if DEBUG
1432 fprintf( stderr, "succeeded\n\n" );
1433 #endif
1434
1435 ijp = shortmtx;
1436 #endif
1437
1438 currentw = w1;
1439 previousw = w2;
1440
1441 match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, doublework, intwork, 1 );
1442
1443 match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, doublework, intwork, 1 );
1444
1445 for( i=1; i<lgth1+1; i++ )
1446 {
1447 initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] );
1448 }
1449 for( j=1; j<lgth2+1; j++ )
1450 {
1451 currentw[j] += ( ogcp2[0] + fgcp2[j-1] );
1452 }
1453
1454 #if STOREWM
1455 WMMTX[0][0] = initverticalw[0];
1456 for( i=1; i<lgth1+1; i++ )
1457 {
1458 WMMTX[i][0] = initverticalw[i];
1459 }
1460 for( j=1; j<lgth2+1; j++ )
1461 {
1462 WMMTX[0][j] = currentw[j];
1463 }
1464 #endif
1465
1466
1467 for( j=1; j<lgth2+1; ++j )
1468 {
1469 m[j] = currentw[j-1] + ogcp1[1];
1470 // m[j] = currentw[j-1];
1471 mp[j] = 0;
1472 }
1473
1474 lastverticalw[0] = currentw[lgth2-1];
1475
1476 imid = lgth1 * 0.5;
1477
1478 jumpi = 0; // atode kawaru.
1479 lasti = lgth1+1;
1480 #if STOREWM
1481 for( i=1; i<lasti; i++ )
1482 #else
1483 for( i=1; i<=imid; i++ )
1484 #endif
1485 {
1486 wtmp = previousw;
1487 previousw = currentw;
1488 currentw = wtmp;
1489
1490 previousw[0] = initverticalw[i-1];
1491
1492 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, doublework, intwork, 0 );
1493 currentw[0] = initverticalw[i];
1494
1495 m[0] = ogcp1[i];
1496 #if STOREM
1497 WMMTX2[i][0] = m[0];
1498 #endif
1499 if( i == imid ) midm[0] = m[0];
1500
1501 mi = previousw[0] + ogcp2[1];
1502 // mi = previousw[0];
1503 mpi = 0;
1504
1505
1506 // ijppt = ijp[i] + 1;
1507 mjpt = m + 1;
1508 prept = previousw;
1509 curpt = currentw + 1;
1510 mpjpt = mp + 1;
1511
1512
1513 lastj = lgth2+1;
1514 for( j=1; j<lastj; j++ )
1515 {
1516
1517 wm = *prept;
1518
1519 #if 0
1520 fprintf( stderr, "%5.0f->", wm );
1521 #endif
1522 g = mi + fgcp2[j-1];
1523 // g = mi + fpenalty;
1524 #if 0
1525 fprintf( stderr, "%5.0f?", g );
1526 #endif
1527 if( g > wm )
1528 {
1529 wm = g;
1530 // *ijppt = -( j - mpi );
1531 }
1532 g = *prept + ogcp2[j];
1533 // g = *prept;
1534 if( g >= mi )
1535 {
1536 mi = g;
1537 mpi = j-1;
1538 }
1539 #if USE_PENALTY_EX
1540 mi += fpenalty_ex;
1541 #endif
1542
1543 g = *mjpt + fgcp1[i-1];
1544 // g = *mjpt + fpenalty;
1545 #if 0
1546 fprintf( stderr, "%5.0f?", g );
1547 #endif
1548 if( g > wm )
1549 {
1550 wm = g;
1551 // *ijppt = +( i - *mpjpt );
1552 }
1553
1554
1555 g = *prept + ogcp1[i];
1556 // g = *prept;
1557 if( g >= *mjpt )
1558 {
1559 *mjpt = g;
1560 *mpjpt = i-1;
1561 }
1562 #if USE_PENALTY_EX
1563 m[j] += fpenalty_ex;
1564 #endif
1565 #if LOCAL
1566 if( wm < localthr )
1567 {
1568 // fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
1569 wm = 0;
1570 }
1571 #endif
1572
1573 #if 0
1574 fprintf( stderr, "%5.0f ", wm );
1575 #endif
1576 *curpt += wm;
1577
1578
1579 #if STOREWM
1580 WMMTX[i][j] = *curpt;
1581 WMMTX2[i][j] = *mjpt;
1582 #endif
1583
1584 if( i == imid ) //muda
1585 {
1586 jumpbackj[j] = *mpjpt; // muda atode matomeru
1587 jumpbacki[j] = mpi; // muda atode matomeru
1588 // fprintf( stderr, "jumpbackj[%d] in forward dp is %d\n", j, *mpjpt );
1589 // fprintf( stderr, "jumpbacki[%d] in forward dp is %d\n", j, mpi );
1590 midw[j] = *curpt;
1591 midm[j] = *mjpt;
1592 midn[j] = mi;
1593 }
1594
1595 // fprintf( stderr, "m[%d] = %f\n", j, m[j] );
1596 mjpt++;
1597 prept++;
1598 mpjpt++;
1599 curpt++;
1600
1601 }
1602 lastverticalw[i] = currentw[lgth2-1];
1603
1604 #if STOREWM
1605 WMMTX2[i][lgth2] = m[lgth2-1];
1606 #endif
1607
1608 #if 0 // ue
1609 if( i == imid )
1610 {
1611 for( j=0; j<lgth2; j++ ) midw[j] = currentw[j];
1612 for( j=0; j<lgth2; j++ ) midm[j] = m[j];
1613 }
1614 #endif
1615 }
1616 // for( j=0; j<lgth2; j++ ) midw[j] = WMMTX[imid][j];
1617 // for( j=0; j<lgth2; j++ ) midm[j] = WMMTX2[imid][j];
1618
1619 #if 0
1620 for( i=0; i<lgth1; i++ )
1621 {
1622 for( j=0; j<lgth2; j++ )
1623 {
1624 fprintf( stderr, "% 10.2f ", WMMTX[i][j] );
1625 }
1626 fprintf( stderr, "\n" );
1627 }
1628 fprintf( stderr, "\n" );
1629 fprintf( stderr, "WMMTX2 = \n" );
1630 for( i=0; i<lgth1; i++ )
1631 {
1632 for( j=0; j<lgth2; j++ )
1633 {
1634 fprintf( stderr, "% 10.2f ", WMMTX2[i][j] );
1635 }
1636 fprintf( stderr, "\n" );
1637 }
1638 fprintf( stderr, "\n" );
1639 #endif
1640
1641 // gyakudp
1642
1643 match_calc( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, doublework, intwork, 1 );
1644 match_calc( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, doublework, intwork, 1 );
1645
1646 for( i=0; i<lgth1-1; i++ )
1647 {
1648 initverticalw[i] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
1649 // initverticalw[i] += fpenalty;
1650 }
1651 for( j=0; j<lgth2-1; j++ )
1652 {
1653 currentw[j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
1654 // currentw[j] += fpenalty;
1655 }
1656
1657 #if STOREWM
1658 for( i=0; i<lgth1-1; i++ )
1659 {
1660 WMMTX[i][lgth2-1] += ( fgcp1[lgth1-1] + ogcp1[i+1] );
1661 // fprintf( stderr, "fgcp1[lgth1-1] + ogcp1[i+1] = %f\n", fgcp1[lgth1-1] + ogcp1[i+1] );
1662 }
1663 for( j=0; j<lgth2-1; j++ )
1664 {
1665 WMMTX[lgth1-1][j] += ( fgcp2[lgth2-1] + ogcp2[j+1] );
1666 // fprintf( stderr, "fgcp2[lgth2-1] + ogcp2[j+1] = %f\n", fgcp2[lgth2-1] + ogcp2[j+1] );
1667 }
1668 #endif
1669
1670
1671
1672
1673
1674
1675 for( j=lgth2-1; j>0; --j )
1676 {
1677 m[j-1] = currentw[j] + fgcp2[lgth2-2];
1678 // m[j-1] = currentw[j];
1679 mp[j] = lgth1-1;
1680 }
1681
1682 // for( j=0; j<lgth2; j++ ) m[j] = 0.0;
1683 // m[lgth2-1] ha irunoka?
1684
1685
1686 // for( i=lgth1-2; i>=imid; i-- )
1687 firstm = -9999999.9;
1688 firstmp = lgth1-1;
1689 for( i=lgth1-2; i>-1; i-- )
1690 {
1691 wtmp = previousw;
1692 previousw = currentw;
1693 currentw = wtmp;
1694 previousw[lgth2-1] = initverticalw[i+1];
1695 // match_calc( currentw, seq1, seq2, i, lgth2 );
1696 match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, doublework, intwork, 0 );
1697
1698 currentw[lgth2-1] = initverticalw[i];
1699
1700 // m[lgth2] = fgcp1[i];
1701 // WMMTX2[i][lgth2] += m[lgth2];
1702 // fprintf( stderr, "m[] = %f\n", m[lgth2] );
1703
1704 mi = previousw[lgth2-1] + fgcp2[lgth2-2];
1705 // mi = previousw[lgth2-1];
1706 mpi = lgth2 - 1;
1707
1708 mjpt = m + lgth2 - 2;
1709 prept = previousw + lgth2 - 1;
1710 curpt = currentw + lgth2 - 2;
1711 mpjpt = mp + lgth2 - 2;
1712
1713
1714 for( j=lgth2-2; j>-1; j-- )
1715 {
1716 wm = *prept;
1717 ijpi = i+1;
1718 ijpj = j+1;
1719
1720 g = mi + ogcp2[j+1];
1721 // g = mi + fpenalty;
1722 if( g > wm )
1723 {
1724 wm = g;
1725 ijpj = mpi;
1726 ijpi = i+1;
1727 }
1728
1729 g = *prept + fgcp2[j];
1730 // g = *prept;
1731 if( g >= mi )
1732 {
1733 // fprintf( stderr, "i,j=%d,%d - renewed! mpi = %d\n", i, j, j+1 );
1734 mi = g;
1735 mpi = j + 1;
1736 }
1737
1738 #if USE_PENALTY_EX
1739 mi += fpenalty_ex;
1740 #endif
1741
1742 // fprintf( stderr, "i,j=%d,%d *mpjpt = %d\n", i, j, *mpjpt );
1743 g = *mjpt + ogcp1[i+1];
1744 // g = *mjpt + fpenalty;
1745 if( g > wm )
1746 {
1747 wm = g;
1748 ijpi = *mpjpt;
1749 ijpj = j+1;
1750 }
1751
1752 // if( i == imid )fprintf( stderr, "i,j=%d,%d \n", i, j );
1753 g = *prept + fgcp1[i];
1754 // g = *prept;
1755 if( g >= *mjpt )
1756 {
1757 *mjpt = g;
1758 *mpjpt = i + 1;
1759 }
1760
1761 #if USE_PENALTY_EX
1762 m[j] += fpenalty_ex;
1763 #endif
1764
1765 if( i == jumpi || i == imid - 1 )
1766 {
1767 jumpforwi[j] = ijpi; //muda
1768 jumpforwj[j] = ijpj; //muda
1769 // fprintf( stderr, "jumpfori[%d] = %d\n", j, ijpi );
1770 // fprintf( stderr, "jumpforj[%d] = %d\n", j, ijpj );
1771 }
1772 if( i == imid ) // muda
1773 {
1774 midw[j] += wm;
1775 // midm[j+1] += *mjpt + fpenalty; //??????
1776 midm[j+1] += *mjpt; //??????
1777 }
1778 if( i == imid - 1 )
1779 {
1780 // midn[j] += mi + fpenalty; //????
1781 midn[j] += mi; //????
1782 }
1783 #if LOCAL
1784 if( wm < localthr )
1785 {
1786 // fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
1787 wm = 0;
1788 }
1789 #endif
1790
1791 #if STOREWM
1792 WMMTX[i][j] += wm;
1793 // WMMTX2[i][j+1] += *mjpt + fpenalty;
1794 WMMTX2[i][j] += *curpt;
1795 #endif
1796 *curpt += wm;
1797
1798 mjpt--;
1799 prept--;
1800 mpjpt--;
1801 curpt--;
1802 }
1803 // fprintf( stderr, "adding *mjpt (=%f) to WMMTX2[%d][%d]\n", *mjpt, i, j+1 );
1804 g = *prept + fgcp1[i];
1805 if( firstm < g )
1806 {
1807 firstm = g;
1808 firstmp = i + 1;
1809 }
1810 #if STOREWM
1811 // WMMTX2[i][j+1] += firstm;
1812 #endif
1813 if( i == imid ) midm[j+1] += firstm;
1814
1815 if( i == imid - 1 )
1816 {
1817 maxwm = midw[1];
1818 jmid = 0;
1819 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1820 for( j=2; j<lgth2-1; j++ )
1821 {
1822 wm = midw[j];
1823 if( wm > maxwm )
1824 {
1825 jmid = j;
1826 maxwm = wm;
1827 }
1828 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1829 }
1830 for( j=0; j<lgth2+1; j++ )
1831 {
1832 wm = midm[j];
1833 if( wm > maxwm )
1834 {
1835 jmid = j;
1836 maxwm = wm;
1837 }
1838 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1839 }
1840
1841 // if( depth == 1 ) fprintf( stderr, "maxwm!! = %f\n", maxwm );
1842
1843
1844 // fprintf( stderr, "### imid=%d, jmid=%d\n", imid, jmid );
1845 wm = midw[jmid];
1846 jumpi = imid-1;
1847 jumpj = jmid-1;
1848 if( jmid > 0 && midn[jmid-1] > wm ) //060413
1849 {
1850 jumpi = imid-1;
1851 jumpj = jumpbacki[jmid];
1852 wm = midn[jmid-1];
1853 // fprintf( stderr, "rejump (n)\n" );
1854 }
1855 if( midm[jmid] > wm )
1856 {
1857 jumpi = jumpbackj[jmid];
1858 jumpj = jmid-1;
1859 wm = midm[jmid];
1860 // fprintf( stderr, "rejump (m) jumpi=%d\n", jumpi );
1861 }
1862
1863
1864 // fprintf( stderr, "--> imid=%d, jmid=%d\n", imid, jmid );
1865 // fprintf( stderr, "--> jumpi=%d, jumpj=%d\n", jumpi, jumpj );
1866 #if 0
1867 fprintf( stderr, "imid = %d\n", imid );
1868 fprintf( stderr, "midn = \n" );
1869 for( j=0; j<lgth2; j++ )
1870 {
1871 fprintf( stderr, "% 7.1f ", midn[j] );
1872 }
1873 fprintf( stderr, "\n" );
1874 fprintf( stderr, "midw = \n" );
1875 for( j=0; j<lgth2; j++ )
1876 {
1877 fprintf( stderr, "% 7.1f ", midw[j] );
1878 }
1879 fprintf( stderr, "\n" );
1880 fprintf( stderr, "midm = \n" );
1881 for( j=0; j<lgth2; j++ )
1882 {
1883 fprintf( stderr, "% 7.1f ", midm[j] );
1884 }
1885 fprintf( stderr, "\n" );
1886 #endif
1887 // fprintf( stderr, "maxwm = %f\n", maxwm );
1888 }
1889 if( i == jumpi ) //saki?
1890 {
1891 // fprintf( stderr, "imid, jumpi = %d,%d\n", imid, jumpi );
1892 // fprintf( stderr, "jmid, jumpj = %d,%d\n", jmid, jumpj );
1893 if( jmid == 0 )
1894 {
1895 // fprintf( stderr, "CHUI2!\n" );
1896 jumpj = 0; jmid = 1;
1897 jumpi = firstmp - 1;
1898 imid = firstmp;
1899 }
1900
1901 #if 0
1902 else if( jmid == lgth2 )
1903 {
1904 fprintf( stderr, "CHUI1!\n" );
1905 jumpi=0; jumpj=0;
1906 imid=jumpforwi[0]; jmid=lgth2-1;
1907 }
1908 #else // 060414
1909 else if( jmid >= lgth2 )
1910 {
1911 // fprintf( stderr, "CHUI1!\n" );
1912 jumpi=imid-1; jmid=lgth2;
1913 jumpj = lgth2-1;
1914 }
1915 #endif
1916 else
1917 {
1918 imid = jumpforwi[jumpj];
1919 jmid = jumpforwj[jumpj];
1920 }
1921 #if 0
1922 fprintf( stderr, "jumpi -> %d\n", jumpi );
1923 fprintf( stderr, "jumpj -> %d\n", jumpj );
1924 fprintf( stderr, "imid -> %d\n", imid );
1925 fprintf( stderr, "jmid -> %d\n", jmid );
1926 #endif
1927
1928 #if STOREWM
1929 // break;
1930 #else
1931 break;
1932 #endif
1933 }
1934 }
1935 #if 0
1936 jumpi=0; jumpj=0;
1937 imid=lgth1-1; jmid=lgth2-1;
1938 }
1939 #endif
1940
1941 // fprintf( stderr, "imid = %d, but jumpi = %d\n", imid, jumpi );
1942 // fprintf( stderr, "jmid = %d, but jumpj = %d\n", jmid, jumpj );
1943
1944 // for( j=0; j<lgth2; j++ ) midw[j] += currentw[j];
1945 // for( j=0; j<lgth2; j++ ) midm[j] += m[j+1];
1946 // for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1947 // for( j=0; j<lgth2; j++ ) midw[j] += WMMTX[imid][j];
1948
1949
1950
1951 for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
1952 map[i][j] = WMMTX[i][j] / maxwm;
1953 // map[i][j] = WMMTX2[i][j] / maxwm;
1954
1955 #if STOREWM
1956
1957 #if 0
1958 for( i=0; i<lgth1; i++ )
1959 {
1960 double maxpairscore = -9999.9;
1961 double tmpscore;
1962
1963 for( j=0; j<lgth2; j++ )
1964 {
1965 if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1966 {
1967 map12[i].pos = j;
1968 map12[i].score = tmpscore;
1969 maxpairscore = tmpscore;
1970 }
1971 }
1972
1973 for( k=0; k<lgth1; k++ )
1974 {
1975 if( i == k ) continue;
1976 if( map12[i].score <= WMMTX[k][map12[i].pos] )
1977 break;
1978 }
1979 if( k != lgth1 )
1980 {
1981 map12[i].pos = -1;
1982 map12[i].score = -1.0;
1983 }
1984 fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", i, map12[i].pos, map12[i].score, seq1[0][i], seq2[0][map12[i].pos] );
1985 }
1986 for( j=0; j<lgth2; j++ )
1987 {
1988 double maxpairscore = -9999.9;
1989 double tmpscore;
1990
1991 for( i=0; i<lgth1; i++ )
1992 {
1993 if( maxpairscore < (tmpscore=WMMTX[i][j]) )
1994 {
1995 map21[j].pos = i;
1996 map21[j].score = tmpscore;
1997 maxpairscore = tmpscore;
1998 }
1999 }
2000
2001 for( k=0; k<lgth2; k++ )
2002 {
2003 if( j == k ) continue;
2004 if( map21[j].score <= WMMTX[map21[j].pos][k] )
2005 break;
2006 }
2007 if( k != lgth2 )
2008 {
2009 map21[j].pos = -1;
2010 map21[j].score = -1.0;
2011 }
2012 fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", j, map21[j].pos, map21[j].score, seq2[0][j], seq1[0][map21[j].pos] );
2013 }
2014
2015 for( i=0; i<lgth1; i++ )
2016 {
2017 if( map12[i].pos != -1 ) if( map21[map12[i].pos].pos != i )
2018 fprintf( stderr, "ERROR i=%d, but map12[i].pos=%d and map21[map12[i].pos]=%d\n", i, map12[i].pos, map21[map12[i].pos].pos );
2019 }
2020 #endif
2021
2022 #if 0
2023 fprintf( stderr, "WMMTX = \n" );
2024 for( i=0; i<lgth1; i++ )
2025 {
2026 fprintf( stderr, "%d ", i );
2027 for( j=0; j<lgth2; j++ )
2028 {
2029 fprintf( stderr, "% 7.2f ", WMMTX[i][j] );
2030 }
2031 fprintf( stderr, "\n" );
2032 }
2033 // fprintf( stderr, "WMMTX2 = (p = %f)\n", fpenalty );
2034 for( i=0; i<lgth1; i++ )
2035 {
2036 fprintf( stderr, "%d ", i );
2037 for( j=0; j<lgth2+1; j++ )
2038 {
2039 fprintf( stderr, "% 7.2f ", WMMTX2[i][j] );
2040 }
2041 fprintf( stderr, "\n" );
2042 }
2043 exit( 1 );
2044 #endif
2045
2046 #if 0
2047 fprintf( stdout, "#WMMTX = \n" );
2048 for( i=0; i<lgth1; i++ )
2049 {
2050 // fprintf( stdout, "%d ", i );
2051 for( j=0; j<lgth2; j++ )
2052 {
2053 // if( WMMTX[i][j] > amino_dis['a']['g'] -1 )
2054 fprintf( stdout, "%d %d %8.1f", i, j, WMMTX[i][j] );
2055 if( WMMTX[i][j] == maxwm )
2056 fprintf( stdout, "selected \n" );
2057 else
2058 fprintf( stdout, "\n" );
2059 }
2060 fprintf( stdout, "\n" );
2061 }
2062 exit( 1 );
2063 #endif
2064
2065 #if 0
2066
2067 fprintf( stderr, "jumpbacki = \n" );
2068 for( j=0; j<lgth2; j++ )
2069 {
2070 fprintf( stderr, "% 7d ", jumpbacki[j] );
2071 }
2072 fprintf( stderr, "\n" );
2073 fprintf( stderr, "jumpbackj = \n" );
2074 for( j=0; j<lgth2; j++ )
2075 {
2076 fprintf( stderr, "% 7d ", jumpbackj[j] );
2077 }
2078 fprintf( stderr, "\n" );
2079 fprintf( stderr, "jumpforwi = \n" );
2080 for( j=0; j<lgth2; j++ )
2081 {
2082 fprintf( stderr, "% 7d ", jumpforwi[j] );
2083 }
2084 fprintf( stderr, "\n" );
2085 fprintf( stderr, "jumpforwj = \n" );
2086 for( j=0; j<lgth2; j++ )
2087 {
2088 fprintf( stderr, "% 7d ", jumpforwj[j] );
2089 }
2090 fprintf( stderr, "\n" );
2091 #endif
2092
2093
2094 #endif
2095
2096
2097 // Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
2098
2099 #if 0 // irukamo
2100 resultlen = strlen( mseq1[0] );
2101 if( alloclen < resultlen || resultlen > N )
2102 {
2103 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
2104 ErrorExit( "LENGTH OVER!\n" );
2105 }
2106 #endif
2107
2108
2109
2110 #if 0
2111 fprintf( stderr, "jumpi = %d, imid = %d\n", jumpi, imid );
2112 fprintf( stderr, "jumpj = %d, jmid = %d\n", jumpj, jmid );
2113
2114 fprintf( stderr, "imid = %d\n", imid );
2115 fprintf( stderr, "jmid = %d\n", jmid );
2116 #endif
2117
2118
2119 FreeFloatVec( w1 );
2120 FreeFloatVec( w2 );
2121 FreeFloatVec( initverticalw );
2122 FreeFloatVec( lastverticalw );
2123 FreeFloatVec( midw );
2124 FreeFloatVec( midm );
2125 FreeFloatVec( midn );
2126
2127 FreeIntVec( jumpbacki );
2128 FreeIntVec( jumpbackj );
2129 FreeIntVec( jumpforwi );
2130 FreeIntVec( jumpforwj );
2131 FreeIntVec( jumpdummi );
2132 FreeIntVec( jumpdummj );
2133
2134 FreeFloatVec( m );
2135 FreeIntVec( mp );
2136
2137 FreeFloatMtx( doublework );
2138 FreeIntMtx( intwork );
2139
2140 #if STOREWM
2141 FreeFloatMtx( WMMTX );
2142 FreeFloatMtx( WMMTX2 );
2143 #endif
2144
2145
2146 free( gaps );
2147 #if MEMSAVE
2148 free( aseq1 );
2149 free( aseq2 );
2150 #else
2151 FreeCharMtx( aseq1 );
2152 FreeCharMtx( aseq2 );
2153 #endif
2154
2155 return( value );
2156 }
2157
2158
2159
Lalignmm_hmout(char ** seq1,char ** seq2,double * eff1,double * eff2,int icyc,int jcyc,int alloclen,char * sgap1,char * sgap2,char * egap1,char * egap2,double ** map)2160 double Lalignmm_hmout( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, char *sgap1, char *sgap2, char *egap1, char *egap2, double **map )
2161 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
2162 {
2163 // int k;
2164 int i, j;
2165 int ll1, ll2;
2166 int lgth1, lgth2;
2167 double wm = 0.0; /* int ?????? */
2168 char **mseq1;
2169 char **mseq2;
2170 // char **mseq;
2171 double *ogcp1;
2172 double *ogcp2;
2173 double *fgcp1;
2174 double *fgcp2;
2175 double **cpmx1;
2176 double **cpmx2;
2177 double **gapinfo;
2178 // double fpenalty;
2179 double fpenalty = (double)RNApenalty;
2180 int nglen1, nglen2;
2181
2182
2183
2184
2185
2186 #if 0
2187 fprintf( stderr, "eff in SA+++align\n" );
2188 for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
2189 #endif
2190
2191 nglen1 = seqlen( seq1[0] );
2192 nglen2 = seqlen( seq2[0] );
2193
2194 lgth1 = strlen( seq1[0] );
2195 lgth2 = strlen( seq2[0] );
2196
2197 ll1 = ( (int)(lgth1) ) + 100;
2198 ll2 = ( (int)(lgth2) ) + 100;
2199
2200 mseq1 = AllocateCharMtx( icyc, ll1+ll2 );
2201 mseq2 = AllocateCharMtx( jcyc, ll1+ll2 );
2202
2203 gapinfo = AllocateFloatMtx( 4, 0 );
2204 ogcp1 = AllocateFloatVec( ll1+2 );
2205 ogcp2 = AllocateFloatVec( ll2+2 );
2206 fgcp1 = AllocateFloatVec( ll1+2 );
2207 fgcp2 = AllocateFloatVec( ll2+2 );
2208
2209
2210 cpmx1 = AllocateFloatMtx( ll1+2, nalphabets+1 );
2211 cpmx2 = AllocateFloatMtx( ll2+2, nalphabets+1 );
2212
2213 for( i=0; i<icyc; i++ )
2214 {
2215 if( strlen( seq1[i] ) != lgth1 )
2216 {
2217 fprintf( stderr, "i = %d / %d\n", i, icyc );
2218 fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2219 exit( 1 );
2220 }
2221 }
2222 for( j=0; j<jcyc; j++ )
2223 {
2224 if( strlen( seq2[j] ) != lgth2 )
2225 {
2226 fprintf( stderr, "j = %d / %d\n", j, jcyc );
2227 fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2228 exit( 1 );
2229 }
2230 }
2231
2232 MScpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
2233 MScpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
2234
2235
2236 #if 1
2237
2238 if( sgap1 )
2239 {
2240 new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 );
2241 new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 );
2242 new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap2 );
2243 new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 );
2244 }
2245 else
2246 {
2247 st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
2248 st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
2249 st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
2250 st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
2251 }
2252
2253 #if 1
2254 for( i=0; i<lgth1; i++ )
2255 {
2256 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty;
2257 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty;
2258 // fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
2259 }
2260 for( i=0; i<lgth2; i++ )
2261 {
2262 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty;
2263 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty;
2264 // fprintf( stderr, "fgcp2[%d] = %f\n", i, fgcp2[i] );
2265 }
2266 #else
2267 for( i=0; i<lgth1; i++ )
2268 {
2269 ogcp1[i] = 0.5 * fpenalty;
2270 fgcp1[i] = 0.5 * fpenalty;
2271 }
2272 for( i=0; i<lgth2; i++ )
2273 {
2274 ogcp2[i] = 0.5 * fpenalty;
2275 fgcp2[i] = 0.5 * fpenalty;
2276 }
2277 #endif
2278
2279 gapinfo[0] = ogcp1;
2280 gapinfo[1] = fgcp1;
2281 gapinfo[2] = ogcp2;
2282 gapinfo[3] = fgcp2;
2283 #endif
2284
2285 #if 0
2286 fprintf( stdout, "in MSalignmm.c\n" );
2287 for( i=0; i<icyc; i++ )
2288 {
2289 fprintf( stdout, ">%d of GROUP1\n", i );
2290 fprintf( stdout, "%s\n", seq1[i] );
2291 }
2292 for( i=0; i<jcyc; i++ )
2293 {
2294 fprintf( stdout, ">%d of GROUP2\n", i );
2295 fprintf( stdout, "%s\n", seq2[i] );
2296 }
2297 fflush( stdout );
2298 #endif
2299
2300 wm = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, 0, lgth1-1, 0, lgth2-1, alloclen, mseq1, mseq2, 0, gapinfo, map );
2301 #if DEBUG
2302 fprintf( stderr, " seq1[0] = %s\n", seq1[0] );
2303 fprintf( stderr, " seq2[0] = %s\n", seq2[0] );
2304 fprintf( stderr, "mseq1[0] = %s\n", mseq1[0] );
2305 fprintf( stderr, "mseq2[0] = %s\n", mseq2[0] );
2306 #endif
2307
2308 // fprintf( stderr, "wm = %f\n", wm );
2309
2310 #if 0
2311
2312 for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
2313 for( i=0; i<jcyc; i++ ) strcpy( seq2[i], mseq2[i] );
2314
2315 if( seqlen( seq1[0] ) != nglen1 )
2316 {
2317 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen1) seqlen(seq1[0])=%d but nglen1=%d\n", seqlen( seq1[0] ), nglen1 );
2318 fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
2319 exit( 1 );
2320 }
2321 if( seqlen( seq2[0] ) != nglen2 )
2322 {
2323 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen2) seqlen(seq2[0])=%d but nglen2=%d\n", seqlen( seq2[0] ), nglen2 );
2324 exit( 1 );
2325 }
2326 #endif
2327
2328 FreeFloatVec( ogcp1 );
2329 FreeFloatVec( ogcp2 );
2330 FreeFloatVec( fgcp1 );
2331 FreeFloatVec( fgcp2 );
2332 FreeFloatMtx( cpmx1 );
2333 FreeFloatMtx( cpmx2 );
2334 free( (void *)gapinfo );
2335
2336 FreeCharMtx( mseq1 );
2337 FreeCharMtx( mseq2 );
2338
2339 lgth1 = strlen( seq1[0] );
2340 lgth2 = strlen( seq2[0] );
2341 for( i=0; i<icyc; i++ )
2342 {
2343 if( strlen( seq1[i] ) != lgth1 )
2344 {
2345 fprintf( stderr, "i = %d / %d\n", i, icyc );
2346 fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2347 exit( 1 );
2348 }
2349 }
2350 for( j=0; j<jcyc; j++ )
2351 {
2352 if( strlen( seq2[j] ) != lgth2 )
2353 {
2354 fprintf( stderr, "j = %d / %d\n", j, jcyc );
2355 fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2356 exit( 1 );
2357 }
2358 }
2359
2360 return( wm );
2361 }
2362
Lalign2m2m_hmout(char ** seq1,char ** seq2,char ** seq1r,char ** seq2r,char * dir1,char * dir2,double * eff1,double * eff2,int icyc,int jcyc,int alloclen,char * sgap1,char * sgap2,char * egap1,char * egap2,double ** map)2363 double Lalign2m2m_hmout( char **seq1, char **seq2, char **seq1r, char **seq2r, char *dir1, char *dir2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, char *sgap1, char *sgap2, char *egap1, char *egap2, double **map )
2364 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
2365 {
2366 // int k;
2367 int i, j;
2368 int ll1, ll2;
2369 int lgth1, lgth2;
2370 double wm = 0.0; /* int ?????? */
2371 char **mseq1;
2372 char **mseq2;
2373 double *ogcp1;
2374 double *ogcp2;
2375 double *fgcp1;
2376 double *fgcp2;
2377 double **cpmx1;
2378 double **cpmx2;
2379 double **gapinfo;
2380 double fpenalty = (double)penalty;
2381 int nglen1, nglen2;
2382
2383 #if 0
2384 fprintf( stderr, "eff in SA+++align\n" );
2385 for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
2386 #endif
2387
2388 nglen1 = seqlen( seq1[0] );
2389 nglen2 = seqlen( seq2[0] );
2390
2391 lgth1 = strlen( seq1[0] );
2392 lgth2 = strlen( seq2[0] );
2393
2394 ll1 = ( (int)(lgth1) ) + 100;
2395 ll2 = ( (int)(lgth2) ) + 100;
2396
2397 mseq1 = AllocateCharMtx( icyc, ll1+ll2 );
2398 mseq2 = AllocateCharMtx( jcyc, ll1+ll2 );
2399
2400 gapinfo = AllocateFloatMtx( 4, 0 );
2401 ogcp1 = AllocateFloatVec( ll1+2 );
2402 ogcp2 = AllocateFloatVec( ll2+2 );
2403 fgcp1 = AllocateFloatVec( ll1+2 );
2404 fgcp2 = AllocateFloatVec( ll2+2 );
2405
2406
2407 cpmx1 = AllocateFloatMtx( ll1+2, 39 );
2408 cpmx2 = AllocateFloatMtx( ll2+2, 39 );
2409
2410 for( i=0; i<icyc; i++ )
2411 {
2412 if( strlen( seq1[i] ) != lgth1 )
2413 {
2414 fprintf( stderr, "i = %d / %d\n", i, icyc );
2415 fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2416 exit( 1 );
2417 }
2418 }
2419 for( j=0; j<jcyc; j++ )
2420 {
2421 if( strlen( seq2[j] ) != lgth2 )
2422 {
2423 fprintf( stderr, "j = %d / %d\n", j, jcyc );
2424 fprintf( stderr, "bug! hairetsu ga kowareta!\n" );
2425 exit( 1 );
2426 }
2427 }
2428
2429 #if 0
2430 MScpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
2431 MScpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
2432 #else
2433 cpmx_ribosum( seq1, seq1r, dir1, cpmx1, eff1, lgth1, icyc );
2434 cpmx_ribosum( seq2, seq2r, dir2, cpmx2, eff2, lgth2, jcyc );
2435 #endif
2436
2437
2438 #if 1
2439
2440 if( sgap1 )
2441 {
2442 new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 );
2443 new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 );
2444 new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap2 );
2445 new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 );
2446 }
2447 else
2448 {
2449 st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
2450 st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
2451 st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
2452 st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
2453 }
2454
2455 #if 1
2456 for( i=0; i<lgth1; i++ )
2457 {
2458 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty;
2459 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty;
2460 // fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
2461 }
2462 for( i=0; i<lgth2; i++ )
2463 {
2464 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty;
2465 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty;
2466 // fprintf( stderr, "fgcp2[%d] = %f\n", i, fgcp2[i] );
2467 }
2468 #else
2469 for( i=0; i<lgth1; i++ )
2470 {
2471 ogcp1[i] = 0.5 * fpenalty;
2472 fgcp1[i] = 0.5 * fpenalty;
2473 }
2474 for( i=0; i<lgth2; i++ )
2475 {
2476 ogcp2[i] = 0.5 * fpenalty;
2477 fgcp2[i] = 0.5 * fpenalty;
2478 }
2479 #endif
2480
2481 gapinfo[0] = ogcp1;
2482 gapinfo[1] = fgcp1;
2483 gapinfo[2] = ogcp2;
2484 gapinfo[3] = fgcp2;
2485 #endif
2486
2487 #if 0
2488 fprintf( stdout, "in MSalignmm.c\n" );
2489 for( i=0; i<icyc; i++ )
2490 {
2491 fprintf( stdout, ">%d of GROUP1\n", i );
2492 fprintf( stdout, "%s\n", seq1[i] );
2493 }
2494 for( i=0; i<jcyc; i++ )
2495 {
2496 fprintf( stdout, ">%d of GROUP2\n", i );
2497 fprintf( stdout, "%s\n", seq2[i] );
2498 }
2499 fflush( stdout );
2500 #endif
2501
2502 wm = MSalign2m2m_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, 0, lgth1-1, 0, lgth2-1, alloclen, mseq1, mseq2, 0, gapinfo, map );
2503 #if DEBUG
2504 fprintf( stderr, " seq1[0] = %s\n", seq1[0] );
2505 fprintf( stderr, " seq2[0] = %s\n", seq2[0] );
2506 fprintf( stderr, "mseq1[0] = %s\n", mseq1[0] );
2507 fprintf( stderr, "mseq2[0] = %s\n", mseq2[0] );
2508 #endif
2509
2510 // fprintf( stderr, "wm = %f\n", wm );
2511
2512 #if 0
2513
2514 for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
2515 for( i=0; i<jcyc; i++ ) strcpy( seq2[i], mseq2[i] );
2516
2517 if( seqlen( seq1[0] ) != nglen1 )
2518 {
2519 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen1) seqlen(seq1[0])=%d but nglen1=%d\n", seqlen( seq1[0] ), nglen1 );
2520 fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
2521 exit( 1 );
2522 }
2523 if( seqlen( seq2[0] ) != nglen2 )
2524 {
2525 fprintf( stderr, "bug! hairetsu ga kowareta! (nglen2) seqlen(seq2[0])=%d but nglen2=%d\n", seqlen( seq2[0] ), nglen2 );
2526 exit( 1 );
2527 }
2528 #endif
2529
2530 FreeFloatVec( ogcp1 );
2531 FreeFloatVec( ogcp2 );
2532 FreeFloatVec( fgcp1 );
2533 FreeFloatVec( fgcp2 );
2534 FreeFloatMtx( cpmx1 );
2535 FreeFloatMtx( cpmx2 );
2536 free( (void *)gapinfo );
2537
2538 FreeCharMtx( mseq1 );
2539 FreeCharMtx( mseq2 );
2540
2541 lgth1 = strlen( seq1[0] );
2542 lgth2 = strlen( seq2[0] );
2543 for( i=0; i<icyc; i++ )
2544 {
2545 if( strlen( seq1[i] ) != lgth1 )
2546 {
2547 fprintf( stderr, "i = %d / %d\n", i, icyc );
2548 fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2549 exit( 1 );
2550 }
2551 }
2552 for( j=0; j<jcyc; j++ )
2553 {
2554 if( strlen( seq2[j] ) != lgth2 )
2555 {
2556 fprintf( stderr, "j = %d / %d\n", j, jcyc );
2557 fprintf( stderr, "hairetsu ga kowareta (end of MSalignmm) !\n" );
2558 exit( 1 );
2559 }
2560 }
2561
2562 return( wm );
2563 }
2564