binaries/src/mafft/core/Falign.c

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