website/archive/binaries/mac/src/mafft/core/tbfast.c

   1 #include "mltaln.h"
   2
   3 #define DEBUG 0
   4 #define IODEBUG 0
   5 #define SCOREOUT 0
   6
   7 static int nadd;
   8 static int treein;
   9 static int topin;
  10 static int treeout;
  11 static int distout;
  12 static int noalign;
  13
  14 typedef struct _jobtable
  15 {
  16     int i;
  17     int j;
  18 } Jobtable;
  19
  20 #ifdef enablemultithread
  21 typedef struct _distancematrixthread_arg
  22 {
  23         int njob;
  24         int thread_no;
  25         float *selfscore;
  26         float **iscore;
  27         char **seq;
  28         Jobtable *jobpospt;
  29         pthread_mutex_t *mutex;
  30 } distancematrixthread_arg_t;
  31
  32 typedef struct _treebasethread_arg
  33 {
  34         int thread_no;
  35         int *nrunpt;
  36         int njob;
  37         int *nlen;
  38         int *jobpospt;
  39         int ***topol;
  40         Treedep *dep;
  41         char **aseq;
  42         double *effarr;
  43         int *alloclenpt;
  44         LocalHom **localhomtable;
  45         RNApair ***singlerna;
  46         double *effarr_kozo;
  47         int *fftlog;
  48         pthread_mutex_t *mutex;
  49         pthread_cond_t *treecond;
  50 } treebasethread_arg_t;
  51 #endif
  52
  53 void arguments( int argc, char *argv[] )
  54 {
  55     int c;
  56
  57         nthread = 1;
  58         outnumber = 0;
  59         scoreout = 0;
  60         treein = 0;
  61         topin = 0;
  62         rnaprediction = 'm';
  63         rnakozo = 0;
  64         nevermemsave = 0;
  65         inputfile = NULL;
  66         addfile = NULL;
  67         addprofile = 1;
  68         fftkeika = 0;
  69         constraint = 0;
  70         nblosum = 62;
  71         fmodel = 0;
  72         calledByXced = 0;
  73         devide = 0;
  74         use_fft = 0; // chuui
  75         force_fft = 0;
  76         fftscore = 1;
  77         fftRepeatStop = 0;
  78         fftNoAnchStop = 0;
  79     weight = 3;
  80     utree = 1;
  81         tbutree = 1;
  82     refine = 0;
  83     check = 1;
  84     cut = 0.0;
  85     disp = 0;
  86     outgap = 1;
  87     alg = 'A';
  88     mix = 0;
  89         tbitr = 0;
  90         scmtd = 5;
  91         tbweight = 0;
  92         tbrweight = 3;
  93         checkC = 0;
  94         treemethod = 'X';
  95         contin = 0;
  96         scoremtx = 1;
  97         kobetsubunkatsu = 0;
  98         dorp = NOTSPECIFIED;
  99         ppenalty = NOTSPECIFIED;
 100         ppenalty_ex = NOTSPECIFIED;
 101         poffset = NOTSPECIFIED;
 102         kimuraR = NOTSPECIFIED;
 103         pamN = NOTSPECIFIED;
 104         geta2 = GETA2;
 105         fftWinSize = NOTSPECIFIED;
 106         fftThreshold = NOTSPECIFIED;
 107         RNAppenalty = NOTSPECIFIED;
 108         RNAppenalty_ex = NOTSPECIFIED;
 109         RNApthr = NOTSPECIFIED;
 110         TMorJTT = JTT;
 111         consweight_multi = 1.0;
 112         consweight_rna = 0.0;
 113
 114     while( --argc > 0 && (*++argv)[0] == '-' )
 115         {
 116         while ( ( c = *++argv[0] ) )
 117                 {
 118             switch( c )
 119             {
 120                                 case 'i':
 121                                         inputfile = *++argv;
 122                                         fprintf( stderr, "inputfile = %s\n", inputfile );
 123                                         --argc;
 124                     goto nextoption;
 125                                 case 'I':
 126                                         nadd = atoi( *++argv );
 127                                         fprintf( stderr, "nadd = %d\n", nadd );
 128                                         --argc;
 129                                         goto nextoption;
 130                                 case 'e':
 131                                         RNApthr = (int)( atof( *++argv ) * 1000 - 0.5 );
 132                                         --argc;
 133                                         goto nextoption;
 134                                 case 'o':
 135                                         RNAppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
 136                                         --argc;
 137                                         goto nextoption;
 138                                 case 'f':
 139                                         ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
 140 //                                      fprintf( stderr, "ppenalty = %d\n", ppenalty );
 141                                         --argc;
 142                                         goto nextoption;
 143                                 case 'g':
 144                                         ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
 145                                         fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
 146                                         --argc;
 147                                         goto nextoption;
 148                                 case 'h':
 149                                         poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
 150 //                                      fprintf( stderr, "poffset = %d\n", poffset );
 151                                         --argc;
 152                                         goto nextoption;
 153                                 case 'k':
 154                                         kimuraR = atoi( *++argv );
 155                                         fprintf( stderr, "kappa = %d\n", kimuraR );
 156                                         --argc;
 157                                         goto nextoption;
 158                                 case 'b':
 159                                         nblosum = atoi( *++argv );
 160                                         scoremtx = 1;
 161                                         fprintf( stderr, "blosum %d / kimura 200\n", nblosum );
 162                                         --argc;
 163                                         goto nextoption;
 164                                 case 'j':
 165                                         pamN = atoi( *++argv );
 166                                         scoremtx = 0;
 167                                         TMorJTT = JTT;
 168                                         fprintf( stderr, "jtt/kimura %d\n", pamN );
 169                                         --argc;
 170                                         goto nextoption;
 171                                 case 'm':
 172                                         pamN = atoi( *++argv );
 173                                         scoremtx = 0;
 174                                         TMorJTT = TM;
 175                                         fprintf( stderr, "tm %d\n", pamN );
 176                                         --argc;
 177                                         goto nextoption;
 178                                 case 'l':
 179                                         fastathreshold = atof( *++argv );
 180                                         constraint = 2;
 181                                         --argc;
 182                                         goto nextoption;
 183                                 case 'r':
 184                                         consweight_rna = atof( *++argv );
 185                                         rnakozo = 1;
 186                                         --argc;
 187                                         goto nextoption;
 188                                 case 'c':
 189                                         consweight_multi = atof( *++argv );
 190                                         --argc;
 191                                         goto nextoption;
 192                                 case 'C':
 193                                         nthread = atoi( *++argv );
 194                                         fprintf( stderr, "nthread = %d\n", nthread );
 195                                         --argc;
 196                                         goto nextoption;
 197                                 case 'R':
 198                                         rnaprediction = 'r';
 199                                         break;
 200                                 case 's':
 201                                         RNAscoremtx = 'r';
 202                                         break;
 203 #if 1
 204                                 case 'a':
 205                                         fmodel = 1;
 206                                         break;
 207 #endif
 208                                 case 'K':
 209                                         addprofile = 0;
 210                                         break;
 211                                 case 'y':
 212                                         distout = 1;
 213                                         break;
 214                                 case 't':
 215                                         treeout = 1;
 216                                         break;
 217                                 case 'T':
 218                                         noalign = 1;
 219                                         break;
 220                                 case 'D':
 221                                         dorp = 'd';
 222                                         break;
 223                                 case 'P':
 224                                         dorp = 'p';
 225                                         break;
 226 #if 1
 227                                 case 'O':
 228                                         outgap = 0;
 229                                         break;
 230 #else
 231                                 case 'O':
 232                                         fftNoAnchStop = 1;
 233                                         break;
 234 #endif
 235                                 case 'S':
 236                                         scoreout = 1;
 237                                         break;
 238 #if 0
 239                                 case 'e':
 240                                         fftscore = 0;
 241                                         break;
 242                                 case 'r':
 243                                         fmodel = -1;
 244                                         break;
 245                                 case 'R':
 246                                         fftRepeatStop = 1;
 247                                         break;
 248                                 case 's':
 249                                         treemethod = 's';
 250                                         break;
 251 #endif
 252                                 case 'X':
 253                                         treemethod = 'X';
 254                                         break;
 255                                 case 'E':
 256                                         treemethod = 'E';
 257                                         break;
 258                                 case 'q':
 259                                         treemethod = 'q';
 260                                         break;
 261                                 case 'n' :
 262                                         outnumber = 1;
 263                                         break;
 264 #if 0
 265                                 case 'a':
 266                                         alg = 'a';
 267                                         break;
 268 #endif
 269                                 case 'Q':
 270                                         alg = 'Q';
 271                                         break;
 272                                 case 'H':
 273                                         alg = 'H';
 274                                         break;
 275                                 case 'A':
 276                                         alg = 'A';
 277                                         break;
 278                                 case 'M':
 279                                         alg = 'M';
 280                                         break;
 281                                 case 'N':
 282                                         nevermemsave = 1;
 283                                         break;
 284                                 case 'B':
 285                                         break;
 286                                 case 'F':
 287                                         use_fft = 1;
 288                                         break;
 289                                 case 'G':
 290                                         force_fft = 1;
 291                                         use_fft = 1;
 292                                         break;
 293                                 case 'U':
 294                                         treein = 1;
 295                                         break;
 296                                 case 'V':
 297                                         topin = 1;
 298                                         break;
 299                                 case 'u':
 300                                         tbrweight = 0;
 301                                         weight = 0;
 302                                         break;
 303                                 case 'v':
 304                                         tbrweight = 3;
 305                                         break;
 306                                 case 'd':
 307                                         disp = 1;
 308                                         break;
 309 /* Modified 01/08/27, default: user tree */
 310                                 case 'J':
 311                                         tbutree = 0;
 312                                         break;
 313 /* modification end. */
 314                                 case 'z':
 315                                         fftThreshold = atoi( *++argv );
 316                                         --argc;
 317                                         goto nextoption;
 318                                 case 'w':
 319                                         fftWinSize = atoi( *++argv );
 320                                         --argc;
 321                                         goto nextoption;
 322                                 case 'Z':
 323                                         checkC = 1;
 324                                         break;
 325                 default:
 326                     fprintf( stderr, "illegal option %c\n", c );
 327                     argc = 0;
 328                     break;
 329             }
 330                 }
 331                 nextoption:
 332                         ;
 333         }
 334     if( argc == 1 )
 335     {
 336         cut = atof( (*argv) );
 337         argc--;
 338     }
 339     if( argc != 0 )
 340     {
 341         fprintf( stderr, "options: Check source file !\n" );
 342         exit( 1 );
 343     }
 344         if( tbitr == 1 && outgap == 0 )
 345         {
 346                 fprintf( stderr, "conflicting options : o, m or u\n" );
 347                 exit( 1 );
 348         }
 349         if( alg == 'C' && outgap == 0 )
 350         {
 351                 fprintf( stderr, "conflicting options : C, o\n" );
 352                 exit( 1 );
 353         }
 354 }
 355
 356 #if 0
 357 static void *distancematrixthread2( void *arg )
 358 {
 359         distancematrixthread_arg_t *targ = (distancematrixthread_arg_t *)arg;
 360         int njob = targ->njob;
 361         int thread_no = targ->thread_no;
 362         float *selfscore = targ->selfscore;
 363         float **iscore = targ->iscore;
 364         char **seq = targ->seq;
 365         Jobtable *jobpospt = targ->jobpospt;
 366
 367         float ssi, ssj, bunbo;
 368         int i, j;
 369
 370         while( 1 )
 371         {
 372                 pthread_mutex_lock( targ->mutex );
 373                 i = jobpospt->i;
 374                 i++;
 375                 if( i == njob-1 )
 376                 {
 377                         pthread_mutex_unlock( targ->mutex );
 378                         return( NULL );
 379                 }
 380                 jobpospt->i = i;
 381                 pthread_mutex_unlock( targ->mutex );
 382
 383                 ssi = selfscore[i];
 384                 if( i % 10 == 0 ) fprintf( stderr, "\r% 5d / %d (thread %4d)", i, njob, thread_no );
 385                 for( j=i+1; j<njob; j++)
 386                 {
 387                         ssj = selfscore[j];
 388                         bunbo = MIN( ssi, ssj );
 389                         if( bunbo == 0.0 )
 390                                 iscore[i][j-i] = 1.0;
 391                         else
 392                                 iscore[i][j-i] = 1.0 - naivepairscore11( seq[i], seq[j], penalty ) / bunbo;
 393                 }
 394         }
 395 }
 396 #endif
 397
 398 #ifdef enablemultithread
 399 static void *distancematrixthread( void *arg )
 400 {
 401         distancematrixthread_arg_t *targ = (distancematrixthread_arg_t *)arg;
 402         int njob = targ->njob;
 403         int thread_no = targ->thread_no;
 404         float *selfscore = targ->selfscore;
 405         float **iscore = targ->iscore;
 406         char **seq = targ->seq;
 407         Jobtable *jobpospt = targ->jobpospt;
 408
 409         float ssi, ssj, bunbo;
 410         int i, j;
 411
 412         while( 1 )
 413         {
 414                 pthread_mutex_lock( targ->mutex );
 415                 j = jobpospt->j;
 416                 i = jobpospt->i;
 417                 j++;
 418                 if( j == njob )
 419                 {
 420                         i++;
 421                         j = i + 1;
 422                         if( i == njob-1 )
 423                         {
 424                                 pthread_mutex_unlock( targ->mutex );
 425                                 return( NULL );
 426                         }
 427                 }
 428                 jobpospt->j = j;
 429                 jobpospt->i = i;
 430                 pthread_mutex_unlock( targ->mutex );
 431
 432
 433                 if( j==i+1 && i % 10 == 0 ) fprintf( stderr, "\r% 5d / %d (thread %4d)", i, njob, thread_no );
 434                 ssi = selfscore[i];
 435                 ssj = selfscore[j];
 436                 bunbo = MIN( ssi, ssj );
 437                 if( bunbo == 0.0 )
 438                         iscore[i][j-i] = 1.0;
 439                 else
 440                         iscore[i][j-i] = 1.0 - naivepairscore11( seq[i], seq[j], penalty ) / bunbo;
 441         }
 442 }
 443
 444 static void *treebasethread( void *arg )
 445 {
 446         treebasethread_arg_t *targ = (treebasethread_arg_t *)arg;
 447         int *nrunpt = targ->nrunpt;
 448         int thread_no = targ->thread_no;
 449         int njob = targ->njob;
 450         int *nlen = targ->nlen;
 451         int *jobpospt = targ->jobpospt;
 452         int ***topol = targ->topol;
 453         Treedep *dep = targ->dep;
 454         char **aseq = targ->aseq;
 455         double *effarr = targ->effarr;
 456         int *alloclen = targ->alloclenpt;
 457         LocalHom **localhomtable = targ->localhomtable;
 458         RNApair ***singlerna = targ->singlerna;
 459         double *effarr_kozo = targ->effarr_kozo;
 460         int *fftlog = targ->fftlog;
 461
 462         char **mseq1, **mseq2;
 463         char **localcopy;
 464         int i, j, l;
 465         int len1, len2;
 466         int clus1, clus2;
 467         float pscore;
 468         char *indication1, *indication2;
 469         double *effarr1 = NULL;
 470         double *effarr2 = NULL;
 471         double *effarr1_kozo = NULL;
 472         double *effarr2_kozo = NULL;
 473         LocalHom ***localhomshrink = NULL;
 474         int m1, m2;
 475         float dumfl = 0.0;
 476         int ffttry;
 477         RNApair ***grouprna1, ***grouprna2;
 478
 479
 480         mseq1 = AllocateCharMtx( njob, 0 );
 481         mseq2 = AllocateCharMtx( njob, 0 );
 482         localcopy = calloc( njob, sizeof( char * ) );
 483
 484         if( rnakozo && rnaprediction == 'm' )
 485         {
 486                 grouprna1 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
 487                 grouprna2 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
 488         }
 489         else
 490         {
 491                 grouprna1 = grouprna2 = NULL;
 492         }
 493
 494         effarr1 = AllocateDoubleVec( njob );
 495         effarr2 = AllocateDoubleVec( njob );
 496         indication1 = AllocateCharVec( 150 );
 497         indication2 = AllocateCharVec( 150 );
 498 #if 0
 499 #else
 500         if( constraint )
 501         {
 502                 localhomshrink = (LocalHom ***)calloc( njob, sizeof( LocalHom ** ) );
 503                 for( i=0; i<njob; i++)
 504                         localhomshrink[i] = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
 505         }
 506 #endif
 507         effarr1_kozo = AllocateDoubleVec( njob ); //tsuneni allocate sareru.
 508         effarr2_kozo = AllocateDoubleVec( njob ); //tsuneni allocate sareru.
 509         for( i=0; i<njob; i++ ) effarr1_kozo[i] = 0.0;
 510         for( i=0; i<njob; i++ ) effarr2_kozo[i] = 0.0;
 511
 512
 513 #if 0
 514 #endif
 515
 516 #if 0
 517         for( i=0; i<njob; i++ )
 518                 fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
 519 #endif
 520
 521 #if 0 //-> main thread
 522         if( constraint )
 523                 calcimportance( njob, effarr, aseq, localhomtable );
 524 #endif
 525
 526
 527 //      writePre( njob, name, nlen, aseq, 0 );
 528
 529
 530 //      for( l=0; l<njob-1; l++ )
 531         while( 1 )
 532         {
 533
 534                 pthread_mutex_lock( targ->mutex );
 535                 l = *jobpospt;
 536                 if( l == njob-1 )
 537                 {
 538                         pthread_mutex_unlock( targ->mutex );
 539                         if( commonIP ) FreeIntMtx( commonIP );
 540                         commonIP = NULL;
 541                         Falign( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
 542                         A__align( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
 543                         free( mseq1 );
 544                         free( mseq2 );
 545                         free( localcopy );
 546                         free( effarr1 );
 547                         free( effarr2 );
 548                         free( effarr1_kozo );
 549                         free( effarr2_kozo );
 550                         free( indication1 );
 551                         free( indication2 );
 552                         if( constraint )
 553                         {
 554                                 for( i=0; i<njob; i++)
 555                                         free( localhomshrink[i] );
 556                                 free( localhomshrink );
 557                         }
 558                         return( NULL );
 559                 }
 560                 *jobpospt = l+1;
 561
 562                 if( dep[l].child0 != -1 )
 563                 {
 564                         while( dep[dep[l].child0].done == 0 )
 565                                 pthread_cond_wait( targ->treecond, targ->mutex );
 566                 }
 567                 if( dep[l].child1 != -1 )
 568                 {
 569                         while( dep[dep[l].child1].done == 0 )
 570                                 pthread_cond_wait( targ->treecond, targ->mutex );
 571                 }
 572 //              while( *nrunpt >= nthread )
 573 //                      pthread_cond_wait( targ->treecond, targ->mutex );
 574                 (*nrunpt)++;
 575
 576 //              pthread_mutex_unlock( targ->mutex );
 577
 578
 579                 m1 = topol[l][0][0];
 580                 m2 = topol[l][1][0];
 581
 582 //              pthread_mutex_lock( targ->mutex );
 583
 584         len1 = strlen( aseq[m1] );
 585         len2 = strlen( aseq[m2] );
 586         if( *alloclen <= len1 + len2 )
 587         {
 588                         fprintf( stderr, "\nReallocating (by thread %d) ..", thread_no );
 589                         *alloclen = ( len1 + len2 ) + 1000;
 590                         ReallocateCharMtx( aseq, njob, *alloclen + 10  );
 591                         fprintf( stderr, "done. *alloclen = %d\n", *alloclen );
 592                 }
 593
 594                 for( i=0; (j=topol[l][0][i])!=-1; i++ )
 595                 {
 596                         localcopy[j] = calloc( *alloclen, sizeof( char ) );
 597                         strcpy( localcopy[j], aseq[j] );
 598                 }
 599                 for( i=0; (j=topol[l][1][i])!=-1; i++ )
 600                 {
 601                         localcopy[j] = calloc( *alloclen, sizeof( char ) );
 602                         strcpy( localcopy[j], aseq[j] );
 603                 }
 604
 605                 pthread_mutex_unlock( targ->mutex );
 606
 607                 if( effarr_kozo )
 608                 {
 609                         clus1 = fastconjuction_noname_kozo( topol[l][0], localcopy, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
 610                         clus2 = fastconjuction_noname_kozo( topol[l][1], localcopy, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
 611                 }
 612                 else
 613                 {
 614                         clus1 = fastconjuction_noname( topol[l][0], localcopy, mseq1, effarr1, effarr, indication1 );
 615                         clus2 = fastconjuction_noname( topol[l][1], localcopy, mseq2, effarr2, effarr, indication2 );
 616                 }
 617
 618
 619
 620 #if 1
 621                 fprintf( stderr, "\rSTEP % 5d /%d (thread %4d) ", l+1, njob-1, thread_no );
 622 #else
 623                 fprintf( stderr, "STEP %d /%d (thread %d) \n", l+1, njob-1, thread_no );
 624                 fprintf( stderr, "group1 = %.66s", indication1 );
 625                 if( strlen( indication1 ) > 66 ) fprintf( stderr, "..." );
 626                 fprintf( stderr, ", child1 = %d\n", dep[l].child0 );
 627                 fprintf( stderr, "group2 = %.66s", indication2 );
 628                 if( strlen( indication2 ) > 66 ) fprintf( stderr, "..." );
 629                 fprintf( stderr, ", child2 = %d\n", dep[l].child1 );
 630
 631                 fprintf( stderr, "Group1's lengths = " );
 632                 for( i=0; i<clus1; i++ ) fprintf( stderr, "%d ", strlen( mseq1[i] ) );
 633                 fprintf( stderr, "\n" );
 634                 fprintf( stderr, "Group2's lengths = " );
 635                 for( i=0; i<clus2; i++ ) fprintf( stderr, "%d ", strlen( mseq2[i] ) );
 636                 fprintf( stderr, "\n" );
 637
 638 #endif
 639
 640
 641
 642 //              for( i=0; i<clus1; i++ ) fprintf( stderr, "## STEP%d-eff for mseq1-%d %f\n", l+1, i, effarr1[i] );
 643
 644                 if( constraint )
 645                 {
 646                         fastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
 647 //                      msfastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
 648 //                      fprintf( stderr, "localhomshrink =\n" );
 649 //                      outlocalhompt( localhomshrink, clus1, clus2 );
 650 //                      weightimportance4( clus1, clus2, effarr1, effarr2, localhomshrink );
 651 //                      fprintf( stderr, "after weight =\n" );
 652 //                      outlocalhompt( localhomshrink, clus1, clus2 );
 653                 }
 654                 if( rnakozo && rnaprediction == 'm' )
 655                 {
 656                         makegrouprna( grouprna1, singlerna, topol[l][0] );
 657                         makegrouprna( grouprna2, singlerna, topol[l][1] );
 658                 }
 659
 660
 661 /*
 662                 fprintf( stderr, "before align all\n" );
 663                 display( localcopy, njob );
 664                 fprintf( stderr, "\n" );
 665                 fprintf( stderr, "before align 1 %s \n", indication1 );
 666                 display( mseq1, clus1 );
 667                 fprintf( stderr, "\n" );
 668                 fprintf( stderr, "before align 2 %s \n", indication2 );
 669                 display( mseq2, clus2 );
 670                 fprintf( stderr, "\n" );
 671 */
 672
 673                 if( !nevermemsave && ( constraint != 2  && alg != 'M'  && ( len1 > 30000 || len2 > 30000 ) ) )
 674                 {
 675                         fprintf( stderr, "\nlen1=%d, len2=%d, Switching to the memsave mode.\n", len1, len2 );
 676                         alg = 'M';
 677                         if( commonIP ) FreeIntMtx( commonIP );
 678                         commonAlloc1 = 0;
 679                         commonAlloc2 = 0;
 680                 }
 681
 682
 683 //              if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 );
 684                 if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 1000 && clus2 < 1000 );
 685                 else                                               ffttry = 0;
 686 //              ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 5000 && clus2 < 5000 ); // v6.708
 687 //              fprintf( stderr, "f=%d, len1/fftlog[m1]=%f, clus1=%d, len2/fftlog[m2]=%f, clus2=%d\n", ffttry, (float)len1/fftlog[m1], clus1, (float)len2/fftlog[m2], clus2 );
 688 //              fprintf( stderr, "f=%d, clus1=%d, fftlog[m1]=%d, clus2=%d, fftlog[m2]=%d\n", ffttry, clus1, fftlog[m1], clus2, fftlog[m2] );
 689                 if( constraint == 2 )
 690                 {
 691                         if( alg == 'M' )
 692                         {
 693                                 fprintf( stderr, "\n\nMemory saving mode is not supported.\n\n" );
 694                                 exit( 1 );
 695                         }
 696                         fprintf( stderr, "c" );
 697                         if( alg == 'A' )
 698                         {
 699                                 imp_match_init_strict( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 700                                 if( rnakozo ) imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
 701                                 pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
 702                         }
 703                         else if( alg == 'H' )
 704                         {
 705                                 imp_match_init_strictH( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 706                                 pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
 707                         }
 708                         else if( alg == 'Q' )
 709                         {
 710                                 imp_match_init_strictQ( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 711                                 if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
 712                                 pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
 713                         }
 714                         else if( alg == 'R' )
 715                         {
 716                                 imp_match_init_strictR( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 717                                 pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
 718                         }
 719                 }
 720                 else if( force_fft || ( use_fft && ffttry ) )
 721                 {
 722                         fprintf( stderr, "f" );
 723                         if( alg == 'M' )
 724                         {
 725                                 fprintf( stderr, "m" );
 726                                 pscore = Falign_udpari_long( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1 );
 727                         }
 728                         else
 729                                 pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL );
 730                 }
 731                 else
 732                 {
 733                         fprintf( stderr, "d" );
 734                         fftlog[m1] = 0;
 735                         switch( alg )
 736                         {
 737                                 case( 'a' ):
 738                                         pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen );
 739                                         break;
 740                                 case( 'M' ):
 741                                         fprintf( stderr, "m" );
 742                                         pscore = MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
 743                                         break;
 744                                 case( 'A' ):
 745                                         pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
 746                                         break;
 747                                 case( 'Q' ):
 748                                         pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
 749                                         break;
 750                                 case( 'R' ):
 751                                         pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
 752                                         break;
 753                                 case( 'H' ):
 754                                         pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
 755                                         break;
 756                                 default:
 757                                         ErrorExit( "ERROR IN SOURCE FILE" );
 758                         }
 759                 }
 760
 761                 nlen[m1] = 0.5 * ( nlen[m1] + nlen[m2] );
 762
 763 #if SCOREOUT
 764                 fprintf( stderr, "score = %10.2f\n", pscore );
 765 #endif
 766
 767 /*
 768                 fprintf( stderr, "after align 1 %s \n", indication1 );
 769                 display( mseq1, clus1 );
 770                 fprintf( stderr, "\n" );
 771                 fprintf( stderr, "after align 2 %s \n", indication2 );
 772                 display( mseq2, clus2 );
 773                 fprintf( stderr, "\n" );
 774 */
 775
 776 //              writePre( njob, name, nlen, localcopy, 0 );
 777
 778                 if( disp ) display( localcopy, njob );
 779
 780                 pthread_mutex_lock( targ->mutex );
 781                 dep[l].done = 1;
 782                 (*nrunpt)--;
 783                 pthread_cond_broadcast( targ->treecond );
 784
 785 //              pthread_mutex_unlock( targ->mutex );
 786 //              pthread_mutex_lock( targ->mutex );
 787
 788                 for( i=0; (j=topol[l][0][i])!=-1; i++ )
 789                         strcpy( aseq[j], localcopy[j] );
 790                 for( i=0; (j=topol[l][1][i])!=-1; i++ )
 791                         strcpy( aseq[j], localcopy[j] );
 792                 pthread_mutex_unlock( targ->mutex );
 793
 794                 for( i=0; (j=topol[l][0][i])!=-1; i++ )
 795                         free( localcopy[j] );
 796                 for( i=0; (j=topol[l][1][i])!=-1; i++ )
 797                         free( localcopy[j] );
 798                 free( topol[l][0] );
 799                 free( topol[l][1] );
 800                 free( topol[l] );
 801
 802         }
 803 }
 804 #endif
 805
 806 void treebase( int *nlen, char **aseq, int nadd, char *mergeoralign, char **mseq1, char **mseq2, int ***topol, double *effarr, int *alloclen, LocalHom **localhomtable, RNApair ***singlerna, double *effarr_kozo )
 807 {
 808         int i, l, m;
 809         int len1nocommongap, len2nocommongap;
 810         int len1, len2;
 811         int clus1, clus2;
 812         float pscore, tscore;
 813         static char *indication1, *indication2;
 814         static double *effarr1 = NULL;
 815         static double *effarr2 = NULL;
 816         static double *effarr1_kozo = NULL;
 817         static double *effarr2_kozo = NULL;
 818         static LocalHom ***localhomshrink = NULL;
 819         static int *fftlog;
 820         int m1, m2;
 821         static int *gaplen;
 822         static int *gapmap;
 823         static int *alreadyaligned;
 824         float dumfl = 0.0;
 825         int ffttry;
 826         RNApair ***grouprna1, ***grouprna2;
 827
 828         if( rnakozo && rnaprediction == 'm' )
 829         {
 830                 grouprna1 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
 831                 grouprna2 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
 832         }
 833         else
 834         {
 835                 grouprna1 = grouprna2 = NULL;
 836         }
 837
 838         if( effarr1 == NULL )
 839         {
 840                 fftlog = AllocateIntVec( njob );
 841                 effarr1 = AllocateDoubleVec( njob );
 842                 effarr2 = AllocateDoubleVec( njob );
 843                 indication1 = AllocateCharVec( 150 );
 844                 indication2 = AllocateCharVec( 150 );
 845                 gaplen = AllocateIntVec( *alloclen+10 );
 846                 gapmap = AllocateIntVec( *alloclen+10 );
 847                 alreadyaligned = AllocateIntVec( njob );
 848 #if 0
 849 #else
 850                 if( constraint )
 851                 {
 852                         localhomshrink = (LocalHom ***)calloc( njob, sizeof( LocalHom ** ) );
 853                         for( i=0; i<njob; i++)
 854                                 localhomshrink[i] = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
 855                 }
 856 #endif
 857                 effarr1_kozo = AllocateDoubleVec( njob ); //tsuneni allocate sareru.
 858                 effarr2_kozo = AllocateDoubleVec( njob ); //tsuneni allocate sareru.
 859                 for( i=0; i<njob; i++ ) effarr1_kozo[i] = 0.0;
 860                 for( i=0; i<njob; i++ ) effarr2_kozo[i] = 0.0;
 861         }
 862
 863         for( i=0; i<njob-nadd; i++ ) alreadyaligned[i] = 1;
 864         for( i=njob-nadd; i<njob; i++ ) alreadyaligned[i] = 0;
 865
 866         for( l=0; l<njob; l++ ) fftlog[l] = 1;
 867
 868 #if 0
 869         fprintf( stderr, "##### fftwinsize = %d, fftthreshold = %d\n", fftWinSize, fftThreshold );
 870 #endif
 871
 872 #if 0
 873         for( i=0; i<njob; i++ )
 874                 fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
 875 #endif
 876
 877
 878         if( constraint )
 879                 calcimportance( njob, effarr, aseq, localhomtable );
 880
 881
 882 //      writePre( njob, name, nlen, aseq, 0 );
 883
 884         tscore = 0.0;
 885         for( l=0; l<njob-1; l++ )
 886         {
 887                 if( mergeoralign[l] == 'n' )
 888                 {
 889 //                      fprintf( stderr, "SKIP!\n" );
 890                         free( topol[l][0] );
 891                         free( topol[l][1] );
 892                         free( topol[l] );
 893                         continue;
 894                 }
 895
 896                 m1 = topol[l][0][0];
 897                 m2 = topol[l][1][0];
 898         len1 = strlen( aseq[m1] );
 899         len2 = strlen( aseq[m2] );
 900         if( *alloclen < len1 + len2 )
 901         {
 902                         fprintf( stderr, "\nReallocating.." );
 903                         *alloclen = ( len1 + len2 ) + 1000;
 904                         ReallocateCharMtx( aseq, njob, *alloclen + 10  );
 905                         gaplen = realloc( gaplen, ( *alloclen + 10 ) * sizeof( int ) );
 906                         if( gaplen == NULL )
 907                         {
 908                                 fprintf( stderr, "Cannot realloc gaplen\n" );
 909                                 exit( 1 );
 910                         }
 911                         gapmap = realloc( gapmap, ( *alloclen + 10 ) * sizeof( int ) );
 912                         if( gapmap == NULL )
 913                         {
 914                                 fprintf( stderr, "Cannot realloc gapmap\n" );
 915                                 exit( 1 );
 916                         }
 917                         fprintf( stderr, "done. *alloclen = %d\n", *alloclen );
 918                 }
 919
 920                 if( effarr_kozo )
 921                 {
 922                         clus1 = fastconjuction_noname_kozo( topol[l][0], aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
 923                         clus2 = fastconjuction_noname_kozo( topol[l][1], aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
 924                 }
 925                 else
 926                 {
 927                         clus1 = fastconjuction_noname( topol[l][0], aseq, mseq1, effarr1, effarr, indication1 );
 928                         clus2 = fastconjuction_noname( topol[l][1], aseq, mseq2, effarr2, effarr, indication2 );
 929                 }
 930
 931                 if( mergeoralign[l] == '1' || mergeoralign[l] == '2' )
 932                 {
 933                         newgapstr = "=";
 934                 }
 935                 else
 936                         newgapstr = "-";
 937
 938
 939                 len1nocommongap = len1;
 940                 len2nocommongap = len2;
 941                 if( mergeoralign[l] == '1' ) // nai
 942                 {
 943                         findcommongaps( clus2, mseq2, gapmap );
 944                         commongappick( clus2, mseq2 );
 945                         len2nocommongap = strlen( mseq2[0] );
 946                 }
 947                 else if( mergeoralign[l] == '2' )
 948                 {
 949                         findcommongaps( clus1, mseq1, gapmap );
 950                         commongappick( clus1, mseq1 );
 951                         len1nocommongap = strlen( mseq1[0] );
 952                 }
 953
 954
 955                 fprintf( trap_g, "\nSTEP-%d\n", l );
 956                 fprintf( trap_g, "group1 = %s\n", indication1 );
 957                 fprintf( trap_g, "group2 = %s\n", indication2 );
 958
 959 #if 1
 960                 fprintf( stderr, "\rSTEP % 5d /%d ", l+1, njob-1 );
 961                 fflush( stderr );
 962 #else
 963                 fprintf( stdout, "STEP %d /%d\n", l+1, njob-1 );
 964                 fprintf( stderr, "STEP %d /%d\n", l+1, njob-1 );
 965                 fprintf( stderr, "group1 = %.66s", indication1 );
 966                 if( strlen( indication1 ) > 66 ) fprintf( stderr, "..." );
 967                 fprintf( stderr, "\n" );
 968                 fprintf( stderr, "group2 = %.66s", indication2 );
 969                 if( strlen( indication2 ) > 66 ) fprintf( stderr, "..." );
 970                 fprintf( stderr, "\n" );
 971 #endif
 972
 973
 974
 975 //              for( i=0; i<clus1; i++ ) fprintf( stderr, "## STEP%d-eff for mseq1-%d %f\n", l+1, i, effarr1[i] );
 976
 977                 if( constraint )
 978                 {
 979                         fastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
 980 //                      msfastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
 981 //                      fprintf( stderr, "localhomshrink =\n" );
 982 //                      outlocalhompt( localhomshrink, clus1, clus2 );
 983 //                      weightimportance4( clus1, clus2, effarr1, effarr2, localhomshrink );
 984 //                      fprintf( stderr, "after weight =\n" );
 985 //                      outlocalhompt( localhomshrink, clus1, clus2 );
 986                 }
 987                 if( rnakozo && rnaprediction == 'm' )
 988                 {
 989                         makegrouprna( grouprna1, singlerna, topol[l][0] );
 990                         makegrouprna( grouprna2, singlerna, topol[l][1] );
 991                 }
 992
 993
 994 /*
 995                 fprintf( stderr, "before align all\n" );
 996                 display( aseq, njob );
 997                 fprintf( stderr, "\n" );
 998                 fprintf( stderr, "before align 1 %s \n", indication1 );
 999                 display( mseq1, clus1 );
1000                 fprintf( stderr, "\n" );
1001                 fprintf( stderr, "before align 2 %s \n", indication2 );
1002                 display( mseq2, clus2 );
1003                 fprintf( stderr, "\n" );
1004 */
1005
1006                 if( !nevermemsave && ( constraint != 2  && alg != 'M'  && ( len1 > 30000 || len2 > 30000 ) ) )
1007                 {
1008                         fprintf( stderr, "\nlen1=%d, len2=%d, Switching to the memsave mode.\n", len1, len2 );
1009                         alg = 'M';
1010                         if( commonIP ) FreeIntMtx( commonIP );
1011                         commonAlloc1 = 0;
1012                         commonAlloc2 = 0;
1013                 }
1014
1015
1016 //              if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 );
1017                 if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 1000 && clus2 < 1000 );
1018                 else                                               ffttry = 0;
1019 //              ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 5000 && clus2 < 5000 ); // v6.708
1020 //              fprintf( stderr, "f=%d, len1/fftlog[m1]=%f, clus1=%d, len2/fftlog[m2]=%f, clus2=%d\n", ffttry, (float)len1/fftlog[m1], clus1, (float)len2/fftlog[m2], clus2 );
1021 //              fprintf( stderr, "f=%d, clus1=%d, fftlog[m1]=%d, clus2=%d, fftlog[m2]=%d\n", ffttry, clus1, fftlog[m1], clus2, fftlog[m2] );
1022                 if( constraint == 2 )
1023                 {
1024                         if( alg == 'M' )
1025                         {
1026                                 fprintf( stderr, "\n\nMemory saving mode is not supported.\n\n" );
1027                                 exit( 1 );
1028                         }
1029                         fprintf( stderr, "c" );
1030                         if( alg == 'A' )
1031                         {
1032                                 imp_match_init_strict( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1033                                 if( rnakozo ) imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
1034                                 pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1035                         }
1036                         else if( alg == 'H' )
1037                         {
1038                                 imp_match_init_strictH( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1039                                 pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
1040                         }
1041                         else if( alg == 'Q' )
1042                         {
1043                                 imp_match_init_strictQ( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1044                                 if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
1045                                 pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
1046                         }
1047                         else if( alg == 'R' )
1048                         {
1049                                 imp_match_init_strictR( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1050                                 pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
1051                         }
1052                 }
1053                 else if( force_fft || ( use_fft && ffttry ) )
1054                 {
1055                         fprintf( stderr, "f" );
1056                         if( alg == 'M' )
1057                         {
1058                                 fprintf( stderr, "m" );
1059                                 pscore = Falign_udpari_long( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1 );
1060                         }
1061                         else
1062                                 pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL );
1063                 }
1064                 else
1065                 {
1066                         fprintf( stderr, "d" );
1067                         fftlog[m1] = 0;
1068                         switch( alg )
1069                         {
1070                                 case( 'a' ):
1071                                         pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen );
1072                                         break;
1073                                 case( 'M' ):
1074                                         fprintf( stderr, "m" );
1075                                         pscore = MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1076                                         break;
1077                                 case( 'A' ):
1078                                         pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1079                                         break;
1080                                 case( 'Q' ):
1081                                         pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
1082                                         break;
1083                                 case( 'R' ):
1084                                         pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
1085                                         break;
1086                                 case( 'H' ):
1087                                         pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
1088                                         break;
1089                                 default:
1090                                         ErrorExit( "ERROR IN SOURCE FILE" );
1091                         }
1092                 }
1093
1094                 nlen[m1] = 0.5 * ( nlen[m1] + nlen[m2] );
1095
1096 #if SCOREOUT
1097                 fprintf( stderr, "score = %10.2f\n", pscore );
1098 #endif
1099                 tscore += pscore;
1100 /*
1101                 fprintf( stderr, "after align 1 %s \n", indication1 );
1102                 display( mseq1, clus1 );
1103                 fprintf( stderr, "\n" );
1104                 fprintf( stderr, "after align 2 %s \n", indication2 );
1105                 display( mseq2, clus2 );
1106                 fprintf( stderr, "\n" );
1107 */
1108
1109 //              writePre( njob, name, nlen, aseq, 0 );
1110
1111                 if( disp ) display( aseq, njob );
1112
1113                 if( mergeoralign[l] == '1' ) // jissainiha nai. atarashii hairetsu ha saigo dakara.
1114                 {
1115                         adjustgapmap( strlen( mseq2[0] )-len2nocommongap+len2, gapmap, mseq2[0] );
1116                         restorecommongaps( njob, aseq, topol[l][0], topol[l][1], gapmap, *alloclen );
1117                         findnewgaps( clus2, mseq2, gaplen );
1118                         insertnewgaps( njob, alreadyaligned, aseq, topol[l][1], topol[l][0], gaplen, gapmap, *alloclen, alg );
1119                         for( i=0; i<njob; i++ ) eq2dash( aseq[i] );
1120                         for( i=0; (m=topol[l][0][i])>-1; i++ ) alreadyaligned[m] = 1;
1121                 }
1122                 if( mergeoralign[l] == '2' )
1123                 {
1124 //                      fprintf( stderr, ">mseq1[0] = \n%s\n", mseq1[0] );
1125 //                      fprintf( stderr, ">mseq2[0] = \n%s\n", mseq2[0] );
1126                         adjustgapmap( strlen( mseq1[0] )-len1nocommongap+len1, gapmap, mseq1[0] );
1127                         restorecommongaps( njob, aseq, topol[l][0], topol[l][1], gapmap, *alloclen );
1128                         findnewgaps( clus1, mseq1, gaplen );
1129                         insertnewgaps( njob, alreadyaligned, aseq, topol[l][0], topol[l][1], gaplen, gapmap, *alloclen, alg );
1130                         for( i=0; i<njob; i++ ) eq2dash( aseq[i] );
1131                         for( i=0; (m=topol[l][1][i])>-1; i++ ) alreadyaligned[m] = 1;
1132                 }
1133
1134                 free( topol[l][0] );
1135                 free( topol[l][1] );
1136                 free( topol[l] );
1137         }
1138 #if SCOREOUT
1139         fprintf( stderr, "totalscore = %10.2f\n\n", tscore );
1140 #endif
1141 }
1142
1143 static void WriteOptions( FILE *fp )
1144 {
1145
1146         if( dorp == 'd' ) fprintf( fp, "DNA\n" );
1147         else
1148         {
1149                 if     ( scoremtx ==  0 ) fprintf( fp, "JTT %dPAM\n", pamN );
1150                 else if( scoremtx ==  1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
1151                 else if( scoremtx ==  2 ) fprintf( fp, "M-Y\n" );
1152         }
1153     fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
1154     if( use_fft ) fprintf( fp, "FFT on\n" );
1155
1156         fprintf( fp, "tree-base method\n" );
1157         if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
1158         else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
1159         if( tbitr || tbweight )
1160         {
1161                 fprintf( fp, "iterate at each step\n" );
1162                 if( tbitr && tbrweight == 0 ) fprintf( fp, "  unweighted\n" );
1163                 if( tbitr && tbrweight == 3 ) fprintf( fp, "  reversely weighted\n" );
1164                 if( tbweight ) fprintf( fp, "  weighted\n" );
1165                 fprintf( fp, "\n" );
1166         }
1167
1168          fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
1169
1170         if( alg == 'a' )
1171                 fprintf( fp, "Algorithm A\n" );
1172         else if( alg == 'A' )
1173                 fprintf( fp, "Algorithm A+\n" );
1174         else if( alg == 'C' )
1175                 fprintf( fp, "Apgorithm A+/C\n" );
1176         else
1177                 fprintf( fp, "Unknown algorithm\n" );
1178
1179         if( treemethod == 'X' )
1180                 fprintf( fp, "Tree = UPGMA (mix).\n" );
1181         else if( treemethod == 'E' )
1182                 fprintf( fp, "Tree = UPGMA (average).\n" );
1183         else if( treemethod == 'q' )
1184                 fprintf( fp, "Tree = Minimum linkage.\n" );
1185         else
1186                 fprintf( fp, "Unknown tree.\n" );
1187
1188     if( use_fft )
1189     {
1190         fprintf( fp, "FFT on\n" );
1191         if( dorp == 'd' )
1192             fprintf( fp, "Basis : 4 nucleotides\n" );
1193         else
1194         {
1195             if( fftscore )
1196                 fprintf( fp, "Basis : Polarity and Volume\n" );
1197             else
1198                 fprintf( fp, "Basis : 20 amino acids\n" );
1199         }
1200         fprintf( fp, "Threshold   of anchors = %d%%\n", fftThreshold );
1201         fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
1202     }
1203         else
1204         fprintf( fp, "FFT off\n" );
1205         fflush( fp );
1206 }
1207
1208
1209 int main( int argc, char *argv[] )
1210 {
1211         static int  *nlen;
1212         static float *selfscore;
1213         int nogaplen;
1214         static char **name, **seq;
1215         static char **mseq1, **mseq2;
1216         static char **bseq;
1217         static float **iscore, **iscore_kozo;
1218         static double *eff, *eff_kozo, *eff_kozo_mapped = NULL;
1219         int i, j, ien, ik, jk;
1220         static int ***topol, ***topol_kozo;
1221         static int *addmem;
1222         static Treedep *dep;
1223         static float **len, **len_kozo;
1224         FILE *prep;
1225         FILE *infp;
1226         FILE *orderfp;
1227         FILE *hat2p;
1228         double unweightedspscore;
1229         int alignmentlength;
1230         char *mergeoralign;
1231         int foundthebranch;
1232
1233         char c;
1234         int alloclen;
1235         LocalHom **localhomtable = NULL;
1236         RNApair ***singlerna;
1237         float ssi, ssj, bunbo;
1238         static char *kozoarivec;
1239         int nkozo;
1240
1241         arguments( argc, argv );
1242 #ifndef enablemultithread
1243         nthread = 0;
1244 #endif
1245
1246         if( inputfile )
1247         {
1248                 infp = fopen( inputfile, "r" );
1249                 if( !infp )
1250                 {
1251                         fprintf( stderr, "Cannot open %s\n", inputfile );
1252                         exit( 1 );
1253                 }
1254         }
1255         else
1256                 infp = stdin;
1257
1258         getnumlen( infp );
1259         rewind( infp );
1260
1261
1262         nkozo = 0;
1263
1264         if( njob < 2 )
1265         {
1266                 fprintf( stderr, "At least 2 sequences should be input!\n"
1267                                                  "Only %d sequence found.\n", njob );
1268                 exit( 1 );
1269         }
1270
1271         seq = AllocateCharMtx( njob, nlenmax+1 );
1272         mseq1 = AllocateCharMtx( njob, 0 );
1273         mseq2 = AllocateCharMtx( njob, 0 );
1274
1275         name = AllocateCharMtx( njob, B+1 );
1276         nlen = AllocateIntVec( njob );
1277         selfscore = AllocateFloatVec( njob );
1278
1279         topol = AllocateIntCub( njob, 2, 0 );
1280         len = AllocateFloatMtx( njob, 2 );
1281         iscore = AllocateFloatHalfMtx( njob );
1282         eff = AllocateDoubleVec( njob );
1283         kozoarivec = AllocateCharVec( njob );
1284
1285         mergeoralign = AllocateCharVec( njob );
1286
1287         dep = (Treedep *)calloc( njob, sizeof( Treedep ) );
1288         if( nadd ) addmem = AllocateIntVec( nadd+1 );
1289
1290         if( constraint )
1291         {
1292                 localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
1293                 for( i=0; i<njob; i++)
1294                 {
1295                         localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
1296                         for( j=0; j<njob; j++)
1297                         {
1298                                 localhomtable[i][j].start1 = -1;
1299                                 localhomtable[i][j].end1 = -1;
1300                                 localhomtable[i][j].start2 = -1;
1301                                 localhomtable[i][j].end2 = -1;
1302                                 localhomtable[i][j].overlapaa = -1.0;
1303                                 localhomtable[i][j].opt = -1.0;
1304                                 localhomtable[i][j].importance = -1.0;
1305                                 localhomtable[i][j].next = NULL;
1306                                 localhomtable[i][j].korh = 'h';
1307                         }
1308                 }
1309
1310                 fprintf( stderr, "Loading 'hat3' ... " );
1311                 prep = fopen( "hat3", "r" );
1312                 if( prep == NULL ) ErrorExit( "Make hat3." );
1313                 readlocalhomtable( prep, njob, localhomtable, kozoarivec );
1314                 fclose( prep );
1315                 fprintf( stderr, "\ndone.\n" );
1316
1317
1318                 nkozo = 0;
1319                 for( i=0; i<njob; i++ )
1320                 {
1321 //                      fprintf( stderr, "kozoarivec[%d] = %d\n", i, kozoarivec[i] );
1322                         if( kozoarivec[i] ) nkozo++;
1323                 }
1324                 if( nkozo )
1325                 {
1326                         topol_kozo = AllocateIntCub( nkozo, 2, 0 );
1327                         len_kozo = AllocateFloatMtx( nkozo, 2 );
1328                         iscore_kozo = AllocateFloatHalfMtx( nkozo );
1329                         eff_kozo = AllocateDoubleVec( nkozo );
1330                         eff_kozo_mapped = AllocateDoubleVec( njob );
1331                 }
1332
1333
1334 //              outlocalhom( localhomtable, njob );
1335
1336 #if 0
1337                 fprintf( stderr, "Extending localhom ... " );
1338                 extendlocalhom2( njob, localhomtable );
1339                 fprintf( stderr, "done.\n" );
1340 #endif
1341         }
1342
1343 #if 0
1344         readData( infp, name, nlen, seq );
1345 #else
1346         readData_pointer( infp, name, nlen, seq );
1347         fclose( infp );
1348 #endif
1349
1350         constants( njob, seq );
1351
1352 #if 0
1353         fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
1354 #endif
1355
1356         initSignalSM();
1357
1358         initFiles();
1359
1360         WriteOptions( trap_g );
1361
1362         c = seqcheck( seq );
1363         if( c )
1364         {
1365                 fprintf( stderr, "Illegal character %c\n", c );
1366                 exit( 1 );
1367         }
1368
1369 //      writePre( njob, name, nlen, seq, 0 );
1370
1371         if( treein )
1372         {
1373 #if 0
1374                 if( nkozo )
1375                 {
1376                         fprintf( stderr, "Both structure and user tree have been given. Not yet supported!\n" );
1377                         exit( 1 );
1378                 }
1379 #endif
1380                 fprintf( stderr, "Loading a tree ... " );
1381                 loadtree( njob, topol, len, name, nlen, dep );
1382                 fprintf( stderr, "\ndone.\n\n" );
1383         }
1384         else
1385         {
1386                 if( tbutree == 0 )
1387                 {
1388                         for( i=1; i<njob; i++ )
1389                         {
1390                                 if( nlen[i] != nlen[0] )
1391                                 {
1392                                         fprintf( stderr, "Input pre-aligned seqences or make hat2.\n" );
1393                                         exit( 1 );
1394                                 }
1395                         }
1396
1397                         fprintf( stderr, "Making a distance matrix .. \n" );
1398                         fflush( stderr );
1399                         ien = njob-1;
1400                         for( i=0; i<njob; i++ )
1401                         {
1402                                 selfscore[i] = naivepairscore11( seq[i], seq[i], penalty );
1403                         }
1404 #ifdef enablemultithread
1405                         if( nthread > 0 )
1406                         {
1407                                 distancematrixthread_arg_t *targ;
1408                                 Jobtable jobpos;
1409                                 pthread_t *handle;
1410                                 pthread_mutex_t mutex;
1411
1412                                 jobpos.i = 0;
1413                                 jobpos.j = 0;
1414
1415                                 targ = calloc( nthread, sizeof( distancematrixthread_arg_t ) );
1416                                 handle = calloc( nthread, sizeof( pthread_t ) );
1417                                 pthread_mutex_init( &mutex, NULL );
1418
1419                                 for( i=0; i<nthread; i++ )
1420                                 {
1421                                         targ[i].thread_no = i;
1422                                         targ[i].njob = njob;
1423                                         targ[i].selfscore = selfscore;
1424                                         targ[i].iscore = iscore;
1425                                         targ[i].seq = seq;
1426                                         targ[i].jobpospt = &jobpos;
1427                                         targ[i].mutex = &mutex;
1428
1429                                         pthread_create( handle+i, NULL, distancematrixthread, (void *)(targ+i) );
1430                                 }
1431
1432                                 for( i=0; i<nthread; i++ )
1433                                 {
1434                                         pthread_join( handle[i], NULL );
1435                                 }
1436                                 pthread_mutex_destroy( &mutex );
1437                                 free( handle );
1438                                 free( targ );
1439                         }
1440                         else
1441 #endif
1442                         {
1443                                 for( i=0; i<ien; i++ )
1444                                 {
1445                                         if( i % 10 == 0 )
1446                                         {
1447                                                 fprintf( stderr, "\r% 5d / %d", i, ien );
1448                                                 fflush( stderr );
1449                                         }
1450                                         ssi = selfscore[i];
1451                                         for( j=i+1; j<njob; j++ )
1452                                         {
1453                                                 ssj = selfscore[j];
1454                                                 bunbo = MIN( ssi, ssj );
1455                                                 if( bunbo == 0.0 )
1456                                                         iscore[i][j-i] = 1.0;
1457                                                 else
1458 //                                                      iscore[i][j-i] = 1.0 - naivepairscore11( seq[i], seq[j], penalty ) / MIN( selfscore[i], selfscore[j] );
1459                                                         iscore[i][j-i] = 1.0 - naivepairscore11( seq[i], seq[j], penalty ) / bunbo;
1460
1461 #if 0
1462                                                 fprintf( stderr, "### ssj = %f\n", ssj );
1463                                                 fprintf( stderr, "### selfscore[i] = %f\n", selfscore[i] );
1464                                                 fprintf( stderr, "### selfscore[j] = %f\n", selfscore[j] );
1465                                                 fprintf( stderr, "### rawscore = %f\n", naivepairscore11( seq[i], seq[j], penalty ) );
1466 #endif
1467                                         }
1468                                 }
1469                         }
1470                         fprintf( stderr, "\ndone.\n\n" );
1471                         fflush( stderr );
1472                 }
1473                 else
1474                 {
1475                         fprintf( stderr, "Loading 'hat2' ... " );
1476                         prep = fopen( "hat2", "r" );
1477                         if( prep == NULL ) ErrorExit( "Make hat2." );
1478                         readhat2_floathalf_pointer( prep, njob, name, iscore );
1479                         fclose( prep );
1480                         fprintf( stderr, "done.\n" );
1481                 }
1482 #if 1
1483                 if( distout )
1484                 {
1485                         hat2p = fopen( "hat2", "w" );
1486                         WriteFloatHat2_pointer_halfmtx( hat2p, njob, name, iscore );
1487                         fclose( hat2p );
1488                 }
1489 #endif
1490                 if( nkozo )
1491                 {
1492                         ien = njob-1;
1493                         ik = 0;
1494                         for( i=0; i<ien; i++ )
1495                         {
1496                                 jk = ik+1;
1497                                 for( j=i+1; j<njob; j++ )
1498                                 {
1499                                         if( kozoarivec[i] && kozoarivec[j] )
1500                                         {
1501                                                 iscore_kozo[ik][jk-ik] = iscore[i][j-i];
1502                                         }
1503                                         if( kozoarivec[j] ) jk++;
1504                                 }
1505                                 if( kozoarivec[i] ) ik++;
1506                         }
1507                 }
1508
1509                 fprintf( stderr, "Constructing a UPGMA tree ... " );
1510                 fflush( stderr );
1511                 if( topin )
1512                 {
1513                         fprintf( stderr, "Loading a topology ... " );
1514                         loadtop( njob, iscore, topol, len );
1515                         fprintf( stderr, "\ndone.\n\n" );
1516                 }
1517                 else if( treeout )
1518                 {
1519                         fixed_musclesupg_float_realloc_nobk_halfmtx_treeout( njob, iscore, topol, len, name, nlen, dep );
1520                 }
1521                 else
1522                 {
1523                         fixed_musclesupg_float_realloc_nobk_halfmtx( njob, iscore, topol, len, dep );
1524                 }
1525 //              else
1526 //                      ErrorExit( "Incorrect tree\n" );
1527
1528                 if( nkozo )
1529                 {
1530 //                      for( i=0; i<nkozo-1; i++ )
1531 //                              for( j=i+1; j<nkozo; j++ )
1532 //                                      fprintf( stderr, "iscore_kozo[%d][%d] =~ %f\n", i, j, iscore_kozo[i][j-i] );
1533                         fixed_musclesupg_float_realloc_nobk_halfmtx( nkozo, iscore_kozo, topol_kozo, len_kozo, NULL );
1534                 }
1535                 fprintf( stderr, "\ndone.\n\n" );
1536                 fflush( stderr );
1537         }
1538
1539
1540         orderfp = fopen( "order", "w" );
1541         if( !orderfp )
1542         {
1543                 fprintf( stderr, "Cannot open 'order'\n" );
1544                 exit( 1 );
1545         }
1546         for( i=0; (j=topol[njob-2][0][i])!=-1; i++ )
1547         {
1548                 fprintf( orderfp, "%d\n", j );
1549         }
1550         for( i=0; (j=topol[njob-2][1][i])!=-1; i++ )
1551         {
1552                 fprintf( orderfp, "%d\n", j );
1553         }
1554         fclose( orderfp );
1555
1556         if( treeout && noalign )
1557         {
1558                 writeData_pointer( prep_g, njob, name, nlen, seq );
1559                 fprintf( stderr, "\n" );
1560                 SHOWVERSION;
1561                 return( 0 );
1562         }
1563
1564 //      countnode( njob, topol, node0 );
1565         if( tbrweight )
1566         {
1567                 weight = 3;
1568 #if 0
1569                 utree = 0; counteff( njob, topol, len, eff ); utree = 1;
1570 #else
1571                 counteff_simple_float( njob, topol, len, eff );
1572
1573                 if( nkozo )
1574                 {
1575 //                      counteff_simple_float( nkozo, topol_kozo, len_kozo, eff_kozo ); // single weight nanode iranai
1576                         for( i=0,j=0; i<njob; i++ )
1577                         {
1578                                 if( kozoarivec[i] )
1579                                 {
1580 //                                      eff_kozo_mapped[i] = eff_kozo[j]; //
1581                                         eff_kozo_mapped[i] = eff[i];      // single weight
1582                                         j++;
1583                                 }
1584                                 else
1585                                         eff_kozo_mapped[i] = 0.0;
1586 //                              fprintf( stderr, "eff_kozo_mapped[%d] = %f\n", i, eff_kozo_mapped[i] );
1587 //                              fprintf( stderr, "            eff[%d] = %f\n", i, eff[i] );
1588                         }
1589                 }
1590
1591
1592 #endif
1593         }
1594         else
1595         {
1596                 for( i=0; i<njob; i++ ) eff[i] = 1.0;
1597                 if( nkozo )
1598                 {
1599                         for( i=0; i<njob; i++ )
1600                         {
1601                                 if( kozoarivec[i] )
1602                                         eff_kozo_mapped[i] = 1.0;
1603                                 else
1604                                         eff_kozo_mapped[i] = 0.0;
1605                         }
1606                 }
1607         }
1608
1609         FreeFloatHalfMtx( iscore, njob );
1610         FreeFloatMtx( len );
1611
1612         alloclen = nlenmax*2+1; //chuui!
1613         bseq = AllocateCharMtx( njob, alloclen );
1614
1615         if( nadd )
1616         {
1617                 alignmentlength = strlen( seq[0] );
1618                 for( i=0; i<njob-nadd; i++ )
1619                 {
1620                         if( alignmentlength != strlen( seq[i] ) )
1621                         {
1622                                 fprintf( stderr, "#################################################################################\n" );
1623                                 fprintf( stderr, "# ERROR!                                                                        #\n" );
1624                                 fprintf( stderr, "# The original%4d sequences must be aligned                                    #\n", njob-nadd );
1625                                 fprintf( stderr, "#################################################################################\n" );
1626                                 exit( 1 );
1627                         }
1628                 }
1629                 if( addprofile )
1630                 {
1631                         alignmentlength = strlen( seq[njob-nadd] );
1632                         for( i=njob-nadd; i<njob; i++ )
1633                         {
1634                                 if( alignmentlength != strlen( seq[i] ) )
1635                                 {
1636                                         fprintf( stderr, "###############################################################################\n" );
1637                                         fprintf( stderr, "# ERROR!                                                                      #\n" );
1638                                         fprintf( stderr, "# The%4d additional sequences must be aligned                                #\n", nadd );
1639                                         fprintf( stderr, "# Otherwise, try the '--add' option, instead of '--addprofile' option.        #\n" );
1640                                         fprintf( stderr, "###############################################################################\n" );
1641                                         exit( 1 );
1642                                 }
1643                         }
1644                         for( i=0; i<nadd; i++ ) addmem[i] = njob-nadd+i;
1645                         addmem[nadd] = -1;
1646                         foundthebranch = 0;
1647                         for( i=0; i<njob-1; i++ )
1648                         {
1649                                 if( samemember( topol[i][0], addmem ) ) // jissainiha nai
1650                                 {
1651                                         mergeoralign[i] = '1';
1652                                         foundthebranch = 1;
1653                                 }
1654                                 else if( samemember( topol[i][1], addmem ) )
1655                                 {
1656                                         mergeoralign[i] = '2';
1657                                         foundthebranch = 1;
1658                                 }
1659                                 else
1660                                 {
1661                                         mergeoralign[i] = 'n';
1662                                 }
1663                         }
1664                         if( !foundthebranch )
1665                         {
1666                                 fprintf( stderr, "###############################################################################\n" );
1667                                 fprintf( stderr, "# ERROR!                                                                      #\n" );
1668                                 fprintf( stderr, "# There is no appropriate position to add the%4d sequences in the guide tree.#\n", nadd );
1669                                 fprintf( stderr, "# Check whether the%4d sequences form a monophyletic cluster.                #\n", nadd );
1670                                 fprintf( stderr, "# If not, try the '--add' option, instead of the '--addprofile' option.       #\n" );
1671                                 fprintf( stderr, "############################################################################### \n" );
1672                                 exit( 1 );
1673                         }
1674                         commongappick( nadd, seq+njob-nadd );
1675                         for( i=njob-nadd; i<njob; i++ ) strcpy( bseq[i], seq[i] );
1676                 }
1677                 else
1678                 {
1679                         for( i=0; i<njob-1; i++ ) mergeoralign[i] = 'n';
1680                         for( j=njob-nadd; j<njob; j++ )
1681                         {
1682                                 addmem[0] = j;
1683                                 addmem[1] = -1;
1684                                 for( i=0; i<njob-1; i++ )
1685                                 {
1686                                         if( samemember( topol[i][0], addmem ) ) // arieru
1687                                         {
1688 //                                              fprintf( stderr, "HIT!\n" );
1689                                                 if( mergeoralign[i] != 'n' ) mergeoralign[i] = 'w';
1690                                                 else mergeoralign[i] = '1';
1691                                         }
1692                                         else if( samemember( topol[i][1], addmem ) )
1693                                         {
1694 //                                              fprintf( stderr, "HIT!\n" );
1695                                                 if( mergeoralign[i] != 'n' ) mergeoralign[i] = 'w';
1696                                                 else mergeoralign[i] = '2';
1697                                         }
1698                                 }
1699                         }
1700
1701                         for( i=0; i<nadd; i++ ) addmem[i] = njob-nadd+i;
1702                         addmem[nadd] = -1;
1703                         for( i=0; i<njob-1; i++ )
1704                         {
1705                                 if( includemember( topol[i][0], addmem ) && includemember( topol[i][1], addmem ) )
1706                                 {
1707                                         mergeoralign[i] = 'w';
1708                                 }
1709                                 else if( includemember( topol[i][0], addmem ) )
1710                                 {
1711                                         mergeoralign[i] = '1';
1712                                 }
1713                                 else if( includemember( topol[i][1], addmem ) )
1714                                 {
1715                                         mergeoralign[i] = '2';
1716                                 }
1717                         }
1718 #if 0
1719                         for( i=0; i<njob-1; i++ )
1720                         {
1721                                 fprintf( stderr, "mem0 = " );
1722                                 for( j=0; topol[i][0][j]>-1; j++ )      fprintf( stderr, "%d ", topol[i][0][j] );
1723                                 fprintf( stderr, "\n" );
1724                                 fprintf( stderr, "mem1 = " );
1725                                 for( j=0; topol[i][1][j]>-1; j++ )      fprintf( stderr, "%d ", topol[i][1][j] );
1726                                 fprintf( stderr, "\n" );
1727                                 fprintf( stderr, "i=%d, mergeoralign[] = %c\n", i, mergeoralign[i] );
1728                         }
1729 #endif
1730                         for( i=njob-nadd; i<njob; i++ ) gappick0( bseq[i], seq[i] );
1731                 }
1732
1733                 commongappick( njob-nadd, seq );
1734                 for( i=0; i<njob-nadd; i++ ) strcpy( bseq[i], seq[i] );
1735         }
1736         else
1737         {
1738                 for( i=0; i<njob; i++ ) gappick0( bseq[i], seq[i] );
1739                 for( i=0; i<njob-1; i++ ) mergeoralign[i] = 'a';
1740         }
1741
1742         if( rnakozo && rnaprediction == 'm' )
1743         {
1744                 singlerna = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
1745                 prep = fopen( "hat4", "r" );
1746                 if( prep == NULL ) ErrorExit( "Make hat4 using mccaskill." );
1747                 fprintf( stderr, "Loading 'hat4' ... " );
1748                 for( i=0; i<njob; i++ )
1749                 {
1750                         nogaplen = strlen( bseq[i] );
1751                         singlerna[i] = (RNApair **)calloc( nogaplen, sizeof( RNApair * ) );
1752                         for( j=0; j<nogaplen; j++ )
1753                         {
1754                                 singlerna[i][j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
1755                                 singlerna[i][j][0].bestpos = -1;
1756                                 singlerna[i][j][0].bestscore = -1.0;
1757                         }
1758                         readmccaskill( prep, singlerna[i], nogaplen );
1759                 }
1760                 fclose( prep );
1761                 fprintf( stderr, "\ndone.\n" );
1762         }
1763         else
1764                 singlerna = NULL;
1765
1766
1767         fprintf( stderr, "Progressive alignment ... \n" );
1768
1769 #ifdef enablemultithread
1770         if( nthread > 0 && nadd == 0 )
1771         {
1772                 treebasethread_arg_t *targ;
1773                 int jobpos;
1774                 pthread_t *handle;
1775                 pthread_mutex_t mutex;
1776                 pthread_cond_t treecond;
1777                 int *fftlog;
1778                 int nrun;
1779                 int nthread_yoyu;
1780
1781                 nthread_yoyu = nthread * 1;
1782                 nrun = 0;
1783                 jobpos = 0;
1784                 targ = calloc( nthread_yoyu, sizeof( treebasethread_arg_t ) );
1785                 fftlog = AllocateIntVec( njob );
1786                 handle = calloc( nthread_yoyu, sizeof( pthread_t ) );
1787                 pthread_mutex_init( &mutex, NULL );
1788                 pthread_cond_init( &treecond, NULL );
1789
1790                 for( i=0; i<njob; i++ ) dep[i].done = 0;
1791                 for( i=0; i<njob; i++ ) fftlog[i] = 1;
1792
1793                 if( constraint )
1794                         calcimportance( njob, eff, bseq, localhomtable );
1795
1796                 for( i=0; i<nthread_yoyu; i++ )
1797                 {
1798                         targ[i].thread_no = i;
1799                         targ[i].nrunpt = &nrun;
1800                         targ[i].njob = njob;
1801                         targ[i].nlen = nlen;
1802                         targ[i].jobpospt = &jobpos;
1803                         targ[i].topol = topol;
1804                         targ[i].dep = dep;
1805                         targ[i].aseq = bseq;
1806                         targ[i].effarr = eff;
1807                         targ[i].alloclenpt = &alloclen;
1808                         targ[i].localhomtable = localhomtable;
1809                         targ[i].singlerna = singlerna;
1810                         targ[i].effarr_kozo = eff_kozo_mapped;
1811                         targ[i].fftlog = fftlog;
1812                         targ[i].mutex = &mutex;
1813                         targ[i].treecond = &treecond;
1814
1815                         pthread_create( handle+i, NULL, treebasethread, (void *)(targ+i) );
1816                 }
1817
1818                 for( i=0; i<nthread_yoyu; i++ )
1819                 {
1820                         pthread_join( handle[i], NULL );
1821                 }
1822                 pthread_mutex_destroy( &mutex );
1823                 pthread_cond_destroy( &treecond );
1824                 free( handle );
1825                 free( targ );
1826                 free( fftlog );
1827         }
1828         else
1829 #endif
1830                 treebase( nlen, bseq, nadd, mergeoralign, mseq1, mseq2, topol, eff, &alloclen, localhomtable, singlerna, eff_kozo_mapped );
1831         fprintf( stderr, "\ndone.\n" );
1832         if( scoreout )
1833         {
1834                 unweightedspscore = plainscore( njob, bseq );
1835                 fprintf( stderr, "\nSCORE %s = %.0f, ", "(treebase)", unweightedspscore );
1836                 fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( njob * strlen( bseq[0] ) ) );
1837                 fprintf( stderr, "\n\n" );
1838         }
1839
1840 #if 0
1841         if( constraint )
1842         {
1843                 LocalHom *tmppt1, *tmppt2;
1844                 for( i=0; i<njob; i++)
1845                 {
1846                         for( j=0; j<njob; j++)
1847                         {
1848                                 tmppt1 = localhomtable[i]+j;
1849                                 while( tmppt2 = tmppt1->next )
1850                                 {
1851                                         free( (void *)tmppt1 );
1852                                         tmppt1 = tmppt2;
1853                                 }
1854                                 free( (void *)tmppt1 );
1855                         }
1856                         free( (void *)(localhomtable[i]+j) );
1857                 }
1858                 free( (void *)localhomtable );
1859         }
1860 #endif
1861
1862         fprintf( trap_g, "done.\n" );
1863         fclose( trap_g );
1864         free( mergeoralign );
1865
1866         writeData_pointer( prep_g, njob, name, nlen, bseq );
1867 #if 0
1868         writeData( stdout, njob, name, nlen, bseq );
1869         writePre( njob, name, nlen, bseq, !contin );
1870         writeData_pointer( prep_g, njob, name, nlen, aseq );
1871 #endif
1872 #if IODEBUG
1873         fprintf( stderr, "OSHIMAI\n" );
1874 #endif
1875
1876         if( constraint ) FreeLocalHomTable( localhomtable, njob );
1877
1878         SHOWVERSION;
1879         return( 0 );
1880 }