Logo Search packages:      
Sourcecode: mafft version File versions  Download package

Falign.c

#include "mltaln.h"

static FILE *fftfp;
static int n20or4or2;

#define KEIKA 0
#define RND   0
#define DEBUG 0

static void generateRndSeq( char *seq, int len )
{
      while( len-- )
#if 1
            *seq++ = (int)( rnd() * n20or4or2 );
#else
            *seq++ = (int)1;
#endif
}

static void vec_init( Fukusosuu *result, int nlen )
{
      while( nlen-- ) result->R = result++->I = 0.0;
}

static void vec_init2( Fukusosuu **result, char *seq, double eff, int st, int ed )
{
      int i;
      for( i=st; i<ed; i++ )
            result[*seq++][i].R += eff;
}

static void seq_vec_2( Fukusosuu *result, double *score, double incr, char *seq )
{
      static int n;
      for( ; *seq; result++ )
      {
            n = amino_n[(int)*seq++];
            if( n < 20 && n >= 0 ) result->R += incr * score[n];
#if 0
            fprintf( stderr, "n=%d, score=%f, inc=%f R=%f\n",n,  score[n], incr * score[n], result->R );
#endif
      }
}

static void seq_vec_3( Fukusosuu **result, double incr, char *seq )
{
      int i;
      int n;
      for( i=0; *seq; i++ )
      {
            n = amino_n[(int)*seq++];
            if( n < n20or4or2 && n >= 0 ) result[n][i].R += incr;
      }
}

static void seq_vec_5( Fukusosuu *result, double *score1, double *score2, double incr, char *seq )
{
      int n;
      for( ; *seq; result++ )
      {
            n = amino_n[(int)*seq++];
            if( n > 20 ) continue;
            result->R += incr * score1[n];
            result->I += incr * score2[n];
#if 0
            fprintf( stderr, "n=%d, score=%f, inc=%f R=%f\n",n,  score[n], incr * score[n], result->R );
#endif
      }
}


static void seq_vec_4( Fukusosuu *result, double incr, char *seq )
{
      char s;
      for( ; *seq; result++ )
      {
            s = *seq++;
            if( s == 'a' )
                  result->R += incr;
            else if( s == 't' )
                  result->R -= incr;
            else if( s == 'g' )
                  result->I += incr;
            else if( s == 'c' )
                  result->I -= incr;
      }
}

static void seq_vec( Fukusosuu *result, char query, double incr, char *seq )
{
#if 0
      int bk = nlen;
#endif
      while( *seq )
      {
            if( *seq++ == query ) result->R += incr;
            result++;
#if 0
fprintf( stderr, "i = %d result->R = %f\n", bk-nlen, (result-1)->R );
#endif
      }
}

static int checkRepeat( int num, int *cutpos )
{
      int tmp, buf;

      buf = *cutpos;
      while( num-- )
      {
            if( ( tmp = *cutpos++ ) < buf ) return( 1 );
            buf = tmp;
      }
      return( 0 );
}

static int segcmp( void *ptr1, void *ptr2 )
{
      int diff;
      Segment **seg1 = (Segment **)ptr1;
      Segment **seg2 = (Segment **)ptr2;
#if 0
      return( (*seg1)->center - (*seg2)->center );
#else
      diff = (*seg1)->center - (*seg2)->center;
      if( diff ) return( diff );

      diff = (*seg1)->start - (*seg2)->start;
      if( diff ) return( diff );

      diff = (*seg1)->end - (*seg2)->end;
      if( diff ) return( diff );

      fprintf( stderr, "USE STABLE SORT !!\n" );
      exit( 1 );
      return( 0 );
#endif
}


static void mymergesort( int first, int last, Segment **seg )
{
      int middle;
      static int i, j, k, p;
      static int allo = 0;
      static Segment **work = NULL;
      if( last > allo )
      {
            allo = last;
            if( work ) free( work );
            work = (Segment **)calloc( allo / 2 + 1, sizeof( Segment *) );
      }

      if( first < last )
      {
            middle = ( first + last ) / 2;
            mymergesort( first, middle, seg );
            mymergesort( middle+1, last, seg );
            p = 0;
            for( i=first; i<=middle; i++ ) work[p++] = seg[i];
            i = middle + 1; j = 0; k = first;
            while( i <= last && j < p )
            {
                  if( work[j]->center <= seg[i]->center ) 
                        seg[k++] = work[j++];
                  else
                        seg[k++] = seg[i++];
            }
            while( j < p ) seg[k++] = work[j++];
      }
}


double Fgetlag( char  **seq1, char  **seq2, 
                      double *eff1, double *eff2, 
                      int    clus1, int    clus2,
                      int alloclen )
{
      int i, j, k, l, m;
      int nlen, nlen2, nlen4;
      static int crossscoresize = 0;
      static char **tmpseq1 = NULL;
      static char **tmpseq2 = NULL;
      static char **tmpptr1 = NULL;
      static char **tmpptr2 = NULL;
      static char **tmpres1 = NULL;
      static char **tmpres2 = NULL;
      static char **result1 = NULL;
      static char **result2 = NULL;
#if RND
      static char **rndseq1 = NULL;
      static char **rndseq2 = NULL;
#endif
      static Fukusosuu **seqVector1 = NULL;
      static Fukusosuu **seqVector2 = NULL;
      static Fukusosuu **naiseki = NULL;   
      static Fukusosuu *naisekiNoWa = NULL; 
      static double *soukan = NULL;
      static double **crossscore = NULL;
      int nlentmp;
      static int *kouho = NULL;
      static Segment *segment = NULL;
      static Segment *segment1 = NULL;
      static Segment *segment2 = NULL;
      static Segment **sortedseg1 = NULL;
      static Segment **sortedseg2 = NULL;
      static int *cut1 = NULL;
      static int *cut2 = NULL;
      static int localalloclen = 0;
      int lag;
      int tmpint;
      int count, count0;
      int len1, len2;
      int totallen;
      double totaleff1, totaleff2;
      float dumfl = 0.0;

      totaleff1 = 0.0; for( i=0; i<clus1; i++ ) totaleff1 += eff1[i];
      totaleff2 = 0.0; for( i=0; i<clus2; i++ ) totaleff2 += eff2[i];

      len1 = strlen( seq1[0] );
      len2 = strlen( seq2[0] );
      nlentmp = MAX( len1, len2 );

      nlen = 1;
      while( nlentmp >= nlen ) nlen <<= 1;
#if 0
      fprintf( stderr, "###   nlen    = %d\n", nlen );
#endif

      nlen2 = nlen/2; nlen4 = nlen2 / 2;

#if DEBUG
      fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
      fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
#endif

      if( !localalloclen )
      {
            kouho = AllocateIntVec( NKOUHO );
            cut1 = AllocateIntVec( MAXSEG );
            cut2 = AllocateIntVec( MAXSEG );
            tmpptr1 = AllocateCharMtx( njob, 0 );
            tmpptr2 = AllocateCharMtx( njob, 0 );
            result1 = AllocateCharMtx( njob, alloclen );
            result2 = AllocateCharMtx( njob, alloclen );
            tmpres1 = AllocateCharMtx( njob, alloclen );
            tmpres2 = AllocateCharMtx( njob, alloclen );
//          crossscore = AllocateDoubleMtx( MAXSEG, MAXSEG );
            segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
            segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
            segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
            sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
            sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
            if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
                  ErrorExit( "Allocation error\n" );

            if     ( scoremtx == -1 ) n20or4or2 = 4;
            else if( fftscore == 1  ) n20or4or2 = 2;
            else                      n20or4or2 = 20;
      }
      if( localalloclen < nlen )
      {
            if( localalloclen )
            {
#if 1
                  FreeFukusosuuMtx ( seqVector1 );
                  FreeFukusosuuMtx ( seqVector2 );
                  FreeFukusosuuVec( naisekiNoWa );
                  FreeFukusosuuMtx( naiseki );
                  FreeDoubleVec( soukan );
                  FreeCharMtx( tmpseq1 );
                  FreeCharMtx( tmpseq2 );
#endif
#if RND
                  FreeCharMtx( rndseq1 );
                  FreeCharMtx( rndseq2 );
#endif
            }


            tmpseq1 = AllocateCharMtx( njob, nlen );
            tmpseq2 = AllocateCharMtx( njob, nlen );
            naisekiNoWa = AllocateFukusosuuVec( nlen );
            naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
            seqVector1 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
            seqVector2 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
            soukan = AllocateDoubleVec( nlen+1 );

#if RND
            rndseq1 = AllocateCharMtx( njob, nlen );
            rndseq2 = AllocateCharMtx( njob, nlen );
            for( i=0; i<njob; i++ )
            {
                  generateRndSeq( rndseq1[i], nlen );
                  generateRndSeq( rndseq2[i], nlen );
            }
#endif
            localalloclen = nlen;
      }
      
      for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
      for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );

