binaries/src/mafft/core/Falign_localhom.c

   1 #include "mltaln.h"
   2
   3 //static FILE *fftfp;
   4 static TLS int n20or4or2;
   5
   6 #define KEIKA 0
   7 #define RND   0
   8 #define DEBUG 0
   9
  10 extern int fft( int, Fukusosuu *, int );
  11
  12
  13 #if 0
  14 static void generateRndSeq( char *seq, int len )
  15 {
  16         while( len-- )
  17 #if 1
  18                 *seq++ = (int)( rnd() * n20or4or2 );
  19 #else
  20                 *seq++ = (int)1;
  21 #endif
  22 }
  23 #endif
  24
  25 static void vec_init( Fukusosuu *result, int nlen )
  26 {
  27         while( nlen-- )
  28         {
  29                 result->R = result->I = 0.0;
  30                 result++;
  31         }
  32 }
  33
  34 #if 0
  35 static void vec_init2( Fukusosuu **result, char *seq, double eff, int st, int ed )
  36 {
  37         int i;
  38         for( i=st; i<ed; i++ )
  39                 result[(int)*seq++][i].R += eff;
  40 }
  41 #endif
  42
  43 static void seq_vec_2( Fukusosuu *result, double *score, double incr, char *seq )
  44 {
  45         static TLS int n;
  46         for( ; *seq; result++ )
  47         {
  48                 n = amino_n[(unsigned char)*seq++];
  49                 if( n < 20 && n >= 0 ) result->R += incr * score[n];
  50 #if 0
  51                 fprintf( stderr, "n=%d, score=%f, inc=%f R=%f\n",n,  score[n], incr * score[n], result->R );
  52 #endif
  53         }
  54 }
  55
  56 static void seq_vec_3( Fukusosuu **result, double incr, char *seq )
  57 {
  58         int i;
  59         int n;
  60         for( i=0; *seq; i++ )
  61         {
  62                 n = amino_n[(unsigned char)*seq++];
  63                 if( n < n20or4or2 && n >= 0 ) result[n][i].R += incr;
  64         }
  65 }
  66
  67
  68 #if 0
  69 static void seq_vec( Fukusosuu *result, char query, double incr, char *seq )
  70 {
  71 #if 0
  72         int bk = nlen;
  73 #endif
  74         while( *seq )
  75         {
  76                 if( *seq++ == query ) result->R += incr;
  77                 result++;
  78 #if 0
  79 fprintf( stderr, "i = %d result->R = %f\n", bk-nlen, (result-1)->R );
  80 #endif
  81         }
  82 }
  83
  84 static int checkRepeat( int num, int *cutpos )
  85 {
  86         int tmp, buf;
  87
  88         buf = *cutpos;
  89         while( num-- )
  90         {
  91                 if( ( tmp = *cutpos++ ) < buf ) return( 1 );
  92                 buf = tmp;
  93         }
  94         return( 0 );
  95 }
  96
  97 static int segcmp( void *ptr1, void *ptr2 )
  98 {
  99         int diff;
 100         Segment **seg1 = (Segment **)ptr1;
 101         Segment **seg2 = (Segment **)ptr2;
 102 #if 0
 103         return( (*seg1)->center - (*seg2)->center );
 104 #else
 105         diff = (*seg1)->center - (*seg2)->center;
 106         if( diff ) return( diff );
 107
 108         diff = (*seg1)->start - (*seg2)->start;
 109         if( diff ) return( diff );
 110
 111         diff = (*seg1)->end - (*seg2)->end;
 112         if( diff ) return( diff );
 113
 114         fprintf( stderr, "USE STABLE SORT !!\n" );
 115         exit( 1 );
 116         return( 0 );
 117 #endif
 118 }
 119
 120 #endif
 121
 122
 123 static void mymergesort( int first, int last, Segment **seg )
 124 {
 125         int middle;
 126         static TLS int i, j, k, p;
 127         static TLS int allo = 0;
 128         static TLS Segment **work = NULL;
 129
 130         if( seg == NULL )
 131         {
 132                 free( work ); work = NULL;
 133                 return;
 134         }
 135
 136         if( last > allo )
 137         {
 138                 allo = last;
 139                 if( work ) free( work );
 140                 work = (Segment **)calloc( allo / 2 + 1, sizeof( Segment *) );
 141         }
 142
 143         if( first < last )
 144         {
 145                 middle = ( first + last ) / 2;
 146                 mymergesort( first, middle, seg );
 147                 mymergesort( middle+1, last, seg );
 148                 p = 0;
 149                 for( i=first; i<=middle; i++ ) work[p++] = seg[i];
 150                 i = middle + 1; j = 0; k = first;
 151                 while( i <= last && j < p )
 152                 {
 153                         if( work[j]->center <= seg[i]->center )
 154                                 seg[k++] = work[j++];
 155                         else
 156                                 seg[k++] = seg[i++];
 157                 }
 158                 while( j < p ) seg[k++] = work[j++];
 159         }
 160 }
 161
 162
 163 double Falign_localhom( int **whichmtx, double ***scoringmatrices, double **n_dynamicmtx,
 164                           char  **seq1, char  **seq2,
 165                           double *eff1, double *eff2,
 166                           double **eff1s, double **eff2s,
 167                           int    clus1, int    clus2,
 168                           int alloclen,
 169                            LocalHom ***localhom, double *totalimpmatch,
 170                            int *gapmap1, int *gapmap2,
 171                                 int *chudanpt, int chudanref, int *chudanres )
 172 {
 173    // tditeration.c deha alloclen ha huhen nanode
 174    // prevalloclen ha iranai.
 175         int i, j, k, l, m, maxk;
 176         int nlen, nlen2, nlen4;
 177         static TLS int crossscoresize = 0;
 178         static TLS char **tmpseq1 = NULL;
 179         static TLS char **tmpseq2 = NULL;
 180         static TLS char **tmpptr1 = NULL;
 181         static TLS char **tmpptr2 = NULL;
 182         static TLS char **tmpres1 = NULL;
 183         static TLS char **tmpres2 = NULL;
 184         static TLS char **result1 = NULL;
 185         static TLS char **result2 = NULL;
 186 #if RND
 187         static TLS char **rndseq1 = NULL;
 188         static TLS char **rndseq2 = NULL;
 189 #endif
 190         static TLS Fukusosuu **seqVector1 = NULL;
 191         static TLS Fukusosuu **seqVector2 = NULL;
 192         static TLS Fukusosuu **naiseki = NULL;
 193         static TLS Fukusosuu *naisekiNoWa = NULL;
 194         static TLS double *soukan = NULL;
 195         static TLS double **crossscore = NULL;
 196         int nlentmp;
 197         static TLS int *kouho = NULL;
 198         static TLS Segment *segment = NULL;
 199         static TLS Segment *segment1 = NULL;
 200         static TLS Segment *segment2 = NULL;
 201         static TLS Segment **sortedseg1 = NULL;
 202         static TLS Segment **sortedseg2 = NULL;
 203         static TLS int *cut1 = NULL;
 204         static TLS int *cut2 = NULL;
 205         static TLS char *sgap1, *egap1, *sgap2, *egap2;
 206         static TLS int localalloclen = 0;
 207         int lag;
 208         int tmpint;
 209         int count, count0;
 210         int len1, len2;
 211         int totallen;
 212         double totalscore;
 213         double impmatch;
 214
 215         extern Fukusosuu   *AllocateFukusosuuVec();
 216         extern Fukusosuu  **AllocateFukusosuuMtx();
 217
 218         if( seq1 == NULL )
 219         {
 220                 if( result1 )
 221                 {
 222 //                      fprintf( stderr, "Freeing localarrays in Falign\n" );
 223                         localalloclen = 0;
 224                         crossscoresize = 0;
 225                         mymergesort( 0, 0, NULL );
 226                         alignableReagion( 0, 0, NULL, NULL, NULL, NULL, NULL );
 227                         fft( 0, NULL, 1 );
 228 //                      A__align( NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0, -1, -1 ); // iru?
 229                         G__align11( NULL, NULL, NULL, 0, 0, 0 );
 230                         partA__align( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL );
 231                         partA__align_variousdist( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL );
 232                         blockAlign2( NULL, NULL, NULL, NULL, NULL, NULL );
 233                         if( crossscore ) FreeDoubleMtx( crossscore );
 234                         FreeCharMtx( result1 );
 235                         FreeCharMtx( result2 );
 236                         FreeCharMtx( tmpres1 );
 237                         FreeCharMtx( tmpres2 );
 238                         FreeCharMtx( tmpseq1 );
 239                         FreeCharMtx( tmpseq2 );
 240                         free( sgap1 );
 241                         free( egap1 );
 242                         free( sgap2 );
 243                         free( egap2 );
 244                         free( kouho );
 245                         free( cut1 );
 246                         free( cut2 );
 247                         free( tmpptr1 );
 248                         free( tmpptr2 );
 249                         free( segment );
 250                         free( segment1 );
 251                         free( segment2 );
 252                         free( sortedseg1 );
 253                         free( sortedseg2 );
 254                         if( !kobetsubunkatsu )
 255                         {
 256                                 FreeFukusosuuMtx ( seqVector1 );
 257                                 FreeFukusosuuMtx ( seqVector2 );
 258                                 FreeFukusosuuVec( naisekiNoWa );
 259                                 FreeFukusosuuMtx( naiseki );
 260                                 FreeDoubleVec( soukan );
 261                         }
 262                 }
 263                 else
 264                 {
 265 //                      fprintf( stderr, "Did not allocate localarrays in Falign\n" );
 266                 }
 267
 268                 return( 0.0 );
 269         }
 270
 271         len1 = strlen( seq1[0] );
 272         len2 = strlen( seq2[0] );
 273         nlentmp = MAX( len1, len2 );
 274
 275         nlen = 1;
 276         while( nlentmp >= nlen ) nlen <<= 1;
 277 #if 0
 278         fprintf( stderr, "###   nlen    = %d\n", nlen );
 279 #endif
 280
 281         nlen2 = nlen/2; nlen4 = nlen2 / 2;
 282
 283 #if DEBUG
 284         fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
 285         fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
 286 #endif
 287
 288         if( !localalloclen )
 289         {
 290                 sgap1 = AllocateCharVec( njob );
 291                 egap1 = AllocateCharVec( njob );
 292                 sgap2 = AllocateCharVec( njob );
 293                 egap2 = AllocateCharVec( njob );
 294                 kouho = AllocateIntVec( NKOUHO );
 295                 cut1 = AllocateIntVec( MAXSEG );
 296                 cut2 = AllocateIntVec( MAXSEG );
 297                 tmpptr1 = AllocateCharMtx( njob, 0 );
 298                 tmpptr2 = AllocateCharMtx( njob, 0 );
 299                 result1 = AllocateCharMtx( njob, alloclen );
 300                 result2 = AllocateCharMtx( njob, alloclen );
 301                 tmpres1 = AllocateCharMtx( njob, alloclen );
 302                 tmpres2 = AllocateCharMtx( njob, alloclen );
 303 //              crossscore = AllocateDoubleMtx( MAXSEG, MAXSEG );
 304                 segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
 305                 segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
 306                 segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
 307                 sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
 308                 sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
 309                 if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
 310                         ErrorExit( "Allocation error\n" );
 311
 312                 if     ( scoremtx == -1 ) n20or4or2 = 4;
 313                 else if( fftscore == 1  ) n20or4or2 = 2;
 314                 else                      n20or4or2 = 20;
 315         }
 316         if( localalloclen < nlen )
 317         {
 318                 if( localalloclen )
 319                 {
 320 #if 1
 321                         if( !kobetsubunkatsu )
 322                         {
 323                                 FreeFukusosuuMtx ( seqVector1 );
 324                                 FreeFukusosuuMtx ( seqVector2 );
 325                                 FreeFukusosuuVec( naisekiNoWa );
 326                                 FreeFukusosuuMtx( naiseki );
 327                                 FreeDoubleVec( soukan );
 328                         }
 329                         FreeCharMtx( tmpseq1 );
 330                         FreeCharMtx( tmpseq2 );
 331 #endif
 332 #if RND
 333                         FreeCharMtx( rndseq1 );
 334                         FreeCharMtx( rndseq2 );
 335 #endif
 336                 }
 337
 338                 tmpseq1 = AllocateCharMtx( njob, nlen );
 339                 tmpseq2 = AllocateCharMtx( njob, nlen );
 340                 if( !kobetsubunkatsu )
 341                 {
 342                         naisekiNoWa = AllocateFukusosuuVec( nlen );
 343                         naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
 344                         seqVector1 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
 345                         seqVector2 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
 346                         soukan = AllocateDoubleVec( nlen+1 );
 347                 }
 348 #if RND
 349                 rndseq1 = AllocateCharMtx( njob, nlen );
 350                 rndseq2 = AllocateCharMtx( njob, nlen );
 351                 for( i=0; i<njob; i++ )
 352                 {
 353                         generateRndSeq( rndseq1[i], nlen );
 354                         generateRndSeq( rndseq2[i], nlen );
 355                 }
 356 #endif
 357                 localalloclen = nlen;
 358         }
 359
 360         for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
 361         for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
 362
 363 #if 0
 364 fftfp = fopen( "input_of_Falign", "w" );
 365 fprintf( fftfp, "nlen = %d\n", nlen );
 366 fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
 367 for( i=0; i<clus1; i++ )
 368         fprintf( fftfp, "%s\n", seq1[i] );
 369 fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
 370 for( i=0; i<clus2; i++ )
 371         fprintf( fftfp, "%s\n", seq2[i] );
 372 fclose( fftfp );
 373 system( "less input_of_Falign < /dev/tty > /dev/tty" );
 374 #endif
 375         if( !kobetsubunkatsu )
 376         {
 377                 fprintf( stderr,  "FFT ... " );
 378
 379                 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
 380                 if( fftscore && scoremtx != -1 )
 381                 {
 382                         for( i=0; i<clus1; i++ )
 383                         {
 384                                 seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
 385                                 seq_vec_2( seqVector1[1], volume,   eff1[i], tmpseq1[i] );
 386                         }
 387                 }
 388                 else
 389                 {
 390 #if 0
 391                         for( i=0; i<clus1; i++ ) for( j=0; j<n20or4or2; j++ )
 392                                 seq_vec( seqVector1[j], amino[j], eff1[i], tmpseq1[i] );
 393 #else
 394                         for( i=0; i<clus1; i++ )
 395                                 seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
 396 #endif
 397                 }
 398 #if RND
 399                 for( i=0; i<clus1; i++ )
 400                 {
 401                         vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
 402                 }
 403 #endif
 404 #if 0
 405 fftfp = fopen( "seqVec", "w" );
 406 fprintf( fftfp, "before transform\n" );
 407 for( k=0; k<n20or4or2; k++ )
 408 {
 409    fprintf( fftfp, "nlen=%d\n", nlen );
 410    fprintf( fftfp, "%c\n", amino[k] );
 411    for( l=0; l<nlen; l++ )
 412    fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
 413 }
 414 fclose( fftfp );
 415 system( "less seqVec < /dev/tty > /dev/tty" );
 416 #endif
 417
 418                 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
 419                 if( fftscore && scoremtx != -1 )
 420                 {
 421                         for( i=0; i<clus2; i++ )
 422                         {
 423                                 seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
 424                                 seq_vec_2( seqVector2[1], volume,   eff2[i], tmpseq2[i] );
 425                         }
 426                 }
 427                 else
 428                 {
 429 #if 0
 430                         for( i=0; i<clus2; i++ ) for( j=0; j<n20or4or2; j++ )
 431                                 seq_vec( seqVector2[j], amino[j], eff2[i], tmpseq2[i] );
 432 #else
 433                         for( i=0; i<clus2; i++ )
 434                                 seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
 435 #endif
 436                 }
 437 #if RND
 438                 for( i=0; i<clus2; i++ )
 439                 {
 440                         vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
 441                 }
 442 #endif
 443
 444 #if 0
 445 fftfp = fopen( "seqVec2", "w" );
 446 fprintf( fftfp, "before fft\n" );
 447 for( k=0; k<n20or4or2; k++ )
 448 {
 449    fprintf( fftfp, "%c\n", amino[k] );
 450    for( l=0; l<nlen; l++ )
 451    fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
 452 }
 453 fclose( fftfp );
 454 system( "less seqVec2 < /dev/tty > /dev/tty" );
 455 #endif
 456
 457                 for( j=0; j<n20or4or2; j++ )
 458                 {
 459                         fft( nlen, seqVector2[j], (j==0) );
 460                         fft( nlen, seqVector1[j], 0 );
 461                 }
 462 #if 0
 463 fftfp = fopen( "seqVec2", "w" );
 464 fprintf( fftfp, "#after fft\n" );
 465 for( k=0; k<n20or4or2; k++ )
 466 {
 467    fprintf( fftfp, "#%c\n", amino[k] );
 468    for( l=0; l<nlen; l++ )
 469            fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
 470 }
 471 fclose( fftfp );
 472 system( "less seqVec2 < /dev/tty > /dev/tty" );
 473 #endif
 474
 475                 for( k=0; k<n20or4or2; k++ )
 476                 {
 477                         for( l=0; l<nlen; l++ )
 478                                 calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
 479                 }
 480                 for( l=0; l<nlen; l++ )
 481                 {
 482                         naisekiNoWa[l].R = 0.0;
 483                         naisekiNoWa[l].I = 0.0;
 484                         for( k=0; k<n20or4or2; k++ )
 485                         {
 486                                 naisekiNoWa[l].R += naiseki[k][l].R;
 487                                 naisekiNoWa[l].I += naiseki[k][l].I;
 488                         }
 489                 }
 490
 491 #if 0
 492         fftfp = fopen( "naisekiNoWa", "w" );
 493         fprintf( fftfp, "#Before fft\n" );
 494         for( l=0; l<nlen; l++ )
 495                 fprintf( fftfp, "%d  %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I );
 496         fclose( fftfp );
 497         system( "less naisekiNoWa < /dev/tty > /dev/tty " );
 498 #endif
 499
 500                 fft( -nlen, naisekiNoWa, 0 );
 501
 502                 for( m=0; m<=nlen2; m++ )
 503                         soukan[m] = naisekiNoWa[nlen2-m].R;
 504                 for( m=nlen2+1; m<nlen; m++ )
 505                         soukan[m] = naisekiNoWa[nlen+nlen2-m].R;
 506
 507 #if 0
 508         fftfp = fopen( "naisekiNoWa", "w" );
 509         fprintf( fftfp, "#After fft\n" );
 510         for( l=0; l<nlen; l++ )
 511                 fprintf( fftfp, "%d  %f\n", l, naisekiNoWa[l].R );
 512         fclose( fftfp );
 513         fftfp = fopen( "list.plot", "w"  );
 514         fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
 515         fclose( fftfp );
 516         system( "/usr/bin/gnuplot list.plot &" );
 517 #endif
 518 #if 0
 519         fprintf( stderr, "frt write start\n" );
 520         fftfp = fopen( "frt", "w" );
 521         for( l=0; l<nlen; l++ )
 522                 fprintf( fftfp, "%d  %f\n", l-nlen2, soukan[l] );
 523         fclose( fftfp );
 524         system( "less frt < /dev/tty > /dev/tty" );
 525 #if 0
 526         fftfp = fopen( "list.plot", "w"  );
 527         fprintf( fftfp, "plot 'frt'\n pause +1" );
 528         fclose( fftfp );
 529         system( "/usr/bin/gnuplot list.plot" );
 530 #endif
 531 #endif
 532
 533
 534                 getKouho( kouho, NKOUHO, soukan, nlen );
 535
 536 #if 0
 537                 for( i=0; i<NKOUHO; i++ )
 538                 {
 539                         fprintf( stderr, "kouho[%d] = %d\n", i, kouho[i] );
 540                 }
 541 #endif
 542         }
 543
 544 #if KEIKA
 545         fprintf( stderr, "Searching anchors ... " );
 546 #endif
 547         count = 0;
 548
 549
 550
 551 #define CAND 0
 552 #if CAND
 553         fftfp = fopen( "cand", "w" );
 554         fclose( fftfp );
 555 #endif
 556         if( kobetsubunkatsu )
 557         {
 558                 maxk = 1;
 559                 kouho[0] = 0;
 560         }
 561         else
 562         {
 563                 maxk = NKOUHO;
 564         }
 565
 566         for( k=0; k<maxk; k++ )
 567         {
 568                 lag = kouho[k];
 569                 zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
 570 #if CAND
 571                 fftfp = fopen( "cand", "a" );
 572                 fprintf( fftfp, "Candidate No.%d lag = %d\n", k+1, lag );
 573                 fprintf( fftfp, "%s\n", tmpptr1[0] );
 574                 fprintf( fftfp, "%s\n", tmpptr2[0] );
 575                 fclose( fftfp );
 576 #endif
 577                 tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
 578
 579                 if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
 580
 581
 582                 while( tmpint-- > 0 )
 583                 {
 584                         if( lag > 0 )
 585                         {
 586                                 segment1[count].start  = segment[count].start ;
 587                                 segment1[count].end    = segment[count].end   ;
 588                                 segment1[count].center = segment[count].center;
 589                                 segment1[count].score  = segment[count].score;
 590
 591                                 segment2[count].start  = segment[count].start  + lag;
 592                                 segment2[count].end    = segment[count].end    + lag;
 593                                 segment2[count].center = segment[count].center + lag;
 594                                 segment2[count].score  = segment[count].score       ;
 595                         }
 596                         else
 597                         {
 598                                 segment1[count].start  = segment[count].start  - lag;
 599                                 segment1[count].end    = segment[count].end    - lag;
 600                                 segment1[count].center = segment[count].center - lag;
 601                                 segment1[count].score  = segment[count].score       ;
 602
 603                                 segment2[count].start  = segment[count].start ;
 604                                 segment2[count].end    = segment[count].end   ;
 605                                 segment2[count].center = segment[count].center;
 606                                 segment2[count].score  = segment[count].score ;
 607                         }
 608 #if 0
 609                         fftfp = fopen( "cand", "a" );
 610                         fprintf( fftfp, "Goukaku=%dko\n", tmpint );
 611                         fprintf( fftfp, "in 1 %d\n", segment1[count].center );
 612                         fprintf( fftfp, "in 2 %d\n", segment2[count].center );
 613                         fclose( fftfp );
 614 #endif
 615                         segment1[count].pair = &segment2[count];
 616                         segment2[count].pair = &segment1[count];
 617                         count++;
 618 #if 0
 619                         fprintf( stderr, "count=%d\n", count );
 620 #endif
 621                 }
 622         }
 623 #if 1
 624         if( !kobetsubunkatsu )
 625                 fprintf( stderr, "%d segments found\n", count );
 626 #endif
 627         if( !count && fftNoAnchStop )
 628                 ErrorExit( "Cannot detect anchor!" );
 629 #if 0
 630         fftfp = fopen( "fft", "a" );
 631         fprintf( fftfp, "RESULT before sort:\n" );
 632         for( l=0; l<count; l++ )
 633         {
 634                 fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
 635                 fprintf( fftfp, "%d score = %f\n", segment2[l].center, segment1[l].score );
 636         }
 637         fclose( fftfp );
 638 #endif
 639
 640 #if KEIKA
 641         fprintf( stderr, "Aligning anchors ... " );
 642 #endif
 643         for( i=0; i<count; i++ )
 644         {
 645                 sortedseg1[i] = &segment1[i];
 646                 sortedseg2[i] = &segment2[i];
 647         }
 648 #if 0
 649         tmpsort( count, sortedseg1 );
 650         tmpsort( count, sortedseg2 );
 651         qsort( sortedseg1, count, sizeof( Segment * ), segcmp );
 652         qsort( sortedseg2, count, sizeof( Segment * ), segcmp );
 653 #else
 654         mymergesort( 0, count-1, sortedseg1 );
 655         mymergesort( 0, count-1, sortedseg2 );
 656 #endif
 657         for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
 658         for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
 659
 660
 661         if( kobetsubunkatsu )
 662         {
 663                 for( i=0; i<count; i++ )
 664             {
 665                         cut1[i+1] = sortedseg1[i]->center;
 666                         cut2[i+1] = sortedseg2[i]->center;
 667                 }
 668                 cut1[0] = 0;
 669                 cut2[0] = 0;
 670                 cut1[count+1] = len1;
 671                 cut2[count+1] = len2;
 672                 count += 2;
 673         }
 674         else
 675         {
 676                 if( crossscoresize < count+2 )
 677                 {
 678                         crossscoresize = count+2;
 679 #if 1
 680                         fprintf( stderr, "######allocating crossscore, size = %d\n", crossscoresize );
 681 #endif
 682                         if( crossscore ) FreeDoubleMtx( crossscore );
 683                         crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
 684                 }
 685                 for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
 686                         crossscore[i][j] = 0.0;
 687                 for( i=0; i<count; i++ )
 688                 {
 689                         crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
 690                         cut1[i+1] = sortedseg1[i]->center;
 691                         cut2[i+1] = sortedseg2[i]->center;
 692                 }
 693
 694 #if DEBUG
 695                 fprintf( stderr, "AFTER SORT\n" );
 696                 for( i=0; i<count; i++ ) fprintf( stderr, "%d, %d\n", segment1[i].start, segment2[i].start );
 697 #endif
 698
 699                 crossscore[0][0] = 10000000.0;
 700                 cut1[0] = 0;
 701                 cut2[0] = 0;
 702                 crossscore[count+1][count+1] = 10000000.0;
 703                 cut1[count+1] = len1;
 704                 cut2[count+1] = len2;
 705                 count += 2;
 706                 count0 = count;
 707
 708                 blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
 709                 if( count0 > count )
 710                 {
 711 #if 0
 712                         fprintf( stderr, "\7 REPEAT!? \n" );
 713 #else
 714                         fprintf( stderr, "REPEAT!? \n" );
 715 #endif
 716                         if( fftRepeatStop ) exit( 1 );
 717                 }
 718 #if KEIKA
 719                 else fprintf( stderr, "done\n" );
 720 #endif
 721         }
 722
 723 #if 0
 724         fftfp = fopen( "fft", "a" );
 725         fprintf( fftfp, "RESULT after sort:\n" );
 726         for( l=0; l<count; l++ )
 727         {
 728                 fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
 729                 fprintf( fftfp, "%d\n", segment2[l].center );
 730         }
 731         fclose( fftfp );
 732 #endif
 733
 734 #if 0
 735         fftfp = fopen( "fft", "a" );
 736         fprintf( fftfp, "RESULT after sort:\n" );
 737         for( l=0; l<count; l++ )
 738         {
 739                 fprintf( fftfp, "cut : %d %d\n", cut1[l], cut2[l] );
 740         }
 741         fclose( fftfp );
 742 #endif
 743
 744 #if KEIKA
 745         fprintf( trap_g, "Devided to %d segments\n", count-1 );
 746         fprintf( trap_g, "%d  %d forg\n", MIN( clus1, clus2 ), count-1 );
 747 #endif
 748
 749         totallen = 0;
 750         for( j=0; j<clus1; j++ ) result1[j][0] = 0;
 751         for( j=0; j<clus2; j++ ) result2[j][0] = 0;
 752         totalscore = 0.0;
 753         *totalimpmatch = 0.0;
 754         for( i=0; i<count-1; i++ )
 755         {
 756 #if DEBUG
 757                 fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
 758 #else
 759 #if KEIKA
 760                 fprintf( stderr, "DP %03d / %03d\r", i+1, count-1 );
 761 #endif
 762 #endif
 763
 764                 if( cut1[i] )
 765                 {
 766                         getkyokaigap( sgap1, seq1, cut1[i]-1, clus1 );
 767                         getkyokaigap( sgap2, seq2, cut2[i]-1, clus2 );
 768                 }
 769                 else
 770                 {
 771                         for( j=0; j<clus1; j++ ) sgap1[j] = 'o';
 772                         for( j=0; j<clus2; j++ ) sgap2[j] = 'o';
 773                 }
 774                 if( cut1[i+1] != len1 )
 775                 {
 776                         getkyokaigap( egap1, seq1, cut1[i+1], clus1 );
 777                         getkyokaigap( egap2, seq2, cut2[i+1], clus2 );
 778                 }
 779                 else
 780                 {
 781                         for( j=0; j<clus1; j++ ) egap1[j] = 'o';
 782                         for( j=0; j<clus2; j++ ) egap2[j] = 'o';
 783                 }
 784
 785                 for( j=0; j<clus1; j++ )
 786                 {
 787                         strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
 788                         tmpres1[j][cut1[i+1]-cut1[i]] = 0;
 789                 }
 790                 if( kobetsubunkatsu ) commongappick_record( clus1, tmpres1, gapmap1 );
 791                 for( j=0; j<clus2; j++ )
 792                 {
 793                         strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
 794                         tmpres2[j][cut2[i+1]-cut2[i]] = 0;
 795                 }
 796                 if( kobetsubunkatsu ) commongappick_record( clus2, tmpres2, gapmap2 );
 797
 798 #if 0
 799                 fprintf( stderr, "count = %d\n", count );
 800                 fprintf( stderr, "### reg1 = %d-%d\n", cut1[i], cut1[i+1]-1 );
 801                 fprintf( stderr, "### reg2 = %d-%d\n", cut2[i], cut2[i+1]-1 );
 802 #endif
 803
 804                 switch( alg )
 805                 {
 806                         case( 'a' ):
 807                                 totalscore += Aalign( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen );
 808                                 break;
 809                         case( 'A' ):
 810                                 if( scoringmatrices ) // called by tditeration.c
 811                                 {
 812                                         totalscore += partA__align_variousdist( whichmtx, scoringmatrices, NULL, tmpres1, tmpres2, eff1, eff2, eff1s, eff2s, clus1, clus2, alloclen, localhom, &impmatch, cut1[i], cut1[i+1]-1, cut2[i], cut2[i+1]-1, gapmap1, gapmap2, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres );
 813                                 }
 814                                 else
 815                                         totalscore += partA__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, localhom, &impmatch, cut1[i], cut1[i+1]-1, cut2[i], cut2[i+1]-1, gapmap1, gapmap2, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres );
 816                                 *totalimpmatch += impmatch;
 817 //                              fprintf( stderr, "*totalimpmatch in Falign_localhom = %f\n", *totalimpmatch );
 818
 819
 820                                 break;
 821                         default:
 822                                 fprintf( stderr, "alg = %c\n", alg );
 823                                 ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
 824                                 break;
 825                 }
 826 #ifdef enablemultithread
 827                 if( chudanres && *chudanres )
 828                 {
 829 //                      fprintf( stderr, "\n\n## CHUUDAN!!! at Falign_localhom\n" );
 830                         return( -1.0 );
 831                 }
 832 #endif
 833
 834                 nlen = strlen( tmpres1[0] );
 835                 if( totallen + nlen > alloclen )
 836                 {
 837                         fprintf( stderr, "totallen=%d +  nlen=%d > alloclen = %d\n", totallen, nlen, alloclen );
 838                         ErrorExit( "LENGTH OVER in Falign\n " );
 839                 }
 840                 for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
 841                 for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
 842                 totallen += nlen;
 843 #if 0
 844                 fprintf( stderr, "%4d\r", totallen );
 845                 fprintf( stderr, "\n\n" );
 846                 for( j=0; j<clus1; j++ )
 847                 {
 848                         fprintf( stderr, "%s\n", tmpres1[j] );
 849                 }
 850                 fprintf( stderr, "-------\n" );
 851                 for( j=0; j<clus2; j++ )
 852                 {
 853                         fprintf( stderr, "%s\n", tmpres2[j] );
 854                 }
 855 #endif
 856         }
 857 #if KEIKA
 858         fprintf( stderr, "DP ... done   \n" );
 859 #endif
 860
 861         for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
 862         for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
 863 #if 0
 864         for( j=0; j<clus1; j++ )
 865         {
 866                 fprintf( stderr, "%s\n", result1[j] );
 867         }
 868         fprintf( stderr, "- - - - - - - - - - -\n" );
 869         for( j=0; j<clus2; j++ )
 870         {
 871                 fprintf( stderr, "%s\n", result2[j] );
 872         }
 873 #endif
 874         return( totalscore );
 875 }