--- /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