#if 0
fftfp = fopen( "input_of_Falign", "w" );
fprintf( fftfp, "nlen = %d\n", nlen );
fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
for( i=0; i<clus1; i++ )
      fprintf( fftfp, "%s\n", seq1[i] );
fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
for( i=0; i<clus2; i++ )
      fprintf( fftfp, "%s\n", seq2[i] );
fclose( fftfp );
system( "less input_of_Falign < /dev/tty > /dev/tty" );
#endif

      if( fftkeika ) fprintf( stderr,  " FFT ... " );

      for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
      if( fftscore && scoremtx != -1 )
      {
            for( i=0; i<clus1; i++ )
            {
                  seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
                  seq_vec_2( seqVector1[1], volume,   eff1[i], tmpseq1[i] );
            }
      }
      else
      {
#if 0
            for( i=0; i<clus1; i++ ) for( j=0; j<n20or4or2; j++ ) 
                  seq_vec( seqVector1[j], amino[j], eff1[i], tmpseq1[i] );
#else
            for( i=0; i<clus1; i++ )
                  seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
#endif
      }
#if RND
      for( i=0; i<clus1; i++ )
      {
            vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
      }
#endif
#if 0
fftfp = fopen( "seqVec", "w" );
fprintf( fftfp, "before transform\n" );
for( k=0; k<n20or4or2; k++ ) 
{
   fprintf( fftfp, "nlen=%d\n", nlen );
   fprintf( fftfp, "%c\n", amino[k] );
   for( l=0; l<nlen; l++ )
   fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
}
fclose( fftfp );
system( "less seqVec < /dev/tty > /dev/tty" );
#endif

      for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
      if( fftscore && scoremtx != -1 )
      {
            for( i=0; i<clus2; i++ )
            {
                  seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
                  seq_vec_2( seqVector2[1], volume,   eff2[i], tmpseq2[i] );
            }
      }
      else
      {
#if 0
            for( i=0; i<clus2; i++ ) for( j=0; j<n20or4or2; j++ ) 
                  seq_vec( seqVector2[j], amino[j], eff2[i], tmpseq2[i] );
#else
            for( i=0; i<clus2; i++ )
                  seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
#endif
      }
#if RND
      for( i=0; i<clus2; i++ )
      {
            vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
      }
#endif

#if 0
fftfp = fopen( "seqVec2", "w" );
fprintf( fftfp, "before fft\n" );
for( k=0; k<n20or4or2; k++ ) 
{
   fprintf( fftfp, "%c\n", amino[k] );
   for( l=0; l<nlen; l++ )
   fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
}
fclose( fftfp );
system( "less seqVec2 < /dev/tty > /dev/tty" );
#endif

      for( j=0; j<n20or4or2; j++ )
      {
            fft( nlen, seqVector2[j], (j==0) );
            fft( nlen, seqVector1[j], 0 );
      }
#if 0
fftfp = fopen( "seqVec2", "w" );
fprintf( fftfp, "#after fft\n" );
for( k=0; k<n20or4or2; k++ ) 
{
   fprintf( fftfp, "#%c\n", amino[k] );
   for( l=0; l<nlen; l++ )
   fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
}
fclose( fftfp );
system( "less seqVec2 < /dev/tty > /dev/tty" );
#endif

      for( k=0; k<n20or4or2; k++ ) 
      {
            for( l=0; l<nlen; l++ ) 
                  calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
      }
      for( l=0; l<nlen; l++ ) 
      {
            naisekiNoWa[l].R = 0.0;
            naisekiNoWa[l].I = 0.0;
            for( k=0; k<n20or4or2; k++ ) 
            {
                  naisekiNoWa[l].R += naiseki[k][l].R;
                  naisekiNoWa[l].I += naiseki[k][l].I;
            }
      }

#if 0
fftfp = fopen( "naisekiNoWa", "w" );
fprintf( fftfp, "#Before fft\n" );
for( l=0; l<nlen; l++ )
      fprintf( fftfp, "%d  %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I ); 
fclose( fftfp );
system( "less naisekiNoWa < /dev/tty > /dev/tty " );
#endif

      fft( -nlen, naisekiNoWa, 0 );

      for( m=0; m<=nlen2; m++ ) 
            soukan[m] = naisekiNoWa[nlen2-m].R;
      for( m=nlen2+1; m<nlen; m++ ) 
            soukan[m] = naisekiNoWa[nlen+nlen2-m].R;

#if 0
fftfp = fopen( "naisekiNoWa", "w" );
fprintf( fftfp, "#After fft\n" );
for( l=0; l<nlen; l++ )
      fprintf( fftfp, "%d  %f\n", l, naisekiNoWa[l].R ); 
fclose( fftfp );
fftfp = fopen( "list.plot", "w"  );
fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
fclose( fftfp );
system( "/usr/bin/gnuplot list.plot &" );
#endif
#if 0
fprintf( stderr, "frt write start\n" );
fftfp = fopen( "frt", "w" );
for( l=0; l<nlen; l++ )
      fprintf( fftfp, "%d  %f\n", l-nlen2, soukan[l] ); 
fclose( fftfp );
system( "less frt < /dev/tty > /dev/tty" );
#if 0
fftfp = fopen( "list.plot", "w"  );
fprintf( fftfp, "plot 'frt'\n pause +1" );
fclose( fftfp );
system( "/usr/bin/gnuplot list.plot" );
#endif
#endif


      getKouho( kouho, NKOUHO, soukan, nlen );

#if 0
      for( i=0; i<NKOUHO; i++ )
      {
            fprintf( stdout, "kouho[%d] = %d\n", i, kouho[i] );
      }
#endif

#if KEIKA
      fprintf( stderr, "Searching anchors ... " );
#endif
      count = 0;



#define CAND 0
#if CAND
      fftfp = fopen( "cand", "w" );
      fclose( fftfp );
#endif

      for( k=0; k<NKOUHO; k++ ) 
      {

            lag = kouho[k];
            zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
#if CAND
            fftfp = fopen( "cand", "a" );
            fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
            fprintf( fftfp, "%s\n", tmpptr1[0] );
            fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
            fprintf( fftfp, "%s\n", tmpptr2[0] );
            fprintf( fftfp, ">\n", k+1, lag );
            fclose( fftfp );
#endif
            tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
            
            if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );


            if( tmpint == 0 ) break; // 060430 iinoka ?
            while( tmpint-- > 0 )
            {
                  if( lag > 0 )
                  {
                        segment1[count].start  = segment[count].start ;
                        segment1[count].end    = segment[count].end   ;
                        segment1[count].center = segment[count].center;
                        segment1[count].score  = segment[count].score;

                        segment2[count].start  = segment[count].start  + lag;
                        segment2[count].end    = segment[count].end    + lag;
                        segment2[count].center = segment[count].center + lag;
                        segment2[count].score  = segment[count].score       ;
                  }
                  else
                  {
                        segment1[count].start  = segment[count].start  - lag;
                        segment1[count].end    = segment[count].end    - lag;
                        segment1[count].center = segment[count].center - lag;
                        segment1[count].score  = segment[count].score       ;

                        segment2[count].start  = segment[count].start ;
                        segment2[count].end    = segment[count].end   ;
                        segment2[count].center = segment[count].center;
                        segment2[count].score  = segment[count].score ;
                  }
#if 0
                  fprintf( stderr, "Goukaku=%dko\n", tmpint ); 
                  fprintf( stderr, "in 1 %d\n", segment1[count].center );
                  fprintf( stderr, "in 2 %d\n", segment2[count].center );
#endif
                  segment1[count].pair = &segment2[count];
                  segment2[count].pair = &segment1[count];
                  count++;
#if 0
                  fprintf( stderr, "count=%d\n", count );
#endif
            }
      }

#if 1
      fprintf( stderr, "done. (%d anchors)\r", count );
#endif
      if( !count && fftNoAnchStop )
            ErrorExit( "Cannot detect anchor!" );
#if 0
      fprintf( stdout, "RESULT before sort:\n" );
      for( l=0; l<count+1; l++ )
      {
            fprintf( stdout, "cut[%d]=%d, ", l, segment1[l].center );
            fprintf( stdout, "%d score = %f\n", segment2[l].center, segment1[l].score );
      }
      exit( 1 );
#endif

#if KEIKA
      fprintf( stderr, "Aligning anchors ... " );
