binaries/src/mafft/core/Falign.c

   1 #include "mltaln.h"
   2
   3 #if 0
   4 static FILE *fftfp;
   5 #endif
   6 static TLS int n20or4or2;
   7
   8 #define KEIKA 0
   9 #define RND   0
  10 #define DEBUG 0
  11
  12 #if RND // by D.Mathog
  13 static void generateRndSeq( char *seq, int len )
  14 {
  15         while( len-- )
  16 #if 1
  17                 *seq++ = (int)( rnd() * n20or4or2 );
  18 #else
  19                 *seq++ = (int)1;
  20 #endif
  21 }
  22 #endif
  23
  24 static void vec_init( Fukusosuu *result, int nlen )
  25 {
  26         while( nlen-- )
  27         {
  28                 result->R = result->I = 0.0;
  29                 result++;
  30         }
  31 }
  32
  33 #if 0 // by D.Mathog
  34 static void vec_init2( Fukusosuu **result, char *seq, double eff, int st, int ed )
  35 {
  36         int i;
  37         for( i=st; i<ed; i++ )
  38                 result[(int)*seq++][i].R += eff;
  39 }
  40 #endif
  41
  42 static void seq_vec_2( Fukusosuu *result, double *score, double incr, char *seq )
  43 {
  44         static TLS int n;
  45         for( ; *seq; result++ )
  46         {
  47                 n = amino_n[(int)*seq++];
  48                 if( n < 20 && n >= 0 ) result->R += incr * score[n];
  49 #if 0
  50                 fprintf( stderr, "n=%d, score=%f, inc=%f R=%f\n",n,  score[n], incr * score[n], result->R );
  51 #endif
  52         }
  53 }
  54
  55 static void seq_vec_3( Fukusosuu **result, double incr, char *seq )
  56 {
  57         int i;
  58         int n;
  59         for( i=0; *seq; i++ )
  60         {
  61                 n = amino_n[(int)*seq++];
  62                 if( n < n20or4or2 && n >= 0 ) result[n][i].R += incr;
  63         }
  64 }
  65
  66 static void seq_vec_5( Fukusosuu *result, double *score1, double *score2, double incr, char *seq )
  67 {
  68         int n;
  69         for( ; *seq; result++ )
  70         {
  71                 n = amino_n[(int)*seq++];
  72                 if( n > 20 ) continue;
  73                 result->R += incr * score1[n];
  74                 result->I += incr * score2[n];
  75 #if 0
  76                 fprintf( stderr, "n=%d, score=%f, inc=%f R=%f\n",n,  score[n], incr * score[n], result->R );
  77 #endif
  78         }
  79 }
  80
  81
  82 static void seq_vec_4( Fukusosuu *result, double incr, char *seq )
  83 {
  84         char s;
  85         for( ; *seq; result++ )
  86         {
  87                 s = *seq++;
  88                 if( s == 'a' )
  89                         result->R += incr;
  90                 else if( s == 't' )
  91                         result->R -= incr;
  92                 else if( s == 'g' )
  93                         result->I += incr;
  94                 else if( s == 'c' )
  95                         result->I -= incr;
  96         }
  97 }
  98
  99 #if 0 // by D.Mathog
 100 static void seq_vec( Fukusosuu *result, char query, double incr, char *seq )
 101 {
 102 #if 0
 103         int bk = nlen;
 104 #endif
 105         while( *seq )
 106         {
 107                 if( *seq++ == query ) result->R += incr;
 108                 result++;
 109 #if 0
 110 fprintf( stderr, "i = %d result->R = %f\n", bk-nlen, (result-1)->R );
 111 #endif
 112         }
 113 }
 114
 115 static int checkRepeat( int num, int *cutpos )
 116 {
 117         int tmp, buf;
 118
 119         buf = *cutpos;
 120         while( num-- )
 121         {
 122                 if( ( tmp = *cutpos++ ) < buf ) return( 1 );
 123                 buf = tmp;
 124         }
 125         return( 0 );
 126 }
 127
 128 static int segcmp( void *ptr1, void *ptr2 )
 129 {
 130         int diff;
 131         Segment **seg1 = (Segment **)ptr1;
 132         Segment **seg2 = (Segment **)ptr2;
 133 #if 0
 134         return( (*seg1)->center - (*seg2)->center );
 135 #else
 136         diff = (*seg1)->center - (*seg2)->center;
 137         if( diff ) return( diff );
 138
 139         diff = (*seg1)->start - (*seg2)->start;
 140         if( diff ) return( diff );
 141
 142         diff = (*seg1)->end - (*seg2)->end;
 143         if( diff ) return( diff );
 144
 145         fprintf( stderr, "USE STABLE SORT !!\n" );
 146         exit( 1 );
 147         return( 0 );
 148 #endif
 149 }
 150 #endif
 151
 152
 153 static void mymergesort( int first, int last, Segment **seg )
 154 {
 155         int middle;
 156         static TLS int i, j, k, p;
 157         static TLS int allo = 0;
 158         static TLS Segment **work = NULL;
 159
 160         if( seg == NULL )
 161         {
 162                 free( work ); work = NULL;
 163                 return;
 164         }
 165
 166         if( last > allo )
 167         {
 168                 allo = last;
 169                 if( work ) free( work );
 170                 work = (Segment **)calloc( allo / 2 + 1, sizeof( Segment *) );
 171         }
 172
 173         if( first < last )
 174         {
 175                 middle = ( first + last ) / 2;
 176                 mymergesort( first, middle, seg );
 177                 mymergesort( middle+1, last, seg );
 178                 p = 0;
 179                 for( i=first; i<=middle; i++ ) work[p++] = seg[i];
 180                 i = middle + 1; j = 0; k = first;
 181                 while( i <= last && j < p )
 182                 {
 183                         if( work[j]->center <= seg[i]->center )
 184                                 seg[k++] = work[j++];
 185                         else
 186                                 seg[k++] = seg[i++];
 187                 }
 188                 while( j < p ) seg[k++] = work[j++];
 189         }
 190 }
 191
 192
 193 double Fgetlag( char  **seq1, char  **seq2,
 194                             double *eff1, double *eff2,
 195                             int    clus1, int    clus2,
 196                             int alloclen )
 197 {
 198         int i, j, k, l, m;
 199         int nlen, nlen2, nlen4;
 200         static TLS int crossscoresize = 0;
 201         static TLS char **tmpseq1 = NULL;
 202         static TLS char **tmpseq2 = NULL;
 203         static TLS char **tmpptr1 = NULL;
 204         static TLS char **tmpptr2 = NULL;
 205         static TLS char **tmpres1 = NULL;
 206         static TLS char **tmpres2 = NULL;
 207         static TLS char **result1 = NULL;
 208         static TLS char **result2 = NULL;
 209 #if RND
 210         static TLS char **rndseq1 = NULL;
 211         static TLS char **rndseq2 = NULL;
 212 #endif
 213         static TLS Fukusosuu **seqVector1 = NULL;
 214         static TLS Fukusosuu **seqVector2 = NULL;
 215         static TLS Fukusosuu **naiseki = NULL;
 216         static TLS Fukusosuu *naisekiNoWa = NULL;
 217         static TLS double *soukan = NULL;
 218         static TLS double **crossscore = NULL;
 219         int nlentmp;
 220         static TLS int *kouho = NULL;
 221         static TLS Segment *segment = NULL;
 222         static TLS Segment *segment1 = NULL;
 223         static TLS Segment *segment2 = NULL;
 224         static TLS Segment **sortedseg1 = NULL;
 225         static TLS Segment **sortedseg2 = NULL;
 226         static TLS int *cut1 = NULL;
 227         static TLS int *cut2 = NULL;
 228         static TLS int localalloclen = 0;
 229         int lag;
 230         int tmpint;
 231         int count, count0;
 232         int len1, len2;
 233         int totallen;
 234         float dumfl = 0.0;
 235         int headgp, tailgp;
 236
 237         len1 = strlen( seq1[0] );
 238         len2 = strlen( seq2[0] );
 239         nlentmp = MAX( len1, len2 );
 240
 241         nlen = 1;
 242         while( nlentmp >= nlen ) nlen <<= 1;
 243 #if 0
 244         fprintf( stderr, "###   nlen    = %d\n", nlen );
 245 #endif
 246
 247         nlen2 = nlen/2; nlen4 = nlen2 / 2;
 248
 249 #if DEBUG
 250         fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
 251         fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
 252 #endif
 253
 254         if( !localalloclen )
 255         {
 256                 kouho = AllocateIntVec( NKOUHO );
 257                 cut1 = AllocateIntVec( MAXSEG );
 258                 cut2 = AllocateIntVec( MAXSEG );
 259                 tmpptr1 = AllocateCharMtx( njob, 0 );
 260                 tmpptr2 = AllocateCharMtx( njob, 0 );
 261                 result1 = AllocateCharMtx( njob, alloclen );
 262                 result2 = AllocateCharMtx( njob, alloclen );
 263                 tmpres1 = AllocateCharMtx( njob, alloclen );
 264                 tmpres2 = AllocateCharMtx( njob, alloclen );
 265 //              crossscore = AllocateDoubleMtx( MAXSEG, MAXSEG );
 266                 segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
 267                 segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
 268                 segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
 269                 sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
 270                 sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
 271                 if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
 272                         ErrorExit( "Allocation error\n" );
 273
 274                 if     ( scoremtx == -1 ) n20or4or2 = 4;
 275                 else if( fftscore == 1  ) n20or4or2 = 2;
 276                 else                      n20or4or2 = 20;
 277         }
 278         if( localalloclen < nlen )
 279         {
 280                 if( localalloclen )
 281                 {
 282 #if 1
 283                         FreeFukusosuuMtx ( seqVector1 );
 284                         FreeFukusosuuMtx ( seqVector2 );
 285                         FreeFukusosuuVec( naisekiNoWa );
 286                         FreeFukusosuuMtx( naiseki );
 287                         FreeDoubleVec( soukan );
 288                         FreeCharMtx( tmpseq1 );
 289                         FreeCharMtx( tmpseq2 );
 290 #endif
 291 #if RND
 292                         FreeCharMtx( rndseq1 );
 293                         FreeCharMtx( rndseq2 );
 294 #endif
 295                 }
 296
 297
 298                 tmpseq1 = AllocateCharMtx( njob, nlen );
 299                 tmpseq2 = AllocateCharMtx( njob, nlen );
 300                 naisekiNoWa = AllocateFukusosuuVec( nlen );
 301                 naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
 302                 seqVector1 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
 303                 seqVector2 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
 304                 soukan = AllocateDoubleVec( nlen+1 );
 305
 306 #if RND
 307                 rndseq1 = AllocateCharMtx( njob, nlen );
 308                 rndseq2 = AllocateCharMtx( njob, nlen );
 309                 for( i=0; i<njob; i++ )
 310                 {
 311                         generateRndSeq( rndseq1[i], nlen );
 312                         generateRndSeq( rndseq2[i], nlen );
 313                 }
 314 #endif
 315                 localalloclen = nlen;
 316         }
 317
 318         for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
 319         for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
 320
 321 #if 0
 322 fftfp = fopen( "input_of_Falign", "w" );
 323 fprintf( fftfp, "nlen = %d\n", nlen );
 324 fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
 325 for( i=0; i<clus1; i++ )
 326         fprintf( fftfp, "%s\n", seq1[i] );
 327 fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
 328 for( i=0; i<clus2; i++ )
 329         fprintf( fftfp, "%s\n", seq2[i] );
 330 fclose( fftfp );
 331 system( "less input_of_Falign < /dev/tty > /dev/tty" );
 332 #endif
 333
 334         if( fftkeika ) fprintf( stderr,  " FFT ... " );
 335
 336         for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
 337         if( fftscore && scoremtx != -1 )
 338         {
 339                 for( i=0; i<clus1; i++ )
 340                 {
 341                         seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
 342                         seq_vec_2( seqVector1[1], volume,   eff1[i], tmpseq1[i] );
 343                 }
 344         }
 345         else
 346         {
 347 #if 0
 348                 for( i=0; i<clus1; i++ ) for( j=0; j<n20or4or2; j++ )
 349                         seq_vec( seqVector1[j], amino[j], eff1[i], tmpseq1[i] );
 350 #else
 351                 for( i=0; i<clus1; i++ )
 352                         seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
 353 #endif
 354         }
 355 #if RND
 356         for( i=0; i<clus1; i++ )
 357         {
 358                 vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
 359         }
 360 #endif
 361 #if 0
 362 fftfp = fopen( "seqVec", "w" );
 363 fprintf( fftfp, "before transform\n" );
 364 for( k=0; k<n20or4or2; k++ )
 365 {
 366    fprintf( fftfp, "nlen=%d\n", nlen );
 367    fprintf( fftfp, "%c\n", amino[k] );
 368    for( l=0; l<nlen; l++ )
 369    fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
 370 }
 371 fclose( fftfp );
 372 system( "less seqVec < /dev/tty > /dev/tty" );
 373 #endif
 374
 375         for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
 376         if( fftscore && scoremtx != -1 )
 377         {
 378                 for( i=0; i<clus2; i++ )
 379                 {
 380                         seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
 381                         seq_vec_2( seqVector2[1], volume,   eff2[i], tmpseq2[i] );
 382                 }
 383         }
 384         else
 385         {
 386 #if 0
 387                 for( i=0; i<clus2; i++ ) for( j=0; j<n20or4or2; j++ )
 388                         seq_vec( seqVector2[j], amino[j], eff2[i], tmpseq2[i] );
 389 #else
 390                 for( i=0; i<clus2; i++ )
 391                         seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
 392 #endif
 393         }
 394 #if RND
 395         for( i=0; i<clus2; i++ )
 396         {
 397                 vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
 398         }
 399 #endif
 400
 401 #if 0
 402 fftfp = fopen( "seqVec2", "w" );
 403 fprintf( fftfp, "before fft\n" );
 404 for( k=0; k<n20or4or2; k++ )
 405 {
 406    fprintf( fftfp, "%c\n", amino[k] );
 407    for( l=0; l<nlen; l++ )
 408    fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
 409 }
 410 fclose( fftfp );
 411 system( "less seqVec2 < /dev/tty > /dev/tty" );
 412 #endif
 413
 414         for( j=0; j<n20or4or2; j++ )
 415         {
 416                 fft( nlen, seqVector2[j], 0 );
 417                 fft( nlen, seqVector1[j], 0 );
 418         }
 419 #if 0
 420 fftfp = fopen( "seqVec2", "w" );
 421 fprintf( fftfp, "#after fft\n" );
 422 for( k=0; k<n20or4or2; k++ )
 423 {
 424    fprintf( fftfp, "#%c\n", amino[k] );
 425    for( l=0; l<nlen; l++ )
 426    fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
 427 }
 428 fclose( fftfp );
 429 system( "less seqVec2 < /dev/tty > /dev/tty" );
 430 #endif
 431
 432         for( k=0; k<n20or4or2; k++ )
 433         {
 434                 for( l=0; l<nlen; l++ )
 435                         calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
 436         }
 437         for( l=0; l<nlen; l++ )
 438         {
 439                 naisekiNoWa[l].R = 0.0;
 440                 naisekiNoWa[l].I = 0.0;
 441                 for( k=0; k<n20or4or2; k++ )
 442                 {
 443                         naisekiNoWa[l].R += naiseki[k][l].R;
 444                         naisekiNoWa[l].I += naiseki[k][l].I;
 445                 }
 446         }
 447
 448 #if 0
 449 fftfp = fopen( "naisekiNoWa", "w" );
 450 fprintf( fftfp, "#Before fft\n" );
 451 for( l=0; l<nlen; l++ )
 452         fprintf( fftfp, "%d  %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I );
 453 fclose( fftfp );
 454 system( "less naisekiNoWa < /dev/tty > /dev/tty " );
 455 #endif
 456
 457         fft( -nlen, naisekiNoWa, 0 );
 458
 459         for( m=0; m<=nlen2; m++ )
 460                 soukan[m] = naisekiNoWa[nlen2-m].R;
 461         for( m=nlen2+1; m<nlen; m++ )
 462                 soukan[m] = naisekiNoWa[nlen+nlen2-m].R;
 463
 464 #if 0
 465 fftfp = fopen( "naisekiNoWa", "w" );
 466 fprintf( fftfp, "#After fft\n" );
 467 for( l=0; l<nlen; l++ )
 468         fprintf( fftfp, "%d  %f\n", l, naisekiNoWa[l].R );
 469 fclose( fftfp );
 470 fftfp = fopen( "list.plot", "w"  );
 471 fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
 472 fclose( fftfp );
 473 system( "/usr/bin/gnuplot list.plot &" );
 474 #endif
 475 #if 0
 476 fprintf( stderr, "frt write start\n" );
 477 fftfp = fopen( "frt", "w" );
 478 for( l=0; l<nlen; l++ )
 479         fprintf( fftfp, "%d  %f\n", l-nlen2, soukan[l] );
 480 fclose( fftfp );
 481 system( "less frt < /dev/tty > /dev/tty" );
 482 #if 0
 483 fftfp = fopen( "list.plot", "w"  );
 484 fprintf( fftfp, "plot 'frt'\n pause +1" );
 485 fclose( fftfp );
 486 system( "/usr/bin/gnuplot list.plot" );
 487 #endif
 488 #endif
 489
 490
 491         getKouho( kouho, NKOUHO, soukan, nlen );
 492
 493 #if 0
 494         for( i=0; i<NKOUHO; i++ )
 495         {
 496                 fprintf( stdout, "kouho[%d] = %d\n", i, kouho[i] );
 497         }
 498 #endif
 499
 500 #if KEIKA
 501         fprintf( stderr, "Searching anchors ... " );
 502 #endif
 503         count = 0;
 504
 505
 506
 507 #define CAND 0
 508 #if CAND
 509         fftfp = fopen( "cand", "w" );
 510         fclose( fftfp );
 511 #endif
 512
 513         for( k=0; k<NKOUHO; k++ )
 514         {
 515
 516                 lag = kouho[k];
 517                 zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
 518 #if CAND
 519                 fftfp = fopen( "cand", "a" );
 520                 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
 521                 fprintf( fftfp, "%s\n", tmpptr1[0] );
 522                 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
 523                 fprintf( fftfp, "%s\n", tmpptr2[0] );
 524                 fprintf( fftfp, ">\n", k+1, lag );
 525                 fclose( fftfp );
 526 #endif
 527                 tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
 528
 529                 if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
 530
 531
 532                 if( tmpint == 0 ) break; // 060430 iinoka ?
 533                 while( tmpint-- > 0 )
 534                 {
 535                         if( lag > 0 )
 536                         {
 537                                 segment1[count].start  = segment[count].start ;
 538                                 segment1[count].end    = segment[count].end   ;
 539                                 segment1[count].center = segment[count].center;
 540                                 segment1[count].score  = segment[count].score;
 541
 542                                 segment2[count].start  = segment[count].start  + lag;
 543                                 segment2[count].end    = segment[count].end    + lag;
 544                                 segment2[count].center = segment[count].center + lag;
 545                                 segment2[count].score  = segment[count].score       ;
 546                         }
 547                         else
 548                         {
 549                                 segment1[count].start  = segment[count].start  - lag;
 550                                 segment1[count].end    = segment[count].end    - lag;
 551                                 segment1[count].center = segment[count].center - lag;
 552                                 segment1[count].score  = segment[count].score       ;
 553
 554                                 segment2[count].start  = segment[count].start ;
 555                                 segment2[count].end    = segment[count].end   ;
 556                                 segment2[count].center = segment[count].center;
 557                                 segment2[count].score  = segment[count].score ;
 558                         }
 559 #if 0
 560                         fprintf( stderr, "Goukaku=%dko\n", tmpint );
 561                         fprintf( stderr, "in 1 %d\n", segment1[count].center );
 562                         fprintf( stderr, "in 2 %d\n", segment2[count].center );
 563 #endif
 564                         segment1[count].pair = &segment2[count];
 565                         segment2[count].pair = &segment1[count];
 566                         count++;
 567 #if 0
 568                         fprintf( stderr, "count=%d\n", count );
 569 #endif
 570                 }
 571         }
 572
 573 #if 1
 574         fprintf( stderr, "done. (%d anchors)\r", count );
 575 #endif
 576         if( !count && fftNoAnchStop )
 577                 ErrorExit( "Cannot detect anchor!" );
 578 #if 0
 579         fprintf( stdout, "RESULT before sort:\n" );
 580         for( l=0; l<count+1; l++ )
 581         {
 582                 fprintf( stdout, "cut[%d]=%d, ", l, segment1[l].center );
 583                 fprintf( stdout, "%d score = %f\n", segment2[l].center, segment1[l].score );
 584         }
 585         exit( 1 );
 586 #endif
 587
 588 #if KEIKA
 589         fprintf( stderr, "Aligning anchors ... " );
 590 #endif
 591         for( i=0; i<count; i++ )
 592         {
 593                 sortedseg1[i] = &segment1[i];
 594                 sortedseg2[i] = &segment2[i];
 595         }
 596
 597         {
 598                 mymergesort( 0, count-1, sortedseg1 );
 599                 mymergesort( 0, count-1, sortedseg2 );
 600                 for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
 601                 for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
 602
 603                 if( crossscoresize < count+2 )
 604                 {
 605                         crossscoresize = count+2;
 606                         fprintf( stderr, "####################################################################################################################################allocating crossscore, size = %d\n", crossscoresize );
 607                         if( crossscore ) FreeDoubleMtx( crossscore );
 608                         crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
 609                 }
 610
 611                 for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
 612                         crossscore[i][j] = 0.0;
 613                 for( i=0; i<count; i++ )
 614                 {
 615                         crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
 616                         cut1[i+1] = sortedseg1[i]->center;
 617                         cut2[i+1] = sortedseg2[i]->center;
 618                 }
 619
 620 #if DEBUG
 621                 fprintf( stderr, "AFTER SORT\n" );
 622                 for( i=0; i<count; i++ ) fprintf( stderr, "%d, %d\n", segment1[i].start, segment2[i].start );
 623 #endif
 624
 625                 crossscore[0][0] = 10000000.0;
 626                 cut1[0] = 0;
 627                 cut2[0] = 0;
 628                 crossscore[count+1][count+1] = 10000000.0;
 629                 cut1[count+1] = len1;
 630                 cut2[count+1] = len2;
 631                 count += 2;
 632                 count0 = count;
 633
 634                 blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
 635         }
 636         if( fftkeika )
 637         {
 638                 if( count0 > count )
 639                 {
 640                         fprintf( stderr, "REPEAT!? \n" );
 641                         if( fftRepeatStop ) exit( 1 );
 642                 }
 643 #if KEIKA
 644                 else
 645                         fprintf( stderr, "done\n" );
 646                         fprintf( stderr, "done. (%d anchors)\n", count );
 647 #endif
 648         }
 649
 650 #if 0
 651         fftfp = fopen( "fft", "a" );
 652         fprintf( fftfp, "RESULT after sort:\n" );
 653         for( l=0; l<count; l++ )
 654         {
 655                 fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
 656                 fprintf( fftfp, "%d\n", segment2[l].center );
 657         }
 658         fclose( fftfp );
 659 #endif
 660
 661 #if 0
 662         fftfp = fopen( "fft", "a" );
 663         fprintf( fftfp, "RESULT after sort:\n" );
 664         for( l=0; l<count; l++ )
 665         {
 666                 fprintf( fftfp, "cut : %d %d\n", cut1[l], cut2[l] );
 667         }
 668         fclose( fftfp );
 669 #endif
 670
 671 #if KEIKA
 672         fprintf( trap_g, "Devided to %d segments\n", count-1 );
 673         fprintf( trap_g, "%d  %d forg\n", MIN( clus1, clus2 ), count-1 );
 674 #endif
 675
 676         totallen = 0;
 677         for( j=0; j<clus1; j++ ) result1[j][0] = 0;
 678         for( j=0; j<clus2; j++ ) result2[j][0] = 0;
 679         for( i=0; i<count-1; i++ )
 680         {
 681                 if( i == 0 ) headgp = outgap; else headgp = 1;
 682                 if( i == count-2 ) tailgp = outgap; else tailgp = 1;
 683
 684 #if DEBUG
 685                 fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
 686 #else
 687 #if KEIKA
 688                 fprintf( stderr, "DP %03d / %03d\r", i+1, count-1 );
 689 #endif
 690 #endif
 691                 for( j=0; j<clus1; j++ )
 692                 {
 693                         strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
 694                         tmpres1[j][cut1[i+1]-cut1[i]] = 0;
 695                 }
 696                 for( j=0; j<clus2; j++ )
 697                 {
 698                         strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
 699                         tmpres2[j][cut2[i+1]-cut2[i]] = 0;
 700                 }
 701                 switch( alg )
 702                 {
 703                         case( 'a' ):
 704                                 Aalign( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen );
 705                                 break;
 706                         case( 'M' ):
 707                                         MSalignmm( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, headgp, tailgp );
 708                                 break;
 709                         case( 'A' ):
 710                                 if( clus1 == 1 && clus2 == 1 )
 711                                         G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
 712                                 else
 713                                         A__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, headgp, tailgp );
 714                                 break;
 715                         case( 'H' ):
 716                                 if( clus1 == 1 && clus2 == 1 )
 717                                         G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
 718                                 else
 719                                         H__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
 720                                 break;
 721                         case( 'Q' ):
 722                                 if( clus1 == 1 && clus2 == 1 )
 723                                         G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
 724                                 else
 725                                         Q__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
 726                                 break;
 727                         default:
 728                                 fprintf( stderr, "alg = %c\n", alg );
 729                                 ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
 730                                 break;
 731                 }
 732
 733                 nlen = strlen( tmpres1[0] );
 734                 if( totallen + nlen > alloclen ) ErrorExit( "LENGTH OVER in Falign\n " );
 735                 for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
 736                 for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
 737                 totallen += nlen;
 738 #if 0
 739                 fprintf( stderr, "%4d\r", totallen );
 740                 fprintf( stderr, "\n\n" );
 741                 for( j=0; j<clus1; j++ )
 742                 {
 743                         fprintf( stderr, "%s\n", tmpres1[j] );
 744                 }
 745                 fprintf( stderr, "-------\n" );
 746                 for( j=0; j<clus2; j++ )
 747                 {
 748                         fprintf( stderr, "%s\n", tmpres2[j] );
 749                 }
 750 #endif
 751         }
 752 #if KEIKA
 753         fprintf( stderr, "DP ... done   \n" );
 754 #endif
 755
 756         for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
 757         for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
 758 #if 0
 759         for( j=0; j<clus1; j++ )
 760         {
 761                 fprintf( stderr, "%s\n", result1[j] );
 762         }
 763         fprintf( stderr, "- - - - - - - - - - -\n" );
 764         for( j=0; j<clus2; j++ )
 765         {
 766                 fprintf( stderr, "%s\n", result2[j] );
 767         }
 768 #endif
 769         return( 0.0 );
 770 }
 771
 772
 773
 774 float Falign( char  **seq1, char  **seq2,
 775                           double *eff1, double *eff2,
 776                           int    clus1, int    clus2,
 777                           int alloclen, int *fftlog,
 778                           int *chudanpt, int chudanref, int *chudanres )
 779 {
 780         int i, j, k, l, m, maxk;
 781         int nlen, nlen2, nlen4;
 782         static TLS int prevalloclen = 0;
 783         static TLS int crossscoresize = 0;
 784         static TLS char **tmpseq1 = NULL;
 785         static TLS char **tmpseq2 = NULL;
 786         static TLS char **tmpptr1 = NULL;
 787         static TLS char **tmpptr2 = NULL;
 788         static TLS char **tmpres1 = NULL;
 789         static TLS char **tmpres2 = NULL;
 790         static TLS char **result1 = NULL;
 791         static TLS char **result2 = NULL;
 792 #if RND
 793         static TLS char **rndseq1 = NULL;
 794         static TLS char **rndseq2 = NULL;
 795 #endif
 796         static TLS Fukusosuu **seqVector1 = NULL;
 797         static TLS Fukusosuu **seqVector2 = NULL;
 798         static TLS Fukusosuu **naiseki = NULL;
 799         static TLS Fukusosuu *naisekiNoWa = NULL;
 800         static TLS double *soukan = NULL;
 801         static TLS double **crossscore = NULL;
 802         int nlentmp;
 803         static TLS int *kouho = NULL;
 804         static TLS Segment *segment = NULL;
 805         static TLS Segment *segment1 = NULL;
 806         static TLS Segment *segment2 = NULL;
 807         static TLS Segment **sortedseg1 = NULL;
 808         static TLS Segment **sortedseg2 = NULL;
 809         static TLS int *cut1 = NULL;
 810         static TLS int *cut2 = NULL;
 811         static TLS char *sgap1, *egap1, *sgap2, *egap2;
 812         static TLS int localalloclen = 0;
 813         int lag;
 814         int tmpint;
 815         int count, count0;
 816         int len1, len2;
 817         int totallen;
 818         float totalscore;
 819         float dumfl = 0.0;
 820         int headgp, tailgp;
 821
 822
 823         if( seq1 == NULL )
 824         {
 825                 if( result1 )
 826                 {
 827 //                      fprintf( stderr, "Freeing localarrays in Falign\n" );
 828                         localalloclen = 0;
 829                         mymergesort( 0, 0, NULL );
 830                         alignableReagion( 0, 0, NULL, NULL, NULL, NULL, NULL );
 831                         fft( 0, NULL, 1 );
 832                         A__align( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
 833                         G__align11( NULL, NULL, 0, 0, 0 );
 834                         blockAlign2( NULL, NULL, NULL, NULL, NULL, NULL );
 835                         if( crossscore ) FreeDoubleMtx( crossscore );
 836                         FreeCharMtx( result1 );
 837                         FreeCharMtx( result2 );
 838                         FreeCharMtx( tmpres1 );
 839                         FreeCharMtx( tmpres2 );
 840                         FreeCharMtx( tmpseq1 );
 841                         FreeCharMtx( tmpseq2 );
 842                         free( sgap1 );
 843                         free( egap1 );
 844                         free( sgap2 );
 845                         free( egap2 );
 846                         free( kouho );
 847                         free( cut1 );
 848                         free( cut2 );
 849                         free( tmpptr1 );
 850                         free( tmpptr2 );
 851                         free( segment );
 852                         free( segment1 );
 853                         free( segment2 );
 854                         free( sortedseg1 );
 855                         free( sortedseg2 );
 856                         if( !kobetsubunkatsu )
 857                         {
 858                                 FreeFukusosuuMtx ( seqVector1 );
 859                                 FreeFukusosuuMtx ( seqVector2 );
 860                                 FreeFukusosuuVec( naisekiNoWa );
 861                                 FreeFukusosuuMtx( naiseki );
 862                                 FreeDoubleVec( soukan );
 863                         }
 864                 }
 865                 else
 866                 {
 867 //                      fprintf( stderr, "Did not allocate localarrays in Falign\n" );
 868                 }
 869
 870                 return( 0.0 );
 871         }
 872
 873         len1 = strlen( seq1[0] );
 874         len2 = strlen( seq2[0] );
 875         nlentmp = MAX( len1, len2 );
 876
 877         nlen = 1;
 878         while( nlentmp >= nlen ) nlen <<= 1;
 879 #if 0
 880         fprintf( stderr, "###   nlen    = %d\n", nlen );
 881 #endif
 882
 883         nlen2 = nlen/2; nlen4 = nlen2 / 2;
 884
 885 #if DEBUG
 886         fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
 887         fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
 888 #endif
 889
 890         if( prevalloclen != alloclen ) // Falign_noudp mo kaeru
 891         {
 892                 if( prevalloclen )
 893                 {
 894                         FreeCharMtx( result1 );
 895                         FreeCharMtx( result2 );
 896                         FreeCharMtx( tmpres1 );
 897                         FreeCharMtx( tmpres2 );
 898                 }
 899 //              fprintf( stderr, "\n\n\nreallocating ...\n" );
 900                 result1 = AllocateCharMtx( njob, alloclen );
 901                 result2 = AllocateCharMtx( njob, alloclen );
 902                 tmpres1 = AllocateCharMtx( njob, alloclen );
 903                 tmpres2 = AllocateCharMtx( njob, alloclen );
 904                 prevalloclen = alloclen;
 905         }
 906         if( !localalloclen )
 907         {
 908                 sgap1 = AllocateCharVec( njob );
 909                 egap1 = AllocateCharVec( njob );
 910                 sgap2 = AllocateCharVec( njob );
 911                 egap2 = AllocateCharVec( njob );
 912                 kouho = AllocateIntVec( NKOUHO );
 913                 cut1 = AllocateIntVec( MAXSEG );
 914                 cut2 = AllocateIntVec( MAXSEG );
 915                 tmpptr1 = AllocateCharMtx( njob, 0 );
 916                 tmpptr2 = AllocateCharMtx( njob, 0 );
 917 //              crossscore = AllocateDoubleMtx( MAXSEG, MAXSEG );
 918                 segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
 919                 segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
 920                 segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
 921                 sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
 922                 sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
 923                 if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
 924                         ErrorExit( "Allocation error\n" );
 925
 926                 if     ( scoremtx == -1 ) n20or4or2 = 1;
 927                 else if( fftscore )       n20or4or2 = 1;
 928                 else                      n20or4or2 = 20;
 929         }
 930         if( localalloclen < nlen )
 931         {
 932                 if( localalloclen )
 933                 {
 934 #if 1
 935                         if( !kobetsubunkatsu )
 936                         {
 937                                 FreeFukusosuuMtx ( seqVector1 );
 938                                 FreeFukusosuuMtx ( seqVector2 );
 939                                 FreeFukusosuuVec( naisekiNoWa );
 940                                 FreeFukusosuuMtx( naiseki );
 941                                 FreeDoubleVec( soukan );
 942                         }
 943                         FreeCharMtx( tmpseq1 );
 944                         FreeCharMtx( tmpseq2 );
 945 #endif
 946 #if RND
 947                         FreeCharMtx( rndseq1 );
 948                         FreeCharMtx( rndseq2 );
 949 #endif
 950                 }
 951
 952                 tmpseq1 = AllocateCharMtx( njob, nlen );
 953                 tmpseq2 = AllocateCharMtx( njob, nlen );
 954                 if( !kobetsubunkatsu )
 955                 {
 956                         naisekiNoWa = AllocateFukusosuuVec( nlen );
 957                         naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
 958                         seqVector1 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
 959                         seqVector2 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
 960                         soukan = AllocateDoubleVec( nlen+1 );
 961                 }
 962 #if RND
 963                 rndseq1 = AllocateCharMtx( njob, nlen );
 964                 rndseq2 = AllocateCharMtx( njob, nlen );
 965                 for( i=0; i<njob; i++ )
 966                 {
 967                         generateRndSeq( rndseq1[i], nlen );
 968                         generateRndSeq( rndseq2[i], nlen );
 969                 }
 970 #endif
 971                 localalloclen = nlen;
 972         }
 973
 974         for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
 975         for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
 976
 977 #if 0
 978 fftfp = fopen( "input_of_Falign", "w" );
 979 fprintf( fftfp, "nlen = %d\n", nlen );
 980 fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
 981 for( i=0; i<clus1; i++ )
 982         fprintf( fftfp, "%s\n", seq1[i] );
 983 fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
 984 for( i=0; i<clus2; i++ )
 985         fprintf( fftfp, "%s\n", seq2[i] );
 986 fclose( fftfp );
 987 system( "less input_of_Falign < /dev/tty > /dev/tty" );
 988 #endif
 989         if( !kobetsubunkatsu )
 990         {
 991                 if( fftkeika ) fprintf( stderr,  " FFT ... " );
 992
 993                 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
 994                 if( fftscore && scoremtx != -1 )
 995                 {
 996                         for( i=0; i<clus1; i++ )
 997                         {
 998 #if 1
 999                                 seq_vec_5( seqVector1[0], polarity, volume, eff1[i], tmpseq1[i] );
1000 #else
1001                                 seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
1002                                 seq_vec_2( seqVector1[1], volume,   eff1[i], tmpseq1[i] );
1003 #endif
1004                         }
1005                 }
1006                 else
1007                 {
1008 #if 0
1009                         for( i=0; i<clus1; i++ ) for( j=0; j<n20or4or2; j++ )
1010                                 seq_vec( seqVector1[j], amino[j], eff1[i], tmpseq1[i] );
1011 #else
1012                         for( i=0; i<clus1; i++ )
1013                                 seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
1014 #endif
1015                 }
1016 #if RND
1017                 for( i=0; i<clus1; i++ )
1018                 {
1019                         vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
1020                 }
1021 #endif
1022 #if 0
1023 fftfp = fopen( "seqVec", "w" );
1024 fprintf( fftfp, "before transform\n" );
1025 for( k=0; k<n20or4or2; k++ )
1026 {
1027    fprintf( fftfp, "nlen=%d\n", nlen );
1028    fprintf( fftfp, "%c\n", amino[k] );
1029    for( l=0; l<nlen; l++ )
1030    fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
1031 }
1032 fclose( fftfp );
1033 system( "less seqVec < /dev/tty > /dev/tty" );
1034 #endif
1035
1036                 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
1037                 if( fftscore && scoremtx != -1 )
1038                 {
1039                         for( i=0; i<clus2; i++ )
1040                         {
1041 #if 1
1042                                 seq_vec_5( seqVector2[0], polarity, volume, eff2[i], tmpseq2[i] );
1043 #else
1044                                 seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
1045                                 seq_vec_2( seqVector2[1], volume,   eff2[i], tmpseq2[i] );
1046 #endif
1047                         }
1048                 }
1049                 else
1050                 {
1051 #if 0
1052                         for( i=0; i<clus2; i++ ) for( j=0; j<n20or4or2; j++ )
1053                                 seq_vec( seqVector2[j], amino[j], eff2[i], tmpseq2[i] );
1054 #else
1055                         for( i=0; i<clus2; i++ )
1056                                 seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
1057 #endif
1058                 }
1059 #if RND
1060                 for( i=0; i<clus2; i++ )
1061                 {
1062                         vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
1063                 }
1064 #endif
1065
1066 #if 0
1067 fftfp = fopen( "seqVec2", "w" );
1068 fprintf( fftfp, "before fft\n" );
1069 for( k=0; k<n20or4or2; k++ )
1070 {
1071    fprintf( fftfp, "%c\n", amino[k] );
1072    for( l=0; l<nlen; l++ )
1073    fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1074 }
1075 fclose( fftfp );
1076 system( "less seqVec2 < /dev/tty > /dev/tty" );
1077 #endif
1078
1079                 for( j=0; j<n20or4or2; j++ )
1080                 {
1081                         fft( nlen, seqVector2[j], 0 );
1082                         fft( nlen, seqVector1[j], 0 );
1083                 }
1084 #if 0
1085 fftfp = fopen( "seqVec2", "w" );
1086 fprintf( fftfp, "#after fft\n" );
1087 for( k=0; k<n20or4or2; k++ )
1088 {
1089    fprintf( fftfp, "#%c\n", amino[k] );
1090    for( l=0; l<nlen; l++ )
1091            fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1092 }
1093 fclose( fftfp );
1094 system( "less seqVec2 < /dev/tty > /dev/tty" );
1095 #endif
1096
1097                 for( k=0; k<n20or4or2; k++ )
1098                 {
1099                         for( l=0; l<nlen; l++ )
1100                                 calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
1101                 }
1102                 for( l=0; l<nlen; l++ )
1103                 {
1104                         naisekiNoWa[l].R = 0.0;
1105                         naisekiNoWa[l].I = 0.0;
1106                         for( k=0; k<n20or4or2; k++ )
1107                         {
1108                                 naisekiNoWa[l].R += naiseki[k][l].R;
1109                                 naisekiNoWa[l].I += naiseki[k][l].I;
1110                         }
1111                 }
1112
1113 #if 0
1114         fftfp = fopen( "naisekiNoWa", "w" );
1115         fprintf( fftfp, "#Before fft\n" );
1116         for( l=0; l<nlen; l++ )
1117                 fprintf( fftfp, "%d  %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I );
1118         fclose( fftfp );
1119         system( "less naisekiNoWa < /dev/tty > /dev/tty " );
1120 #endif
1121
1122                 fft( -nlen, naisekiNoWa, 0 );
1123
1124                 for( m=0; m<=nlen2; m++ )
1125                         soukan[m] = naisekiNoWa[nlen2-m].R;
1126                 for( m=nlen2+1; m<nlen; m++ )
1127                         soukan[m] = naisekiNoWa[nlen+nlen2-m].R;
1128
1129 #if 0
1130         fftfp = fopen( "naisekiNoWa", "w" );
1131         fprintf( fftfp, "#After fft\n" );
1132         for( l=0; l<nlen; l++ )
1133                 fprintf( fftfp, "%d  %f\n", l, naisekiNoWa[l].R );
1134         fclose( fftfp );
1135         fftfp = fopen( "list.plot", "w"  );
1136         fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
1137         fclose( fftfp );
1138         system( "/usr/bin/gnuplot list.plot &" );
1139 #endif
1140 #if 0
1141         fprintf( stderr, "soukan\n" );
1142         for( l=0; l<nlen; l++ )
1143                 fprintf( stderr, "%d  %f\n", l-nlen2, soukan[l] );
1144 #if 0
1145         fftfp = fopen( "list.plot", "w"  );
1146         fprintf( fftfp, "plot 'frt'\n pause +1" );
1147         fclose( fftfp );
1148         system( "/usr/bin/gnuplot list.plot" );
1149 #endif
1150 #endif
1151
1152
1153                 getKouho( kouho, NKOUHO, soukan, nlen );
1154
1155 #if 0
1156                 for( i=0; i<NKOUHO; i++ )
1157                 {
1158                         fprintf( stderr, "kouho[%d] = %d\n", i, kouho[i] );
1159                 }
1160 #endif
1161         }
1162
1163 #if KEIKA
1164         fprintf( stderr, "Searching anchors ... " );
1165 #endif
1166         count = 0;
1167
1168
1169
1170 #define CAND 0
1171 #if CAND
1172         fftfp = fopen( "cand", "w" );
1173         fclose( fftfp );
1174 #endif
1175         if( kobetsubunkatsu )
1176         {
1177                 maxk = 1;
1178                 kouho[0] = 0;
1179         }
1180         else
1181         {
1182                 maxk = NKOUHO;
1183         }
1184
1185         for( k=0; k<maxk; k++ )
1186         {
1187                 lag = kouho[k];
1188                 if( lag <= -len1 || len2 <= lag ) continue;
1189                 zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
1190 #if CAND
1191                 fftfp = fopen( "cand", "a" );
1192                 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
1193                 fprintf( fftfp, "%s\n", tmpptr1[0] );
1194                 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
1195                 fprintf( fftfp, "%s\n", tmpptr2[0] );
1196                 fprintf( fftfp, ">\n", k+1, lag );
1197                 fclose( fftfp );
1198 #endif
1199
1200 //              fprintf( stderr, "lag = %d\n", lag );
1201                 tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
1202
1203 //              if( lag == -50 ) exit( 1 );
1204
1205                 if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
1206
1207
1208                 if( tmpint == 0 ) break; // 060430 iinoka ?
1209                 while( tmpint-- > 0 )
1210                 {
1211 #if 0
1212                         if( segment[count].end - segment[count].start < fftWinSize )
1213                         {
1214                                 count++;
1215                                 continue;
1216                         }
1217 #endif
1218                         if( lag > 0 )
1219                         {
1220                                 segment1[count].start  = segment[count].start ;
1221                                 segment1[count].end    = segment[count].end   ;
1222                                 segment1[count].center = segment[count].center;
1223                                 segment1[count].score  = segment[count].score;
1224
1225                                 segment2[count].start  = segment[count].start  + lag;
1226                                 segment2[count].end    = segment[count].end    + lag;
1227                                 segment2[count].center = segment[count].center + lag;
1228                                 segment2[count].score  = segment[count].score       ;
1229                         }
1230                         else
1231                         {
1232                                 segment1[count].start  = segment[count].start  - lag;
1233                                 segment1[count].end    = segment[count].end    - lag;
1234                                 segment1[count].center = segment[count].center - lag;
1235                                 segment1[count].score  = segment[count].score       ;
1236
1237                                 segment2[count].start  = segment[count].start ;
1238                                 segment2[count].end    = segment[count].end   ;
1239                                 segment2[count].center = segment[count].center;
1240                                 segment2[count].score  = segment[count].score ;
1241                         }
1242 #if 0
1243                         fprintf( stderr, "in 1 %d\n", segment1[count].center );
1244                         fprintf( stderr, "in 2 %d\n", segment2[count].center );
1245 #endif
1246                         segment1[count].pair = &segment2[count];
1247                         segment2[count].pair = &segment1[count];
1248                         count++;
1249                 }
1250         }
1251 #if 0
1252         if( !kobetsubunkatsu && fftkeika )
1253                 fprintf( stderr, "%d anchors found\r", count );
1254 #endif
1255         if( !count && fftNoAnchStop )
1256                 ErrorExit( "Cannot detect anchor!" );
1257 #if 0
1258         fprintf( stderr, "RESULT before sort:\n" );
1259         for( l=0; l<count+1; l++ )
1260         {
1261                 fprintf( stderr, "cut[%d]=%d, ", l, segment1[l].center );
1262                 fprintf( stderr, "%d score = %f\n", segment2[l].center, segment1[l].score );
1263         }
1264 #endif
1265
1266 #if KEIKA
1267         fprintf( stderr, "done. (%d anchors)\n", count );
1268         fprintf( stderr, "Aligning anchors ... " );
1269 #endif
1270         for( i=0; i<count; i++ )
1271         {
1272                 sortedseg1[i] = &segment1[i];
1273                 sortedseg2[i] = &segment2[i];
1274         }
1275 #if 0
1276         tmpsort( count, sortedseg1 );
1277         tmpsort( count, sortedseg2 );
1278         qsort( sortedseg1, count, sizeof( Segment * ), segcmp );
1279         qsort( sortedseg2, count, sizeof( Segment * ), segcmp );
1280 #else
1281         mymergesort( 0, count-1, sortedseg1 );
1282         mymergesort( 0, count-1, sortedseg2 );
1283 #endif
1284         for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
1285         for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
1286
1287
1288         if( kobetsubunkatsu )
1289         {
1290                 for( i=0; i<count; i++ )
1291             {
1292                         cut1[i+1] = sortedseg1[i]->center;
1293                         cut2[i+1] = sortedseg2[i]->center;
1294                 }
1295                 cut1[0] = 0;
1296                 cut2[0] = 0;
1297                 cut1[count+1] = len1;
1298                 cut2[count+1] = len2;
1299                 count += 2;
1300         }
1301         else
1302         {
1303                 if( crossscoresize < count+2 )
1304                 {
1305                         crossscoresize = count+2;
1306 #if 1
1307                         if( fftkeika ) fprintf( stderr, "######allocating crossscore, size = %d\n", crossscoresize );
1308 #endif
1309                         if( crossscore ) FreeDoubleMtx( crossscore );
1310                         crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
1311                 }
1312                 for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
1313                         crossscore[i][j] = 0.0;
1314                 for( i=0; i<count; i++ )
1315                 {
1316                         crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
1317                         cut1[i+1] = sortedseg1[i]->center;
1318                         cut2[i+1] = sortedseg2[i]->center;
1319                 }
1320
1321 #if 0
1322                 fprintf( stderr, "AFTER SORT\n" );
1323                 for( i=0; i<count+1; i++ ) fprintf( stderr, "%d, %d\n", cut1[i], cut2[i] );
1324                 fprintf( stderr, "crossscore = \n" );
1325                 for( i=0; i<count+1; i++ )
1326                 {
1327                         for( j=0; j<count+1; j++ )
1328                                 fprintf( stderr, "%.0f ", crossscore[i][j] );
1329                         fprintf( stderr, "\n" );
1330                 }
1331 #endif
1332
1333                 crossscore[0][0] = 10000000.0;
1334                 cut1[0] = 0;
1335                 cut2[0] = 0;
1336                 crossscore[count+1][count+1] = 10000000.0;
1337                 cut1[count+1] = len1;
1338                 cut2[count+1] = len2;
1339                 count += 2;
1340                 count0 = count;
1341
1342                 blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
1343
1344 //              if( count-count0 )
1345 //                      fprintf( stderr, "%d unused anchors\n", count0-count );
1346
1347                 if( !kobetsubunkatsu && fftkeika )
1348                         fprintf( stderr, "%d anchors found\n", count );
1349                 if( fftkeika )
1350                 {
1351                         if( count0 > count )
1352                         {
1353 #if 0
1354                                 fprintf( stderr, "\7 REPEAT!? \n" );
1355 #else
1356                                 fprintf( stderr, "REPEAT!? \n" );
1357 #endif
1358                                 if( fftRepeatStop ) exit( 1 );
1359                         }
1360 #if KEIKA
1361                         else fprintf( stderr, "done\n" );
1362 #endif
1363                 }
1364         }
1365
1366 #if 0
1367         fftfp = fopen( "fft", "a" );
1368         fprintf( fftfp, "RESULT after sort:\n" );
1369         for( l=0; l<count; l++ )
1370         {
1371                 fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
1372                 fprintf( fftfp, "%d\n", segment2[l].center );
1373         }
1374         fclose( fftfp );
1375 #endif
1376
1377 #if 0
1378         fprintf( stderr, "RESULT after blckalign:\n" );
1379         for( l=0; l<count+1; l++ )
1380         {
1381                 fprintf( stderr, "cut : %d %d\n", cut1[l], cut2[l] );
1382         }
1383 #endif
1384
1385 #if 0
1386         fprintf( trap_g, "Devided to %d segments\n", count-1 );
1387         fprintf( trap_g, "%d  %d forg\n", MIN( clus1, clus2 ), count-1 );
1388 #endif
1389
1390         totallen = 0;
1391         for( j=0; j<clus1; j++ ) result1[j][0] = 0;
1392         for( j=0; j<clus2; j++ ) result2[j][0] = 0;
1393         totalscore = 0.0;
1394         *fftlog = -1;
1395         for( i=0; i<count-1; i++ )
1396         {
1397                 *fftlog += 1;
1398                 if( i == 0 ) headgp = outgap; else headgp = 1;
1399                 if( i == count-2 ) tailgp = outgap; else tailgp = 1;
1400
1401
1402                 if( cut1[i] ) // chuui
1403                 {
1404 //                      getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
1405 //                      getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
1406                         getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
1407                         getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
1408                 }
1409                 else
1410                 {
1411                         for( j=0; j<clus1; j++ ) sgap1[j] = 'o';
1412                         for( j=0; j<clus2; j++ ) sgap2[j] = 'o';
1413                 }
1414                 if( cut1[i+1] != len1 ) // chuui
1415                 {
1416                         getkyokaigap( egap1, seq1, cut1[i+1], clus1 );
1417                         getkyokaigap( egap2, seq2, cut2[i+1], clus2 );
1418                 }
1419                 else
1420                 {
1421                         for( j=0; j<clus1; j++ ) egap1[j] = 'o';
1422                         for( j=0; j<clus2; j++ ) egap2[j] = 'o';
1423                 }
1424 #if 0
1425                 {
1426                         fprintf( stderr, "kyokkaigap1(%d)=", cut1[i]-1 );
1427                         for( j=0; j<clus1; j++ )
1428                                 fprintf( stderr, "%c", sgap1[j] );
1429                         fprintf( stderr, "=kyokkaigap1-start\n" );
1430                 }
1431                 {
1432                         fprintf( stderr, "kyokkaigap2(%d)=", cut2[i]-1 );
1433                         for( j=0; j<clus2; j++ )
1434                                 fprintf( stderr, "%c", sgap2[j] );
1435                         fprintf( stderr, "=kyokkaigap2-start\n" );
1436                 }
1437                 {
1438                         fprintf( stderr, "kyokkaigap1(%d)=", cut1[i]-1 );
1439                         for( j=0; j<clus1; j++ )
1440                                 fprintf( stderr, "%c", egap1[j] );
1441                         fprintf( stderr, "=kyokkaigap1-end\n" );
1442                 }
1443                 {
1444                         fprintf( stderr, "kyokkaigap2(%d)=", cut2[i]-1 );
1445                         for( j=0; j<clus2; j++ )
1446                                 fprintf( stderr, "%c", egap2[j] );
1447                         fprintf( stderr, "=kyokkaigap2-end\n" );
1448                 }
1449 #endif
1450
1451 #if DEBUG
1452                 fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
1453 #else
1454 #if KEIKA
1455                 fprintf( stderr, "DP %03d / %03d\r", i+1, count-1 );
1456 #endif
1457 #endif
1458                 for( j=0; j<clus1; j++ )
1459                 {
1460                         strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
1461                         tmpres1[j][cut1[i+1]-cut1[i]] = 0;
1462                 }
1463                 if( kobetsubunkatsu && fftkeika ) commongappick( clus1, tmpres1 ); //dvtditr \e$B$K8F$P$l$?$H$-\e(B fftkeika=1
1464 //              if( kobetsubunkatsu ) commongappick( clus1, tmpres1 );
1465                 for( j=0; j<clus2; j++ )
1466                 {
1467                         strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
1468                         tmpres2[j][cut2[i+1]-cut2[i]] = 0;
1469                 }
1470                 if( kobetsubunkatsu && fftkeika ) commongappick( clus2, tmpres2 ); //dvtditr \e$B$K8F$P$l$?$H$-\e(B fftkeika=1
1471 //              if( kobetsubunkatsu ) commongappick( clus2, tmpres2 );
1472
1473                 if( constraint )
1474                 {
1475                         fprintf( stderr, "Not supported\n" );
1476                         exit( 1 );
1477                 }
1478 #if 0
1479                 fprintf( stderr, "i=%d, before alignment", i );
1480                 fprintf( stderr, "%4d\n", totallen );
1481                 fprintf( stderr, "\n\n" );
1482                 for( j=0; j<clus1; j++ )
1483                 {
1484                         fprintf( stderr, "%s\n", tmpres1[j] );
1485                 }
1486                 fprintf( stderr, "-------\n" );
1487                 for( j=0; j<clus2; j++ )
1488                 {
1489                         fprintf( stderr, "%s\n", tmpres2[j] );
1490                 }
1491 #endif
1492
1493 #if 0
1494                 fprintf( stdout, "writing input\n" );
1495                 for( j=0; j<clus1; j++ )
1496                 {
1497                         fprintf( stdout, ">%d of GROUP1\n", j );
1498                         fprintf( stdout, "%s\n", tmpres1[j] );
1499                 }
1500                 for( j=0; j<clus2; j++ )
1501                 {
1502                         fprintf( stdout, ">%d of GROUP2\n", j );
1503                         fprintf( stdout, "%s\n", tmpres2[j] );
1504                 }
1505                 fflush( stdout );
1506 #endif
1507                 switch( alg )
1508                 {
1509                         case( 'a' ):
1510                                 totalscore += Aalign( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen );
1511                                 break;
1512                         case( 'M' ):
1513                                         totalscore += MSalignmm( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres, headgp, tailgp );
1514                                 break;
1515                         case( 'A' ):
1516                                 if( clus1 == 1 && clus2 == 1 )
1517                                 {
1518                                         totalscore += G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
1519                                 }
1520                                 else
1521                                         totalscore += A__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres, headgp, tailgp );
1522                                 break;
1523                         case( 'H' ):
1524                                 if( clus1 == 1 && clus2 == 1 )
1525                                 {
1526                                         totalscore += G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
1527                                 }
1528                                 else
1529                                         totalscore += H__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, sgap1, sgap2, egap1, egap2 );
1530                                 break;
1531                         case( 'Q' ):
1532                                 if( clus1 == 1 && clus2 == 1 )
1533                                 {
1534                                         totalscore += G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
1535                                 }
1536                                 else
1537                                         totalscore += Q__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, sgap1, sgap2, egap1, egap2 );
1538                                 break;
1539                         default:
1540                                 fprintf( stderr, "alg = %c\n", alg );
1541                                 ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
1542                                 break;
1543                 }
1544
1545 #ifdef enablemultithread
1546                 if( chudanres && *chudanres )
1547                 {
1548 //                      fprintf( stderr, "\n\n## CHUUDAN!!! at Falign_localhom\n" );
1549                         return( -1.0 );
1550                 }
1551 #endif
1552
1553                 nlen = strlen( tmpres1[0] );
1554                 if( totallen + nlen > alloclen )
1555                 {
1556                         fprintf( stderr, "totallen=%d +  nlen=%d > alloclen = %d\n", totallen, nlen, alloclen );
1557                         ErrorExit( "LENGTH OVER in Falign\n " );
1558                 }
1559                 for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
1560                 for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
1561                 totallen += nlen;
1562 #if 0
1563                 fprintf( stderr, "i=%d", i );
1564                 fprintf( stderr, "%4d\n", totallen );
1565                 fprintf( stderr, "\n\n" );
1566                 for( j=0; j<clus1; j++ )
1567                 {
1568                         fprintf( stderr, "%s\n", tmpres1[j] );
1569                 }
1570                 fprintf( stderr, "-------\n" );
1571                 for( j=0; j<clus2; j++ )
1572                 {
1573                         fprintf( stderr, "%s\n", tmpres2[j] );
1574                 }
1575 #endif
1576         }
1577 #if KEIKA
1578         fprintf( stderr, "DP ... done   \n" );
1579 #endif
1580
1581         for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
1582         for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
1583 #if 0
1584         for( j=0; j<clus1; j++ )
1585         {
1586                 fprintf( stderr, "%s\n", result1[j] );
1587         }
1588         fprintf( stderr, "- - - - - - - - - - -\n" );
1589         for( j=0; j<clus2; j++ )
1590         {
1591                 fprintf( stderr, "%s\n", result2[j] );
1592         }
1593 #endif
1594         return( totalscore );
1595 }
1596
1597
1598
1599
1600
1601
1602
1603
1604 /*
1605 sakujo wo kentou (2010/10/05)
1606 */
1607 float Falign_udpari_long( char  **seq1, char  **seq2,
1608                           double *eff1, double *eff2,
1609                           int    clus1, int    clus2,
1610                           int alloclen, int *fftlog )
1611 {
1612         int i, j, k, l, m, maxk;
1613         int nlen, nlen2, nlen4;
1614         static TLS int prevalloclen = 0;
1615         static TLS int crossscoresize = 0;
1616         static TLS char **tmpseq1 = NULL;
1617         static TLS char **tmpseq2 = NULL;
1618         static TLS char **tmpptr1 = NULL;
1619         static TLS char **tmpptr2 = NULL;
1620         static TLS char **tmpres1 = NULL;
1621         static TLS char **tmpres2 = NULL;
1622         static TLS char **result1 = NULL;
1623         static TLS char **result2 = NULL;
1624 #if RND
1625         static TLS char **rndseq1 = NULL;
1626         static TLS char **rndseq2 = NULL;
1627 #endif
1628         static TLS Fukusosuu **seqVector1 = NULL;
1629         static TLS Fukusosuu **seqVector2 = NULL;
1630         static TLS Fukusosuu **naiseki = NULL;
1631         static TLS Fukusosuu *naisekiNoWa = NULL;
1632         static TLS double *soukan = NULL;
1633         static TLS double **crossscore = NULL;
1634         int nlentmp;
1635         static TLS int *kouho = NULL;
1636         static TLS Segment *segment = NULL;
1637         static TLS Segment *segment1 = NULL;
1638         static TLS Segment *segment2 = NULL;
1639         static TLS Segment **sortedseg1 = NULL;
1640         static TLS Segment **sortedseg2 = NULL;
1641         static TLS int *cut1 = NULL;
1642         static TLS int *cut2 = NULL;
1643         static TLS char *sgap1, *egap1, *sgap2, *egap2;
1644         static TLS int localalloclen = 0;
1645         int lag;
1646         int tmpint;
1647         int count, count0;
1648         int len1, len2;
1649         int totallen;
1650         float totalscore;
1651         int nkouho = 0;
1652         int headgp, tailgp;
1653 //      float dumfl = 0.0;
1654
1655
1656
1657         len1 = strlen( seq1[0] );
1658         len2 = strlen( seq2[0] );
1659         nlentmp = MAX( len1, len2 );
1660
1661         nlen = 1;
1662         while( nlentmp >= nlen ) nlen <<= 1;
1663 #if 0
1664         fprintf( stderr, "###   nlen    = %d\n", nlen );
1665 #endif
1666
1667         nlen2 = nlen/2; nlen4 = nlen2 / 2;
1668
1669 #if 0
1670         fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
1671         fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
1672 #endif
1673
1674         if( prevalloclen != alloclen ) // Falign_noudp mo kaeru
1675         {
1676                 if( prevalloclen )
1677                 {
1678                         FreeCharMtx( result1 );
1679                         FreeCharMtx( result2 );
1680                         FreeCharMtx( tmpres1 );
1681                         FreeCharMtx( tmpres2 );
1682                 }
1683 //              fprintf( stderr, "\n\n\nreallocating ...\n" );
1684                 result1 = AllocateCharMtx( njob, alloclen );
1685                 result2 = AllocateCharMtx( njob, alloclen );
1686                 tmpres1 = AllocateCharMtx( njob, alloclen );
1687                 tmpres2 = AllocateCharMtx( njob, alloclen );
1688                 prevalloclen = alloclen;
1689         }
1690
1691         if( !localalloclen )
1692         {
1693                 sgap1 = AllocateCharVec( njob );
1694                 egap1 = AllocateCharVec( njob );
1695                 sgap2 = AllocateCharVec( njob );
1696                 egap2 = AllocateCharVec( njob );
1697                 kouho = AllocateIntVec( NKOUHO_LONG );
1698                 cut1 = AllocateIntVec( MAXSEG );
1699                 cut2 = AllocateIntVec( MAXSEG );
1700                 tmpptr1 = AllocateCharMtx( njob, 0 );
1701                 tmpptr2 = AllocateCharMtx( njob, 0 );
1702                 segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
1703                 segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
1704                 segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
1705                 sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
1706                 sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
1707                 if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
1708                         ErrorExit( "Allocation error\n" );
1709
1710                 if     ( scoremtx == -1 ) n20or4or2 = 1;
1711                 else if( fftscore )       n20or4or2 = 1;
1712                 else                      n20or4or2 = 20;
1713         }
1714         if( localalloclen < nlen )
1715         {
1716                 if( localalloclen )
1717                 {
1718 #if 1
1719                         if( !kobetsubunkatsu )
1720                         {
1721                                 FreeFukusosuuMtx ( seqVector1 );
1722                                 FreeFukusosuuMtx ( seqVector2 );
1723                                 FreeFukusosuuVec( naisekiNoWa );
1724                                 FreeFukusosuuMtx( naiseki );
1725                                 FreeDoubleVec( soukan );
1726                         }
1727                         FreeCharMtx( tmpseq1 );
1728                         FreeCharMtx( tmpseq2 );
1729 #endif
1730 #if RND
1731                         FreeCharMtx( rndseq1 );
1732                         FreeCharMtx( rndseq2 );
1733 #endif
1734                 }
1735
1736
1737                 tmpseq1 = AllocateCharMtx( njob, nlen );
1738                 tmpseq2 = AllocateCharMtx( njob, nlen );
1739                 if( !kobetsubunkatsu )
1740                 {
1741                         naisekiNoWa = AllocateFukusosuuVec( nlen );
1742                         naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
1743                         seqVector1 = AllocateFukusosuuMtx( n20or4or2, nlen+1 );
1744                         seqVector2 = AllocateFukusosuuMtx( n20or4or2, nlen+1 );
1745                         soukan = AllocateDoubleVec( nlen+1 );
1746                 }
1747 #if RND
1748                 rndseq1 = AllocateCharMtx( njob, nlen );
1749                 rndseq2 = AllocateCharMtx( njob, nlen );
1750                 for( i=0; i<njob; i++ )
1751                 {
1752                         generateRndSeq( rndseq1[i], nlen );
1753                         generateRndSeq( rndseq2[i], nlen );
1754                 }
1755 #endif
1756                 localalloclen = nlen;
1757         }
1758
1759         for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
1760         for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
1761
1762 #if 0
1763 fftfp = fopen( "input_of_Falign", "w" );
1764 fprintf( fftfp, "nlen = %d\n", nlen );
1765 fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
1766 for( i=0; i<clus1; i++ )
1767         fprintf( fftfp, "%s\n", seq1[i] );
1768 fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
1769 for( i=0; i<clus2; i++ )
1770         fprintf( fftfp, "%s\n", seq2[i] );
1771 fclose( fftfp );
1772 system( "less input_of_Falign < /dev/tty > /dev/tty" );
1773 #endif
1774         if( !kobetsubunkatsu )
1775         {
1776                 fprintf( stderr,  " FFT ... " );
1777
1778                 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
1779                 if( scoremtx == -1 )
1780                 {
1781                         for( i=0; i<clus1; i++ )
1782                                 seq_vec_4( seqVector1[0], eff1[i], tmpseq1[i] );
1783                 }
1784                 else if( fftscore )
1785                 {
1786                         for( i=0; i<clus1; i++ )
1787                         {
1788 #if 0
1789                                 seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
1790                                 seq_vec_2( seqVector1[1], volume,   eff1[i], tmpseq1[i] );
1791 #else
1792                                 seq_vec_5( seqVector1[0], polarity, volume, eff1[i], tmpseq1[i] );
1793 #endif
1794                         }
1795                 }
1796                 else
1797                 {
1798                         for( i=0; i<clus1; i++ )
1799                                 seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
1800                 }
1801 #if RND
1802                 for( i=0; i<clus1; i++ )
1803                 {
1804                         vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
1805                 }
1806 #endif
1807 #if 0
1808 fftfp = fopen( "seqVec", "w" );
1809 fprintf( fftfp, "before transform\n" );
1810 for( k=0; k<n20or4or2; k++ )
1811 {
1812    fprintf( fftfp, "nlen=%d\n", nlen );
1813    fprintf( fftfp, "%c\n", amino[k] );
1814    for( l=0; l<nlen; l++ )
1815    fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
1816 }
1817 fclose( fftfp );
1818 system( "less seqVec < /dev/tty > /dev/tty" );
1819 #endif
1820
1821                 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
1822                 if( scoremtx == -1 )
1823                 {
1824                         for( i=0; i<clus2; i++ )
1825                                 seq_vec_4( seqVector2[0], eff2[i], tmpseq2[i] );
1826                 }
1827                 else if( fftscore )
1828                 {
1829                         for( i=0; i<clus2; i++ )
1830                         {
1831 #if 0
1832                                 seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
1833                                 seq_vec_2( seqVector2[1], volume,   eff2[i], tmpseq2[i] );
1834 #else
1835                                 seq_vec_5( seqVector2[0], polarity, volume, eff2[i], tmpseq2[i] );
1836 #endif
1837                         }
1838                 }
1839                 else
1840                 {
1841                         for( i=0; i<clus2; i++ )
1842                                 seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
1843                 }
1844 #if RND
1845                 for( i=0; i<clus2; i++ )
1846                 {
1847                         vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
1848                 }
1849 #endif
1850
1851 #if 0
1852 fftfp = fopen( "seqVec2", "w" );
1853 fprintf( fftfp, "before fft\n" );
1854 for( k=0; k<n20or4or2; k++ )
1855 {
1856    fprintf( fftfp, "%c\n", amino[k] );
1857    for( l=0; l<nlen; l++ )
1858    fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1859 }
1860 fclose( fftfp );
1861 system( "less seqVec2 < /dev/tty > /dev/tty" );
1862 #endif
1863
1864                 for( j=0; j<n20or4or2; j++ )
1865                 {
1866                         fft( nlen, seqVector2[j], 0 );
1867                         fft( nlen, seqVector1[j], 0 );
1868                 }
1869 #if 0
1870 fftfp = fopen( "seqVec2", "w" );
1871 fprintf( fftfp, "#after fft\n" );
1872 for( k=0; k<n20or4or2; k++ )
1873 {
1874    fprintf( fftfp, "#%c\n", amino[k] );
1875    for( l=0; l<nlen; l++ )
1876            fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1877 }
1878 fclose( fftfp );
1879 system( "less seqVec2 < /dev/tty > /dev/tty" );
1880 #endif
1881
1882                 for( k=0; k<n20or4or2; k++ )
1883                 {
1884                         for( l=0; l<nlen; l++ )
1885                                 calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
1886                 }
1887                 for( l=0; l<nlen; l++ )
1888                 {
1889                         naisekiNoWa[l].R = 0.0;
1890                         naisekiNoWa[l].I = 0.0;
1891                         for( k=0; k<n20or4or2; k++ )
1892                         {
1893                                 naisekiNoWa[l].R += naiseki[k][l].R;
1894                                 naisekiNoWa[l].I += naiseki[k][l].I;
1895                         }
1896                 }
1897
1898 #if 0
1899         fftfp = fopen( "naisekiNoWa", "w" );
1900         fprintf( fftfp, "#Before fft\n" );
1901         for( l=0; l<nlen; l++ )
1902                 fprintf( fftfp, "%d  %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I );
1903         fclose( fftfp );
1904         system( "less naisekiNoWa < /dev/tty > /dev/tty " );
1905 #endif
1906
1907                 fft( -nlen, naisekiNoWa, 0 );
1908
1909                 for( m=0; m<=nlen2; m++ )
1910                         soukan[m] = naisekiNoWa[nlen2-m].R;
1911                 for( m=nlen2+1; m<nlen; m++ )
1912                         soukan[m] = naisekiNoWa[nlen+nlen2-m].R;
1913
1914 #if 0
1915         fftfp = fopen( "naisekiNoWa", "w" );
1916         fprintf( fftfp, "#After fft\n" );
1917         for( l=0; l<nlen; l++ )
1918                 fprintf( fftfp, "%d  %f\n", l, naisekiNoWa[l].R );
1919         fclose( fftfp );
1920         fftfp = fopen( "list.plot", "w"  );
1921         fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
1922         fclose( fftfp );
1923         system( "/usr/bin/gnuplot list.plot &" );
1924 #endif
1925 #if 0
1926         fprintf( stderr, "soukan\n" );
1927         for( l=0; l<nlen; l++ )
1928                 fprintf( stderr, "%d  %f\n", l-nlen2, soukan[l] );
1929 #if 0
1930         fftfp = fopen( "list.plot", "w"  );
1931         fprintf( fftfp, "plot 'frt'\n pause +1" );
1932         fclose( fftfp );
1933         system( "/usr/bin/gnuplot list.plot" );
1934 #endif
1935 #endif
1936
1937
1938                 nkouho = getKouho( kouho, NKOUHO_LONG, soukan, nlen );
1939
1940 #if 0
1941                 for( i=0; i<nkouho; i++ )
1942                 {
1943                         fprintf( stderr, "kouho[%d] = %d\n", i, kouho[i] );
1944                 }
1945 #endif
1946         }
1947
1948 #if KEIKA
1949         fprintf( stderr, "Searching anchors ... " );
1950 #endif
1951         count = 0;
1952
1953
1954
1955 #define CAND 0
1956 #if CAND
1957         fftfp = fopen( "cand", "w" );
1958         fclose( fftfp );
1959 #endif
1960         if( kobetsubunkatsu )
1961         {
1962                 maxk = 1;
1963                 kouho[0] = 0;
1964         }
1965         else
1966         {
1967                 maxk = nkouho;
1968         }
1969
1970         for( k=0; k<maxk; k++ )
1971         {
1972                 lag = kouho[k];
1973                 if( lag <= -len1 || len2 <= lag ) continue;
1974 //              fprintf( stderr, "k=%d, lag=%d\n", k, lag );
1975                 zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
1976 #if CAND
1977                 fftfp = fopen( "cand", "a" );
1978                 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
1979                 fprintf( fftfp, "%s\n", tmpptr1[0] );
1980                 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
1981                 fprintf( fftfp, "%s\n", tmpptr2[0] );
1982                 fprintf( fftfp, ">\n", k+1, lag );
1983                 fclose( fftfp );
1984 #endif
1985
1986 //              fprintf( stderr, "lag = %d\n", lag );
1987                 tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
1988 //              fprintf( stderr, "lag = %d, %d found\n", lag, tmpint );
1989
1990 //              if( lag == -50 ) exit( 1 );
1991
1992                 if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
1993
1994 //              fprintf( stderr, "##### k=%d / %d\n", k, maxk );
1995 //              if( tmpint == 0 ) break; // 060430 iinoka ? // 090530 yameta
1996                 while( tmpint-- > 0 )
1997                 {
1998 #if 0
1999                         if( segment[count].end - segment[count].start < fftWinSize )
2000                         {
2001                                 count++;
2002                                 continue;
2003                         }
2004 #endif
2005                         if( lag > 0 )
2006                         {
2007                                 segment1[count].start  = segment[count].start ;
2008                                 segment1[count].end    = segment[count].end   ;
2009                                 segment1[count].center = segment[count].center;
2010                                 segment1[count].score  = segment[count].score;
2011
2012                                 segment2[count].start  = segment[count].start  + lag;
2013                                 segment2[count].end    = segment[count].end    + lag;
2014                                 segment2[count].center = segment[count].center + lag;
2015                                 segment2[count].score  = segment[count].score       ;
2016                         }
2017                         else
2018                         {
2019                                 segment1[count].start  = segment[count].start  - lag;
2020                                 segment1[count].end    = segment[count].end    - lag;
2021                                 segment1[count].center = segment[count].center - lag;
2022                                 segment1[count].score  = segment[count].score       ;
2023
2024                                 segment2[count].start  = segment[count].start ;
2025                                 segment2[count].end    = segment[count].end   ;
2026                                 segment2[count].center = segment[count].center;
2027                                 segment2[count].score  = segment[count].score ;
2028                         }
2029 #if 0
2030                         fprintf( stderr, "##### k=%d / %d\n", k, maxk );
2031                         fprintf( stderr, "anchor %d, score = %f\n", count, segment1[count].score );
2032                         fprintf( stderr, "in 1 %d\n", segment1[count].center );
2033                         fprintf( stderr, "in 2 %d\n", segment2[count].center );
2034 #endif
2035                         segment1[count].pair = &segment2[count];
2036                         segment2[count].pair = &segment1[count];
2037                         count++;
2038 #if 0
2039                         fprintf( stderr, "count=%d\n", count );
2040 #endif
2041                 }
2042         }
2043 #if 1
2044         if( !kobetsubunkatsu )
2045                 fprintf( stderr, "done. (%d anchors) ", count );
2046 #endif
2047         if( !count && fftNoAnchStop )
2048                 ErrorExit( "Cannot detect anchor!" );
2049 #if 0
2050         fprintf( stderr, "RESULT before sort:\n" );
2051         for( l=0; l<count+1; l++ )
2052         {
2053                 fprintf( stderr, "cut[%d]=%d, ", l, segment1[l].center );
2054                 fprintf( stderr, "%d score = %f\n", segment2[l].center, segment1[l].score );
2055         }
2056 #endif
2057
2058         for( i=0; i<count; i++ )
2059         {
2060                 sortedseg1[i] = &segment1[i];
2061                 sortedseg2[i] = &segment2[i];
2062         }
2063 #if 0
2064         tmpsort( count, sortedseg1 );
2065         tmpsort( count, sortedseg2 );
2066         qsort( sortedseg1, count, sizeof( Segment * ), segcmp );
2067         qsort( sortedseg2, count, sizeof( Segment * ), segcmp );
2068 #else
2069         mymergesort( 0, count-1, sortedseg1 );
2070         mymergesort( 0, count-1, sortedseg2 );
2071 #endif
2072         for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
2073         for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
2074
2075
2076
2077         if( kobetsubunkatsu )
2078         {
2079                 for( i=0; i<count; i++ )
2080             {
2081                         cut1[i+1] = sortedseg1[i]->center;
2082                         cut2[i+1] = sortedseg2[i]->center;
2083                 }
2084                 cut1[0] = 0;
2085                 cut2[0] = 0;
2086                 cut1[count+1] = len1;
2087                 cut2[count+1] = len2;
2088                 count += 2;
2089         }
2090
2091         else
2092         {
2093                 if( count < 5000 )
2094                 {
2095                         if( crossscoresize < count+2 )
2096                         {
2097                                 crossscoresize = count+2;
2098 #if 1
2099                                 if( fftkeika ) fprintf( stderr, "######allocating crossscore, size = %d\n", crossscoresize );
2100 #endif
2101                                 if( crossscore ) FreeDoubleMtx( crossscore );
2102                                 crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
2103                         }
2104                         for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
2105                                 crossscore[i][j] = 0.0;
2106                         for( i=0; i<count; i++ )
2107                         {
2108                                 crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
2109                                 cut1[i+1] = sortedseg1[i]->center;
2110                                 cut2[i+1] = sortedseg2[i]->center;
2111                         }
2112
2113 #if 0
2114                         fprintf( stderr, "AFTER SORT\n" );
2115                         for( i=0; i<count+1; i++ ) fprintf( stderr, "%d, %d\n", cut1[i], cut2[i] );
2116                         fprintf( stderr, "crossscore = \n" );
2117                         for( i=0; i<count+1; i++ )
2118                         {
2119                                 for( j=0; j<count+1; j++ )
2120                                         fprintf( stderr, "%.0f ", crossscore[i][j] );
2121                                 fprintf( stderr, "\n" );
2122                         }
2123 #endif
2124
2125                         crossscore[0][0] = 10000000.0;
2126                         cut1[0] = 0;
2127                         cut2[0] = 0;
2128                         crossscore[count+1][count+1] = 10000000.0;
2129                         cut1[count+1] = len1;
2130                         cut2[count+1] = len2;
2131                         count += 2;
2132                         count0 = count;
2133
2134 //                      fprintf( stderr, "\n\n\ncalling blockAlign2\n\n\n\n" );
2135                         blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
2136
2137 //                      if( count-count0 )
2138 //                              fprintf( stderr, "%d unused anchors\n", count0-count );
2139
2140                         if( !kobetsubunkatsu && fftkeika )
2141                                 fprintf( stderr, "%d anchors found\n", count );
2142                         if( fftkeika )
2143                         {
2144                                 if( count0 > count )
2145                                 {
2146 #if 0
2147                                         fprintf( stderr, "\7 REPEAT!? \n" );
2148 #else
2149                                         fprintf( stderr, "REPEAT!? \n" );
2150 #endif
2151                                         if( fftRepeatStop ) exit( 1 );
2152                                 }
2153 #if KEIKA
2154                                 else fprintf( stderr, "done\n" );
2155 #endif
2156                         }
2157                 }
2158
2159
2160                 else
2161                 {
2162                         fprintf( stderr, "\nMany anchors were found. The upper-level DP is skipped.\n\n" );
2163
2164                         cut1[0] = 0;
2165                         cut2[0] = 0;
2166                         count0 = 0;
2167                         for( i=0; i<count; i++ )
2168                         {
2169 //                              fprintf( stderr, "i=%d, %d-%d ?\n", i, sortedseg1[i]->center, sortedseg1[i]->pair->center );
2170                                 if( sortedseg1[i]->center > cut1[count0]
2171                                  && sortedseg1[i]->pair->center > cut2[count0] )
2172                                 {
2173                                         count0++;
2174                                         cut1[count0] = sortedseg1[i]->center;
2175                                         cut2[count0] = sortedseg1[i]->pair->center;
2176                                 }
2177                                 else
2178                                 {
2179                                         if( i && sortedseg1[i]->score > sortedseg1[i-1]->score )
2180                                         {
2181                                                 if( sortedseg1[i]->center > cut1[count0-1]
2182                                                  && sortedseg1[i]->pair->center > cut2[count0-1] )
2183                                                 {
2184                                                         cut1[count0] = sortedseg1[i]->center;
2185                                                         cut2[count0] = sortedseg1[i]->pair->center;
2186                                                 }
2187                                                 else
2188                                                 {
2189 //                                                      count0--;
2190                                                 }
2191                                         }
2192                                 }
2193                         }
2194 //                      if( count-count0 )
2195 //                              fprintf( stderr, "%d anchors unused\n", count-count0 );
2196                         cut1[count0+1] = len1;
2197                         cut2[count0+1] = len2;
2198                         count = count0 + 2;
2199                         count0 = count;
2200
2201                 }
2202         }
2203
2204 //      exit( 0 );
2205
2206 #if 0
2207         fftfp = fopen( "fft", "a" );
2208         fprintf( fftfp, "RESULT after sort:\n" );
2209         for( l=0; l<count; l++ )
2210         {
2211                 fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
2212                 fprintf( fftfp, "%d\n", segment2[l].center );
2213         }
2214         fclose( fftfp );
2215 #endif
2216
2217 #if 0
2218         fprintf( stderr, "RESULT after blckalign:\n" );
2219         for( l=0; l<count+1; l++ )
2220         {
2221                 fprintf( stderr, "cut : %d %d\n", cut1[l], cut2[l] );
2222         }
2223 #endif
2224
2225 #if 0
2226         fprintf( trap_g, "Devided to %d segments\n", count-1 );
2227         fprintf( trap_g, "%d  %d forg\n", MIN( clus1, clus2 ), count-1 );
2228 #endif
2229
2230         totallen = 0;
2231         for( j=0; j<clus1; j++ ) result1[j][0] = 0;
2232         for( j=0; j<clus2; j++ ) result2[j][0] = 0;
2233         totalscore = 0.0;
2234         *fftlog = -1;
2235         for( i=0; i<count-1; i++ )
2236         {
2237                 *fftlog += 1;
2238                 if( i == 0 ) headgp = outgap; else headgp = 1;
2239                 if( i == count-2 ) tailgp = outgap; else tailgp = 1;
2240
2241                 if( cut1[i] )
2242                 {
2243 //                      getkyokaigap( sgap1, seq1, cut1[i]-1, clus1 );
2244 //                      getkyokaigap( sgap2, seq2, cut2[i]-1, clus2 );
2245                         getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
2246                         getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
2247                 }
2248                 else
2249                 {
2250                         for( j=0; j<clus1; j++ ) sgap1[j] = 'o';
2251                         for( j=0; j<clus2; j++ ) sgap2[j] = 'o';
2252                 }
2253                 if( cut1[i+1] != len1 )
2254                 {
2255                         getkyokaigap( egap1, seq1, cut1[i+1], clus1 );
2256                         getkyokaigap( egap2, seq2, cut2[i+1], clus2 );
2257                 }
2258                 else
2259                 {
2260                         for( j=0; j<clus1; j++ ) egap1[j] = 'o';
2261                         for( j=0; j<clus2; j++ ) egap2[j] = 'o';
2262                 }
2263 #if DEBUG
2264                 fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
2265 #else
2266 #if 1
2267                 fprintf( stderr, "DP %05d / %05d \b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b", i+1, count-1 );
2268 #endif
2269 #endif
2270                 for( j=0; j<clus1; j++ )
2271                 {
2272                         strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
2273                         tmpres1[j][cut1[i+1]-cut1[i]] = 0;
2274                 }
2275                 if( kobetsubunkatsu && fftkeika ) commongappick( clus1, tmpres1 ); //dvtditr \e$B$K8F$P$l$?$H$-\e(B fftkeika=1
2276 //              if( kobetsubunkatsu ) commongappick( clus1, tmpres1 );
2277                 for( j=0; j<clus2; j++ )
2278                 {
2279 //                      fprintf( stderr, "### cut2[i+1]-cut2[i] = %d\n", cut2[i+1]-cut2[i] );
2280                         if( cut2[i+1]-cut2[i] <= 0 )
2281                                 fprintf( stderr, "### cut2[i+1]=%d, cut2[i]=%d\n", cut2[i+1], cut2[i] );
2282                         strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
2283                         tmpres2[j][cut2[i+1]-cut2[i]] = 0;
2284                 }
2285                 if( kobetsubunkatsu && fftkeika ) commongappick( clus2, tmpres2 ); //dvtditr \e$B$K8F$P$l$?$H$-\e(B fftkeika=1
2286 //              if( kobetsubunkatsu ) commongappick( clus2, tmpres2 );
2287
2288                 if( constraint )
2289                 {
2290                         fprintf( stderr, "Not supported\n" );
2291                         exit( 1 );
2292                 }
2293 #if 0
2294                 fprintf( stderr, "i=%d, before alignment", i );
2295                 fprintf( stderr, "%4d\n", totallen );
2296                 fprintf( stderr, "\n\n" );
2297                 for( j=0; j<clus1; j++ )
2298                 {
2299                         fprintf( stderr, "%s\n", tmpres1[j] );
2300                 }
2301                 fprintf( stderr, "-------\n" );
2302                 for( j=0; j<clus2; j++ )
2303                 {
2304                         fprintf( stderr, "%s\n", tmpres2[j] );
2305                 }
2306 #endif
2307
2308 #if 0
2309                 fprintf( stdout, "writing input\n" );
2310                 for( j=0; j<clus1; j++ )
2311                 {
2312                         fprintf( stdout, ">%d of GROUP1\n", j );
2313                         fprintf( stdout, "%s\n", tmpres1[j] );
2314                 }
2315                 for( j=0; j<clus2; j++ )
2316                 {
2317                         fprintf( stdout, ">%d of GROUP2\n", j );
2318                         fprintf( stdout, "%s\n", tmpres2[j] );
2319                 }
2320                 fflush( stdout );
2321 #endif
2322                 switch( alg )
2323                 {
2324                         case( 'M' ):
2325                                         totalscore += MSalignmm( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, sgap1, sgap2, egap1, egap2, NULL, 0, NULL, headgp, tailgp );
2326                                 break;
2327                         default:
2328                                 fprintf( stderr, "alg = %c\n", alg );
2329                                 ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
2330                                 break;
2331                 }
2332
2333                 nlen = strlen( tmpres1[0] );
2334                 if( totallen + nlen > alloclen )
2335                 {
2336                         fprintf( stderr, "totallen=%d +  nlen=%d > alloclen = %d\n", totallen, nlen, alloclen );
2337                         ErrorExit( "LENGTH OVER in Falign\n " );
2338                 }
2339                 for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
2340                 for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
2341                 totallen += nlen;
2342 #if 0
2343                 fprintf( stderr, "i=%d", i );
2344                 fprintf( stderr, "%4d\n", totallen );
2345                 fprintf( stderr, "\n\n" );
2346                 for( j=0; j<clus1; j++ )
2347                 {
2348                         fprintf( stderr, "%s\n", tmpres1[j] );
2349                 }
2350                 fprintf( stderr, "-------\n" );
2351                 for( j=0; j<clus2; j++ )
2352                 {
2353                         fprintf( stderr, "%s\n", tmpres2[j] );
2354                 }
2355 #endif
2356         }
2357 #if KEIKA
2358         fprintf( stderr, "DP ... done   \n" );
2359 #endif
2360
2361         for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
2362         for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
2363 #if 0
2364         for( j=0; j<clus1; j++ )
2365         {
2366                 fprintf( stderr, "%s\n", result1[j] );
2367         }
2368         fprintf( stderr, "- - - - - - - - - - -\n" );
2369         for( j=0; j<clus2; j++ )
2370         {
2371                 fprintf( stderr, "%s\n", result2[j] );
2372         }
2373 #endif
2374         return( totalscore );
2375 }