+++ /dev/null
-#include "mltaln.h"
-
-#define DEBUG 0
-
-static char *whereismccaskillmea;
-
-#ifdef enablemultithread
-typedef struct _thread_arg
-{
- int thread_no;
- int njob;
- int *jobpospt;
- int **gapmap;
- char **nogap;
- int nlenmax;
- RNApair ***pairprob;
- pthread_mutex_t *mutex;
-} thread_arg_t;
-#endif
-
-void outmccaskill( FILE *fp, RNApair **pairprob, int length )
-{
- int i;
- RNApair *pt;
- for( i=0; i<length; i++ ) for( pt=pairprob[i]; pt->bestpos!=-1; pt++ )
- {
- if( pt->bestpos > i )
- fprintf( fp, "%d %d %50.40f\n", i, pt->bestpos, pt->bestscore );
- }
-}
-
-#if 1
-static void readrawmccaskill( FILE *fp, RNApair **pairprob, int length )
-{
- char gett[1000];
- int *pairnum;
- int i;
- int left, right;
- float prob;
-
- pairnum = (int *)calloc( length, sizeof( int ) );
- for( i=0; i<length; i++ ) pairnum[i] = 0;
-
- while( 1 )
- {
- fgets( gett, 999, fp );
- if( feof( fp ) ) break;
- if( gett[0] == '>' ) continue;
- sscanf( gett, "%d %d %f", &left, &right, &prob );
- if( prob < 0.01 ) continue; // mxscarna to mafft ryoho ni eikyou
-//fprintf( stderr, "gett = %s\n", gett );
-
- if( left != right && prob > 0.0 )
- {
- pairprob[left] = (RNApair *)realloc( pairprob[left], (pairnum[left]+2) * sizeof( RNApair ) );
- pairprob[left][pairnum[left]].bestscore = prob;
- pairprob[left][pairnum[left]].bestpos = right;
- pairnum[left]++;
- pairprob[left][pairnum[left]].bestscore = -1.0;
- pairprob[left][pairnum[left]].bestpos = -1;
-// fprintf( stderr, "%d-%d, %f\n", left, right, prob );
-
- pairprob[right] = (RNApair *)realloc( pairprob[right], (pairnum[right]+2) * sizeof( RNApair ) );
- pairprob[right][pairnum[right]].bestscore = prob;
- pairprob[right][pairnum[right]].bestpos = left;
- pairnum[right]++;
- pairprob[right][pairnum[right]].bestscore = -1.0;
- pairprob[right][pairnum[right]].bestpos = -1;
-// fprintf( stderr, "%d-%d, %f\n", right, left, prob );
- }
- }
- free( pairnum );
-}
-#endif
-
-#ifdef enablemultithread
-static void *athread( void *arg )
-{
- thread_arg_t *targ = (thread_arg_t *)arg;
- int thread_no = targ->thread_no;
- int njob = targ->njob;
- int *jobpospt = targ->jobpospt;
- int **gapmap = targ->gapmap;
- char **nogap = targ->nogap;
- int nlenmax = targ->nlenmax;
- RNApair ***pairprob = targ->pairprob;
-
- int i, res;
- FILE *infp;
- char *com;
- char *dirname;
-
- dirname = calloc( 100, sizeof( char ) );
- com = calloc( 1000, sizeof( char ) );
-
-
- while( 1 )
- {
- pthread_mutex_lock( targ->mutex );
- i = *jobpospt;
- if( i == njob )
- {
- pthread_mutex_unlock( targ->mutex );
- return( NULL );
- }
- *jobpospt = i+1;
- pthread_mutex_unlock( targ->mutex );
-
-
- sprintf( dirname, "_%d", i );
- sprintf( com, "rm -rf %s", dirname );
- system( com );
- sprintf( com, "mkdir %s", dirname );
- 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] );
- fprintf( infp, ">in\n" );
- write1seq( infp, nogap[i] );
- fclose( infp );
-
- 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 );
- res = system( com );
-
- if( res )
- {
- fprintf( stderr, "ERROR IN mccaskill_mea\n" );
- exit( 1 );
- }
-
- sprintf( com, "%s/_mccaskillout", dirname );
- infp = fopen( com, "r" );
- readrawmccaskill( infp, pairprob[i], nlenmax );
- fclose( infp );
-
- sprintf( com, "rm -rf %s > /dev/null 2>&1", dirname );
- if( system( com ) )
- {
- fprintf( stderr, "retrying to rmdir\n" );
-// nanosleep( 100000 );
- sleep( 1 );
- system( com );
- }
- }
- free( dirname );
- free( com );
-}
-#endif
-
-void arguments( int argc, char *argv[] )
-{
- int c;
- nthread = 1;
- inputfile = NULL;
- dorp = NOTSPECIFIED;
- kimuraR = NOTSPECIFIED;
- pamN = NOTSPECIFIED;
- whereismccaskillmea = NULL;
-
- while( --argc > 0 && (*++argv)[0] == '-' )
- {
- while ( (c = *++argv[0]) )
- {
- switch( c )
- {
- case 'i':
- inputfile = *++argv;
- fprintf( stderr, "inputfile = %s\n", inputfile );
- --argc;
- goto nextoption;
- case 'd':
- whereismccaskillmea = *++argv;
- fprintf( stderr, "whereismccaskillmea = %s\n", whereismccaskillmea );
- --argc;
- goto nextoption;
- case 'C':
- nthread = atoi( *++argv );
- fprintf( stderr, "nthread = %d\n", nthread );
- --argc;
- goto nextoption;
- default:
- fprintf( stderr, "illegal option %c\n", c );
- argc = 0;
- break;
- }
- }
- nextoption:
- ;
- }
- if( argc != 0 )
- {
- fprintf( stderr, "options: Check source file !\n" );
- exit( 1 );
- }
-}
-
-
-int main( int argc, char *argv[] )
-{
- static char com[10000];
- static int *nlen;
- int left, right;
- int res;
- static char **name, **seq, **nogap;
- static int **gapmap;
- static int *order;
- int i, j;
- FILE *infp;
- RNApair ***pairprob;
- RNApair **alnpairprob;
- RNApair *pairprobpt;
- RNApair *pt;
- int *alnpairnum;
- float prob;
- int adpos;
-
- arguments( argc, argv );
-#ifndef enablemultithread
- nthread = 0;
-#endif
-
- if( inputfile )
- {
- infp = fopen( inputfile, "r" );
- if( !infp )
- {
- fprintf( stderr, "Cannot open %s\n", inputfile );
- exit( 1 );
- }
- }
- else
- infp = stdin;
-
- if( !whereismccaskillmea )
- whereismccaskillmea = "";
-
- getnumlen( infp );
- rewind( infp );
-
- if( dorp != 'd' )
- {
- fprintf( stderr, "nuc only\n" );
- exit( 1 );
- }
-
- seq = AllocateCharMtx( njob, nlenmax*2+1 );
- nogap = AllocateCharMtx( njob, nlenmax*2+1 );
- gapmap = AllocateIntMtx( njob, nlenmax*2+1 );
- order = AllocateIntVec( njob );
- name = AllocateCharMtx( njob, B+1 );
- nlen = AllocateIntVec( njob );
- pairprob = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
- alnpairprob = (RNApair **)calloc( nlenmax, sizeof( RNApair * ) );
- alnpairnum = AllocateIntVec( nlenmax );
-
- for( i=0; i<nlenmax; i++ ) alnpairnum[i] = 0;
-
- readData_pointer( infp, name, nlen, seq );
- fclose( infp );
-
- for( i=0; i<njob; i++ )
- {
- pairprob[i] = (RNApair **)calloc( nlenmax, sizeof( RNApair * ) );
- for( j=0; j<nlenmax; j++ )
- {
- pairprob[i][j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
- pairprob[i][j][0].bestpos = -1;
- pairprob[i][j][0].bestscore = -1.0;
- }
- strcpy( nogap[i], seq[i] );
- order[i] = i;
- }
- for( j=0; j<nlenmax; j++ )
- {
- alnpairprob[j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
- alnpairprob[j][0].bestpos = -1;
- alnpairprob[j][0].bestscore = -1.0;
- }
-
-
- constants( njob, seq );
-
- fprintf( stderr, "Running mxscarna with the mccaskill_mea mode.\n" );
-#ifdef enablemultithread
- if( nthread > 0 )
- {
- int jobpos;
- pthread_t *handle;
- pthread_mutex_t mutex;
- thread_arg_t *targ;
- jobpos = 0;
-
- targ = calloc( nthread, sizeof( thread_arg_t ) );
- handle = calloc( nthread, sizeof( pthread_t ) );
- pthread_mutex_init( &mutex, NULL );
-
- for( i=0; i<nthread; i++ )
- {
- targ[i].thread_no = i;
- targ[i].njob = njob;
- targ[i].jobpospt = &jobpos;
- targ[i].gapmap = gapmap;
- targ[i].nogap = nogap;
- targ[i].nlenmax = nlenmax;
- targ[i].pairprob = pairprob;
- targ[i].mutex = &mutex;
-
-// athread( targ );
- pthread_create( handle+i, NULL, athread, (void *)(targ+i) );
-
- }
-
- for( i=0; i<nthread; i++ )
- {
- pthread_join( handle[i], NULL );
- }
- pthread_mutex_destroy( &mutex );
-
-
- for( i=0; i<njob; i++ )
- {
- fprintf( stdout, ">%d\n", i );
- outmccaskill( stdout, pairprob[i], nlenmax );
- }
- }
- else
-#endif
- {
- for( i=0; i<njob; i++ )
- {
- fprintf( stderr, "%d / %d\n", i+1, njob );
- commongappick_record( 1, nogap+i, gapmap[i] );
- infp = fopen( "_mccaskillinorg", "w" );
-// fprintf( infp, ">in\n%s\n", nogap[i] );
- fprintf( infp, ">in\n" );
- write1seq( infp, nogap[i] );
- fclose( infp );
-
- system( "tr -d '\\r' < _mccaskillinorg > _mccaskillin" ); // for cygwin, wakaran
- sprintf( com, "env PATH=%s mxscarnamod -m -writebpp _mccaskillin > _mccaskillout 2>_dum", whereismccaskillmea );
- res = system( com );
-
- if( res )
- {
- fprintf( stderr, "ERROR IN mccaskill_mea\n" );
- exit( 1 );
- }
-
- infp = fopen( "_mccaskillout", "r" );
- readrawmccaskill( infp, pairprob[i], nlenmax );
- fclose( infp );
- fprintf( stdout, ">%d\n", i );
- outmccaskill( stdout, pairprob[i], nlenmax );
- }
- }
-
- for( i=0; i<njob; i++ )
- {
- for( j=0; j<nlen[i]; j++ ) for( pairprobpt=pairprob[i][j]; pairprobpt->bestpos!=-1; pairprobpt++ )
- {
- left = gapmap[i][j];
- right = gapmap[i][pairprobpt->bestpos];
- prob = pairprobpt->bestscore;
-
- for( pt=alnpairprob[left]; pt->bestpos!=-1; pt++ )
- if( pt->bestpos == right ) break;
-
- if( pt->bestpos == -1 )
- {
- alnpairprob[left] = (RNApair *)realloc( alnpairprob[left], (alnpairnum[left]+2) * sizeof( RNApair ) );
- adpos = alnpairnum[left];
- alnpairnum[left]++;
- alnpairprob[left][adpos].bestscore = 0.0;
- alnpairprob[left][adpos].bestpos = right;
- alnpairprob[left][adpos+1].bestscore = -1.0;
- alnpairprob[left][adpos+1].bestpos = -1;
- pt = alnpairprob[left]+adpos;
- }
- else
- adpos = pt-alnpairprob[left];
-
- pt->bestscore += prob;
- if( pt->bestpos != right )
- {
- fprintf( stderr, "okashii!\n" );
- exit( 1 );
- }
-// fprintf( stderr, "adding %d-%d, %f\n", left, right, prob );
- }
- }
-
- for( i=0; i<njob; i++ )
- {
- for( j=0; j<nlenmax; j++ ) free( pairprob[i][j] );
- free( pairprob[i] );
- }
- free( pairprob );
- for( j=0; j<nlenmax; j++ ) free( alnpairprob[j] );
- free( alnpairprob );
- free( alnpairnum );
- fprintf( stderr, "%d thread(s)\n", nthread );
- return( 0 );
-
-#if 0
- fprintf( stdout, "result=\n" );
-
- for( i=0; i<nlenmax; i++ ) for( pairprobpt=alnpairprob[i]; pairprobpt->bestpos!=-1; pairprobpt++ )
- {
- pairprobpt->bestscore /= (float)njob;
- left = i;
- right = pairprobpt->bestpos;
- prob = pairprobpt->bestscore;
- fprintf( stdout, "%d-%d, %f\n", left, right, prob );
- }
-
- return( 0 );
-#endif
-}
git://source.jalview.org / jabaws.git / blobdiff