#endif
      for( i=0; i<count; i++ )
      {
            sortedseg1[i] = &segment1[i];
            sortedseg2[i] = &segment2[i];
      }

      {
            mymergesort( 0, count-1, sortedseg1 ); 
            mymergesort( 0, count-1, sortedseg2 ); 
            for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
            for( i=0; i<count; i++ ) sortedseg2[i]->number = i;

            if( crossscoresize < count+2 )
            {
                  crossscoresize = count+2;
                  fprintf( stderr, "####################################################################################################################################allocating crossscore, size = %d\n", crossscoresize );
                  if( crossscore ) FreeDoubleMtx( crossscore );
                  crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
            }

            for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
                  crossscore[i][j] = 0.0;
            for( i=0; i<count; i++ )
            {
                  crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
                  cut1[i+1] = sortedseg1[i]->center;
                  cut2[i+1] = sortedseg2[i]->center;
            }

#if DEBUG
            fprintf( stderr, "AFTER SORT\n" );
            for( i=0; i<count; i++ ) fprintf( stderr, "%d, %d\n", segment1[i].start, segment2[i].start );
#endif

            crossscore[0][0] = 10000000.0;
            cut1[0] = 0; 
            cut2[0] = 0;
            crossscore[count+1][count+1] = 10000000.0;
            cut1[count+1] = len1;
            cut2[count+1] = len2;
            count += 2;
            count0 = count;

            blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
      }
      if( fftkeika )
      {
            if( count0 > count )
            {
                  fprintf( stderr, "REPEAT!? \n" ); 
                  if( fftRepeatStop ) exit( 1 );
            }
#if KEIKA
            else 
                  fprintf( stderr, "done\n" );
                  fprintf( stderr, "done. (%d anchors)\n", count );
#endif
      }

#if 0
      fftfp = fopen( "fft", "a" );
      fprintf( fftfp, "RESULT after sort:\n" );
      for( l=0; l<count; l++ )
      {
            fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
            fprintf( fftfp, "%d\n", segment2[l].center );
      }
      fclose( fftfp );
#endif

#if 0
      fftfp = fopen( "fft", "a" );
      fprintf( fftfp, "RESULT after sort:\n" );
      for( l=0; l<count; l++ )
      {
            fprintf( fftfp, "cut : %d %d\n", cut1[l], cut2[l] );
      }
      fclose( fftfp );
#endif

#if KEIKA
      fprintf( trap_g, "Devided to %d segments\n", count-1 );
      fprintf( trap_g, "%d  %d forg\n", MIN( clus1, clus2 ), count-1 );
#endif

      totallen = 0;
      for( j=0; j<clus1; j++ ) result1[j][0] = 0;
      for( j=0; j<clus2; j++ ) result2[j][0] = 0;
      for( i=0; i<count-1; i++ )
      {
#if DEBUG
            fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
#else
#if KEIKA
            fprintf( stderr, "DP %03d / %03d\r", i+1, count-1 );
#endif
#endif
            for( j=0; j<clus1; j++ )
            {
                  strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
                  tmpres1[j][cut1[i+1]-cut1[i]] = 0;
            }
            for( j=0; j<clus2; j++ )
            {
                  strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
                  tmpres2[j][cut2[i+1]-cut2[i]] = 0;
            }
            switch( alg )
            {
                  case( 'a' ):
                        Aalign( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen );
                        break;
                  case( 'M' ):
                              MSalignmm( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL );
                        break;
                  case( 'A' ):
                        if( clus1 == 1 && clus2 == 1 )
                              G__align11( tmpres1, tmpres2, alloclen );
                        else
                              A__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
                        break;
                  case( 'H' ):
                        if( clus1 == 1 && clus2 == 1 )
                              G__align11( tmpres1, tmpres2, alloclen );
                        else
                              H__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
                        break;
                  case( 'Q' ):
                        if( clus1 == 1 && clus2 == 1 )
                              G__align11( tmpres1, tmpres2, alloclen );
                        else
                              Q__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
                        break;
                  default:
                        fprintf( stderr, "alg = %c\n", alg );
                        ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
                        break;
            }

            nlen = strlen( tmpres1[0] );
            if( totallen + nlen > alloclen ) ErrorExit( "LENGTH OVER in Falign\n " );
            for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
            for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
            totallen += nlen;
#if 0
            fprintf( stderr, "%4d\r", totallen );
            fprintf( stderr, "\n\n" );
            for( j=0; j<clus1; j++ ) 
            {
                  fprintf( stderr, "%s\n", tmpres1[j] );
            }
            fprintf( stderr, "-------\n" );
            for( j=0; j<clus2; j++ ) 
            {
                  fprintf( stderr, "%s\n", tmpres2[j] );
            }
#endif
      }
#if KEIKA
      fprintf( stderr, "DP ... done   \n" );
#endif

      for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
      for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
#if 0
      for( j=0; j<clus1; j++ ) 
      {
            fprintf( stderr, "%s\n", result1[j] );
      }
      fprintf( stderr, "- - - - - - - - - - -\n" );
      for( j=0; j<clus2; j++ ) 
      {
            fprintf( stderr, "%s\n", result2[j] );
      }
#endif
      return( 0.0 );
}
float Falign( char  **seq1, char  **seq2, 
                    double *eff1, double *eff2, 
                    int    clus1, int    clus2,
                    int alloclen, int *fftlog )
{
      int i, j, k, l, m, maxk;
      int nlen, nlen2, nlen4;
      static int prevalloclen = 0;
      static int crossscoresize = 0;
      static char **tmpseq1 = NULL;
      static char **tmpseq2 = NULL;
      static char **tmpptr1 = NULL;
      static char **tmpptr2 = NULL;
      static char **tmpres1 = NULL;
      static char **tmpres2 = NULL;
      static char **result1 = NULL;
      static char **result2 = NULL;
#if RND
      static char **rndseq1 = NULL;
      static char **rndseq2 = NULL;
#endif
      static Fukusosuu **seqVector1 = NULL;
      static Fukusosuu **seqVector2 = NULL;
      static Fukusosuu **naiseki = NULL;   
      static Fukusosuu *naisekiNoWa = NULL; 
      static double *soukan = NULL;
      static double **crossscore = NULL;
      int nlentmp;
      static int *kouho = NULL;
      static Segment *segment = NULL;
      static Segment *segment1 = NULL;
      static Segment *segment2 = NULL;
      static Segment **sortedseg1 = NULL;
      static Segment **sortedseg2 = NULL;
      static int *cut1 = NULL;
      static int *cut2 = NULL;
      static char *sgap1, *egap1, *sgap2, *egap2;
      static int localalloclen = 0;
      int lag;
      int tmpint;
      int count, count0;
      int len1, len2;
      int totallen;
      double totaleff1, totaleff2;
      float totalscore;
      float dumfl = 0.0;


      totaleff1 = 0.0; for( i=0; i<clus1; i++ ) totaleff1 += eff1[i];
      totaleff2 = 0.0; for( i=0; i<clus2; i++ ) totaleff2 += eff2[i];

      len1 = strlen( seq1[0] );
      len2 = strlen( seq2[0] );
      nlentmp = MAX( len1, len2 );

      nlen = 1;
      while( nlentmp >= nlen ) nlen <<= 1;
#if 0
      fprintf( stderr, "###   nlen    = %d\n", nlen );
#endif

      nlen2 = nlen/2; nlen4 = nlen2 / 2;

#if DEBUG
      fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
      fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
#endif

      if( prevalloclen != alloclen ) // Falign_noudp mo kaeru
      {
            if( prevalloclen )
            {
                  FreeCharMtx( result1 );
                  FreeCharMtx( result2 );
                  FreeCharMtx( tmpres1 );
                  FreeCharMtx( tmpres2 );
            }
//          fprintf( stderr, "\n\n\nreallocating ...\n" ); 
            result1 = AllocateCharMtx( njob, alloclen );
            result2 = AllocateCharMtx( njob, alloclen );
            tmpres1 = AllocateCharMtx( njob, alloclen );
            tmpres2 = AllocateCharMtx( njob, alloclen );
            prevalloclen = alloclen;
      }
      if( !localalloclen )
      {
            sgap1 = AllocateCharVec( njob );
            egap1 = AllocateCharVec( njob );
            sgap2 = AllocateCharVec( njob );
            egap2 = AllocateCharVec( njob );
            kouho = AllocateIntVec( NKOUHO );
            cut1 = AllocateIntVec( MAXSEG );
            cut2 = AllocateIntVec( MAXSEG );
            tmpptr1 = AllocateCharMtx( njob, 0 );
            tmpptr2 = AllocateCharMtx( njob, 0 );
//          crossscore = AllocateDoubleMtx( MAXSEG, MAXSEG );
            segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
            segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
            segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
            sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
            sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
            if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
                  ErrorExit( "Allocation error\n" );

            if     ( scoremtx == -1 ) n20or4or2 = 1;
            else if( fftscore )       n20or4or2 = 1;
            else                      n20or4or2 = 20;
      }
      if( localalloclen < nlen )
      {
            if( localalloclen )
            {
#if 1
                  if( !kobetsubunkatsu )
                  {
                        FreeFukusosuuMtx ( seqVector1 );
                        FreeFukusosuuMtx ( seqVector2 );
                        FreeFukusosuuVec( naisekiNoWa );
                        FreeFukusosuuMtx( naiseki );
                        FreeDoubleVec( soukan );
                  }
                  FreeCharMtx( tmpseq1 );
                  FreeCharMtx( tmpseq2 );
#endif
#if RND
                  FreeCharMtx( rndseq1 );
                  FreeCharMtx( rndseq2 );
#endif
            }

            tmpseq1 = AllocateCharMtx( njob, nlen );
            tmpseq2 = AllocateCharMtx( njob, nlen );
            if( !kobetsubunkatsu )
            {
                  naisekiNoWa = AllocateFukusosuuVec( nlen );
                  naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
                  seqVector1 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
                  seqVector2 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
                  soukan = AllocateDoubleVec( nlen+1 );
            }
#if RND
            rndseq1 = AllocateCharMtx( njob, nlen );
            rndseq2 = AllocateCharMtx( njob, nlen );
            for( i=0; i<njob; i++ )
            {
                  generateRndSeq( rndseq1[i], nlen );
                  generateRndSeq( rndseq2[i], nlen );
            }
#endif
            localalloclen = nlen;
      }
      
      for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
      for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );

