8 #define USE_PENALTY_EX 0
9 #define FASTMATCHCALC 1
11 static void OpeningGapCount( float *ogcp, int clus, char **seq, double *eff, int len )
17 for( i=0; i<len; i++ ) ogcp[i] = 0.0;
18 for( j=0; j<clus; j++ )
22 for( i=0; i<len; i++ )
25 gc = ( seq[j][i] == '-' );
27 if( !gb * gc ) ogcp[i] += feff;
32 // for( i=0; i<len; i++ )
34 // fprintf( stderr, "ogcp[%d] = %30.20f\n", i, ogcp[i] );
35 // ogcp[i] /= totaleff;
39 static void FinalGapCount( float *fgcp, int clus, char **seq, double *eff, int len )
44 for( i=0; i<len+1; i++ ) fgcp[i] = 0.0;
45 for( j=0; j<clus; j++ )
48 gc = ( seq[j][0] == '-' );
49 for( i=1; i<len+1; i++ )
52 gc = ( seq[j][i] == '-' );
54 if( gb * !gc ) fgcp[i-1] += feff;
57 // totaleff += eff[j];
58 // fprintf( stderr, "eff[%d] = %30.20f\n", j, eff[j] );
60 // for( i=0; i<len; i++ )
62 // fgcp[i] /= totaleff;
63 // fprintf( stderr, "fgcp[%d] = %30.20f\n", i, fgcp[i] );
66 static float **impmtx = NULL;
67 float imp_match_out_sc( int i1, int j1 )
69 // fprintf( stderr, "imp+match = %f\n", impmtx[i1][j1] * fastathreshold );
70 // fprintf( stderr, "val = %f\n", impmtx[i1][j1] );
71 return( impmtx[i1][j1] );
73 static void imp_match_out_vead_gapmap( float *imp, int i1, int lgth2, int *gapmap2 )
76 float *pt = impmtx[i1];
77 int *gapmappt = gapmap2;
79 *imp++ += pt[*gapmappt++];
82 float *pt = impmtx[i1];
83 for( j=0; j<lgth2; j++ )
84 *imp++ += pt[gapmap2[j]];
88 static void imp_match_out_vead_tate_gapmap( float *imp, int j1, int lgth1, int *gapmap1 )
91 int *gapmappt = gapmap1;
93 *imp++ += impmtx[*gapmappt++][j1];
96 for( i=0; i<lgth1; i++ )
97 *imp++ += impmtx[gapmap1[i]][j1];
101 static void imp_match_out_vead( float *imp, int i1, int lgth2 )
104 float *pt = impmtx[i1];
109 float *pt = impmtx[i1];
110 for( j=0; j<lgth2; j++ )
114 static void imp_match_out_vead_tate( float *imp, int j1, int lgth1 )
117 for( i=0; i<lgth1; i++ )
118 *imp++ += impmtx[i][j1];
121 void imp_match_init_strict( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, LocalHom ***localhom, int forscore )
123 int i, j, k1, k2, tmpint, start1, start2, end1, end2;
124 static int impalloclen = 0;
127 char *pt, *pt1, *pt2;
128 static char *nocount1 = NULL;
129 static char *nocount2 = NULL;
132 if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
134 if( impmtx ) FreeFloatMtx( impmtx );
135 if( nocount1 ) free( nocount1 );
136 if( nocount2 ) free( nocount2 );
137 impalloclen = MAX( lgth1, lgth2 ) + 2;
138 impmtx = AllocateFloatMtx( impalloclen, impalloclen );
139 nocount1 = AllocateCharVec( impalloclen );
140 nocount2 = AllocateCharVec( impalloclen );
143 for( i=0; i<lgth1; i++ )
145 for( j=0; j<clus1; j++ )
146 if( seq1[j][i] == '-' ) break;
147 if( j != clus1 ) nocount1[i] = 1;
148 else nocount1[i] = 0;
150 for( i=0; i<lgth2; i++ )
152 for( j=0; j<clus2; j++ )
153 if( seq2[j][i] == '-' ) break;
154 if( j != clus2 ) nocount2[i] = 1;
155 else nocount2[i] = 0;
159 fprintf( stderr, "nocount2 =\n" );
160 for( i = 0; i<impalloclen; i++ )
162 fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );
169 fprintf( stderr, "eff1 in _init_strict = \n" );
170 for( i=0; i<clus1; i++ )
171 fprintf( stderr, "eff1[] = %f\n", eff1[i] );
172 for( i=0; i<clus2; i++ )
173 fprintf( stderr, "eff2[] = %f\n", eff2[i] );
176 for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
178 effijx = fastathreshold;
179 for( i=0; i<clus1; i++ )
181 for( j=0; j<clus2; j++ )
183 effij = (float)( eff1[i] * eff2[j] * effijx );
184 tmpptr = localhom[i][j];
187 // fprintf( stderr, "start1 = %d\n", tmpptr->start1 );
188 // fprintf( stderr, "end1 = %d\n", tmpptr->end1 );
189 // fprintf( stderr, "i = %d, seq1 = \n%s\n", i, seq1[i] );
190 // fprintf( stderr, "j = %d, seq2 = \n%s\n", j, seq2[j] );
195 if( *pt++ != '-' ) tmpint++;
196 if( tmpint == tmpptr->start1 ) break;
198 start1 = pt - seq1[i] - 1;
200 if( tmpptr->start1 == tmpptr->end1 ) end1 = start1;
206 // fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );
207 if( tmpint == tmpptr->end1 ) break;
208 if( *pt++ != '-' ) tmpint++;
210 end1 = pt - seq1[i] - 0;
214 // fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );
215 if( *pt++ != '-' ) tmpint++;
216 if( tmpint == tmpptr->end1 ) break;
218 end1 = pt - seq1[i] - 1;
226 if( *pt++ != '-' ) tmpint++;
227 if( tmpint == tmpptr->start2 ) break;
229 start2 = pt - seq2[j] - 1;
230 if( tmpptr->start2 == tmpptr->end2 ) end2 = start2;
236 if( tmpint == tmpptr->end2 ) break;
237 if( *pt++ != '-' ) tmpint++;
239 end2 = pt - seq2[j] - 0;
243 if( *pt++ != '-' ) tmpint++;
244 if( tmpint == tmpptr->end2 ) break;
246 end2 = pt - seq2[j] - 1;
249 // fprintf( stderr, "start1 = %d (%c), end1 = %d (%c), start2 = %d (%c), end2 = %d (%c)\n", start1, seq1[i][start1], end1, seq1[i][end1], start2, seq2[j][start2], end2, seq2[j][end2] );
250 // fprintf( stderr, "step 0\n" );
251 if( end1 - start1 != end2 - start2 )
253 // fprintf( stderr, "CHUUI!!, start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
257 k1 = start1; k2 = start2;
260 while( *pt1 && *pt2 )
262 if( *pt1 != '-' && *pt2 != '-' )
264 // ½Å¤ß¤òÆó½Å¤Ë¤«¤±¤Ê¤¤¤è¤¦¤ËÃðդ·¤Æ²¼¤µ¤¤¡£
265 // impmtx[k1][k2] += tmpptr->wimportance * fastathreshold;
266 // impmtx[k1][k2] += tmpptr->importance * effij;
267 impmtx[k1][k2] += tmpptr->fimportance * effij;
268 // fprintf( stderr, "#### impmtx[k1][k2] = %f, tmpptr->fimportance=%f, effij=%f\n", impmtx[k1][k2], tmpptr->fimportance, effij );
269 // fprintf( stderr, "mark, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
270 // fprintf( stderr, "%d (%c) - %d (%c) - %f\n", k1, *pt1, k2, *pt2, tmpptr->fimportance * effij );
274 else if( *pt1 != '-' && *pt2 == '-' )
276 // fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
279 else if( *pt1 == '-' && *pt2 != '-' )
281 // fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
284 else if( *pt1 == '-' && *pt2 == '-' )
286 // fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
290 if( k1 > end1 || k2 > end2 ) break;
293 while( k1 <= end1 && k2 <= end2 )
295 fprintf( stderr, "k1,k2=%d,%d - ", k1, k2 );
296 if( !nocount1[k1] && !nocount2[k2] )
298 impmtx[k1][k2] += tmpptr->wimportance * eff1[i] * eff2[j] * fastathreshold;
299 fprintf( stderr, "marked\n" );
302 fprintf( stderr, "no count\n" );
306 tmpptr = tmpptr->next;
311 if( clus1 == 1 && clus2 == 6 )
313 fprintf( stderr, "\n" );
314 fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
315 fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
316 fprintf( stderr, "impmtx = \n" );
317 for( k2=0; k2<lgth2; k2++ )
318 fprintf( stderr, "%6.3f ", (double)k2 );
319 fprintf( stderr, "\n" );
320 for( k1=0; k1<lgth1; k1++ )
322 fprintf( stderr, "%d ", k1 );
323 for( k2=0; k2<3; k2++ )
324 fprintf( stderr, "%2.1f ", impmtx[k1][k2] );
325 fprintf( stderr, "\n" );
333 void imp_match_init( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, LocalHom ***localhom )
335 int dif, i, j, k1, k2, tmpint, start1, start2, end1, end2;
336 static int impalloclen = 0;
339 static char *nocount1 = NULL;
340 static char *nocount2 = NULL;
342 if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
344 if( impmtx ) FreeFloatMtx( impmtx );
345 if( nocount1 ) free( nocount1 );
346 if( nocount2 ) free( nocount2 );
347 impalloclen = MAX( lgth1, lgth2 ) + 2;
348 impmtx = AllocateFloatMtx( impalloclen, impalloclen );
349 nocount1 = AllocateCharVec( impalloclen );
350 nocount2 = AllocateCharVec( impalloclen );
353 for( i=0; i<lgth1; i++ )
355 for( j=0; j<clus1; j++ )
356 if( seq1[j][i] == '-' ) break;
357 if( j != clus1 ) nocount1[i] = 1;
358 else nocount1[i] = 0;
360 for( i=0; i<lgth2; i++ )
362 for( j=0; j<clus2; j++ )
363 if( seq2[j][i] == '-' ) break;
364 if( j != clus2 ) nocount2[i] = 1;
365 else nocount2[i] = 0;
369 fprintf( stderr, "nocount2 =\n" );
370 for( i = 0; i<impalloclen; i++ )
372 fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );
376 for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
378 for( i=0; i<clus1; i++ )
380 fprintf( stderr, "i = %d, seq1 = %s\n", i, seq1[i] );
381 for( j=0; j<clus2; j++ )
383 fprintf( stderr, "start1 = %d\n", localhom[i][j]->start1 );
384 fprintf( stderr, "end1 = %d\n", localhom[i][j]->end1 );
385 fprintf( stderr, "j = %d, seq2 = %s\n", j, seq2[j] );
390 if( *pt++ != '-' ) tmpint++;
391 if( tmpint == localhom[i][j]->start1 ) break;
393 start1 = pt - seq1[i] - 1;
397 // fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, localhom[i][j].end1, pt-seq1[i] );
398 if( *pt++ != '-' ) tmpint++;
399 if( tmpint == localhom[i][j]->end1 ) break;
401 end1 = pt - seq1[i] - 1;
407 if( *pt++ != '-' ) tmpint++;
408 if( tmpint == localhom[i][j]->start2 ) break;
410 start2 = pt - seq2[j] - 1;
413 if( *pt++ != '-' ) tmpint++;
414 if( tmpint == localhom[i][j]->end2 ) break;
416 end2 = pt - seq2[j] - 1;
417 // fprintf( stderr, "start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
420 fprintf( stderr, "step 0\n" );
421 while( k1 <= end1 && k2 <= end2 )
424 if( !nocount1[k1] && !nocount2[k2] )
425 impmtx[k1][k2] += localhom[i][j].wimportance * eff1[i] * eff2[j];
428 if( !nocount1[k1] && !nocount2[k2] )
429 impmtx[k1][k2] += localhom[i][j]->wimportance * eff1[i] * eff2[j];
434 dif = ( end1 - start1 ) - ( end2 - start2 );
435 fprintf( stderr, "dif = %d\n", dif );
440 fprintf( stderr, "dif = %d\n", dif );
443 while( k1 <= end1 && k2 <= end2 )
445 if( 0 <= k2 && start2 <= k2 && !nocount1[k1] && !nocount2[k2] )
446 impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
460 if( k1 >= 0 && k1 >= start1 && !nocount1[k1] && !nocount2[k2] )
461 impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
470 fprintf( stderr, "impmtx = \n" );
471 for( k2=0; k2<lgth2; k2++ )
472 fprintf( stderr, "%6.3f ", (double)k2 );
473 fprintf( stderr, "\n" );
474 for( k1=0; k1<lgth1; k1++ )
476 fprintf( stderr, "%d", k1 );
477 for( k2=0; k2<lgth2; k2++ )
478 fprintf( stderr, "%6.3f ", impmtx[k1][k2] );
479 fprintf( stderr, "\n" );
485 static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
490 float **cpmxpd = floatwork;
491 int **cpmxpdn = intwork;
492 float *matchpt, *cpmxpdpt, **cpmxpdptpt;
493 int *cpmxpdnpt, **cpmxpdnptpt;
497 for( j=0; j<lgth2; j++ )
500 for( l=0; l<26; l++ )
504 cpmxpd[j][count] = cpmx2[l][j];
505 cpmxpdn[j][count] = l;
509 cpmxpdn[j][count] = -1;
514 for( l=0; l<26; l++ )
517 for( j=0; j<26; j++ )
518 scarr[l] += n_dis[j][l] * cpmx1[j][i1];
521 cpmxpdnptpt = cpmxpdn;
526 cpmxpdnpt = *cpmxpdnptpt++;
527 cpmxpdpt = *cpmxpdptpt++;
528 while( *cpmxpdnpt>-1 )
529 *matchpt += scarr[*cpmxpdnpt++] * *cpmxpdpt++;
536 float **cpmxpd = floatwork;
537 int **cpmxpdn = intwork;
542 for( j=0; j<lgth2; j++ )
545 for( l=0; l<26; l++ )
549 cpmxpd[count][j] = cpmx2[l][j];
550 cpmxpdn[count][j] = l;
554 cpmxpdn[count][j] = -1;
557 for( l=0; l<26; l++ )
560 for( k=0; k<26; k++ )
561 scarr[l] += n_dis[k][l] * cpmx1[k][i1];
563 for( j=0; j<lgth2; j++ )
566 for( k=0; cpmxpdn[k][j]>-1; k++ )
567 match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
572 static void Atracking_localhom( float *impwmpt, float *lasthorizontalw, float *lastverticalw,
573 char **seq1, char **seq2,
574 char **mseq1, char **mseq2,
575 float **cpmx1, float **cpmx2,
576 short **ijp, int icyc, int jcyc )
578 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
581 lgth1 = strlen( seq1[0] );
582 lgth2 = strlen( seq2[0] );
585 for( i=0; i<lgth1; i++ )
587 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
595 wm = lastverticalw[0];
596 for( i=0; i<lgth1; i++ )
598 if( lastverticalw[i] >= wm )
600 wm = lastverticalw[i];
601 iin = i; jin = lgth2-1;
602 ijp[lgth1][lgth2] = +( lgth1 - i );
605 for( j=0; j<lgth2; j++ )
607 if( lasthorizontalw[j] >= wm )
609 wm = lasthorizontalw[j];
610 iin = lgth1-1; jin = j;
611 ijp[lgth1][lgth2] = -( lgth2 - j );
616 for( i=0; i<lgth1+1; i++ )
620 for( j=0; j<lgth2+1; j++ )
622 ijp[0][j] = -( j + 1 );
625 for( i=0; i<icyc; i++ )
627 mseq1[i] += lgth1+lgth2;
630 for( j=0; j<jcyc; j++ )
632 mseq2[j] += lgth1+lgth2;
635 iin = lgth1; jin = lgth2;
637 for( k=0; k<=lgth1+lgth2; k++ )
639 if( ijp[iin][jin] < 0 )
641 ifi = iin-1; jfi = jin+ijp[iin][jin];
643 else if( ijp[iin][jin] > 0 )
645 ifi = iin-ijp[iin][jin]; jfi = jin-1;
649 ifi = iin-1; jfi = jin-1;
654 for( i=0; i<icyc; i++ )
655 *--mseq1[i] = seq1[i][ifi+l];
656 for( j=0; j<jcyc; j++ )
663 for( i=0; i<icyc; i++ )
665 for( j=0; j<jcyc; j++ )
666 *--mseq2[j] = seq2[j][jfi+l];
669 if( iin == lgth1 || jin == lgth2 )
673 *impwmpt += imp_match_out_sc( iin, jin );
675 // fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
677 if( iin <= 0 || jin <= 0 ) break;
678 for( i=0; i<icyc; i++ )
679 *--mseq1[i] = seq1[i][ifi];
680 for( j=0; j<jcyc; j++ )
681 *--mseq2[j] = seq2[j][jfi];
683 iin = ifi; jin = jfi;
686 static void Atracking_localhom_gapmap( float *impwmpt, float *lasthorizontalw, float *lastverticalw,
687 char **seq1, char **seq2,
688 char **mseq1, char **mseq2,
689 float **cpmx1, float **cpmx2,
690 short **ijp, int icyc, int jcyc,
691 int *gapmap1, int *gapmap2 )
693 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
696 lgth1 = strlen( seq1[0] );
697 lgth2 = strlen( seq2[0] );
700 for( i=0; i<lgth1; i++ )
702 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
710 wm = lastverticalw[0];
711 for( i=0; i<lgth1; i++ )
713 if( lastverticalw[i] >= wm )
715 wm = lastverticalw[i];
716 iin = i; jin = lgth2-1;
717 ijp[lgth1][lgth2] = +( lgth1 - i );
720 for( j=0; j<lgth2; j++ )
722 if( lasthorizontalw[j] >= wm )
724 wm = lasthorizontalw[j];
725 iin = lgth1-1; jin = j;
726 ijp[lgth1][lgth2] = -( lgth2 - j );
731 for( i=0; i<lgth1+1; i++ )
735 for( j=0; j<lgth2+1; j++ )
737 ijp[0][j] = -( j + 1 );
740 for( i=0; i<icyc; i++ )
742 mseq1[i] += lgth1+lgth2;
745 for( j=0; j<jcyc; j++ )
747 mseq2[j] += lgth1+lgth2;
750 iin = lgth1; jin = lgth2;
752 for( k=0; k<=lgth1+lgth2; k++ )
754 if( ijp[iin][jin] < 0 )
756 ifi = iin-1; jfi = jin+ijp[iin][jin];
758 else if( ijp[iin][jin] > 0 )
760 ifi = iin-ijp[iin][jin]; jfi = jin-1;
764 ifi = iin-1; jfi = jin-1;
769 for( i=0; i<icyc; i++ )
770 *--mseq1[i] = seq1[i][ifi+l];
771 for( j=0; j<jcyc; j++ )
778 for( i=0; i<icyc; i++ )
780 for( j=0; j<jcyc; j++ )
781 *--mseq2[j] = seq2[j][jfi+l];
784 if( iin == lgth1 || jin == lgth2 )
788 *impwmpt += imp_match_out_sc( gapmap1[iin], gapmap2[jin] );
790 // fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
792 if( iin <= 0 || jin <= 0 ) break;
793 for( i=0; i<icyc; i++ )
794 *--mseq1[i] = seq1[i][ifi];
795 for( j=0; j<jcyc; j++ )
796 *--mseq2[j] = seq2[j][jfi];
798 iin = ifi; jin = jfi;
801 static float Atracking( float *lasthorizontalw, float *lastverticalw,
802 char **seq1, char **seq2,
803 char **mseq1, char **mseq2,
804 float **cpmx1, float **cpmx2,
805 short **ijp, int icyc, int jcyc )
807 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
810 lgth1 = strlen( seq1[0] );
811 lgth2 = strlen( seq2[0] );
814 for( i=0; i<lgth1; i++ )
816 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
824 wm = lastverticalw[0];
825 for( i=0; i<lgth1; i++ )
827 if( lastverticalw[i] >= wm )
829 wm = lastverticalw[i];
830 iin = i; jin = lgth2-1;
831 ijp[lgth1][lgth2] = +( lgth1 - i );
834 for( j=0; j<lgth2; j++ )
836 if( lasthorizontalw[j] >= wm )
838 wm = lasthorizontalw[j];
839 iin = lgth1-1; jin = j;
840 ijp[lgth1][lgth2] = -( lgth2 - j );
845 for( i=0; i<lgth1+1; i++ )
849 for( j=0; j<lgth2+1; j++ )
851 ijp[0][j] = -( j + 1 );
854 for( i=0; i<icyc; i++ )
856 mseq1[i] += lgth1+lgth2;
859 for( j=0; j<jcyc; j++ )
861 mseq2[j] += lgth1+lgth2;
864 iin = lgth1; jin = lgth2;
865 for( k=0; k<=lgth1+lgth2; k++ )
867 if( ijp[iin][jin] < 0 )
869 ifi = iin-1; jfi = jin+ijp[iin][jin];
871 else if( ijp[iin][jin] > 0 )
873 ifi = iin-ijp[iin][jin]; jfi = jin-1;
877 ifi = iin-1; jfi = jin-1;
882 for( i=0; i<icyc; i++ )
883 *--mseq1[i] = seq1[i][ifi+l];
884 for( j=0; j<jcyc; j++ )
891 for( i=0; i<icyc; i++ )
893 for( j=0; j<jcyc; j++ )
894 *--mseq2[j] = seq2[j][jfi+l];
897 if( iin <= 0 || jin <= 0 ) break;
898 for( i=0; i<icyc; i++ )
899 *--mseq1[i] = seq1[i][ifi];
900 for( j=0; j<jcyc; j++ )
901 *--mseq2[j] = seq2[j][jfi];
903 iin = ifi; jin = jfi;
908 float A__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch )
909 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
913 int lasti; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
916 float wm = 0.0; /* int ?????? */
918 float *currentw, *previousw;
919 // float fpenalty = (float)penalty;
921 float fpenalty_ex = (float)penalty_ex;
926 float *mjpt, *prept, *curpt;
932 static float *w1, *w2;
934 static float *initverticalw; /* kufuu sureba iranai */
935 static float *lastverticalw; /* kufuu sureba iranai */
943 static float **cpmx1;
944 static float **cpmx2;
945 static int **intwork;
946 static float **floatwork;
947 static int orlgth1 = 0, orlgth2 = 0;
948 float fpenalty = (float)penalty;
952 fprintf( stderr, "eff in SA+++align\n" );
953 for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
957 mseq1 = AllocateCharMtx( njob, 0 );
958 mseq2 = AllocateCharMtx( njob, 0 );
962 lgth1 = strlen( seq1[0] );
963 lgth2 = strlen( seq2[0] );
965 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
969 if( orlgth1 > 0 && orlgth2 > 0 )
973 FreeFloatVec( match );
974 FreeFloatVec( initverticalw );
975 FreeFloatVec( lastverticalw );
982 FreeFloatVec( ogcp1 );
983 FreeFloatVec( ogcp2 );
984 FreeFloatVec( fgcp1 );
985 FreeFloatVec( fgcp2 );
988 FreeFloatMtx( cpmx1 );
989 FreeFloatMtx( cpmx2 );
991 FreeFloatMtx( floatwork );
992 FreeIntMtx( intwork );
995 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
996 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
999 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
1002 w1 = AllocateFloatVec( ll2+2 );
1003 w2 = AllocateFloatVec( ll2+2 );
1004 match = AllocateFloatVec( ll2+2 );
1006 initverticalw = AllocateFloatVec( ll1+2 );
1007 lastverticalw = AllocateFloatVec( ll1+2 );
1009 m = AllocateFloatVec( ll2+2 );
1010 mp = AllocateIntVec( ll2+2 );
1012 mseq = AllocateCharMtx( njob, ll1+ll2 );
1014 ogcp1 = AllocateFloatVec( ll1+2 );
1015 ogcp2 = AllocateFloatVec( ll2+2 );
1016 fgcp1 = AllocateFloatVec( ll1+2 );
1017 fgcp2 = AllocateFloatVec( ll2+2 );
1019 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
1020 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
1023 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
1024 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
1026 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
1027 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
1031 fprintf( stderr, "succeeded\n" );
1034 orlgth1 = ll1 - 100;
1035 orlgth2 = ll2 - 100;
1039 for( i=0; i<icyc; i++ )
1044 for( j=0; j<jcyc; j++ )
1046 mseq2[j] = mseq[icyc+j];
1051 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
1055 if( commonAlloc1 && commonAlloc2 )
1057 FreeShortMtx( commonIP );
1060 ll1 = MAX( orlgth1, commonAlloc1 );
1061 ll2 = MAX( orlgth2, commonAlloc2 );
1064 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
1067 commonIP = AllocateShortMtx( ll1+10, ll2+10 );
1070 fprintf( stderr, "succeeded\n\n" );
1081 for( i=0; i<icyc; i++ )
1083 fprintf( stderr, "## totaleff = %f\n", t );
1086 cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
1087 cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
1089 OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
1090 OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1091 FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1092 FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1094 for( i=0; i<lgth1; i++ )
1096 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty;
1097 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty;
1099 for( i=0; i<lgth2; i++ )
1101 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty;
1102 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty;
1105 for( i=0; i<lgth1; i++ )
1106 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
1112 match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
1114 imp_match_out_vead_tate( initverticalw, 0, lgth1 ); // 060306
1116 match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
1118 imp_match_out_vead( currentw, 0, lgth2 ); // 060306
1119 #if 0 // -> tbfast.c
1121 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
1127 for( i=1; i<lgth1+1; i++ )
1129 initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
1131 for( j=1; j<lgth2+1; j++ )
1133 currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
1139 for( j=1; j<lgth2+1; j++ )
1140 currentw[j] -= offset * j / 2.0;
1141 for( i=1; i<lgth1+1; i++ )
1142 initverticalw[i] -= offset * i / 2.0;
1146 for( j=1; j<lgth2+1; ++j )
1148 m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0;
1151 lastverticalw[0] = currentw[lgth2-1];
1153 if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1156 fprintf( stderr, "currentw = \n" );
1157 for( i=0; i<lgth1+1; i++ )
1159 fprintf( stderr, "%5.2f ", currentw[i] );
1161 fprintf( stderr, "\n" );
1162 fprintf( stderr, "initverticalw = \n" );
1163 for( i=0; i<lgth2+1; i++ )
1165 fprintf( stderr, "%5.2f ", initverticalw[i] );
1167 fprintf( stderr, "\n" );
1168 fprintf( stderr, "fcgp\n" );
1169 for( i=0; i<lgth1; i++ )
1170 fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1171 for( i=0; i<lgth2; i++ )
1172 fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1175 for( i=1; i<lasti; i++ )
1178 previousw = currentw;
1181 previousw[0] = initverticalw[i-1];
1183 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1185 fprintf( stderr, "\n" );
1186 fprintf( stderr, "i=%d\n", i );
1187 fprintf( stderr, "currentw = \n" );
1188 for( j=0; j<lgth2; j++ )
1190 fprintf( stderr, "%5.2f ", currentw[j] );
1192 fprintf( stderr, "\n" );
1196 // fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1198 imp_match_out_vead( currentw, i, lgth2 );
1200 imp_match_out_vead( currentw, i, lgth2 );
1204 fprintf( stderr, "\n" );
1205 fprintf( stderr, "i=%d\n", i );
1206 fprintf( stderr, "currentw = \n" );
1207 for( j=0; j<lgth2; j++ )
1209 fprintf( stderr, "%5.2f ", currentw[j] );
1211 fprintf( stderr, "\n" );
1213 currentw[0] = initverticalw[i];
1217 mi = previousw[0] + ogcp2[1]; mpi = 0;
1222 curpt = currentw + 1;
1224 for( j=1; j<lgth2+1; j++ )
1230 fprintf( stderr, "%5.0f->", wm );
1232 g = mi + fgcp2[j-1];
1234 fprintf( stderr, "%5.0f?", g );
1239 *ijppt = -( j - mpi );
1241 g = *prept + ogcp2[j];
1251 g = *mjpt + fgcp1[i-1];
1253 fprintf( stderr, "%5.0f?", g );
1258 *ijppt = +( i - *mpjpt );
1260 g = *prept + ogcp1[i];
1267 m[j] += fpenalty_ex;
1271 fprintf( stderr, "%5.0f ", wm );
1280 lastverticalw[i] = currentw[lgth2-1];
1282 floatncpy( previousw, currentw, lgth2+1 );
1283 previousw[0] = initverticalw[i-1];
1286 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1287 currentw[0] = initverticalw[i];
1289 mi = previousw[0] + ogcp2[1]; mpi = 0;
1290 for( j=1; j<lgth2+1; j++ )
1292 wm = previousw[j-1];
1299 ijp[i][j] = -( j - mpi );
1302 if( mi <= previousw[j-1] + g )
1304 mi = previousw[j-1] + g;
1312 ijp[i][j] = +( i - mp[j] );
1315 if( m[j] <= previousw[j-1] + g )
1317 m[j] = previousw[j-1] + g ;
1322 lastverticalw[i] = currentw[lgth2-1];
1327 // fprintf( stderr, "wm = %f\n", wm );
1332 for( j=1; j<lgth2+1; j++ )
1333 currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1334 for( i=1; i<lgth1+1; i++ )
1335 lastverticalw[i] -= offset * ( lgth1 - i / 2.0);
1340 fprintf( stderr, "\n" );
1341 for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1342 fprintf( stderr, "#####\n" );
1343 for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1344 fprintf( stderr, "====>" );
1345 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1346 for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1350 Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1353 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1355 // fprintf( stderr, "### impmatch = %f\n", *impmatch );
1357 resultlen = strlen( mseq1[0] );
1358 if( alloclen < resultlen || resultlen > N )
1360 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1361 ErrorExit( "LENGTH OVER!\n" );
1365 for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1366 for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1368 fprintf( stderr, "\n" );
1369 for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1370 fprintf( stderr, "#####\n" );
1371 for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1374 // fprintf( stderr, "wm = %f\n", wm );
1379 float A__align_gapmap( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch, int *gapmap1, int *gapmap2 )
1380 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1384 int lasti; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
1387 float wm = 0.0; /* int ?????? */
1389 float *currentw, *previousw;
1390 // float fpenalty = (float)penalty;
1392 float fpenalty_ex = (float)penalty_ex;
1397 float *mjpt, *prept, *curpt;
1400 static float mi, *m;
1402 static int mpi, *mp;
1403 static float *w1, *w2;
1404 static float *match;
1405 static float *initverticalw; /* kufuu sureba iranai */
1406 static float *lastverticalw; /* kufuu sureba iranai */
1407 static char **mseq1;
1408 static char **mseq2;
1410 static float *ogcp1;
1411 static float *ogcp2;
1412 static float *fgcp1;
1413 static float *fgcp2;
1414 static float **cpmx1;
1415 static float **cpmx2;
1416 static int **intwork;
1417 static float **floatwork;
1418 static int orlgth1 = 0, orlgth2 = 0;
1426 fprintf( stderr, "eff in SA+++align\n" );
1427 for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
1431 mseq1 = AllocateCharMtx( njob, 0 );
1432 mseq2 = AllocateCharMtx( njob, 0 );
1436 lgth1 = strlen( seq1[0] );
1437 lgth2 = strlen( seq2[0] );
1439 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
1443 if( orlgth1 > 0 && orlgth2 > 0 )
1447 FreeFloatVec( match );
1448 FreeFloatVec( initverticalw );
1449 FreeFloatVec( lastverticalw );
1454 FreeCharMtx( mseq );
1456 FreeFloatVec( ogcp1 );
1457 FreeFloatVec( ogcp2 );
1458 FreeFloatVec( fgcp1 );
1459 FreeFloatVec( fgcp2 );
1462 FreeFloatMtx( cpmx1 );
1463 FreeFloatMtx( cpmx2 );
1465 FreeFloatMtx( floatwork );
1466 FreeIntMtx( intwork );
1469 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
1470 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
1473 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
1476 w1 = AllocateFloatVec( ll2+2 );
1477 w2 = AllocateFloatVec( ll2+2 );
1478 match = AllocateFloatVec( ll2+2 );
1480 initverticalw = AllocateFloatVec( ll1+2 );
1481 lastverticalw = AllocateFloatVec( ll1+2 );
1483 m = AllocateFloatVec( ll2+2 );
1484 mp = AllocateIntVec( ll2+2 );
1486 mseq = AllocateCharMtx( njob, ll1+ll2 );
1488 ogcp1 = AllocateFloatVec( ll1+2 );
1489 ogcp2 = AllocateFloatVec( ll2+2 );
1490 fgcp1 = AllocateFloatVec( ll1+2 );
1491 fgcp2 = AllocateFloatVec( ll2+2 );
1493 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
1494 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
1497 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
1498 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
1500 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
1501 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
1505 fprintf( stderr, "succeeded\n" );
1508 orlgth1 = ll1 - 100;
1509 orlgth2 = ll2 - 100;
1513 for( i=0; i<icyc; i++ )
1518 for( j=0; j<jcyc; j++ )
1520 mseq2[j] = mseq[icyc+j];
1525 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
1529 if( commonAlloc1 && commonAlloc2 )
1531 FreeShortMtx( commonIP );
1534 ll1 = MAX( orlgth1, commonAlloc1 );
1535 ll2 = MAX( orlgth2, commonAlloc2 );
1538 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
1541 commonIP = AllocateShortMtx( ll1+10, ll2+10 );
1544 fprintf( stderr, "succeeded\n\n" );
1552 cpmx_calc_new( seq1, cpmx1, eff1, strlen( seq1[0] ), icyc );
1553 cpmx_calc_new( seq2, cpmx2, eff2, strlen( seq2[0] ), jcyc );
1555 OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
1556 OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1557 FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1558 FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1560 for( i=0; i<lgth1; i++ )
1562 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] );
1563 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] );
1565 for( i=0; i<lgth2; i++ )
1567 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] );
1568 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] );
1571 for( i=0; i<lgth1; i++ )
1572 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
1579 match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
1581 imp_match_out_vead_tate_gapmap( initverticalw, gapmap2[0], lgth1, gapmap1 ); // 060306
1584 match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
1586 imp_match_out_vead_gapmap( currentw, gapmap1[0], lgth2, gapmap2 ); // 060306
1587 #if 0 // -> tbfast.c
1589 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
1595 for( i=1; i<lgth1+1; i++ )
1597 initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
1599 for( j=1; j<lgth2+1; j++ )
1601 currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
1607 for( j=1; j<lgth2+1; j++ )
1608 currentw[j] -= offset * j / 2.0;
1609 for( i=1; i<lgth1+1; i++ )
1610 initverticalw[i] -= offset * i / 2.0;
1614 for( j=1; j<lgth2+1; ++j )
1616 m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0;
1619 lastverticalw[0] = currentw[lgth2-1];
1621 if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1624 fprintf( stderr, "currentw = \n" );
1625 for( i=0; i<lgth1+1; i++ )
1627 fprintf( stderr, "%5.2f ", currentw[i] );
1629 fprintf( stderr, "\n" );
1630 fprintf( stderr, "initverticalw = \n" );
1631 for( i=0; i<lgth2+1; i++ )
1633 fprintf( stderr, "%5.2f ", initverticalw[i] );
1635 fprintf( stderr, "\n" );
1636 fprintf( stderr, "fcgp\n" );
1637 for( i=0; i<lgth1; i++ )
1638 fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1639 for( i=0; i<lgth2; i++ )
1640 fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1643 for( i=1; i<lasti; i++ )
1646 previousw = currentw;
1649 previousw[0] = initverticalw[i-1];
1651 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1653 fprintf( stderr, "\n" );
1654 fprintf( stderr, "i=%d\n", i );
1655 fprintf( stderr, "currentw = \n" );
1656 for( j=0; j<lgth2; j++ )
1658 fprintf( stderr, "%5.2f ", currentw[j] );
1660 fprintf( stderr, "\n" );
1664 // fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1666 imp_match_out_vead( currentw, i, lgth2 );
1668 imp_match_out_vead_gapmap( currentw, gapmap1[i], lgth2, gapmap2 );
1672 fprintf( stderr, "\n" );
1673 fprintf( stderr, "i=%d\n", i );
1674 fprintf( stderr, "currentw = \n" );
1675 for( j=0; j<lgth2; j++ )
1677 fprintf( stderr, "%5.2f ", currentw[j] );
1679 fprintf( stderr, "\n" );
1681 currentw[0] = initverticalw[i];
1684 mi = previousw[0] + ogcp2[1]; mpi = 0;
1689 curpt = currentw + 1;
1692 ogcp2pt = ogcp2 + 1;
1693 fgcp1va = fgcp1[i-1];
1695 for( j=1; j<lgth2+1; j++ )
1701 fprintf( stderr, "%5.0f->", wm );
1705 fprintf( stderr, "%5.0f?", g );
1710 *ijppt = -( j - mpi );
1712 g = *prept + *ogcp2pt;
1722 g = *mjpt + fgcp1va;
1724 fprintf( stderr, "%5.0f?", g );
1729 *ijppt = +( i - *mpjpt );
1731 g = *prept + ogcp1va;
1738 m[j] += fpenalty_ex;
1742 fprintf( stderr, "%5.0f ", wm );
1753 lastverticalw[i] = currentw[lgth2-1];
1759 for( j=1; j<lgth2+1; j++ )
1760 currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1761 for( i=1; i<lgth1+1; i++ )
1762 lastverticalw[i] -= offset * ( lgth1 - i / 2.0);
1767 fprintf( stderr, "\n" );
1768 for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1769 fprintf( stderr, "#####\n" );
1770 for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1771 fprintf( stderr, "====>" );
1772 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1773 for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1777 Atracking_localhom_gapmap( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, gapmap1, gapmap2 );
1780 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1782 // fprintf( stderr, "### impmatch = %f\n", *impmatch );
1784 resultlen = strlen( mseq1[0] );
1785 if( alloclen < resultlen || resultlen > N )
1787 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1788 ErrorExit( "LENGTH OVER!\n" );
1792 for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1793 for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1795 fprintf( stderr, "\n" );
1796 for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1797 fprintf( stderr, "#####\n" );
1798 for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1801 // fprintf( stderr, "wm = %f\n", wm );
1806 float translate_and_Calign( char **mseq1, char **mseq2, double *effarr1, double *effarr2, int clus1, int clus2, int alloclen )
1816 if ( clus1 == 1 && clus2 != 1 )
1818 seq = mseq1[0]; aseq = mseq2; effarr = effarr2; nseq = clus2+1;
1820 printf( "effarr in transl... = \n" );
1821 for( i=0; i<clus2; i++ ) printf( "%f ", effarr2[i] );
1824 else if( clus1 != 1 && clus2 == 1 )
1826 seq = mseq2[0]; aseq = mseq1; effarr = effarr1; nseq = clus1+1;
1828 else ErrorExit( "ERROR in translate_and_Calign" );
1830 result = Calignm1( &wm, aseq, seq, effarr, nseq-2, 0 );
1832 resultlen = strlen( result[0] );
1833 if( alloclen < resultlen || resultlen > N )
1835 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1836 ErrorExit( "LENGTH OVER!\n" );
1838 for( i=0; i<nseq-1; i++ ) strcpy( aseq[i], result[i] );
1839 strcpy( seq, result[nseq-1] );