JWS-112 Bumping version of Mafft to version 7.310.

[jabaws.git] / binaries / src / mafft / core / mccaskillwrap.c

diff --git a/binaries/src/mafft/core/mccaskillwrap.c b/binaries/src/mafft/core/mccaskillwrap.c
index a7ea4ea..0f86a50 100644 (file)
--- a/binaries/src/mafft/core/mccaskillwrap.c
+++ b/binaries/src/mafft/core/mccaskillwrap.c
@@ -36,7 +36,7 @@ static void readrawmccaskill( FILE *fp, RNApair **pairprob, int length )
        int *pairnum;
        int i;
        int left, right;
-       float prob;
+       double prob;
 
        pairnum = (int *)calloc( length, sizeof( int ) );
        for( i=0; i<length; i++ ) pairnum[i] = 0;
@@ -46,7 +46,7 @@ static void readrawmccaskill( FILE *fp, RNApair **pairprob, int length )
                fgets( gett, 999, fp );
                if( feof( fp ) ) break;
                if( gett[0] == '>' ) continue;
-               sscanf( gett, "%d %d %f", &left, &right, &prob );
+               sscanf( gett, "%d %d %lf", &left, &right, &prob );
                if( prob < 0.01 ) continue; // mxscarna to mafft ryoho ni eikyou
 //fprintf( stderr, "gett = %s\n", gett );
 
@@ -101,11 +101,14 @@ static void *athread( void *arg )
                if( i == njob )
                {
                        pthread_mutex_unlock( targ->mutex );
-                       return( NULL );
+//                     return( NULL );
+                       break;
                }
                *jobpospt = i+1;
                pthread_mutex_unlock( targ->mutex );
 
+               commongappick_record( 1, nogap+i, gapmap[i] );
+               if( strlen( nogap[i] ) == 0 ) continue;
 
                sprintf( dirname, "_%d", i );
                sprintf( com, "rm -rf %s", dirname );
@@ -114,7 +117,6 @@ static void *athread( void *arg )
                system( com );
 
                fprintf( stderr, "%d / %d (by thread %4d)\n", i+1, njob, thread_no );
-               commongappick_record( 1, nogap+i, gapmap[i] );
                sprintf( com, "%s/_mccaskillinorg", dirname );
                infp = fopen( com, "w" );
 //             fprintf( infp, ">in\n%s\n", nogap[i] );
@@ -124,7 +126,10 @@ static void *athread( void *arg )
 
                sprintf( com, "tr -d '\\r' < %s/_mccaskillinorg > %s/_mccaskillin", dirname, dirname );
                system( com ); // for cygwin, wakaran
-               sprintf( com, "cd %s; %s/mxscarnamod -m -writebpp  _mccaskillin > _mccaskillout 2>_dum", dirname, whereismccaskillmea );
+               if( alg == 'G' )
+                       sprintf( com, "cd %s; %s/dafs --mafft-out _mccaskillout _mccaskillin > _dum1 2>_dum", dirname, whereismccaskillmea );
+               else
+                       sprintf( com, "cd %s; %s/mxscarnamod -m -writebpp  _mccaskillin > _mccaskillout 2>_dum", dirname, whereismccaskillmea );
                res = system( com );
 
                if( res )
@@ -149,6 +154,7 @@ static void *athread( void *arg )
        }
        free( dirname );
        free( com );
+       return( NULL );
 }
 #endif
 
@@ -161,6 +167,7 @@ void arguments( int argc, char *argv[] )
        kimuraR = NOTSPECIFIED;
        pamN = NOTSPECIFIED;
        whereismccaskillmea = NULL;
