JWS-112 Bumping version of Mafft to version 7.310.
[jabaws.git] / binaries / src / mafft / core / splittbfast.c
index a810285..2cbcdfb 100644 (file)
@@ -26,14 +26,11 @@ static int uselongest;
 static int treeout;
 static int classsize;
 static int picksize;
-static int maxl;
-static int tsize;
 static int reorder;
 static int pid;
 static int maxdepth = 0;
 static double tokyoripara;
 
-static double lenfaca, lenfacb, lenfacc, lenfacd;
 #define PLENFACA 0.01
 #define PLENFACB 10000
 #define PLENFACC 10000
@@ -89,73 +86,6 @@ int dcompare( const Scores *a, const Scores *b )
        }
 }
 
-#if 0
-static void    gappickandx0( char *out, char *in )
-{
-       char c;
-       if( scoremtx == -1 )
-       {
-               while( *in )
-               {
-                       if( (c=*in++) == '-' )
-                               ;       
-                       else if( c == 'u' )
-                               *out++ = 't';
-                       else if( amino_n[c] < 4 && amino_n[c] > -1 )
-                               *out++ = c;
-                       else
-                               *out++ = 'n';
-               }
-       }
-       else
-       {
-               while( *in )
-               {
-                       if( (c=*in++) == '-' )
-                               ;       
-                       else if( amino_n[c] < 20 && amino_n[c] > -1 )
-                               *out++ = c;
-                       else
-                               *out++ = 'X';
-               }
-       }
-       *out = 0;
-}      
-
-static int getkouho( int *pickkouho, double prob, int nin, Scores *scores, char **seq ) // 0 < prob < 1
-{
-       int nkouho = 0;
-       int i, j;
-       int *iptr = pickkouho;
-       for( i=1; i<nin; i++ )
-       {
-               if( ( nkouho==0 || rnd() < prob ) && ( scores[i].shimon != scores->shimon || strcmp( seq[scores->numinseq], seq[scores[i].numinseq] ) ) )
-               {
-#if 0
-                       for( j=0; j<nkouho; j++ )
-                       {
-                               if( scores[i].shimon == scores[pickkouho[j]].shimon || !strcmp( seq[scores[pickkouho[j]].numinseq], seq[scores[i].numinseq] ) ) 
-                                       break;
-                       }
-                       if( j == nkouho )
-#endif
-                       {
-                               *iptr++ = i;
-                               nkouho++;
-//                             fprintf( stderr, "ok! nkouho=%d\n", nkouho );
-                       }
-               }
-               else
-               {
-                       ;
-//                     fprintf( stderr, "no! %d-%d\n", 0, scores[i].numinseq );
-               }
-       }
-       fprintf( stderr, "\ndone\n\n"  );
-       return nkouho;
-}
-
-#endif
 
 static void getfastascoremtx( int **tmpaminodis )
 {
@@ -636,12 +566,14 @@ void arguments( int argc, char *argv[] )
        tbrweight = 3;
        checkC = 0;
        treemethod = 'X';
+       sueff_global = 0.1;
        contin = 0;
        scoremtx = 1;
        kobetsubunkatsu = 0;
        dorp = NOTSPECIFIED;
        ppenalty = -1530;
        ppenalty_ex = NOTSPECIFIED;
+       penalty_shift_factor = 1000.0;
        poffset = -123;
        kimuraR = NOTSPECIFIED;
        pamN = NOTSPECIFIED;
@@ -652,6 +584,9 @@ void arguments( int argc, char *argv[] )
        classsize = NOTSPECIFIED;
        picksize = NOTSPECIFIED;
        tokyoripara = NOTSPECIFIED;
+       legacygapcost = 0;
+       nwildcard = 0;
+       outnumber = 0;
 
     while( --argc > 0 && (*++argv)[0] == '-' )
        {
@@ -660,12 +595,12 @@ void arguments( int argc, char *argv[] )
             switch( c )
             {
                                case 'p':
-                                       picksize = atoi( *++argv );
+                                       picksize = myatoi( *++argv );
                                        fprintf( stderr, "picksize = %d\n", picksize );
                                        --argc;
                                        goto nextoption;
                                case 's':
-                                       classsize = atoi( *++argv );
+                                       classsize = myatoi( *++argv );
                                        fprintf( stderr, "groupsize = %d\n", classsize );
                                        --argc;
                                        goto nextoption;
@@ -679,6 +614,10 @@ void arguments( int argc, char *argv[] )
 //                                     fprintf( stderr, "ppenalty = %d\n", ppenalty );
                                        --argc;
                                        goto nextoption;
+                               case 'Q':
+                                       penalty_shift_factor = atof( *++argv );
+                                       --argc;
+                                       goto nextoption;
                                case 'g':
                                        ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
                                        fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
@@ -690,25 +629,25 @@ void arguments( int argc, char *argv[] )
                                        --argc;
                                        goto nextoption;
                                case 'k':
-                                       kimuraR = atoi( *++argv );
+                                       kimuraR = myatoi( *++argv );
                                        fprintf( stderr, "kimuraR = %d\n", kimuraR );
                                        --argc;
                                        goto nextoption;
                                case 'b':
-                                       nblosum = atoi( *++argv );
+                                       nblosum = myatoi( *++argv );
                                        scoremtx = 1;
 //                                     fprintf( stderr, "blosum %d\n", nblosum );
                                        --argc;
                                        goto nextoption;
                                case 'j':
-                                       pamN = atoi( *++argv );
+                                       pamN = myatoi( *++argv );
                                        scoremtx = 0;
                                        TMorJTT = JTT;
                                        fprintf( stderr, "jtt %d\n", pamN );
                                        --argc;
                                        goto nextoption;
                                case 'm':
-                                       pamN = atoi( *++argv );
+                                       pamN = myatoi( *++argv );
                                        scoremtx = 0;
                                        TMorJTT = TM;
                                        fprintf( stderr, "tm %d\n", pamN );
@@ -721,6 +660,9 @@ void arguments( int argc, char *argv[] )
                                case 'l':
                                        uselongest = 0;
                                        break;
+                               case 'n' :
+                                       outnumber = 1;
+                                       break;
 #if 1
                                case 'a':
                                        fmodel = 1;
@@ -732,7 +674,7 @@ void arguments( int argc, char *argv[] )
                                case 'Z':
                                        fromaln = 1;
                                        break;
-                               case 'L':
+                               case 'U':
                                        doalign = 1;
                                        break;
                                case 'x':
@@ -756,6 +698,9 @@ void arguments( int argc, char *argv[] )
                                case 'O':
                                        fftNoAnchStop = 1;
                                        break;
+                               case 'L':
+                                       legacygapcost = 1;
+                                       break;
 #if 0
                                case 'R':
                                        fftRepeatStop = 1;
@@ -766,13 +711,13 @@ void arguments( int argc, char *argv[] )
                                case 'a':
                                        alg = 'a';
                                        break;
-#endif
                                case 'R':
                                        alg = 'R';
                                        break;
                                case 'Q':
                                        alg = 'Q';
                                        break;
+#endif
                                case 'A':
                                        alg = 'A';
                                        break;
@@ -805,8 +750,11 @@ void arguments( int argc, char *argv[] )
                                        tbutree = 0;
                                        break;
                                case 'X':
-                                       treemethod = 'X'; // mix
-                                       break;
+                                       treemethod = 'X'; // tsukawareteiru ????
+                                       sueff_global = atof( *++argv );
+                                       fprintf( stderr, "sueff_global = %f\n", sueff_global );
+                                       --argc;
+                                       goto nextoption;
                                case 'E':
                                        treemethod = 'E'; // upg (average)
                                        break;
@@ -814,13 +762,16 @@ void arguments( int argc, char *argv[] )
                                        treemethod = 'q'; // minimum
                                        break;
                                case 'z':
-                                       fftThreshold = atoi( *++argv );
+                                       fftThreshold = myatoi( *++argv );
                                        --argc; 
                                        goto nextoption;
                                case 'w':
-                                       fftWinSize = atoi( *++argv );
+                                       fftWinSize = myatoi( *++argv );
                                        --argc;
                                        goto nextoption;
+                               case ':':
+                                       nwildcard = 1;
+                                       break;
                 default:
                     fprintf( stderr, "illegal option %c\n", c );
                     argc = 0;
@@ -852,8 +803,7 @@ void arguments( int argc, char *argv[] )
        }
 }
 
-static int maxl;
-static int tsize;
+static int nunknown = 0;
 
 int seq_grp_nuc( int *grp, char *seq )
 {
@@ -865,7 +815,7 @@ int seq_grp_nuc( int *grp, char *seq )
                if( tmp < 4 )
                        *grp++ = tmp;
                else
-                       fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
+                       nunknown++;
        }
        *grp = END_OF_VEC;
        return( grp-grpbk );
@@ -881,7 +831,7 @@ int seq_grp( int *grp, char *seq )
                if( tmp < 6 )
                        *grp++ = tmp;
                else
-                       fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
+                       nunknown++;
        }
        *grp = END_OF_VEC;
        return( grp-grpbk );