#if 0
fftfp = fopen( "input_of_Falign", "w" );
fprintf( fftfp, "nlen = %d\n", nlen );
fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
for( i=0; i<clus1; i++ )
      fprintf( fftfp, "%s\n", seq1[i] );
fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
for( i=0; i<clus2; i++ )
      fprintf( fftfp, "%s\n", seq2[i] );
fclose( fftfp );
system( "less input_of_Falign < /dev/tty > /dev/tty" );
#endif
      if( !kobetsubunkatsu )
      {
            if( fftkeika ) fprintf( stderr,  " FFT ... " );

            for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
            if( fftscore && scoremtx != -1 )
            {
                  for( i=0; i<clus1; i++ )
                  {
#if 1
                        seq_vec_5( seqVector1[0], polarity, volume, eff1[i], tmpseq1[i] );
#else
                        seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
                        seq_vec_2( seqVector1[1], volume,   eff1[i], tmpseq1[i] );
#endif
                  }
            }
            else
            {
#if 0
                  for( i=0; i<clus1; i++ ) for( j=0; j<n20or4or2; j++ ) 
                        seq_vec( seqVector1[j], amino[j], eff1[i], tmpseq1[i] );
#else
                  for( i=0; i<clus1; i++ )
                        seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
#endif
            }
#if RND
            for( i=0; i<clus1; i++ )
            {
                  vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
            }
#endif
#if 0
fftfp = fopen( "seqVec", "w" );
fprintf( fftfp, "before transform\n" );
for( k=0; k<n20or4or2; k++ ) 
{
   fprintf( fftfp, "nlen=%d\n", nlen );
   fprintf( fftfp, "%c\n", amino[k] );
   for( l=0; l<nlen; l++ )
   fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
}
fclose( fftfp );
system( "less seqVec < /dev/tty > /dev/tty" );
#endif

            for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
            if( fftscore && scoremtx != -1 )
            {
                  for( i=0; i<clus2; i++ )
                  {
#if 1
                        seq_vec_5( seqVector2[0], polarity, volume, eff2[i], tmpseq2[i] );
#else
                        seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
                        seq_vec_2( seqVector2[1], volume,   eff2[i], tmpseq2[i] );
#endif
                  }
            }
            else
            {
#if 0
                  for( i=0; i<clus2; i++ ) for( j=0; j<n20or4or2; j++ ) 
                        seq_vec( seqVector2[j], amino[j], eff2[i], tmpseq2[i] );
#else
                  for( i=0; i<clus2; i++ )
                        seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
#endif
            }
#if RND
            for( i=0; i<clus2; i++ )
            {
                  vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
            }
#endif

#if 0
fftfp = fopen( "seqVec2", "w" );
fprintf( fftfp, "before fft\n" );
for( k=0; k<n20or4or2; k++ ) 
{
   fprintf( fftfp, "%c\n", amino[k] );
   for( l=0; l<nlen; l++ )
   fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
}
fclose( fftfp );
system( "less seqVec2 < /dev/tty > /dev/tty" );
#endif

            for( j=0; j<n20or4or2; j++ )
            {
                  fft( nlen, seqVector2[j], (j==0) );
                  fft( nlen, seqVector1[j], 0 );
            }
#if 0
fftfp = fopen( "seqVec2", "w" );
fprintf( fftfp, "#after fft\n" );
for( k=0; k<n20or4or2; k++ ) 
{
   fprintf( fftfp, "#%c\n", amino[k] );
   for( l=0; l<nlen; l++ )
         fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
}
fclose( fftfp );
system( "less seqVec2 < /dev/tty > /dev/tty" );
#endif

            for( k=0; k<n20or4or2; k++ ) 
            {
                  for( l=0; l<nlen; l++ ) 
                        calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
            }
            for( l=0; l<nlen; l++ ) 
            {
                  naisekiNoWa[l].R = 0.0;
                  naisekiNoWa[l].I = 0.0;
                  for( k=0; k<n20or4or2; k++ ) 
                  {
                        naisekiNoWa[l].R += naiseki[k][l].R;
                        naisekiNoWa[l].I += naiseki[k][l].I;
                  }
            }
      
#if 0
      fftfp = fopen( "naisekiNoWa", "w" );
      fprintf( fftfp, "#Before fft\n" );
      for( l=0; l<nlen; l++ )
            fprintf( fftfp, "%d  %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I ); 
      fclose( fftfp );
      system( "less naisekiNoWa < /dev/tty > /dev/tty " );
#endif

            fft( -nlen, naisekiNoWa, 0 );
      
            for( m=0; m<=nlen2; m++ ) 
                  soukan[m] = naisekiNoWa[nlen2-m].R;
            for( m=nlen2+1; m<nlen; m++ ) 
                  soukan[m] = naisekiNoWa[nlen+nlen2-m].R;

#if 0
      fftfp = fopen( "naisekiNoWa", "w" );
      fprintf( fftfp, "#After fft\n" );
      for( l=0; l<nlen; l++ )
            fprintf( fftfp, "%d  %f\n", l, naisekiNoWa[l].R ); 
      fclose( fftfp );
      fftfp = fopen( "list.plot", "w"  );
      fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
      fclose( fftfp );
      system( "/usr/bin/gnuplot list.plot &" );
#endif
#if 0
      fprintf( stderr, "soukan\n" );
      for( l=0; l<nlen; l++ )
            fprintf( stderr, "%d  %f\n", l-nlen2, soukan[l] ); 
#if 0
      fftfp = fopen( "list.plot", "w"  );
      fprintf( fftfp, "plot 'frt'\n pause +1" );
      fclose( fftfp );
      system( "/usr/bin/gnuplot list.plot" );
#endif
#endif


            getKouho( kouho, NKOUHO, soukan, nlen );

#if 0
            for( i=0; i<NKOUHO; i++ )
            {
                  fprintf( stderr, "kouho[%d] = %d\n", i, kouho[i] );
            }
#endif
      }

#if KEIKA
      fprintf( stderr, "Searching anchors ... " );
#endif
      count = 0;



#define CAND 0
#if CAND
      fftfp = fopen( "cand", "w" );
      fclose( fftfp );
#endif
      if( kobetsubunkatsu )
      {
            maxk = 1;
            kouho[0] = 0;
      }
      else
      {
            maxk = NKOUHO;
      }

      for( k=0; k<maxk; k++ ) 
      {
            lag = kouho[k];
            if( lag <= -len1 || len2 <= lag ) continue;
            zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
#if CAND
            fftfp = fopen( "cand", "a" );
            fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
            fprintf( fftfp, "%s\n", tmpptr1[0] );
            fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
            fprintf( fftfp, "%s\n", tmpptr2[0] );
            fprintf( fftfp, ">\n", k+1, lag );
            fclose( fftfp );
#endif

//          fprintf( stderr, "lag = %d\n", lag );
            tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );

//          if( lag == -50 ) exit( 1 );
            
            if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );


            if( tmpint == 0 ) break; // 060430 iinoka ?
            while( tmpint-- > 0 )
            {
#if 0
                  if( segment[count].end - segment[count].start < fftWinSize )
                  {
                        count++;
                        continue;
                  }
#endif
                  if( lag > 0 )
                  {
                        segment1[count].start  = segment[count].start ;
                        segment1[count].end    = segment[count].end   ;
                        segment1[count].center = segment[count].center;
                        segment1[count].score  = segment[count].score;

                        segment2[count].start  = segment[count].start  + lag;
                        segment2[count].end    = segment[count].end    + lag;
                        segment2[count].center = segment[count].center + lag;
                        segment2[count].score  = segment[count].score       ;
                  }
                  else
                  {
                        segment1[count].start  = segment[count].start  - lag;
                        segment1[count].end    = segment[count].end    - lag;
                        segment1[count].center = segment[count].center - lag;
                        segment1[count].score  = segment[count].score       ;

                        segment2[count].start  = segment[count].start ;
                        segment2[count].end    = segment[count].end   ;
                        segment2[count].center = segment[count].center;
                        segment2[count].score  = segment[count].score ;
                  }
#if 0
                  fprintf( stderr, "in 1 %d\n", segment1[count].center );
                  fprintf( stderr, "in 2 %d\n", segment2[count].center );
#endif
                  segment1[count].pair = &segment2[count];
                  segment2[count].pair = &segment1[count];
                  count++;
            }
      }
