binaries/src/mafft/core/dndpre.c

   1 #include "mltaln.h"
   2
   3 #define TEST 0
   4
   5 static int treeout = 0;
   6
   7 #ifdef enablemultithread
   8 typedef struct _jobtable
   9 {
  10     int i;
  11     int j;
  12 } Jobtable;
  13
  14 typedef struct _thread_arg
  15 {
  16         int njob;
  17         int thread_no;
  18         float *selfscore;
  19         double **mtx;
  20         char **seq;
  21         Jobtable *jobpospt;
  22         pthread_mutex_t *mutex;
  23 } thread_arg_t;
  24
  25 void *athread( void *arg )
  26 {
  27         thread_arg_t *targ = (thread_arg_t *)arg;
  28         int njob = targ->njob;
  29         int thread_no = targ->thread_no;
  30         float *selfscore = targ->selfscore;
  31         double **mtx = targ->mtx;
  32         char **seq = targ->seq;
  33         Jobtable *jobpospt = targ->jobpospt;
  34
  35         int i, j;
  36         float ssi, ssj, bunbo;
  37
  38         while( 1 )
  39         {
  40                 pthread_mutex_lock( targ->mutex );
  41                 j = jobpospt->j;
  42                 i = jobpospt->i;
  43                 j++;
  44                 if( j == njob )
  45                 {
  46                         fprintf( stderr, "%4d/%4d (thread %4d), dndpre\r", i+1, njob, thread_no );
  47                         i++;
  48                         j = i + 1;
  49                         if( i == njob-1 )
  50                         {
  51                                 pthread_mutex_unlock( targ->mutex );
  52                                 return( NULL );
  53                         }
  54                 }
  55                 jobpospt->j = j;
  56                 jobpospt->i = i;
  57                 pthread_mutex_unlock( targ->mutex );
  58
  59                 ssi = selfscore[i];
  60                 ssj = selfscore[j];
  61
  62                 bunbo = MIN( ssi, ssj );
  63                 if( bunbo == 0.0 )
  64                         mtx[i][j] = 1.0;
  65                 else
  66                         mtx[i][j] = 1.0 - (double)naivepairscore11( seq[i], seq[j], penalty ) / bunbo;
  67         }
  68 }
  69
  70 #endif
  71
  72 void arguments( int argc, char *argv[] )
  73 {
  74     int c;
  75
  76         nthread = 1;
  77         alg = 'X';
  78         fmodel = 0;
  79         treeout = 0;
  80         scoremtx = 1;
  81         nblosum = 62;
  82         dorp = NOTSPECIFIED;
  83         inputfile = NULL;
  84         ppenalty = NOTSPECIFIED; //?
  85         ppenalty_ex = NOTSPECIFIED; //?
  86         poffset = NOTSPECIFIED; //?
  87         kimuraR = NOTSPECIFIED;
  88         pamN = NOTSPECIFIED;
  89
  90     while( --argc > 0 && (*++argv)[0] == '-' )
  91         {
  92         while ( (c = *++argv[0]) )
  93                 {
  94             switch( c )
  95             {
  96                                 case 't':
  97                                         treeout = '1';
  98                                         break;
  99                                 case 'D':
 100                                         dorp = 'd';
 101                                         break;
 102                                 case 'a':
 103                                         fmodel = 1;
 104                                         break;
 105                                 case 'P':
 106                                         dorp = 'p';
 107                                         break;
 108                                 case 'i':
 109                                         inputfile = *++argv;
 110                                         fprintf( stderr, "inputfile = %s\n", inputfile );
 111                                         --argc;
 112                                         goto nextoption;
 113                                 case 'C':
 114                                         nthread = atoi( *++argv );
 115                                         fprintf( stderr, "nthread = %d\n", nthread );
 116                                         --argc;
 117                                         goto nextoption;
 118             }
 119                 }
 120                 nextoption:
 121                         ;
 122         }
 123     if( argc != 0 )
 124         {
 125                 fprintf( stderr, "options: Check source file !\n" );
 126                 exit( 1 );
 127         }
 128 }
 129
 130 int main( int argc, char **argv )
 131 {
 132         int i, j;
 133         char **seq;
 134         static char **name;
 135         static int nlen[M];
 136         float *selfscore;
 137         double **mtx;
 138         FILE *fp;
 139         FILE *infp;
 140         float ssi, ssj, bunbo;
 141
 142
 143         arguments( argc, argv );
 144 #ifndef enablemultithread
 145         nthread = 0;
 146 #endif
 147
 148         if( inputfile )
 149         {
 150                 infp = fopen( inputfile, "r" );
 151                 if( !infp )
 152                 {
 153                         fprintf( stderr, "Cannot open %s\n", inputfile );
 154                         exit( 1 );
 155                 }
 156         }
 157         else
 158                 infp = stdin;
 159
 160 #if 0
 161         PreRead( stdin, &njob, &nlenmax );
 162 #else
 163         getnumlen( infp );
 164 #endif
 165         rewind( infp );
 166
 167         seq = AllocateCharMtx( njob, nlenmax+1 );
 168         name = AllocateCharMtx( njob, B+1 );
 169         mtx = AllocateDoubleMtx( njob, njob );
 170         selfscore = AllocateFloatVec( njob );
 171
 172 #if 0
 173         FRead( stdin, name, nlen, seq );
 174 #else
 175         readData_pointer( infp, name, nlen, seq );
 176 #endif
 177         fclose( infp );
 178
 179         constants( njob, seq );
 180
 181 #if 0
 182         for( i=0; i<njob-1; i++ )
 183         {
 184                 fprintf( stderr, "%4d/%4d\r", i+1, njob );
 185                 for( j=i+1; j<njob; j++ )
 186                         mtx[i][j] = (double)substitution_hosei( seq[i], seq[j] );
 187 //                      fprintf( stderr, "i=%d,j=%d, l=%d &&&  %f\n", i, j, nlen[0], mtx[i][j] );
 188         }
 189 #else // 061003
 190         for( i=0; i<njob; i++ )
 191         {
 192                 selfscore[i] = (float)naivepairscore11( seq[i], seq[i], penalty );
 193
 194         }
 195 #ifdef enablemultithread
 196         if( nthread > 0 )
 197         {
 198                 thread_arg_t *targ;
 199                 Jobtable jobpos;
 200                 pthread_t *handle;
 201                 pthread_mutex_t mutex;
 202
 203                 jobpos.i = 0;
 204                 jobpos.j = 0;
 205
 206                 targ = calloc( nthread, sizeof( thread_arg_t ) );
 207                 handle = calloc( nthread, sizeof( pthread_t ) );
 208                 pthread_mutex_init( &mutex, NULL );
 209
 210                 for( i=0; i<nthread; i++ )
 211                 {
 212                         targ[i].thread_no = i;
 213                         targ[i].njob = njob;
 214                         targ[i].selfscore = selfscore;
 215                         targ[i].mtx = mtx;
 216                         targ[i].seq = seq;
 217                         targ[i].jobpospt = &jobpos;
 218                         targ[i].mutex = &mutex;
 219
 220                         pthread_create( handle+i, NULL, athread, (void *)(targ+i) );
 221                 }
 222
 223                 for( i=0; i<nthread; i++ )
 224                 {
 225                         pthread_join( handle[i], NULL );
 226                 }
 227                 pthread_mutex_destroy( &mutex );
 228         }
 229         else
 230 #endif
 231         {
 232                 for( i=0; i<njob-1; i++ )
 233                 {
 234                         ssi = selfscore[i];
 235                         fprintf( stderr, "%4d/%4d\r", i+1, njob );
 236                         for( j=i+1; j<njob; j++ )
 237                         {
 238                                 ssj = selfscore[j];
 239                                 bunbo = MIN( ssi, ssj );
 240                                 if( bunbo == 0.0 )
 241                                         mtx[i][j] = 1.0;
 242                                 else
 243                                         mtx[i][j] = 1.0 - (double)naivepairscore11( seq[i], seq[j], penalty ) / bunbo;
 244 //                                      mtx[i][j] = 1.0 - (double)naivepairscore11( seq[i], seq[j], penalty ) / MIN( selfscore[i], selfscore[j] );
 245 //                              fprintf( stderr, "i=%d,j=%d, l=%d### %f, score = %d\n", i, j, nlen[0], mtx[i][j], naivepairscore11( seq[i], seq[j], penalty )  );
 246                         }
 247                 }
 248         }
 249 #endif
 250
 251 #if TEST
 252         for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
 253                 fprintf( stdout, "i=%d, j=%d, mtx[][] = %f\n", i, j, mtx[i][j] );
 254 #endif
 255
 256         fp = fopen( "hat2", "w" );
 257         WriteHat2_pointer( fp, njob, name, mtx );
 258         fclose( fp );
 259 #if 0
 260         if( treeout )
 261         {
 262                 int ***topol;
 263                 double **len;
 264
 265                 topol = AllocateIntCub( njob, 2, njob );
 266                 len = AllocateDoubleMtx( njob, njob );
 267                 veryfastsupg_double_outtree( njob, mtx, topol, len );
 268         }
 269 #endif
 270         SHOWVERSION;
 271         exit( 0 );
 272 /*
 273         res = system( ALNDIR "/spgsdl < hat2"  );
 274         if( res ) exit( 1 );
 275         else exit( 0 );
 276 */
 277 }