@@ -895,7 +845,7 @@ void makecompositiontable_p( short *table, int *pointt )
                table[point]++;
 }
 
-int commonsextet_p( short *table, int *pointt )
+static int localcommonsextet_p( short *table, int *pointt )
 {
        int value = 0;
        short tmp;
@@ -982,13 +932,14 @@ static void pairalign( int nseq, int *nlen, char **seq, int *mem1, int *mem2, do
 {
        int l, len1, len2;
        int clus1, clus2;
-       float pscore, tscore;
+       double pscore, tscore;
        static int *fftlog;
        static char *indication1, *indication2;
        static double *effarr1 = NULL;
        static double *effarr2 = NULL;
        static char **mseq1, **mseq2;
-       float dumfl = 0.0;
+//     double dumfl = 0.0;
+       double dumdb = 0.0;
        int ffttry;
        int m1, m2;
 #if 0
@@ -1022,8 +973,8 @@ static void pairalign( int nseq, int *nlen, char **seq, int *mem1, int *mem2, do
        }
 
 #if WEIGHT
-       clus1 = fastconjuction_noname( mem1, seq, mseq1, effarr1, weight, indication1 );
-       clus2 = fastconjuction_noname( mem2, seq, mseq2, effarr2, weight, indication2 );
+       clus1 = fastconjuction_noname( mem1, seq, mseq1, effarr1, weight, indication1, 0.0 );
+       clus2 = fastconjuction_noname( mem2, seq, mseq2, effarr2, weight, indication2, 0.0 );
 #else
        clus1 = fastconjuction_noweight( mem1, seq, mseq1, effarr1, indication1 );
        clus2 = fastconjuction_noweight( mem2, seq, mseq2, effarr2, indication2 );
@@ -1053,6 +1004,7 @@ static void pairalign( int nseq, int *nlen, char **seq, int *mem1, int *mem2, do
                fprintf( stderr, "\nlen1=%d, len2=%d, Switching to the memsave mode\n", len1, len2 );
                alg = 'M';
                if( commonIP ) FreeIntMtx( commonIP );
+               commonIP = 0;
                commonAlloc1 = 0;
                commonAlloc2 = 0;
        }
@@ -1067,12 +1019,12 @@ static void pairalign( int nseq, int *nlen, char **seq, int *mem1, int *mem2, do
                {
                        fprintf( stderr, "\bm" );
 //                     fprintf( stderr, "%d-%d", clus1, clus2 );
-                       pscore = Falign_udpari_long( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1 );
+                       pscore = Falign_udpari_long( NULL, NULL, n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, NULL, NULL, clus1, clus2, *alloclen, fftlog+m1 );
                }
                else
                {
 //                     fprintf( stderr, "%d-%d", clus1, clus2 );
-                       pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL );
+                       pscore = Falign( NULL, NULL, n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, NULL, NULL, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL );
                }
        }
        else
@@ -1087,30 +1039,18 @@ static void pairalign( int nseq, int *nlen, char **seq, int *mem1, int *mem2, do
                        case( 'M' ):
                                fprintf( stderr, "\bm" );
 //                             fprintf( stderr, "%d-%d", clus1, clus2 );
-                               pscore = MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
-                               break;
-                       case( 'Q' ):
-                               if( clus1 == 1 && clus2 == 1 )
-                               {
-//                                     fprintf( stderr, "%d-%d", clus1, clus2 );
-                                       pscore = G__align11( mseq1, mseq2, *alloclen, outgap, outgap );
-                               }
-                               else
-                               {
-//                                     fprintf( stderr, "%d-%d", clus1, clus2 );
-                                       pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
-                               }
+                               pscore = MSalignmm( n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
                                break;
                        case( 'A' ):
                                if( clus1 == 1 && clus2 == 1 )
                                {
 //                                     fprintf( stderr, "%d-%d", clus1, clus2 );
-                                       pscore = G__align11( mseq1, mseq2, *alloclen, outgap, outgap );
+                                       pscore = G__align11( n_dis_consweight_multi, mseq1, mseq2, *alloclen, outgap, outgap );
                                }
                                else
                                {
 //                                     fprintf( stderr, "%d-%d", clus1, clus2 );
-                                       pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
+                                       pscore = A__align( n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumdb, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap, -1, -1 );
                                }
                                break;
                        default:
@@ -1225,8 +1165,8 @@ static int splitseq_mq( Scores *scores, int nin, int *nlen, char **seq, char **o
        int selfscore0;
        double **dfromc;
        double **dfromcp;
-       float **pickmtx;
-       float **yukomtx;
+       double **pickmtx;
+       double **yukomtx;
        static short *table1;
        Scores **outs, *ptr;
        int *numin;
@@ -1253,10 +1193,10 @@ static int splitseq_mq( Scores *scores, int nin, int *nlen, char **seq, char **o
        int ***topol;
        int *treeorder;
        int picktmp;
-       float **len;
+       double **len;
        double minscore;
 //     double *minscoreinpick;
-       float *hanni;
+       double *hanni;
        double lenfac;
        double longer;
        double shorter;
@@ -1264,7 +1204,7 @@ static int splitseq_mq( Scores *scores, int nin, int *nlen, char **seq, char **o
        static char **mseq2 = NULL;
        double *blastresults = NULL; // by Mathog, a guess
        static int palloclen = 0;
-       float maxdist;
+       double maxdist;
 
        if( orderpos == NULL )
                orderpos = order;
@@ -1322,16 +1262,24 @@ static int splitseq_mq( Scores *scores, int nin, int *nlen, char **seq, char **o
                        }
                        free( tmptree );
        
-                       *tree = (char *)calloc( treelen + nin + 5, sizeof( char ) );
-                       if( nin > 1 ) **tree = '(';
-                       else              **tree = '\0';
-//                     **tree = '\0';
+                       *tree = (char *)calloc( treelen + nin + 15, sizeof( char ) );
+                       **tree = '\n';
+                       if( nin > 1 ) 
+                       {
+                               *(*tree+1) = '(';
+                               *(*tree+2) = '\0';
+                       }
+                       else
+                       {
+                               *(*tree+1) = '\0';
+                       }
                        for( j=0; j<nin-1; j++ )
                        {
                                sprintf( *tree+strlen( *tree ), "%d,", scores[j].numinseq+1 );
                        }
                        sprintf( *tree+strlen( *tree ), "%d", scores[j].numinseq+1 );
-                       if( nin > 1 ) strcat( *tree, ")" );
+                       if( nin > 1 ) strcat( *tree, ")\n" );
+                       else strcat( *tree, "\n" );
 //                     fprintf( stdout, "*tree = %s\n", *tree );
                }
 
@@ -1458,7 +1406,7 @@ static int splitseq_mq( Scores *scores, int nin, int *nlen, char **seq, char **o
                                {
                                        if( fromaln )
                                        {
-//                                             scores[i].score = ( 1.0 - (double)G__align11_noalign( amino_disLN, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
+//                                             scores[i].score = ( 1.0 - (double)G__align11_noalign( n_disLN, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
                                                scores[i].score = ( 1.0 - (double)naivepairscore11( orialn[scores[i].numinseq], orialn[scores->numinseq], penalty ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
                                        }
                                        else
@@ -1466,7 +1414,7 @@ static int splitseq_mq( Scores *scores, int nin, int *nlen, char **seq, char **o
                                                if( *depthpt == 0 ) fprintf( stderr, "\r%d / %d   ", i, nin );
                                                gappick0( mseq2[0], seq[scores[i].numinseq] );
 //                                             fprintf( stdout, "### before calc scores[%d] = %f (%c)\n", i, scores[i].score, qinoya == scores->numinseq?'o':'x' );
-                                               scores[i].score = ( 1.0 - (double)G__align11_noalign( amino_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
+                                               scores[i].score = ( 1.0 - (double)G__align11_noalign( n_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
 //                                             fprintf( stderr, "scores[i] = %f\n", scores[i].score );
 //                                             fprintf( stderr, "m1=%s\n", seq[scores[0].numinseq] );
 //                                             fprintf( stderr, "m2=%s\n", seq[scores[i].numinseq] );
@@ -1476,7 +1424,7 @@ static int splitseq_mq( Scores *scores, int nin, int *nlen, char **seq, char **o
                        }
                        else
                        {
-                               scores[i].score = ( 1.0 - (double)commonsextet_p( table1, scores[i].pointt ) / MIN( selfscore0, scores[i].selfscore ) ) * lenfac;
+                               scores[i].score = ( 1.0 - (double)localcommonsextet_p( table1, scores[i].pointt ) / MIN( selfscore0, scores[i].selfscore ) ) * lenfac;
                                if( scores[i].score > MAX6DIST ) scores[i].score = MAX6DIST;
                        }
 //                     if( i ) fprintf( stderr, "%d-%d d %4.2f len %d %d\n", 1, i+1, scores[i].score, scores->orilen, scores[i].orilen );
@@ -1655,7 +1603,7 @@ exit( 1 );
        {
                if( s_p_map[j] != -1 )
                {
-                       pickmtx[0][s_p_map[j]] = (float)scores[j].score;
+                       pickmtx[0][s_p_map[j]] = (double)scores[j].score;
 //                     fprintf( stderr, "pickmtx[0][%d] = %f\n", s_p_map[j], pickmtx[0][s_p_map[j]] );
                }
        }
@@ -1730,14 +1678,14 @@ exit( 1 );
                                        {
 //                                             fprintf( stderr, "\r%d / %d   ", i, nin );
                                                gappick0( mseq2[0], seq[scores[picks[i]].numinseq] );
-                                               pickmtx[j][i-j] = ( 1.0 - (double)G__align11_noalign( amino_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[picks[i]].selfscore ) ) * 1;
+                                               pickmtx[j][i-j] = ( 1.0 - (double)G__align11_noalign( n_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[picks[i]].selfscore ) ) * 1;
        //                                      fprintf( stderr, "scores[picks[i]] = %f\n", scores[picks[i]].score );
                                        }
                                }
                        }
                        else
                        {
-                               pickmtx[j][i-j] = ( 1.0 - (double)commonsextet_p( table1, scores[picks[i]].pointt ) / MIN( selfscore0, scores[picks[i]].selfscore ) ) * lenfac;
+                               pickmtx[j][i-j] = ( 1.0 - (double)localcommonsextet_p( table1, scores[picks[i]].pointt ) / MIN( selfscore0, scores[picks[i]].selfscore ) ) * lenfac;
                                if( pickmtx[j][i-j] > MAX6DIST ) pickmtx[j][i-j] = MAX6DIST;
                        }
 
@@ -1800,10 +1748,10 @@ exit( 1 );
        fprintf( stderr, "DIANA!!\n" );
        if( npick > 2 )
        {
-               float avdist;
-               float avdist1;
-               float avdist2;
-               float maxavdist;
+               double avdist;
+               double avdist1;
+               double avdist2;
+               double maxavdist;
                int splinter;
                int count;
                int dochokoho;
@@ -1861,7 +1809,7 @@ exit( 1 );
                                        }
                                }
                                if( count < 1 ) avdist1 = 0.0;
-                               else avdist1 /= (float)count;
+                               else avdist1 /= (double)count;
                                fprintf( stderr, "docho %d (%dinori), avdist1 = %f\n", dochokoho, p_o_map[dochokoho] + 1, avdist1 );
 
                                count = 0;
@@ -1882,7 +1830,7 @@ exit( 1 );
                                        }
                                }
                                if( count < 1 ) avdist2 = 0.0;
-                               else avdist2 /= (float)count;
+                               else avdist2 /= (double)count;
                                fprintf( stderr, "docho %d (%dinori), avdist2 = %f\n", dochokoho, p_o_map[dochokoho] + 1, avdist2 );
 
                                if( avdist2 < avdist1 ) 
@@ -1941,8 +1889,8 @@ exit( 1 );
        if( npick > 2 )
        {
 #if 0
-               float avdist;
-               float maxavdist;
+               double avdist;
+               double maxavdist;
                int count;
                int splinter;
                maxavdist = 0.0;
@@ -1978,9 +1926,9 @@ exit( 1 );
                        if( tsukau[i] == 0 ) continue;
                        for( j=i+1; j<npick; j++ )
                        {
-//                             float kijun = maxdist *  1/(npick-2);
-//                             float kijun = maxavdist * tokyoripara;
-                               float kijun;
+//                             double kijun = maxdist *  1/(npick-2);
+//                             double kijun = maxavdist * tokyoripara;
+                               double kijun;
                                kijun = maxdist * tokyoripara;  // atode kakunin
 //                             fprintf( stderr, "%d-%d\n", i, j );
 //                             fprintf( stderr, "maxdist = %f\n", maxdist );
@@ -2249,13 +2197,13 @@ exit( 1 );
                                                        else
                                                        {
                                                                gappick0( mseq2[0], seq[scores[j].numinseq] );
-                                                               dfromc[i][j] = ( 1.0 - (double)G__align11_noalign( amino_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[j].selfscore ) ) * 1;
+                                                               dfromc[i][j] = ( 1.0 - (double)G__align11_noalign( n_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[j].selfscore ) ) * 1;
                                                        }
                                                }
                                        }
                                        else
                                        {
-                                               dfromc[i][j] = ( 1.0 - (double)commonsextet_p( table1, scores[j].pointt ) / MIN( selfscore0, scores[j].selfscore ) ) * lenfac;
+                                               dfromc[i][j] = ( 1.0 - (double)localcommonsextet_p( table1, scores[j].pointt ) / MIN( selfscore0, scores[j].selfscore ) ) * lenfac;
                                                if( dfromc[i][j] > MAX6DIST ) dfromc[i][j] = MAX6DIST;
                                        }
                                }
@@ -2345,8 +2293,8 @@ exit( 1 );
                if( nyuko > 2 )
                {
                        fprintf( stderr, "upgma       " );
-//                     veryfastsupg_float_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len );
-                       fixed_musclesupg_float_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len, NULL );
+//                     veryfastsupg_double_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len );
+                       fixed_musclesupg_double_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len, NULL, 1, 1 );
                        fprintf( stderr, "\r                      \r" );
                }
                else
@@ -2470,9 +2418,9 @@ exit( 1 );
                        mem2 = AllocateIntVec( njob+1 );
                }
 
-//             veryfastsupg_float_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len );
+//             veryfastsupg_double_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len );
        
-//             counteff_simple_float( nyuko, topol, len, eff );
+//             counteff_simple_double( nyuko, topol, len, eff );
 
 
                nlim = nyuko-1;
@@ -2898,6 +2846,7 @@ int main( int argc, char *argv[] )
 //             fprintf( stdout, ">%s\n", name[i] );
 //             fprintf( stdout, "%s\n", seq[i] );
        }
+       if( nunknown ) fprintf( stderr, "\nThere are %d ambiguous characters\n", nunknown );
 //     exit( 1 );
 
 #if 0
@@ -3033,6 +2982,7 @@ int main( int argc, char *argv[] )
                                pscore = 0;
                                for( pt=seq[i]; *pt; pt++ )
                                {
+//                                     pscore += amino_dis[(int)*pt][(int)*pt];
                                        pscore += amino_dis[(int)*pt][(int)*pt];
                                }
                                scores[i].selfscore = pscore; 
@@ -3044,7 +2994,7 @@ int main( int argc, char *argv[] )
                        table1 = (short *)calloc( tsize, sizeof( short ) );
                        if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
                        makecompositiontable_p( table1, pointt[i] );
-                       scores[i].selfscore = commonsextet_p( table1, pointt[i] );
+                       scores[i].selfscore = localcommonsextet_p( table1, pointt[i] );
                        free( table1 );
                }
        }