#if 0
      if( !kobetsubunkatsu && fftkeika )
            fprintf( stderr, "%d anchors found\r", count );
#endif
      if( !count && fftNoAnchStop )
            ErrorExit( "Cannot detect anchor!" );
#if 0
      fprintf( stderr, "RESULT before sort:\n" );
      for( l=0; l<count+1; l++ )
      {
            fprintf( stderr, "cut[%d]=%d, ", l, segment1[l].center );
            fprintf( stderr, "%d score = %f\n", segment2[l].center, segment1[l].score );
      }
#endif

#if KEIKA
      fprintf( stderr, "done. (%d anchors)\n", count );
      fprintf( stderr, "Aligning anchors ... " );
#endif
      for( i=0; i<count; i++ )
      {
            sortedseg1[i] = &segment1[i];
            sortedseg2[i] = &segment2[i];
      }
#if 0
      tmpsort( count, sortedseg1 ); 
      tmpsort( count, sortedseg2 ); 
      qsort( sortedseg1, count, sizeof( Segment * ), segcmp );
      qsort( sortedseg2, count, sizeof( Segment * ), segcmp );
#else
      mymergesort( 0, count-1, sortedseg1 ); 
      mymergesort( 0, count-1, sortedseg2 ); 
#endif
      for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
      for( i=0; i<count; i++ ) sortedseg2[i]->number = i;


      if( kobetsubunkatsu )
      {
            for( i=0; i<count; i++ )
          {
                  cut1[i+1] = sortedseg1[i]->center;
                  cut2[i+1] = sortedseg2[i]->center;
            }
            cut1[0] = 0;
            cut2[0] = 0;
            cut1[count+1] = len1;
            cut2[count+1] = len2;
            count += 2;
      }
      else
      {
            if( crossscoresize < count+2 )
            {
                  crossscoresize = count+2;
#if 1
                  if( fftkeika ) fprintf( stderr, "######allocating crossscore, size = %d\n", crossscoresize );
#endif
                  if( crossscore ) FreeDoubleMtx( crossscore );
                  crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
            }
            for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
                  crossscore[i][j] = 0.0;
            for( i=0; i<count; i++ )
            {
                  crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
                  cut1[i+1] = sortedseg1[i]->center;
                  cut2[i+1] = sortedseg2[i]->center;
            }

#if 0
            fprintf( stderr, "AFTER SORT\n" );
            for( i=0; i<count+1; i++ ) fprintf( stderr, "%d, %d\n", cut1[i], cut2[i] );
            fprintf( stderr, "crossscore = \n" );
            for( i=0; i<count+1; i++ )
            {
                  for( j=0; j<count+1; j++ )
                        fprintf( stderr, "%.0f ", crossscore[i][j] );
                  fprintf( stderr, "\n" );
            }
#endif

            crossscore[0][0] = 10000000.0;
            cut1[0] = 0; 
            cut2[0] = 0;
            crossscore[count+1][count+1] = 10000000.0;
            cut1[count+1] = len1;
            cut2[count+1] = len2;
            count += 2;
            count0 = count;
      
            blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );

//          if( count-count0 )
//                fprintf( stderr, "%d unused anchors\n", count0-count );

            if( !kobetsubunkatsu && fftkeika )
                  fprintf( stderr, "%d anchors found\n", count );
            if( fftkeika )
            {
                  if( count0 > count )
                  {
#if 0
                        fprintf( stderr, "\7 REPEAT!? \n" ); 
#else
                        fprintf( stderr, "REPEAT!? \n" ); 
#endif
                        if( fftRepeatStop ) exit( 1 );
                  }
#if KEIKA
                  else fprintf( stderr, "done\n" );
#endif
            }
      }

#if 0
      fftfp = fopen( "fft", "a" );
      fprintf( fftfp, "RESULT after sort:\n" );
      for( l=0; l<count; l++ )
      {
            fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
            fprintf( fftfp, "%d\n", segment2[l].center );
      }
      fclose( fftfp );
#endif

#if 0
      fprintf( stderr, "RESULT after blckalign:\n" );
      for( l=0; l<count+1; l++ )
      {
            fprintf( stderr, "cut : %d %d\n", cut1[l], cut2[l] );
      }
#endif

#if 0
      fprintf( trap_g, "Devided to %d segments\n", count-1 );
      fprintf( trap_g, "%d  %d forg\n", MIN( clus1, clus2 ), count-1 );
#endif

      totallen = 0;
      for( j=0; j<clus1; j++ ) result1[j][0] = 0;
      for( j=0; j<clus2; j++ ) result2[j][0] = 0;
      totalscore = 0.0;
      *fftlog = -1;
      for( i=0; i<count-1; i++ )
      {
            *fftlog += 1;

            if( cut1[i] ) // chuui
            {
//                getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
//                getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
                  getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
                  getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
            }
            else
            {
                  for( j=0; j<clus1; j++ ) sgap1[j] = 'o';
                  for( j=0; j<clus2; j++ ) sgap2[j] = 'o';
            }
            if( cut1[i+1] != len1 ) // chuui
            {       
                  getkyokaigap( egap1, seq1, cut1[i+1], clus1 );
                  getkyokaigap( egap2, seq2, cut2[i+1], clus2 );
            }       
            else    
            {       
                  for( j=0; j<clus1; j++ ) egap1[j] = 'o';
                  for( j=0; j<clus2; j++ ) egap2[j] = 'o';
            }
#if 0
            {
                  fprintf( stderr, "kyokkaigap1(%d)=", cut1[i]-1 );
                  for( j=0; j<clus1; j++ )
                        fprintf( stderr, "%c", sgap1[j] );
                  fprintf( stderr, "=kyokkaigap1-start\n" );
            }
            {
                  fprintf( stderr, "kyokkaigap2(%d)=", cut2[i]-1 );
                  for( j=0; j<clus2; j++ )
                        fprintf( stderr, "%c", sgap2[j] );
                  fprintf( stderr, "=kyokkaigap2-start\n" );
            }
            {
                  fprintf( stderr, "kyokkaigap1(%d)=", cut1[i]-1 );
                  for( j=0; j<clus1; j++ )
                        fprintf( stderr, "%c", egap1[j] );
                  fprintf( stderr, "=kyokkaigap1-end\n" );
            }
            {
                  fprintf( stderr, "kyokkaigap2(%d)=", cut2[i]-1 );
                  for( j=0; j<clus2; j++ )
                        fprintf( stderr, "%c", egap2[j] );
                  fprintf( stderr, "=kyokkaigap2-end\n" );
            }
#endif

#if DEBUG
            fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
#else
#if KEIKA
            fprintf( stderr, "DP %03d / %03d\r", i+1, count-1 );
#endif
#endif
            for( j=0; j<clus1; j++ )
            {
                  strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
                  tmpres1[j][cut1[i+1]-cut1[i]] = 0;
            }
            if( kobetsubunkatsu && fftkeika ) commongappick( clus1, tmpres1 ); //dvtditr に呼ばれたとき fftkeika=1
//          if( kobetsubunkatsu ) commongappick( clus1, tmpres1 );
            for( j=0; j<clus2; j++ )
            {
                  strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
                  tmpres2[j][cut2[i+1]-cut2[i]] = 0;
            }
            if( kobetsubunkatsu && fftkeika ) commongappick( clus2, tmpres2 ); //dvtditr に呼ばれたとき fftkeika=1
//          if( kobetsubunkatsu ) commongappick( clus2, tmpres2 );

            if( constraint )
            {
                  fprintf( stderr, "Not supported\n" );
                  exit( 1 );
            }
#if 0
            fprintf( stderr, "i=%d, before alignment", i );
            fprintf( stderr, "%4d\n", totallen );
            fprintf( stderr, "\n\n" );
            for( j=0; j<clus1; j++ ) 
            {
                  fprintf( stderr, "%s\n", tmpres1[j] );
            }
            fprintf( stderr, "-------\n" );
            for( j=0; j<clus2; j++ ) 
            {
                  fprintf( stderr, "%s\n", tmpres2[j] );
            }
#endif

#if 0
            fprintf( stdout, "writing input\n" );
            for( j=0; j<clus1; j++ )
            {
                  fprintf( stdout, ">%d of GROUP1\n", j );
                  fprintf( stdout, "%s\n", tmpres1[j] );
            }
            for( j=0; j<clus2; j++ )
            {
                  fprintf( stdout, ">%d of GROUP2\n", j );
                  fprintf( stdout, "%s\n", tmpres2[j] );
            }
            fflush( stdout );
#endif
            switch( alg )
            {
                  case( 'a' ):
                        totalscore += Aalign( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen );
                        break;
                  case( 'M' ):
                              totalscore += MSalignmm( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, sgap1, sgap2, egap1, egap2 );
                        break;
                  case( 'A' ):
                        if( clus1 == 1 && clus2 == 1 )
                        {
                              totalscore += G__align11( tmpres1, tmpres2, alloclen );
                        }
                        else
                              totalscore += A__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, sgap1, sgap2, egap1, egap2 );
                        break;
                  case( 'H' ):
                        if( clus1 == 1 && clus2 == 1 )
                        {
                              totalscore += G__align11( tmpres1, tmpres2, alloclen );
                        }
                        else
                              totalscore += H__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, sgap1, sgap2, egap1, egap2 );
                        break;
                  case( 'Q' ):
                        if( clus1 == 1 && clus2 == 1 )
                        {
                              totalscore += G__align11( tmpres1, tmpres2, alloclen );
                        }
                        else
                              totalscore += Q__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, sgap1, sgap2, egap1, egap2 );
                        break;
                  default:
                        fprintf( stderr, "alg = %c\n", alg );
                        ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
                        break;
            }

            nlen = strlen( tmpres1[0] );
            if( totallen + nlen > alloclen )
            {
                  fprintf( stderr, "totallen=%d +  nlen=%d > alloclen = %d\n", totallen, nlen, alloclen );
                  ErrorExit( "LENGTH OVER in Falign\n " );
            }
            for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
            for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
            totallen += nlen;