+       alg = 's';
 
     while( --argc > 0 && (*++argv)[0] == '-' )
        {
@@ -179,10 +186,16 @@ void arguments( int argc, char *argv[] )
                                        --argc;
                                        goto nextoption;
                                case 'C':
-                                       nthread = atoi( *++argv );
+                                       nthread = myatoi( *++argv );
                                        fprintf( stderr, "nthread = %d\n", nthread );
                                        --argc; 
                                        goto nextoption;
+                               case 's':
+                                       alg = 's'; // use scarna; default
+                                       break;
+                               case 'G':
+                                       alg = 'G'; // use dafs, instead of scarna
+                                       break;
                 default:
                     fprintf( stderr, "illegal option %c\n", c );
                     argc = 0;
@@ -216,7 +229,7 @@ int main( int argc, char *argv[] )
        RNApair *pairprobpt;
        RNApair *pt;
        int *alnpairnum;
-       float prob;
+       double prob;
        int adpos;
 
        arguments( argc, argv );
@@ -285,7 +298,10 @@ int main( int argc, char *argv[] )
 
        constants( njob, seq );
 
-       fprintf( stderr, "Running mxscarna with the mccaskill_mea mode.\n" );
+       if( alg == 'G' )
+               fprintf( stderr, "Running DAFS (Sato et al. 2012; http://www.ncrna.org/).\n" );
+       else
+               fprintf( stderr, "Running mxscarna with the mccaskill_mea mode.\n" );
 #ifdef enablemultithread
        if( nthread > 0 )
        {
@@ -321,6 +337,9 @@ int main( int argc, char *argv[] )
                }
                pthread_mutex_destroy( &mutex );
 
+               free( handle );
+               free( targ );
+
 
                for( i=0; i<njob; i++ )
                {
@@ -335,6 +354,12 @@ int main( int argc, char *argv[] )
                {
                        fprintf( stderr, "%d / %d\n", i+1, njob );
                        commongappick_record( 1, nogap+i, gapmap[i] );
+                       if( strlen( nogap[i] ) == 0 ) 
+                       {
+                               fprintf( stdout, ">%d\n", i );
+                               continue;
+                       }
+
                        infp = fopen( "_mccaskillinorg", "w" );
 //                     fprintf( infp, ">in\n%s\n", nogap[i] );
                        fprintf( infp, ">in\n" );
@@ -342,7 +367,10 @@ int main( int argc, char *argv[] )
                        fclose( infp );
        
                        system( "tr -d '\\r' < _mccaskillinorg > _mccaskillin" ); // for cygwin, wakaran
-                       sprintf( com, "env PATH=%s mxscarnamod -m -writebpp  _mccaskillin > _mccaskillout 2>_dum", whereismccaskillmea );
+                       if( alg == 'G' )
+                               sprintf( com, "env PATH=%s dafs --mafft-out _mccaskillout _mccaskillin > _dum1 2>_dum", whereismccaskillmea );
+                       else
+                               sprintf( com, "env PATH=%s mxscarnamod -m -writebpp  _mccaskillin > _mccaskillout 2>_dum", whereismccaskillmea );
                        res = system( com );
        
                        if( res )
@@ -403,6 +431,13 @@ int main( int argc, char *argv[] )
        for( j=0; j<nlenmax; j++ ) free( alnpairprob[j] );
        free( alnpairprob );
        free( alnpairnum );
+       free( order );
+       free( nlen );
+       FreeCharMtx( seq );
+       FreeCharMtx( nogap );
+       FreeCharMtx( name );
+       FreeIntMtx( gapmap );
+       freeconstants();
        fprintf( stderr, "%d thread(s)\n", nthread );
        return( 0 );
 
@@ -411,7 +446,7 @@ int main( int argc, char *argv[] )
 
        for( i=0; i<nlenmax; i++ ) for( pairprobpt=alnpairprob[i]; pairprobpt->bestpos!=-1; pairprobpt++ )
        {
-               pairprobpt->bestscore /= (float)njob;
+               pairprobpt->bestscore /= (double)njob;
                left = i;
                right = pairprobpt->bestpos;
                prob = pairprobpt->bestscore;