#if 0
            fprintf( stderr, "i=%d", i );
            fprintf( stderr, "%4d\n", totallen );
            fprintf( stderr, "\n\n" );
            for( j=0; j<clus1; j++ ) 
            {
                  fprintf( stderr, "%s\n", tmpres1[j] );
            }
            fprintf( stderr, "-------\n" );
            for( j=0; j<clus2; j++ ) 
            {
                  fprintf( stderr, "%s\n", tmpres2[j] );
            }
#endif
      }
#if KEIKA
      fprintf( stderr, "DP ... done   \n" );
#endif

      for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
      for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
#if 0
      for( j=0; j<clus1; j++ ) 
      {
            fprintf( stderr, "%s\n", result1[j] );
      }
      fprintf( stderr, "- - - - - - - - - - -\n" );
      for( j=0; j<clus2; j++ ) 
      {
            fprintf( stderr, "%s\n", result2[j] );
      }
#endif
      return( totalscore );
}








float Falign_noudp( char  **seq1, char  **seq2, 
                    double *eff1, double *eff2, 
                    int    clus1, int    clus2,
                    int alloclen, int *fftlog )
{
      int i, j, k, l, m, maxk;
      int nlen, nlen2, nlen4;
      static int prevalloclen = 0;
      static int crossscoresize = 0;
      static char **tmpseq1 = NULL;
      static char **tmpseq2 = NULL;
      static char **tmpptr1 = NULL;
      static char **tmpptr2 = NULL;
      static char **tmpres1 = NULL;
      static char **tmpres2 = NULL;
      static char **result1 = NULL;
      static char **result2 = NULL;
#if RND
      static char **rndseq1 = NULL;
      static char **rndseq2 = NULL;
#endif
      static Fukusosuu **seqVector1 = NULL;
      static Fukusosuu **seqVector2 = NULL;
      static Fukusosuu **naiseki = NULL;   
      static Fukusosuu *naisekiNoWa = NULL; 
      static double *soukan = NULL;
      int nlentmp;
      static int *kouho = NULL;
      static Segment *segment = NULL;
      static Segment *segment1 = NULL;
      static Segment *segment2 = NULL;
      static Segment **sortedseg1 = NULL;
      static Segment **sortedseg2 = NULL;
      static int *cut1 = NULL;
      static int *cut2 = NULL;
      static char *sgap1, *egap1, *sgap2, *egap2;
      static int localalloclen = 0;
      int lag;
      int tmpint;
      int count, count0;
      int len1, len2;
      int totallen;
      double totaleff1, totaleff2;
      float totalscore;
      int nkouho;
      float dumfl = 0.0;

      totaleff1 = 0.0; for( i=0; i<clus1; i++ ) totaleff1 += eff1[i];
      totaleff2 = 0.0; for( i=0; i<clus2; i++ ) totaleff2 += eff2[i];

      len1 = strlen( seq1[0] );
      len2 = strlen( seq2[0] );
      nlentmp = MAX( len1, len2 );

      nlen = 1;
      while( nlentmp >= nlen ) nlen <<= 1;
#if 0
      fprintf( stderr, "###   nlen    = %d\n", nlen );
#endif

      nlen2 = nlen/2; nlen4 = nlen2 / 2;

#if 0
      fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
      fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
#endif

      if( prevalloclen != alloclen ) // Falign_noudp mo kaeru
      {
            if( prevalloclen )
            {
                  FreeCharMtx( result1 );
                  FreeCharMtx( result2 );
                  FreeCharMtx( tmpres1 );
                  FreeCharMtx( tmpres2 );
            }
//          fprintf( stderr, "\n\n\nreallocating ...\n" ); 
            result1 = AllocateCharMtx( njob, alloclen );
            result2 = AllocateCharMtx( njob, alloclen );
            tmpres1 = AllocateCharMtx( njob, alloclen );
            tmpres2 = AllocateCharMtx( njob, alloclen );
            prevalloclen = alloclen;
      }

      if( !localalloclen )
      {
            sgap1 = AllocateCharVec( njob );
            egap1 = AllocateCharVec( njob );
            sgap2 = AllocateCharVec( njob );
            egap2 = AllocateCharVec( njob );
            kouho = AllocateIntVec( NKOUHO_LONG );
            cut1 = AllocateIntVec( MAXSEG );
            cut2 = AllocateIntVec( MAXSEG );
            tmpptr1 = AllocateCharMtx( njob, 0 );
            tmpptr2 = AllocateCharMtx( njob, 0 );
            segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
            segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
            segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
            sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
            sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
            if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
                  ErrorExit( "Allocation error\n" );

            if     ( scoremtx == -1 ) n20or4or2 = 1;
            else if( fftscore )       n20or4or2 = 1;
            else                      n20or4or2 = 20;
      }
      if( localalloclen < nlen )
      {
            if( localalloclen )
            {
#if 1
                  if( !kobetsubunkatsu )
                  {
                        FreeFukusosuuMtx ( seqVector1 );
                        FreeFukusosuuMtx ( seqVector2 );
                        FreeFukusosuuVec( naisekiNoWa );
                        FreeFukusosuuMtx( naiseki );
                        FreeDoubleVec( soukan );
                  }
                  FreeCharMtx( tmpseq1 );
                  FreeCharMtx( tmpseq2 );
#endif
#if RND
                  FreeCharMtx( rndseq1 );
                  FreeCharMtx( rndseq2 );
#endif
            }


            tmpseq1 = AllocateCharMtx( njob, nlen );
            tmpseq2 = AllocateCharMtx( njob, nlen );
            if( !kobetsubunkatsu )
            {
                  naisekiNoWa = AllocateFukusosuuVec( nlen );
                  naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
                  seqVector1 = AllocateFukusosuuMtx( n20or4or2, nlen+1 );
                  seqVector2 = AllocateFukusosuuMtx( n20or4or2, nlen+1 );
                  soukan = AllocateDoubleVec( nlen+1 );
            }
#if RND
            rndseq1 = AllocateCharMtx( njob, nlen );
            rndseq2 = AllocateCharMtx( njob, nlen );
            for( i=0; i<njob; i++ )
            {
                  generateRndSeq( rndseq1[i], nlen );
                  generateRndSeq( rndseq2[i], nlen );
            }
#endif
            localalloclen = nlen;
      }
      
      for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
      for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );

#if 0
fftfp = fopen( "input_of_Falign", "w" );
fprintf( fftfp, "nlen = %d\n", nlen );
fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
for( i=0; i<clus1; i++ )
      fprintf( fftfp, "%s\n", seq1[i] );
fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
for( i=0; i<clus2; i++ )
      fprintf( fftfp, "%s\n", seq2[i] );
fclose( fftfp );
system( "less input_of_Falign < /dev/tty > /dev/tty" );
#endif
      if( !kobetsubunkatsu )
      {
            fprintf( stderr,  " FFT ... " );

            for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
            if( scoremtx == -1 )
            {
                  for( i=0; i<clus1; i++ )
                        seq_vec_4( seqVector1[0], eff1[i], tmpseq1[i] );
            }
            else if( fftscore )
            {
                  for( i=0; i<clus1; i++ )
                  {
#if 0
                        seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
                        seq_vec_2( seqVector1[1], volume,   eff1[i], tmpseq1[i] );
#else
                        seq_vec_5( seqVector1[0], polarity, volume, eff1[i], tmpseq1[i] );
#endif
                  }
            }
            else
            {
                  for( i=0; i<clus1; i++ )
                        seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
            }
#if RND
            for( i=0; i<clus1; i++ )
            {
                  vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
            }
#endif
#if 0
fftfp = fopen( "seqVec", "w" );
fprintf( fftfp, "before transform\n" );
for( k=0; k<n20or4or2; k++ ) 
{
   fprintf( fftfp, "nlen=%d\n", nlen );
   fprintf( fftfp, "%c\n", amino[k] );
   for( l=0; l<nlen; l++ )
   fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
}
fclose( fftfp );
system( "less seqVec < /dev/tty > /dev/tty" );
#endif

            for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
            if( scoremtx == -1 )
            {
                  for( i=0; i<clus2; i++ )
                        seq_vec_4( seqVector2[0], eff2[i], tmpseq2[i] );
            }
            else if( fftscore )
            {
                  for( i=0; i<clus2; i++ )
                  {
#if 0
                        seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
                        seq_vec_2( seqVector2[1], volume,   eff2[i], tmpseq2[i] );
#else
                        seq_vec_5( seqVector2[0], polarity, volume, eff2[i], tmpseq2[i] );
#endif
                  }
            }
            else
            {
                  for( i=0; i<clus2; i++ )
                        seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
            }
#if RND
            for( i=0; i<clus2; i++ )
            {
                  vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
            }
#endif

#if 0
fftfp = fopen( "seqVec2", "w" );
fprintf( fftfp, "before fft\n" );
for( k=0; k<n20or4or2; k++ ) 
{
   fprintf( fftfp, "%c\n", amino[k] );
   for( l=0; l<nlen; l++ )
   fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
}
fclose( fftfp );
system( "less seqVec2 < /dev/tty > /dev/tty" );
#endif

            for( j=0; j<n20or4or2; j++ )
            {
                  fft( nlen, seqVector2[j], (j==0) );
                  fft( nlen, seqVector1[j], 0 );
            }
#if 0
fftfp = fopen( "seqVec2", "w" );
fprintf( fftfp, "#after fft\n" );
for( k=0; k<n20or4or2; k++ ) 
{
   fprintf( fftfp, "#%c\n", amino[k] );
   for( l=0; l<nlen; l++ )
         fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
}
fclose( fftfp );
system( "less seqVec2 < /dev/tty > /dev/tty" );
#endif

            for( k=0; k<n20or4or2; k++ ) 
            {
                  for( l=0; l<nlen; l++ ) 
                        calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
            }
            for( l=0; l<nlen; l++ ) 
            {
                  naisekiNoWa[l].R = 0.0;
                  naisekiNoWa[l].I = 0.0;
                  for( k=0; k<n20or4or2; k++ ) 
                  {
                        naisekiNoWa[l].R += naiseki[k][l].R;
                        naisekiNoWa[l].I += naiseki[k][l].I;
                  }
            }
      
#if 0
      fftfp = fopen( "naisekiNoWa", "w" );
      fprintf( fftfp, "#Before fft\n" );
      for( l=0; l<nlen; l++ )
            fprintf( fftfp, "%d  %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I ); 
      fclose( fftfp );
      system( "less naisekiNoWa < /dev/tty > /dev/tty " );
#endif

            fft( -nlen, naisekiNoWa, 0 );
      
            for( m=0; m<=nlen2; m++ ) 
                  soukan[m] = naisekiNoWa[nlen2-m].R;
            for( m=nlen2+1; m<nlen; m++ ) 
                  soukan[m] = naisekiNoWa[nlen+nlen2-m].R;

#if 0
      fftfp = fopen( "naisekiNoWa", "w" );
      fprintf( fftfp, "#After fft\n" );
      for( l=0; l<nlen; l++ )
            fprintf( fftfp, "%d  %f\n", l, naisekiNoWa[l].R ); 
      fclose( fftfp );
      fftfp = fopen( "list.plot", "w"  );
      fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
      fclose( fftfp );
      system( "/usr/bin/gnuplot list.plot &" );
#endif
#if 0
      fprintf( stderr, "soukan\n" );
      for( l=0; l<nlen; l++ )
            fprintf( stderr, "%d  %f\n", l-nlen2, soukan[l] ); 
#if 0
      fftfp = fopen( "list.plot", "w"  );
      fprintf( fftfp, "plot 'frt'\n pause +1" );
      fclose( fftfp );
      system( "/usr/bin/gnuplot list.plot" );
#endif
#endif


            nkouho = getKouho( kouho, NKOUHO_LONG, soukan, nlen );

#if 0
            for( i=0; i<nkouho; i++ )
            {
                  fprintf( stderr, "kouho[%d] = %d\n", i, kouho[i] );
            }
#endif
      }

#if KEIKA
      fprintf( stderr, "Searching anchors ... " );
#endif
      count = 0;



#define CAND 0
#if CAND
      fftfp = fopen( "cand", "w" );
      fclose( fftfp );
#endif
      if( kobetsubunkatsu )
      {
            maxk = 1;
            kouho[0] = 0;
      }
      else
      {
            maxk = nkouho;
      }

      for( k=0; k<maxk; k++ ) 
      {
            lag = kouho[k];
            if( lag <= -len1 || len2 <= lag ) continue;
//          fprintf( stderr, "k=%d, lag=%d\n", k, lag );
            zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
#if CAND
            fftfp = fopen( "cand", "a" );
            fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
            fprintf( fftfp, "%s\n", tmpptr1[0] );
            fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
            fprintf( fftfp, "%s\n", tmpptr2[0] );
            fprintf( fftfp, ">\n", k+1, lag );
            fclose( fftfp );
#endif

//          fprintf( stderr, "lag = %d\n", lag );
            tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
//          fprintf( stderr, "lag = %d, %d found\n", lag, tmpint );

//          if( lag == -50 ) exit( 1 );
            
            if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );

//          fprintf( stderr, "##### k=%d / %d\n", k, maxk );
            if( tmpint == 0 ) break; // 060430 iinoka ?
            while( tmpint-- > 0 )
            {
#if 0
                  if( segment[count].end - segment[count].start < fftWinSize )
                  {
                        count++;
                        continue;
                  }
#endif
                  if( lag > 0 )
                  {
                        segment1[count].start  = segment[count].start ;
                        segment1[count].end    = segment[count].end   ;
                        segment1[count].center = segment[count].center;
                        segment1[count].score  = segment[count].score;

                        segment2[count].start  = segment[count].start  + lag;
                        segment2[count].end    = segment[count].end    + lag;
                        segment2[count].center = segment[count].center + lag;
                        segment2[count].score  = segment[count].score       ;
                  }
                  else
                  {
                        segment1[count].start  = segment[count].start  - lag;
                        segment1[count].end    = segment[count].end    - lag;
                        segment1[count].center = segment[count].center - lag;
                        segment1[count].score  = segment[count].score       ;

                        segment2[count].start  = segment[count].start ;
                        segment2[count].end    = segment[count].end   ;
                        segment2[count].center = segment[count].center;
                        segment2[count].score  = segment[count].score ;
                  }
#if 0
                  fprintf( stderr, "##### k=%d / %d\n", k, maxk );
                  fprintf( stderr, "anchor %d, score = %f\n", count, segment1[count].score );
                  fprintf( stderr, "in 1 %d\n", segment1[count].center );
                  fprintf( stderr, "in 2 %d\n", segment2[count].center );
#endif
                  segment1[count].pair = &segment2[count];
                  segment2[count].pair = &segment1[count];
                  count++;
#if 0
                  fprintf( stderr, "count=%d\n", count );
#endif
            }
      }
#if 1
      if( !kobetsubunkatsu )
            fprintf( stderr, "done. (%d anchors) ", count );
#endif
      if( !count && fftNoAnchStop )
            ErrorExit( "Cannot detect anchor!" );
#if 0
      fprintf( stderr, "RESULT before sort:\n" );
      for( l=0; l<count+1; l++ )
      {
            fprintf( stderr, "cut[%d]=%d, ", l, segment1[l].center );
            fprintf( stderr, "%d score = %f\n", segment2[l].center, segment1[l].score );
      }
#endif

      for( i=0; i<count; i++ )
      {
            sortedseg1[i] = &segment1[i];
            sortedseg2[i] = &segment2[i];
      }
#if 0
      tmpsort( count, sortedseg1 ); 
      tmpsort( count, sortedseg2 ); 
      qsort( sortedseg1, count, sizeof( Segment * ), segcmp );
      qsort( sortedseg2, count, sizeof( Segment * ), segcmp );
#else
      mymergesort( 0, count-1, sortedseg1 ); 
      mymergesort( 0, count-1, sortedseg2 ); 
#endif
      for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
      for( i=0; i<count; i++ ) sortedseg2[i]->number = i;



      if( kobetsubunkatsu )
      {
            for( i=0; i<count; i++ )
          {
                  cut1[i+1] = sortedseg1[i]->center;
                  cut2[i+1] = sortedseg2[i]->center;
            }
            cut1[0] = 0;
            cut2[0] = 0;
            cut1[count+1] = len1;
            cut2[count+1] = len2;
            count += 2;
      }
      else
      {
            cut1[0] = 0; 
            cut2[0] = 0;
            count0 = 0;
            for( i=0; i<count; i++ )
            {
//                fprintf( stderr, "i=%d, %d-%d ?\n", i, sortedseg1[i]->center, sortedseg1[i]->pair->center );
                  if( sortedseg1[i]->center > cut1[count0]
                   && sortedseg1[i]->pair->center > cut2[count0] )
                  {
                        count0++;
                        cut1[count0] = sortedseg1[i]->center;
                        cut2[count0] = sortedseg1[i]->pair->center;
                  }
                  else
                  {
                        if( i && sortedseg1[i]->score > sortedseg1[i-1]->score )
                        {
                              if( sortedseg1[i]->center > cut1[count0-1]
                               && sortedseg1[i]->pair->center > cut2[count0-1] )
                              {
                                    cut1[count0] = sortedseg1[i]->center;
                                    cut2[count0] = sortedseg1[i]->pair->center;
                              }
                              else
                              {
//                                  count0--;
                              }
                        }
                  }
            }
//          if( count-count0 )
//                fprintf( stderr, "%d anchors unused\n", count-count0 );
            cut1[count0+1] = len1;
            cut2[count0+1] = len2;
            count = count0 + 2;
            count0 = count;

      }

//    exit( 0 );

#if 0
      fftfp = fopen( "fft", "a" );
      fprintf( fftfp, "RESULT after sort:\n" );
      for( l=0; l<count; l++ )
      {
            fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
            fprintf( fftfp, "%d\n", segment2[l].center );
      }
      fclose( fftfp );
#endif

#if 0
      fprintf( stderr, "RESULT after blckalign:\n" );
      for( l=0; l<count+1; l++ )
      {
            fprintf( stderr, "cut : %d %d\n", cut1[l], cut2[l] );
      }
#endif

#if 0
      fprintf( trap_g, "Devided to %d segments\n", count-1 );
      fprintf( trap_g, "%d  %d forg\n", MIN( clus1, clus2 ), count-1 );
#endif

      totallen = 0;
      for( j=0; j<clus1; j++ ) result1[j][0] = 0;
      for( j=0; j<clus2; j++ ) result2[j][0] = 0;
      totalscore = 0.0;
      *fftlog = -1;
      for( i=0; i<count-1; i++ )
      {
            *fftlog += 1;

            if( cut1[i] )
            {
//                getkyokaigap( sgap1, seq1, cut1[i]-1, clus1 );
//                getkyokaigap( sgap2, seq2, cut2[i]-1, clus2 );
                  getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
                  getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
            }
            else
            {
                  for( j=0; j<clus1; j++ ) sgap1[j] = 'o';
                  for( j=0; j<clus2; j++ ) sgap2[j] = 'o';
            }
            if( cut1[i+1] != len1 )
            {       
                  getkyokaigap( egap1, seq1, cut1[i+1], clus1 );
                  getkyokaigap( egap2, seq2, cut2[i+1], clus2 );
            }       
            else    
            {       
                  for( j=0; j<clus1; j++ ) egap1[j] = 'o';
                  for( j=0; j<clus2; j++ ) egap2[j] = 'o';
            }
#if DEBUG
            fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
#else
#if 1
            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 );
#endif
#endif
            for( j=0; j<clus1; j++ )
            {
                  strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
                  tmpres1[j][cut1[i+1]-cut1[i]] = 0;
            }
            if( kobetsubunkatsu && fftkeika ) commongappick( clus1, tmpres1 ); //dvtditr に呼ばれたとき fftkeika=1
//          if( kobetsubunkatsu ) commongappick( clus1, tmpres1 );
            for( j=0; j<clus2; j++ )
            {
//                fprintf( stderr, "### cut2[i+1]-cut2[i] = %d\n", cut2[i+1]-cut2[i] );
                  if( cut2[i+1]-cut2[i] <= 0 )
                        fprintf( stderr, "### cut2[i+1]=%d, cut2[i]=%d\n", cut2[i+1], cut2[i] );
                  strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
                  tmpres2[j][cut2[i+1]-cut2[i]] = 0;
            }
            if( kobetsubunkatsu && fftkeika ) commongappick( clus2, tmpres2 ); //dvtditr に呼ばれたとき fftkeika=1
//          if( kobetsubunkatsu ) commongappick( clus2, tmpres2 );

            if( constraint )
            {
                  fprintf( stderr, "Not supported\n" );
                  exit( 1 );
            }
#if 0
            fprintf( stderr, "i=%d, before alignment", i );
            fprintf( stderr, "%4d\n", totallen );
            fprintf( stderr, "\n\n" );
            for( j=0; j<clus1; j++ ) 
            {
                  fprintf( stderr, "%s\n", tmpres1[j] );
            }
            fprintf( stderr, "-------\n" );
            for( j=0; j<clus2; j++ ) 
            {
                  fprintf( stderr, "%s\n", tmpres2[j] );
            }
#endif

#if 0
            fprintf( stdout, "writing input\n" );
            for( j=0; j<clus1; j++ )
            {
                  fprintf( stdout, ">%d of GROUP1\n", j );
                  fprintf( stdout, "%s\n", tmpres1[j] );
            }
            for( j=0; j<clus2; j++ )
            {
                  fprintf( stdout, ">%d of GROUP2\n", j );
                  fprintf( stdout, "%s\n", tmpres2[j] );
            }
            fflush( stdout );
#endif
            switch( alg )
            {
                  case( 'M' ):
                              totalscore += MSalignmm( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, sgap1, sgap2, egap1, egap2 );
                        break;
                  default:
                        fprintf( stderr, "alg = %c\n", alg );
                        ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
                        break;
            }

            nlen = strlen( tmpres1[0] );
            if( totallen + nlen > alloclen )
            {
                  fprintf( stderr, "totallen=%d +  nlen=%d > alloclen = %d\n", totallen, nlen, alloclen );
                  ErrorExit( "LENGTH OVER in Falign\n " );
            }
            for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
            for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
            totallen += nlen;
#if 0
            fprintf( stderr, "i=%d", i );
            fprintf( stderr, "%4d\n", totallen );
            fprintf( stderr, "\n\n" );
            for( j=0; j<clus1; j++ ) 
            {
                  fprintf( stderr, "%s\n", tmpres1[j] );
            }
            fprintf( stderr, "-------\n" );
            for( j=0; j<clus2; j++ ) 
            {
                  fprintf( stderr, "%s\n", tmpres2[j] );
            }
#endif
      }
#if KEIKA
      fprintf( stderr, "DP ... done   \n" );
#endif

      for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
      for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
#if 0
      for( j=0; j<clus1; j++ ) 
      {
            fprintf( stderr, "%s\n", result1[j] );
      }
      fprintf( stderr, "- - - - - - - - - - -\n" );
      for( j=0; j<clus2; j++ ) 
      {
            fprintf( stderr, "%s\n", result2[j] );
      }
#endif
      return( totalscore );
}

Generated by  Doxygen 1.6.0   Back to index