binaries/src/mafft/core/tditeration.c

   1
   2 /*
   3         tree-dependent iteration
   4     algorithm A+ when group-to-group, C when group-to-singleSeqence
   5                          OR
   6     algorithm A+
   7 */
   8
   9 #include "mltaln.h"
  10
  11
  12 #define DEBUG 0
  13 #define RECORD 0
  14
  15 extern char **seq_g;
  16 extern char **res_g;
  17
  18 static int nwa;
  19
  20 #ifdef enablemultithread
  21 typedef struct _threadarg
  22 {
  23         int thread_no;
  24         int *jobposintpt;
  25         int *ndonept;
  26         int *ntrypt;
  27         int *collectingpt;
  28         int njob;
  29         int nbranch;
  30         int maxiter;
  31         int nkozo;
  32         int *subgenerationpt;
  33         float *basegainpt;
  34         float *gainlist;
  35         float *tscorelist;
  36         int *generationofinput;
  37         char *kozoarivec;
  38         char **mastercopy;
  39         char ***candidates;
  40         int *generationofmastercopypt;
  41         int *branchtable;
  42         RNApair ***singlerna;
  43         LocalHom **localhomtable;
  44         int alloclen;
  45         Node *stopol;
  46         int ***topol;
  47 //      double **len;
  48         float **tscorehistory_detail;
  49         int *finishpt;
  50         pthread_mutex_t *mutex;
  51         pthread_cond_t *collection_end;
  52         pthread_cond_t *collection_start;
  53 } threadarg_t;
  54 #endif
  55
  56 #if 1
  57 static void shuffle( int *arr, int n )
  58 {
  59         int i;
  60         int x;
  61         int b;
  62
  63         for( i=1; i<n; i++ )
  64         {
  65                 x = rand() % (i+1);
  66                 if( x != i )
  67                 {
  68                         b = arr[i];
  69                         arr[i] = arr[x];
  70                         arr[x] = b;
  71                 }
  72         }
  73 }
  74 #endif
  75
  76 static void Writeoption2( FILE *fp, int cycle, double cut )
  77 {
  78         fprintf( fp, "%dth cycle\n", cycle );
  79     fprintf( fp, "marginal score to search : current score * (100-%d) / 100\n", (int)cut );
  80 }
  81
  82 static void Writeoptions( FILE *fp )
  83 {
  84         fprintf( fp, "Tree-dependent-iteration\n" );
  85     if( scoremtx == 1 )
  86         fprintf( fp, "Blosum %d\n", nblosum );
  87     else if( scoremtx == -1 )
  88         fprintf( fp, "DNA\n" );
  89     else if( scoremtx == 2 )
  90         fprintf( fp, "Miyata-Yasunaga\n" );
  91         else
  92                 fprintf( fp, "JTT %dPAM\n", pamN );
  93
  94         if( scoremtx == 0 || scoremtx == 1 )
  95         fprintf( fp, "Gap Penalty = %+5.3f, %5.2f, %+5.3f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
  96         else
  97                 fprintf( fp, "Gap Penalty = %+5.3f\n", (double)penalty/1000 );
  98
  99
 100     if( scmtd == 3 )
 101         fprintf( fp, "score of rnd or sco\n" );
 102     else if( scmtd == 4 )
 103         fprintf( fp, "score = sigma( score for a pair of homologous amino acids ) / ( number of amino acids pairs )\n" );
 104     else if( scmtd == 5 )
 105         fprintf( fp, "score : SP\n" );
 106     if( mix )
 107         fprintf( fp, "?\n" );
 108     else
 109     {
 110         if( weight == 2 )
 111             fprintf( fp, "weighted rationale-1,  geta2 = %f\n", geta2 );
 112         else if( weight == 3 )
 113             fprintf( fp, "weighted like ClustalW," );
 114         else if( weight == 4 )
 115             fprintf( fp, "weighted rationale-2,  geta2 = %f\n", geta2 );
 116         else
 117             fprintf( fp, "unweighted\n" );
 118     }
 119     if( weight && utree )
 120         fprintf( fp, "using tree defined by the file hat2.\n" );
 121     if( weight && !utree )
 122         fprintf( fp, "using temporary tree.\n" );
 123
 124         if( treemethod == 'n' )
 125                 fprintf( fp, "Tree is calculated with Neighbor-Joining Method.\n" );
 126         else if( treemethod == 'q' )
 127                 fprintf( fp, "Tree is calculated with Minimum linkage.\n" );
 128         else if( treemethod == 'X' )
 129                 fprintf( fp, "Tree is calculated with simplified UPG Method and UPG Method.\n" );
 130         else if( treemethod == 'E' )
 131                 fprintf( fp, "Tree is calculated with UPG Method.\n" );
 132         else
 133                 fprintf( fp, "Tree is calculated with unknown Method.\n" );
 134
 135         if( alg == 'C' )
 136                 fprintf( fp, "Algorithm A+ / C\n" );
 137         else if( alg == 'A' )
 138                 fprintf( fp, "Algorithm A+ \n" );
 139         else if( alg == 'a' )
 140                 fprintf( fp, "Algorithm A \n" );
 141         else
 142                 fprintf( fp, "Algorithm ? \n" );
 143
 144         if( use_fft )
 145         {
 146                 if( scoremtx == -1 )
 147                 {
 148                         fprintf( fp, "Basis : 4 nucleotides\n" );
 149                 }
 150                 else
 151                 {
 152                         if( fftscore )
 153                                 fprintf( fp, "Basis : Polarity and Volume\n" );
 154                         else
 155                                 fprintf( fp, "Basis : 20 amino acids\n" );
 156                 }
 157                 fprintf( fp, "Threshold   of anchors = %d%%\n", fftThreshold );
 158                 fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
 159         }
 160 }
 161
 162 #ifdef enablemultithread
 163
 164 static void freelocalarrays(
 165         float *tscorehistory,
 166         RNApair ***grouprna1, RNApair ***grouprna2,
 167         RNApair *rnapairboth,
 168         char *indication1, char *indication2,
 169         double *effarr, double *effarrforlocalhom, double *effarr1, double *effarr2,
 170         char **mseq1, char **mseq2,
 171         char **localcopy,
 172         int *gapmap1, int *gapmap2,
 173         double *effarr1_kozo, double *effarr2_kozo, double *effarr_kozo,
 174         char **pair,
 175         LocalHom *** localhomshrink
 176 )
 177 {
 178 //      fprintf( stderr, "Skipping freelocalarrays\n" );
 179 //      return;
 180         int i;
 181         if( commonIP ) FreeIntMtx( commonIP );
 182         commonIP = NULL;
 183         Falign( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
 184         Falign_localhom( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL,NULL, 0, NULL );
 185         if( rnakozo && rnaprediction == 'm' )
 186         {
 187                 free( grouprna1 ); // nakamimo?
 188                 free( grouprna2 ); // nakamimo?
 189         }
 190
 191         free( tscorehistory );
 192         free( indication1 );
 193         free( indication2 );
 194         free( effarr );
 195         free( effarrforlocalhom );
 196         free( effarr1 );
 197         free( effarr2 );
 198         free( mseq1 );
 199         free( mseq2 );
 200         FreeCharMtx( localcopy );
 201         free( gapmap1 );
 202         free( gapmap2 );
 203
 204         free( effarr1_kozo );
 205         free( effarr2_kozo );
 206         free( effarr_kozo );
 207
 208         FreeCharMtx( pair );
 209
 210         if( rnakozo ) free( rnapairboth );
 211
 212         if( constraint )
 213         {
 214                 for( i=0; i<njob; i++)
 215                 {
 216                         free( localhomshrink[i] ); // nakamimo??
 217                 }
 218                 free( localhomshrink );
 219         }
 220 }
 221
 222
 223 static void *athread( void *arg )
 224 {
 225
 226         threadarg_t *targ = (threadarg_t *)arg;
 227         int thread_no = targ->thread_no;
 228         int njob = targ->njob;
 229         int nbranch = targ->nbranch;
 230         int maxiter = targ->maxiter;
 231         int *ndonept = targ->ndonept;
 232         int *ntrypt = targ->ntrypt;
 233         int *collectingpt = targ->collectingpt;
 234         int *jobposintpt = targ->jobposintpt;
 235         int nkozo = targ->nkozo;
 236         float *gainlist = targ->gainlist;
 237         float *tscorelist = targ->tscorelist;
 238         int *generationofinput = targ->generationofinput;
 239         int *subgenerationpt = targ->subgenerationpt;
 240         float *basegainpt = targ->basegainpt;
 241         char *kozoarivec = targ->kozoarivec;
 242         char **mastercopy = targ->mastercopy;
 243         char ***candidates = targ->candidates;
 244         int *generationofmastercopypt = targ->generationofmastercopypt;
 245         int *branchtable = targ->branchtable;
 246         RNApair ***singlerna = targ->singlerna;
 247         LocalHom **localhomtable = targ->localhomtable;
 248         int alloclen = targ->alloclen;
 249         Node * stopol = targ->stopol;
 250         int ***topol = targ->topol;
 251 //      double **len = targ->len;
 252         float **tscorehistory_detail = targ->tscorehistory_detail;
 253         int *finishpt = targ->finishpt;
 254
 255         int i, j, k, l, ii;
 256         float gain;
 257         int iterate;
 258         char **pair;
 259         int locnjob;
 260         int s1, s2;
 261         int clus1, clus2;
 262         char **localcopy;
 263         char **mseq1, **mseq2;
 264         double *effarr, *effarr_kozo; // re-calc
 265         double *effarr1, *effarr2, *effarr1_kozo, *effarr2_kozo;
 266         char *indication1, *indication2;
 267         int length;
 268         RNApair ***grouprna1, ***grouprna2;
 269         RNApair *rnapairboth;
 270         LocalHom ***localhomshrink;
 271         int *gapmap1, *gapmap2;
 272         float tscore, mscore, oimpmatch, impmatch;
 273         int identity;
 274         double tmpdouble;
 275         float naivescore0 = 0, naivescore1;
 276         double *effarrforlocalhom;
 277         float *tscorehistory;
 278         int intdum;
 279 #if 0
 280         int oscillating;
 281         int lin, ldf;
 282 #endif
 283         float maxgain;
 284         int bestthread;
 285         int branchpos;
 286         int subgenerationatfirst;
 287         double unweightedspscore;
 288         int myjob;
 289         int converged2 = 0;
 290         int chudanres;
 291
 292
 293         locnjob = njob;
 294
 295         if( utree == 0 )
 296         {
 297                 fprintf( stderr, "Dynamic tree is not supported in the multithread version.\n" );
 298                 exit( 1 );
 299         }
 300         if( score_check == 2 )
 301         {
 302                 fprintf( stderr, "Score_check 2 is not supported in the multithread version.\n" );
 303                 exit( 1 );
 304         }
 305
 306         if( weight == 2 )
 307         {
 308                 fprintf( stderr, "Weight 2 is not supported in the multithread version.\n" );
 309                 exit( 1 );
 310         }
 311         if( cooling &&  cut > 0.0 )
 312         {
 313                 fprintf( stderr, "Cooling is not supported in the multithread version.\n" );
 314                 exit( 1 );
 315         }
 316
 317         tscorehistory = calloc( maxiter, sizeof( float ) );
 318
 319         if( rnakozo && rnaprediction == 'm' )
 320         {
 321                 grouprna1 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
 322                 grouprna2 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
 323         }
 324         else
 325         {
 326                 grouprna1 = grouprna2 = NULL;
 327         }
 328
 329         indication1 = AllocateCharVec( njob*3+50 );
 330         indication2 = AllocateCharVec( njob*3+50 );
 331         effarr = AllocateDoubleVec( locnjob );
 332         effarrforlocalhom = AllocateDoubleVec( locnjob );
 333         effarr1 = AllocateDoubleVec( locnjob );
 334         effarr2 = AllocateDoubleVec( locnjob );
 335         mseq1 = AllocateCharMtx( locnjob, 0 );
 336         mseq2 = AllocateCharMtx( locnjob, 0 );
 337         localcopy = AllocateCharMtx( locnjob, alloclen );
 338         gapmap1 = AllocateIntVec( alloclen );
 339         gapmap2 = AllocateIntVec( alloclen );
 340
 341         effarr1_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
 342         effarr2_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
 343         effarr_kozo = AllocateDoubleVec( locnjob );
 344         for( i=0; i<locnjob; i++ )
 345                 effarr_kozo[i] = effarr1_kozo[i] = effarr2_kozo[i] = 0.0;
 346
 347         pair = AllocateCharMtx( locnjob, locnjob );
 348
 349
 350         if( rnakozo ) rnapairboth = (RNApair *)calloc( alloclen, sizeof( RNApair ) );
 351
 352         if( constraint )
 353         {
 354                 localhomshrink = (LocalHom ***)calloc( njob, sizeof( LocalHom ** ) );
 355                 for( i=0; i<njob; i++)
 356                 {
 357                         localhomshrink[i] = (LocalHom **)calloc( njob, sizeof( LocalHom * ) );
 358                 }
 359         }
 360
 361
 362         if( thread_no == 0 )
 363         {
 364                 *ntrypt = 0;
 365                 srand( randomseed );
 366                 *finishpt = 0;
 367                 for( iterate=0; iterate<maxiter; iterate++ )
 368                 {
 369                         pthread_mutex_lock( targ->mutex );
 370
 371                         if( *collectingpt == 1 )
 372                         {
 373                                 *collectingpt = 0;
 374                                 *generationofmastercopypt = iterate;
 375                                 *subgenerationpt = 0;
 376                                 *basegainpt = 0.0;
 377                                 *ndonept = 0;
 378                                 *jobposintpt = 0;
 379                                 for( i=0; i<nwa; i++ ) gainlist[i] = 0;
 380                                 for( i=0; i<nwa; i++ ) tscorelist[i] = 0.0;
 381                                 for( i=0; i<nbranch; i++ ) generationofinput[i] = -1;
 382                                 if( parallelizationstrategy != BESTFIRST && randomseed != 0 ) shuffle( branchtable, nbranch );
 383                                 pthread_cond_broadcast( targ->collection_end );
 384                                 pthread_mutex_unlock( targ->mutex );
 385                         }
 386                         else
 387                         {
 388                                 pthread_cond_broadcast( targ->collection_end );
 389                                 pthread_mutex_unlock( targ->mutex );
 390                                 freelocalarrays
 391                                 (
 392                                         tscorehistory,
 393                                         grouprna1, grouprna2,
 394                                         rnapairboth,
 395                                         indication1, indication2,
 396                                         effarr, effarrforlocalhom, effarr1, effarr2,
 397                                         mseq1, mseq2,
 398                                         localcopy,
 399                                         gapmap1, gapmap2,
 400                                         effarr1_kozo, effarr2_kozo, effarr_kozo,
 401                                         pair,
 402                                         localhomshrink
 403                                 );
 404 //                              return( NULL );
 405                                 pthread_exit( NULL );
 406                         }
 407
 408                         pthread_mutex_lock( targ->mutex );
 409                         while( *ndonept < nbranch )
 410                                 pthread_cond_wait( targ->collection_start, targ->mutex );
 411                         pthread_mutex_unlock( targ->mutex );
 412 //                      fprintf( stderr, "Thread 0 got a signal, *collectionpt = %d\n", *collectingpt );
 413
 414 /*
 415                         Hoka no thread ga keisan
 416 */
 417
 418                         pthread_mutex_lock( targ->mutex );
 419                         *collectingpt = 1; // chofuku
 420
 421 #if 0
 422                         for( i=0; i<nbranch; i++ )
 423                         {
 424                                 if( generationofinput[i] != iterate )
 425                                 {
 426                                         fprintf( stderr, "Error! generationofinput[%d] = %d, but iterate=%d\n", i, generationofinput[i], iterate );
 427                                         exit( 1 );
 428
 429                                 }
 430                         }
 431 #endif
 432
 433                         maxgain = gainlist[1];
 434                         bestthread = 1;
 435                         for( i=2; i<nwa; i++ )
 436                         {
 437                                 if( gainlist[i] > maxgain )
 438                                 {
 439                                         maxgain = gainlist[i];
 440                                         bestthread = i;
 441                                 }
 442                         }
 443
 444                         if( maxgain > 0.0 )
 445                         {
 446 //                              fprintf( stderr, "\nGain = %f\n", maxgain );
 447 //                              fprintf( stderr, "best gain = %f by thread %d\n", gainlist[bestthread], bestthread );
 448 //                              fprintf( stderr, "tscorelist[best] = %f by thread %d\n", tscorelist[bestthread], bestthread );
 449                                 if( parallelizationstrategy == BESTFIRST )
 450                                 {
 451                                         for( i=0; i<locnjob; i++ ) strcpy( mastercopy[i], candidates[bestthread][i] );
 452                                         if( scoreout )
 453                                         {
 454                                                 unweightedspscore = plainscore( locnjob, mastercopy );
 455                                                 fprintf( stderr, "\nSCORE %d = %.0f, ", iterate * nbranch, unweightedspscore );
 456                                                 fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( locnjob * strlen( mastercopy[0] ) ) );
 457                                                 if( weight || constraint ) fprintf( stderr, " (differs from the objective score)" );
 458                                                 fprintf( stderr, "\n" );
 459                                         }
 460                                 }
 461 #if 1
 462 //                              fprintf( stderr, "gain(%d, by %d) = %f\n", iterate, bestthread, maxgain );
 463                                 for( i=iterate-1; i>0; i-- )
 464                                 {
 465 //                                      if( iterate-i < 15 ) fprintf( stderr, "hist[%d] = %f\n", i, tscorehistory[i] );
 466                                         if( tscorehistory[i] == tscorelist[bestthread] )
 467                                         {
 468                                                 fprintf( stderr, "\nOscillating? %f == %f\n", tscorehistory[i], tscorelist[bestthread] );
 469                                                 *collectingpt = -1;
 470                                                 break;
 471                                         }
 472                                 }
 473                                 tscorehistory[iterate] = tscorelist[bestthread];
 474 #endif
 475                         }
 476                         else
 477                         {
 478                                 fprintf( stderr, "\nConverged.\n" );
 479                                 *collectingpt = -1;
 480 //                              pthread_cond_broadcast( targ->collection_end );
 481 //                              pthread_mutex_unlock( targ->mutex );
 482 //                              freelocalarrays();
 483 //                              return( NULL );
 484 //                              pthread_exit( NULL );
 485                         }
 486
 487 #if 1
 488                         if( *finishpt )
 489                         {
 490                                 fprintf( stderr, "\nConverged2.\n" );
 491                                 *collectingpt = -1;
 492                         }
 493 #endif
 494
 495                         pthread_mutex_unlock( targ->mutex );
 496                 }
 497                 pthread_mutex_lock( targ->mutex );
 498                 fprintf( stderr, "\nReached %d\n", maxiter );
 499                 *collectingpt = -1;
 500                 pthread_cond_broadcast( targ->collection_end );
 501                 pthread_mutex_unlock( targ->mutex );
 502                 freelocalarrays
 503                 (
 504                         tscorehistory,
 505                         grouprna1, grouprna2,
 506                         rnapairboth,
 507                         indication1, indication2,
 508                         effarr, effarrforlocalhom, effarr1, effarr2,
 509                         mseq1, mseq2,
 510                         localcopy,
 511                         gapmap1, gapmap2,
 512                         effarr1_kozo, effarr2_kozo, effarr_kozo,
 513                         pair,
 514                         localhomshrink
 515                 );
 516                 return( NULL );
 517                 pthread_exit( NULL );
 518         }
 519         else
 520         {
 521                 while( 1 )
 522                 {
 523 #if 0
 524                         if( iterate % 2 == 0 )
 525                         {
 526                                 lin = 0; ldf = +1;
 527                         }
 528                         else
 529                         {
 530                                 lin = locnjob - 2; ldf = -1;
 531                         }
 532                         for( l=lin; l < locnjob-1 && l >= 0 ; l+=ldf )
 533                                 for( k=0; k<2; k++ )
 534 #endif
 535
 536                         pthread_mutex_lock( targ->mutex );
 537                         while( *collectingpt > 0 )
 538                                 pthread_cond_wait( targ->collection_end, targ->mutex );
 539                         if( *collectingpt == -1 )
 540                         {
 541                                 pthread_mutex_unlock( targ->mutex );
 542                                 freelocalarrays
 543                                 (
 544                                         tscorehistory,
 545                                         grouprna1, grouprna2,
 546                                         rnapairboth,
 547                                         indication1, indication2,
 548                                         effarr, effarrforlocalhom, effarr1, effarr2,
 549                                         mseq1, mseq2,
 550                                         localcopy,
 551                                         gapmap1, gapmap2,
 552                                         effarr1_kozo, effarr2_kozo, effarr_kozo,
 553                                         pair,
 554                                         localhomshrink
 555                                 );
 556                                 return( NULL );
 557                                 pthread_exit( NULL );
 558                         }
 559 //                      pthread_mutex_unlock( targ->mutex );
 560
 561
 562 //                      pthread_mutex_lock( targ->mutex );
 563                         if( *jobposintpt == nbranch )
 564                         {
 565                                 if( *collectingpt != -1 ) *collectingpt = 1; // chofuku
 566                                 pthread_mutex_unlock( targ->mutex );
 567                                 continue;
 568                         }
 569 //                      fprintf( stderr, "JOB jobposintpt=%d\n", *jobposintpt );
 570                         myjob = branchtable[*jobposintpt];
 571                         l = myjob / 2;
 572                         if( l == locnjob-2 ) k = 1;
 573                         else k = myjob - l * 2;
 574 //                      fprintf( stderr, "JOB l=%d, k=%d\n", l, k );
 575                         branchpos = myjob;
 576                         (*jobposintpt)++;
 577                         iterate = *generationofmastercopypt;
 578                         (*ntrypt)++;
 579                         pthread_mutex_unlock( targ->mutex );
 580
 581
 582
 583 //                      fprintf( stderr, "branchpos = %d (thread %d)\n", branchpos, thread_no );
 584
 585 //                      fprintf( stderr, "iterate=%d, l=%d, k=%d (thread %d)\n", iterate, l, k, thread_no );
 586
 587 #if 0
 588                         fprintf( stderr, "STEP %03d-%03d-%d (Thread %d)    ", iterate+1, l+1, k, thread_no );
 589                         fprintf( stderr, "STEP %03d-%03d-%d (thread %d) %s       ", iterate+1, l+1, k, thread_no, use_fft?"\n":"\n" );
 590 #endif
 591                         for( i=0; i<locnjob; i++ ) for( j=0; j<locnjob; j++ ) pair[i][j] = 0;
 592
 593                         OneClusterAndTheOther( locnjob, pair, &s1, &s2, topol, l, k );
 594 #if 0
 595                         fprintf( stderr, "STEP%d-%d\n", l, k );
 596                         for( i=0; i<locnjob; i++ )
 597                         {
 598                                 for( j=0; j<locnjob; j++ )
 599                                 {
 600                                         fprintf( stderr, "%#3d", pair[i][j] );
 601                                 }
 602                                 fprintf( stderr, "\n" );
 603                         }
 604 #endif
 605                         if( !weight )
 606                         {
 607                                 for( i=0; i<locnjob; i++ ) effarr[i] = 1.0;
 608                                 if( nkozo )
 609                                 {
 610                                         for( i=0; i<locnjob; i++ )
 611                                         {
 612                                                 if( kozoarivec[i] )
 613                                                         effarr_kozo[i] = 1.0;
 614                                                 else
 615                                                         effarr_kozo[i] = 0.0;
 616                                         }
 617                                 }
 618                         }
 619                         else if( weight == 4 )
 620                         {
 621                                 weightFromABranch( locnjob, effarr, stopol, topol, l, k );
 622                                 if( nkozo ) // hitomadu single weight.
 623                                 {
 624                                         for( i=0; i<locnjob; i++ )
 625                                         {
 626                                                 if( kozoarivec[i] ) effarr_kozo[i] = effarr[i];
 627                                                 else effarr_kozo[i] = 0.0;
 628                                         }
 629                                 }
 630                         }
 631                         else
 632                         {
 633                                 fprintf( stderr, "weight error!\n" );
 634                                 exit( 1 );
 635                         }
 636
 637                         yarinaoshi:
 638
 639                         pthread_mutex_lock( targ->mutex );
 640                         for( i=0; i<locnjob; i++ ) strcpy( localcopy[i], mastercopy[i] );
 641                         subgenerationatfirst = *subgenerationpt;
 642                         pthread_mutex_unlock( targ->mutex );
 643                         length = strlen( localcopy[0] );
 644
 645                         if( nkozo )
 646                         {
 647                                 double tmptmptmp;
 648                                 tmptmptmp = 0.0;
 649                                 clus1 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s1, localcopy, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
 650                                 for( i=0; i<clus1; i++ ) effarr1_kozo[i] *= 1.0; // 0.5 ga sairyo ?
 651                                 tmptmptmp = 0.0;
 652                                 clus2 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s2, localcopy, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
 653                                 for( i=0; i<clus2; i++ ) effarr2_kozo[i] *= 1.0; // 0.5 ga sairyo ?
 654
 655 #if 0
 656                                 fprintf( stderr, "\ngroup1 = %s\n", indication1 );
 657                                 for( i=0; i<clus1; i++ ) fprintf( stderr, "effarr1_kozo[%d], effarr1[]  = %f, %f\n", i, effarr1_kozo[i], effarr1[i] );
 658                                 fprintf( stderr, "\ngroup2 = %s\n", indication2 );
 659                                 for( i=0; i<clus2; i++ ) fprintf( stderr, "effarr2_kozo[%d], effarr2[]  = %f, %f\n", i, effarr2_kozo[i], effarr2[i] );
 660 #endif
 661                         }
 662                         else
 663                         {
 664                                 clus1 = conjuctionfortbfast( pair, s1, localcopy, mseq1, effarr1, effarr, indication1 );
 665                                 clus2 = conjuctionfortbfast( pair, s2, localcopy, mseq2, effarr2, effarr, indication2 );
 666                         }
 667
 668                 if( rnakozo && rnaprediction == 'm' )
 669                 {
 670                                 makegrouprnait( grouprna1, singlerna, pair, s1 );
 671                                 makegrouprnait( grouprna2, singlerna, pair, s2 );
 672                 }
 673
 674                         if( score_check == 2 )
 675                         {
 676                                 fprintf( stderr, "Score_check 2 is not supported in the multithread version.\n" );
 677                                 exit( 1 );
 678                         }
 679                         else if( score_check )
 680                         {
 681                                 if( RNAscoremtx == 'r' )
 682                                         intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
 683                                 else
 684                                         intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
 685
 686                                 if( constraint )
 687                                 {
 688                                         shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
 689 //                                      weightimportance4( clus1, clus2,  effarr1, effarr2, localhomshrink ); // >>>
 690                                         oimpmatch = 0.0;
 691                                         if( use_fft )
 692                                         {
 693                                                 if( alg == 'Q' )
 694                                                 {
 695                                                         part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 696                                                         if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
 697                                                         for(  i=length-1; i>=0; i-- )
 698                                                         {
 699                                                                 oimpmatch += part_imp_match_out_scQ( i, i );
 700 //                                                              fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
 701                                                         }
 702                                                 }
 703                                                 else
 704                                                 {
 705                                                         part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 706                                                         if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
 707                                                         for(  i=length-1; i>=0; i-- )
 708                                                         {
 709                                                                 oimpmatch += part_imp_match_out_sc( i, i );
 710 //                                                              fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
 711                                                         }
 712                                                 }
 713 //                                              fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
 714                                         }
 715                                         else
 716                                         {
 717                                                 if( alg == 'Q' )
 718                                                 {
 719                                                         imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 720                                                         if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
 721
 722                                                         for(  i=length-1; i>=0; i-- )
 723                                                         {
 724                                                                 oimpmatch += imp_match_out_scQ( i, i );
 725 //                                                              fprintf( stderr, "#### i=%d, initial impmatch = %f\n", i, oimpmatch );
 726                                                         }
 727                                                 }
 728                                                 else
 729                                                 {
 730                                                         imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 731
 732                                                         fprintf( stderr, "not supported\n" );
 733                                                         exit( 1 );
 734
 735                                                         for(  i=length-1; i>=0; i-- )
 736                                                         {
 737                                                                 oimpmatch += imp_match_out_sc( i, i );
 738 //                                                              fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
 739                                                         }
 740                                                 }
 741 //                                              fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
 742                                         }
 743 //                                      fprintf( stderr, "#### initial impmatch = %f\n", oimpmatch );
 744                                 }
 745                                 else
 746                                 {
 747                                         oimpmatch = 0.0;
 748                                 }
 749
 750
 751 //                              fprintf( stderr, "#### tmpdouble = %f\n", tmpdouble );
 752                                 mscore = (double)oimpmatch + tmpdouble;
 753                         }
 754                         else
 755                         {
 756                                 fprintf( stderr, "score_check = %d\n", score_check );
 757                                 fprintf( stderr, "Not supported\n" );
 758                                 exit( 1 );
 759                         }
 760
 761
 762 //                      if( rnakozo ) foldalignedrna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, rnapairboth );
 763
 764 //                      if( !use_fft && !rnakozo )
 765                         if( !use_fft )
 766                         {
 767                                 commongappick_record( clus1, mseq1, gapmap1 );
 768                                 commongappick_record( clus2, mseq2, gapmap2 );
 769                         }
 770
 771 #if 0
 772                         fprintf( stderr, "##### mscore = %f\n", mscore );
 773 #endif
 774
 775 #if DEBUG
 776                         if( !devide )
 777                         {
 778                         fprintf( trap_g, "\n%d-%d-%d\n", iterate+1, l+1, k );
 779                         fprintf( trap_g, "group1 = %s\n", indication1 );
 780                         fprintf( trap_g, "group2 = %s\n", indication2 );
 781                                 fflush( trap_g );
 782                         }
 783
 784 #endif
 785 #if 0
 786                         printf( "STEP %d-%d-%d\n", iterate, l, k );
 787                         for( i=0; i<clus2; i++ ) printf( "%f ", effarr2[i] );
 788                         printf( "\n" );
 789 #endif
 790                         if( constraint == 2 )
 791                         {
 792                                 if( use_fft )
 793                                 {
 794 //                                      if( alg == 'Q' )
 795 //                                              part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 796 //                                      else
 797 //                                              part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 798                                         chudanres = 0;
 799                                         Falign_localhom( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2, subgenerationpt, subgenerationatfirst, &chudanres );
 800 //                                      fprintf( stderr, "##### impmatch = %f\n", impmatch );
 801                                         if( chudanres && parallelizationstrategy == BAATARI2 )
 802                                         {
 803 //                                              fprintf( stderr, "#### yarinaoshi!!! INS-i\n" );
 804                                                 goto yarinaoshi;
 805                                         }
 806                                 }
 807                                 else
 808                                 {
 809                                         if( alg == 'Q' )
 810                                         {
 811                                                 float wm;
 812 //                                              imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap  wo tsukuttakara iranai.
 813 //                                              if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, gapmap1, gapmap2, rnapairboth );
 814
 815                                                 wm = Q__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL, gapmap1, gapmap2 );
 816                                                 fprintf( stderr, "wm = %f\n", wm );
 817 #if 0
 818                                                 fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
 819                                                 naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
 820                                                 fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
 821
 822                                                 if( naivescore1 > naivescore0 )
 823                                                         fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
 824                                                 else if( naivescore1 < naivescore0 )
 825                                                         fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
 826                                                 else
 827                                                         fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 828 #if 0 // chuui
 829                                                 if( abs( wm - naivescore1 ) > 100 )
 830                                                 {
 831 //                                                      fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
 832 //                                                      rewind( stdout );
 833 //                                                      for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 834 //                                                      for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 835 //                                                      exit( 1 );
 836                                                 }
 837 #endif
 838 #endif
 839                                         }
 840                                         else if( alg == 'R' )
 841                                         {
 842                                                 float wm;
 843                                                 imp_match_init_strictR( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
 844                                                 wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
 845 //                                              fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
 846                                                 naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
 847 //                                              fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
 848
 849                                                 if( naivescore1 > naivescore0 )
 850                                                         fprintf( stderr, "%d-%d, ns: %f->%f UP!\n", clus1, clus2, naivescore0, naivescore1 );
 851                                                 else if( naivescore1 < naivescore0 )
 852                                                         fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
 853                                                 else
 854                                                         fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 855 #if 0 // chuui
 856                                                 if( abs( wm - naivescore1 ) > 100 )
 857                                                 {
 858 //                                                      fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
 859                                                         rewind( stdout );
 860                                                         for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 861                                                         for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 862                                                         exit( 1 );
 863                                                 }
 864 #endif
 865                                         }
 866                                         else if( alg == 'H' )
 867                                         {
 868                                                         float wm;
 869                                                 imp_match_init_strictH( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
 870                                                 wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
 871                                                 fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
 872                                                 naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
 873                                                 fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
 874
 875                                                 if( naivescore1 > naivescore0 )
 876                                                         fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
 877                                                 else if( naivescore1 < naivescore0 )
 878                                                         fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
 879                                                 else
 880                                                         fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 881 #if 0 // chuui
 882                                                 if( abs( wm - naivescore1 ) > 100 )
 883                                                 {
 884 //                                                      fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
 885 //                                                      for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 886 //                                                      for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 887 //                                                      exit( 1 );
 888                                                 }
 889 #endif
 890                                         }
 891                                         else
 892                                         {
 893 //                                              imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 894                                                 A__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2 );
 895                                                 fprintf( stderr, "A__align_gapmap\n" );
 896 //                                              fprintf( stderr, "##### impmatch = %f\n", impmatch );
 897                                         }
 898                                 }
 899                         }
 900                         else if( use_fft )
 901                         {
 902                                 float totalwm;
 903                                 chudanres = 0;
 904                                 totalwm = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, subgenerationpt, subgenerationatfirst, &chudanres );
 905                                 if( chudanres && parallelizationstrategy == BAATARI2 )
 906                                 {
 907 //                                      fprintf( stderr, "#### yarinaoshi!!! FFT-NS-i\n" );
 908                                         goto yarinaoshi;
 909                                 }
 910
 911 //                              fprintf( stderr, "totalwm = %f\n", totalwm );
 912 #if 0
 913                                 if( alg == 'Q' )
 914                                 {
 915                                         fprintf( stderr, "totalwm = %f\n", totalwm );
 916                                         naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
 917
 918                                         if( naivescore1 > naivescore0 )
 919                                                 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
 920                                         else if( naivescore1 < naivescore0 )
 921                                                 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
 922                                         else
 923                                                 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 924 #if 1 // chuui
 925                                         if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
 926                                         {
 927 //                                              fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
 928 //                                              for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 929 //                                              for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 930 //                                              exit( 1 );
 931                                         }
 932 #endif
 933                                 }
 934 #endif
 935                                 if( alg == 'R' )
 936                                 {
 937                                         fprintf( stderr, "totalwm = %f\n", totalwm );
 938                                         naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
 939
 940                                         if( naivescore1 > naivescore0 )
 941                                                 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
 942                                         else if( naivescore1 < naivescore0 )
 943                                                 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
 944                                         else
 945                                                 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 946 #if 1 // chuui
 947                                         if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
 948                                         {
 949 //                                              fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
 950 //                                              for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 951 //                                              for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 952 //                                              exit( 1 );
 953                                         }
 954                                 }
 955 #endif
 956                         }
 957                         else
 958                         {
 959                                 if( alg == 'M' )
 960                                 {
 961                                         chudanres = 0;
 962                                         MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL, subgenerationpt, subgenerationatfirst, &chudanres, outgap, outgap );
 963                                         if( chudanres && parallelizationstrategy == BAATARI2 )
 964                                         {
 965 //                                              fprintf( stderr, "#### yarinaoshi!!! NW-NS-i\n" );
 966                                                 goto yarinaoshi;
 967                                         }
 968                                 }
 969                                 else if( alg == 'A' )
 970                                 {
 971                                         A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL, NULL, 0, NULL, 1, 1 ); //outgap==1
 972                                 }
 973                                 else if( alg == 'Q' )
 974                                 {
 975                                         float wm;
 976                                         wm = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
 977                                         fprintf( stderr, "wm = %f\n", wm );
 978                                         fprintf( stderr, "impmatch = %f\n", impmatch );
 979                                         naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
 980
 981                                         if( naivescore1 > naivescore0 )
 982                                                 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
 983                                         else if( naivescore1 < naivescore0 )
 984                                                 fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
 985                                         else
 986                                                 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 987 #if 1 // chuui
 988                                         if( abs( wm - naivescore1 ) > 100 )
 989                                         {
 990 //                                              fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
 991 //                                              rewind( stderr );
 992 //                                              rewind( stdout );
 993 //                                              for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 994 //                                              for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 995 //                                              exit( 1 );
 996                                         }
 997 #endif
 998                                 }
 999                                 else if( alg == 'R' )
1000                                 {
1001                                         float wm;
1002                                         wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
1003                                         naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
1004
1005                                         if( naivescore1 > naivescore0 )
1006                                                 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
1007                                         else if( naivescore1 < naivescore0 )
1008                                                 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
1009                                         else
1010                                                 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
1011 #if 1 // chuui
1012                                         if( abs( wm - naivescore1 ) > 100 )
1013                                         {
1014 //                                              fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
1015 //                                              rewind( stderr );
1016 //                                              rewind( stdout );
1017 //                                              for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
1018 //                                              for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
1019 //                                              exit( 1 );
1020                                         }
1021 #endif
1022                                 }
1023                                 else if( alg == 'H' )
1024                                 {
1025                                         float wm;
1026                                         wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
1027                                         naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
1028
1029                                         if( naivescore1 > naivescore0 )
1030                                                 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
1031                                         else if( naivescore1 < naivescore0 )
1032                                         {
1033                                                 fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
1034                                         }
1035                                         else
1036                                                 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
1037
1038 #if 0 // chuui
1039                                         if( abs( wm - naivescore1 ) > 100 )
1040                                         {
1041                                                 rewind( stdout );
1042                                                 for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
1043                                                 for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
1044                                                 exit( 1 );
1045                                         }
1046 #endif
1047                                 }
1048                                 else if( alg == 'a' )
1049                                 {
1050                                         Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
1051                                 }
1052                                 else ErrorExit( "Sorry!" );
1053                         }
1054 //                      fprintf( stderr, "## impmatch = %f\n", impmatch );
1055
1056 #if 1
1057                         if( parallelizationstrategy == BAATARI2 && *subgenerationpt != subgenerationatfirst )
1058                         {
1059 //                              fprintf( stderr, "\nYarinaoshi2!! (Thread %d)\n", thread_no );
1060                                 goto yarinaoshi;
1061                         }
1062 #endif
1063
1064                         identity = !strcmp( localcopy[s1], mastercopy[s1] );
1065                         identity *= !strcmp( localcopy[s2], mastercopy[s2] );
1066
1067
1068 /* Bug?  : idnetitcal but score change when scoreing mtx != JTT  */
1069
1070                         length = strlen( mseq1[0] );
1071
1072                         if( identity )
1073                         {
1074                                 tscore = mscore;
1075 //                              if( !devide ) fprintf( trap_g, "tscore =  %f   identical.\n", tscore );
1076 //                              fprintf( stderr, " identical." );
1077                                 fprintf( stderr, "%03d-%04d-%d (thread %4d) identical     \r", iterate+1, *ndonept, k, thread_no );
1078                         }
1079                         else
1080                         {
1081                                 if( score_check )
1082                                 {
1083                                         if( constraint == 2 )
1084                                         {
1085 #if 1
1086                                                 if( RNAscoremtx == 'r' )
1087                                                         intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1088                                                 else
1089                                                         intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1090 #else
1091                                                 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1092 #endif
1093
1094                                                 tscore = impmatch + tmpdouble;
1095
1096 //                                              fprintf( stderr, "tmpdouble=%f, impmatch = %f -> %f, tscore = %f\n", tmpdouble, oimpmatch, impmatch, tscore );
1097                                         }
1098                                         else
1099                                         {
1100                                                 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1101                                                 tscore = tmpdouble;
1102                                         }
1103 //                                      fprintf( stderr, "#######ii=%d, iterate=%d score = %f -> %f \n", ii, iterate , mscore, tscore );
1104         #if 0
1105                                         for( i=0; i<1; i++ )
1106                                                 fprintf( stderr, "%s\n", mseq1[i] );
1107                                         fprintf( stderr, "+++++++\n" );
1108                                         for( i=0; i<1; i++ )
1109                                                 fprintf( stderr, "%s\n", mseq2[i] );
1110         #endif
1111
1112                                 }
1113                                 else
1114                                 {
1115                                         tscore = mscore + 1.0;
1116 //                                      tscore = 0.0;
1117 //                                      fprintf( stderr, "in line 705, tscore=%f\n", tscore );
1118 //                                      for( i=0; i<length; i++ )
1119 //                                              tscore = tscore + (double)mseq1[0][i];
1120 //                                      mscore = tscore - 1.0;
1121                                 }
1122
1123                                 if( isnan( mscore ) )
1124                                 {
1125                                         fprintf( stderr, "\n\nmscore became NaN\n" );
1126                                         exit( 1 );
1127                                 }
1128                                 if( isnan( tscore ) )
1129                                 {
1130                                         fprintf( stderr, "\n\ntscore became NaN\n" );
1131                                         exit( 1 );
1132                                 }
1133
1134
1135
1136 //                              fprintf( stderr, "@@@@@ mscore,tscore = %f,%f\n", mscore, tscore );
1137
1138 #if 1
1139                                 if( parallelizationstrategy == BAATARI1 && *subgenerationpt != subgenerationatfirst )
1140                                 {
1141 //                                      fprintf( stderr, "\nYarinaoshi1!! (Thread %d)\n", thread_no );
1142                                         goto yarinaoshi;
1143                                 }
1144 #endif
1145                                 gain = tscore - ( mscore - cut/100.0*mscore );
1146                                 if( gain > 0 )
1147                                 {
1148                                         if( parallelizationstrategy == BESTFIRST )
1149                                         {
1150                                                 if( gain > gainlist[thread_no] )
1151                                                 {
1152                                                         gainlist[thread_no] = gain;
1153                                                         for( i=0; i<locnjob; i++ ) strcpy( candidates[thread_no][i], localcopy[i] );
1154                                                         tscorelist[thread_no] = tscore;
1155 //                                              if( iterate == 0 ) fprintf( stderr, "hist %d-%d-%d, gain=%f (Thread %d)\n", iterate, l, k, gain, thread_no );
1156                                                 }
1157                                         }
1158                                         else
1159                                         {
1160                                                 pthread_mutex_lock( targ->mutex );
1161                                                 for( i=0; i<locnjob; i++ ) strcpy( mastercopy[i], localcopy[i] );
1162                                                 *subgenerationpt += 1;
1163                                                 gainlist[thread_no] = *basegainpt + gain;
1164                                                 *basegainpt += gain;
1165
1166                                                 if( scoreout )
1167                                                 {
1168                                                         unweightedspscore = plainscore( locnjob, localcopy );
1169                                                         fprintf( stderr, "\nSCORE %d = %.0f, ", *ntrypt, unweightedspscore );
1170                                                         fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( locnjob * strlen( mastercopy[0] ) ) );
1171                                                         if( weight || constraint ) fprintf( stderr, " (differs from the objective score)" );
1172                                                         fprintf( stderr, "\n" );
1173                                                 }
1174
1175                                                 pthread_mutex_unlock( targ->mutex );
1176                                                 tscorelist[thread_no] = tscore;
1177                                         }
1178 #if 0
1179                                         fprintf( stderr, "tscore =  %f   mscore = %f  accepted.\n", tscore, mscore );
1180                                         fprintf( stderr, "\nbetter! gain = %f (thread %d)\r", gain, thread_no );
1181 #else
1182                                         fprintf( stderr, "%03d-%04d-%d (thread %4d) better     \r", iterate+1, *ndonept, k, thread_no );
1183 #endif
1184
1185                                 }
1186                                 else
1187                                 {
1188 #if 0
1189                                         fprintf( stderr, "tscore =  %f   mscore = %f  rejected.\r", tscore, mscore );
1190                                         fprintf( stderr, "worse! gain = %f", gain );
1191 #else
1192                                         fprintf( stderr, "%03d-%04d-%d (thread %4d) worse      \r", iterate+1, *ndonept, k, thread_no );
1193 #endif
1194                                         tscore = mscore;
1195                                 }
1196                         }
1197                         converged2 = 0;
1198                         for( ii=iterate-2; ii>=0; ii-=1 )
1199                         {
1200 //                              fprintf( stderr, "Checking tscorehistory %f ?= %f\n", tscore, tscorehistory_detail[ii][branchpos] );
1201                                 if( tscore == tscorehistory_detail[ii][branchpos] )
1202                                 {
1203                                         converged2 = 1;
1204                                         break;
1205                                 }
1206                         }
1207                         if( parallelizationstrategy != BESTFIRST && converged2 )
1208                         {
1209 //                              fprintf( stderr, "\nFINISH!\n" );
1210                                 pthread_mutex_lock( targ->mutex );
1211                                 *finishpt = 1;
1212                                 pthread_mutex_unlock( targ->mutex );
1213                         }
1214
1215                         tscorehistory_detail[iterate][branchpos] = tscore;
1216                         fprintf( stderr, "\r" );
1217
1218                         pthread_mutex_lock( targ->mutex );
1219                         (*ndonept)++;
1220 //                      fprintf( stderr, "*ndonept = %d, nbranch = %d (thread %d) iterate=%d\n", *ndonept, nbranch, thread_no, iterate );
1221                         generationofinput[branchpos] = iterate;
1222                         if( *ndonept == nbranch )
1223                         {
1224                                 if( *collectingpt != -1 ) *collectingpt = 1; // chofuku
1225 //                              fprintf( stderr, "Thread %d sends a signal, *ndonept = %d\n", thread_no, *ndonept );
1226                                 pthread_cond_signal( targ->collection_start );
1227                         }
1228                         pthread_mutex_unlock( targ->mutex );
1229                 }            /* while( 1 ) */
1230
1231         }                  /* for( iterate ) */
1232 //      return( NULL );
1233 }
1234 #endif
1235
1236
1237 int TreeDependentIteration( int locnjob, char **name, int nlen[M],
1238                              char **aseq, char **bseq, int ***topol, double **len,
1239                              int alloclen, LocalHom **localhomtable,
1240                                                          RNApair ***singlerna,
1241                                                          int nkozo, char *kozoarivec )
1242 {
1243         int i, j, k, l, iterate, ii, iu, ju;
1244         int lin, ldf, length;
1245         int clus1, clus2;
1246         int s1, s2;
1247         static double **imanoten;
1248         static Node *stopol;
1249         static double *effarrforlocalhom = NULL;
1250         static double *effarr = NULL;
1251         static double *effarr1 = NULL;
1252         static double *effarr2 = NULL;
1253         static double *effarr_kozo = NULL;
1254         static double *effarr1_kozo = NULL;
1255         static double *effarr2_kozo = NULL;
1256         static double **mtx = NULL;
1257         static int **node = NULL;
1258         static int *branchnode = NULL;
1259         static double **branchWeight = NULL;
1260         static char **mseq1, **mseq2;
1261         static float ***history;
1262         FILE *trap;
1263         double tscore, mscore;
1264         int identity;
1265         int converged;
1266         int oscillating;
1267         float naivescore0 = 0.0; // by D.Mathog, a guess
1268         float naivescore1;
1269 #if 0
1270         char pair[njob][njob];
1271 #else
1272         static char **pair;
1273 #endif
1274 #if DEBUG + RECORD
1275         double score_for_check0, score_for_check1;
1276         static double **effmtx = NULL;
1277         extern double score_calc0();
1278 #endif
1279         static char *indication1, *indication2;
1280         static LocalHom ***localhomshrink = NULL;
1281         float impmatch = 0.0, oimpmatch = 0.0;
1282         static int *gapmap1;
1283         static int *gapmap2;
1284         double tmpdouble;
1285         int intdum;
1286         static RNApair *rnapairboth;
1287         RNApair ***grouprna1, ***grouprna2;
1288         double unweightedspscore;
1289
1290         if( rnakozo && rnaprediction == 'm' )
1291         {
1292                 grouprna1 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
1293                 grouprna2 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
1294         }
1295         else
1296         {
1297                 grouprna1 = grouprna2 = NULL;
1298         }
1299
1300         Writeoptions( trap_g );
1301         fflush( trap_g );
1302
1303         if( effarr == NULL ) /* locnjob == njob ni kagiru */
1304         {
1305                 indication1 = AllocateCharVec( njob*3+50 );
1306                 indication2 = AllocateCharVec( njob*3+50 );
1307                 effarr = AllocateDoubleVec( locnjob );
1308                 effarrforlocalhom = AllocateDoubleVec( locnjob );
1309                 effarr1 = AllocateDoubleVec( locnjob );
1310                 effarr2 = AllocateDoubleVec( locnjob );
1311                 mseq1 = AllocateCharMtx( locnjob, 0 );
1312                 mseq2 = AllocateCharMtx( locnjob, 0 );
1313                 mtx = AllocateDoubleMtx( locnjob, locnjob );
1314                 node = AllocateIntMtx( locnjob, locnjob );
1315                 branchnode = AllocateIntVec( locnjob );
1316                 branchWeight = AllocateDoubleMtx( locnjob, 2 );
1317                 history = AllocateFloatCub( niter, locnjob, 2 );
1318                 stopol = (Node *)calloc( locnjob * 2, sizeof( Node ) );
1319                 gapmap1 = AllocateIntVec( alloclen );
1320                 gapmap2 = AllocateIntVec( alloclen );
1321                 if( score_check == 2 ) imanoten = AllocateDoubleMtx( njob, njob );
1322
1323                 effarr1_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
1324                 effarr2_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
1325                 effarr_kozo = AllocateDoubleVec( locnjob );
1326                 for( i=0; i<locnjob; i++ )
1327                         effarr_kozo[i] = effarr1_kozo[i] = effarr2_kozo[i] = 0.0;
1328
1329 #if 0
1330 #else
1331                 pair = AllocateCharMtx( locnjob, locnjob );
1332                 if( rnakozo ) rnapairboth = (RNApair *)calloc( alloclen, sizeof( RNApair ) );
1333
1334                 if( constraint )
1335                 {
1336                         localhomshrink = (LocalHom ***)calloc( njob, sizeof( LocalHom ** ) );
1337                         for( i=0; i<njob; i++)
1338                         {
1339                                 localhomshrink[i] = (LocalHom **)calloc( njob, sizeof( LocalHom * ) );
1340                         }
1341                 }
1342 #endif
1343         }
1344 #if DEBUG + RECORD
1345         if( !effmtx ) effmtx = AllocateDoubleMtx( locnjob, locnjob );
1346         for( i=0; i<locnjob; i++ ) for( j=0; j<locnjob; j++ ) effmtx[i][j] = 1.0;
1347 #endif
1348
1349         for( i=0; i<locnjob; i++ ) strcpy( bseq[i], aseq[i] );
1350
1351         writePre( locnjob, name, nlen, aseq, 0 );
1352
1353         if( utree )
1354         {
1355                 if( constraint )
1356                 {
1357                         counteff_simple( locnjob, topol, len, effarrforlocalhom );
1358                         calcimportance( locnjob, effarrforlocalhom, aseq, localhomtable );
1359                 }
1360
1361                 if( weight == 2 )
1362                 {
1363                         countnode_int( locnjob, topol, node );
1364                         if( nkozo )
1365                         {
1366                                 fprintf( stderr, "Not supported, weight=%d nkozo=%d.\n", weight, nkozo );
1367                         }
1368                 }
1369                 else if( weight == 4 )
1370                 {
1371                         treeCnv( stopol, locnjob, topol, len, branchWeight );
1372                         calcBranchWeight( branchWeight, locnjob, stopol, topol, len );
1373                 }
1374         }
1375
1376 #ifdef enablemultithread
1377         if( nthread > 0 )
1378         {
1379                 threadarg_t *targ;
1380                 pthread_t *handle;
1381                 pthread_mutex_t mutex;
1382                 pthread_cond_t collection_end;
1383                 pthread_cond_t collection_start;
1384                 int jobposint;
1385                 int generationofmastercopy;
1386                 int subgeneration;
1387                 float basegain;
1388                 int *generationofinput;
1389                 float *gainlist;
1390                 float *tscorelist;
1391                 int ndone;
1392                 int ntry;
1393                 int collecting;
1394                 int nbranch;
1395                 int maxiter;
1396                 char ***candidates;
1397                 int *branchtable;
1398                 float **tscorehistory_detail;
1399                 int finish;
1400
1401                 nwa = nthread + 1;
1402                 nbranch = (njob-1) * 2 - 1;
1403                 maxiter = niter;
1404
1405                 targ = calloc( nwa, sizeof( threadarg_t ) );
1406                 handle = calloc( nwa, sizeof( pthread_t ) );
1407                 pthread_mutex_init( &mutex, NULL );
1408                 pthread_cond_init( &collection_end, NULL );
1409                 pthread_cond_init( &collection_start, NULL );
1410
1411                 gainlist = calloc( nwa, sizeof( float ) );
1412                 tscorelist = calloc( nwa, sizeof( float ) );
1413                 branchtable = calloc( nbranch, sizeof( int ) );
1414                 generationofinput = calloc( nbranch, sizeof( int ) );
1415                 if( parallelizationstrategy == BESTFIRST )
1416                         candidates = AllocateCharCub( nwa, locnjob, alloclen );
1417                 for( i=0; i<nbranch; i++ ) branchtable[i] = i;
1418                 tscorehistory_detail = AllocateFloatMtx( maxiter, nbranch );
1419
1420                 collecting = 1;
1421
1422                 for( i=0; i<nwa; i++ )
1423                 {
1424                         targ[i].thread_no = i;
1425                         targ[i].njob = njob;
1426                         targ[i].nbranch = nbranch;
1427                         targ[i].maxiter = maxiter;
1428                         targ[i].ndonept = &ndone;
1429                         targ[i].ntrypt = &ntry;
1430                         targ[i].collectingpt = &collecting;
1431                         targ[i].jobposintpt = &jobposint;
1432                         targ[i].gainlist = gainlist;
1433                         targ[i].tscorelist = tscorelist;
1434                         targ[i].nkozo = nkozo;
1435                         targ[i].kozoarivec = kozoarivec;
1436                         targ[i].mastercopy = bseq;
1437                         targ[i].candidates = candidates;
1438                         targ[i].subgenerationpt = &subgeneration;
1439                         targ[i].basegainpt = &basegain;
1440                         targ[i].generationofmastercopypt = &generationofmastercopy;
1441                         targ[i].generationofinput = generationofinput;
1442                         targ[i].branchtable = branchtable;
1443                         targ[i].singlerna = singlerna;
1444                         targ[i].localhomtable = localhomtable;
1445                         targ[i].alloclen = alloclen;
1446                         targ[i].stopol = stopol;
1447                         targ[i].topol = topol;
1448 //                      targ[i].len = len;
1449                         targ[i].mutex = &mutex;
1450                         targ[i].collection_end = &collection_end;
1451                         targ[i].collection_start = &collection_start;
1452                         targ[i].tscorehistory_detail = tscorehistory_detail;
1453                         targ[i].finishpt = &finish;
1454
1455                         pthread_create( handle+i, NULL, athread, (void *)(targ+i) );
1456                 }
1457
1458                 for( i=0; i<nwa; i++ )
1459                 {
1460                         pthread_join( handle[i], NULL );
1461                 }
1462
1463                 pthread_mutex_destroy( &mutex );
1464                 pthread_cond_destroy( &collection_end );
1465                 pthread_cond_destroy( &collection_start );
1466
1467                 free( targ );
1468                 free( handle );
1469                 free( gainlist );
1470                 free( tscorelist );
1471                 free( branchtable );
1472                 free( generationofinput );
1473                 if( parallelizationstrategy == BESTFIRST )
1474                         FreeCharCub( candidates );
1475                 FreeFloatMtx( tscorehistory_detail );
1476         }
1477         else
1478 #endif
1479         {
1480 #if 0
1481                 int *branchtable;
1482                 int jobpos;
1483                 int myjob;
1484
1485                 int nbranch;
1486                 nbranch = (njob-1) * 2 - 1;
1487
1488                 branchtable = calloc( nbranch, sizeof( int ) );
1489                 for( i=0; i<nbranch; i++ ) branchtable[i] = i;
1490
1491                 srand( randomseed );
1492 #endif
1493
1494                 if( parallelizationstrategy == BESTFIRST )
1495                 {
1496                         fprintf( stderr, "Not implemented.  Try --thread 1 --bestfirst\n" );
1497                         exit( 1 );
1498                 }
1499                 converged = 0;
1500                 if( cooling ) cut *= 2.0;
1501                 for( iterate = 0; iterate<niter; iterate++ )
1502                 {
1503                         if( cooling ) cut *= 0.5; /* ... */
1504
1505 #if 0
1506                         if( randomseed != 0 ) shuffle( branchtable, nbranch );
1507 #endif
1508
1509                         fprintf( trap_g, "\n" );
1510                         Writeoption2( trap_g, iterate, cut );
1511                         fprintf( trap_g, "\n" );
1512
1513
1514                         if( utree == 0 )
1515                         {
1516                                 if( nkozo )
1517                                 {
1518                                         fprintf( stderr, "The combination of dynamic tree and kozo is not supported.\n" );
1519                                         exit( 1 );
1520                                 }
1521                                 if( devide )
1522                                 {
1523                                         static char *buff1 = NULL;
1524                                         static char *buff2 = NULL;
1525                                         if( !buff1 )
1526                                         {
1527                                                 buff1 = AllocateCharVec( alloclen );
1528                                                 buff2 = AllocateCharVec( alloclen );
1529                                         }
1530
1531                                         for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ )
1532                                         {
1533                                                 buff1[0] = buff2[0] = 0;
1534                                                 strcat( buff1, res_g[i] ); strcat( buff2, res_g[j] );
1535                                                 strcat( buff1,  bseq[i] ); strcat( buff2,  bseq[j] );
1536                                                 strcat( buff1, seq_g[i] ); strcat( buff2, seq_g[j] );
1537
1538                                                 mtx[i][j] = (double)substitution_hosei( buff1, buff2 );
1539                                         }
1540                                 }
1541                                 else
1542                                 {
1543                                         for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ )
1544                                                 mtx[i][j] = (double)substitution_hosei( bseq[i], bseq[j] );
1545                                 }
1546
1547                                 if     ( treemethod == 'n' )
1548                                         nj( locnjob, mtx, topol, len );
1549                                 else if( treemethod == 's' )
1550                                         spg( locnjob, mtx, topol, len );
1551                                 else if( treemethod == 'X' )
1552                                         supg( locnjob, mtx, topol, len );
1553                                 else if( treemethod == 'p' )
1554                                         upg2( locnjob, mtx, topol, len );
1555                                 /* veryfastsupg\e$B$O!":#$N$H$3$m;H$($^$;$s!#\e(B*/
1556                                 /* \e$B=gHV$NLdBj$,$"$k$N$G\e(B                  */
1557
1558                                 if( weight == 2 )
1559                                         countnode_int( locnjob, topol, node );
1560                                 else if( weight == 4 )
1561                                 {
1562                                         treeCnv( stopol, locnjob, topol, len, branchWeight );
1563                                         calcBranchWeight( branchWeight, locnjob, stopol, topol, len );
1564                                 }
1565                                 trap = fopen( "hat2", "w" );
1566                                 if( !trap ) ErrorExit( "Cannot open hat2." );
1567                                 WriteHat2_pointer( trap, locnjob, name, mtx );
1568                                 fclose( trap );
1569                                 if( constraint )
1570                                 {
1571                                         counteff_simple( locnjob, topol, len, effarrforlocalhom );
1572                                         calcimportance( locnjob, effarrforlocalhom, aseq, localhomtable );
1573                                 }
1574                         }
1575
1576                         if( iterate % 2 == 0 )
1577                         {
1578                                 lin = 0; ldf = +1;
1579                         }
1580                         else
1581                         {
1582                                 lin = locnjob - 2; ldf = -1;
1583                         }
1584
1585                         if( score_check == 2 )
1586                         {
1587                                 effarr1[0] = 1.0;
1588                                 effarr2[0] = 1.0;
1589                                 length = strlen( bseq[0] );
1590                                 for( i=0; i<locnjob-1; i++ )
1591                                         for( j=i+1; j<locnjob; j++ )
1592                                                 intergroup_score( bseq+i, bseq+j, effarr1, effarr2, 1, 1, length, imanoten[i]+j );
1593                         }
1594
1595 #if 1
1596                         for( l=lin; l < locnjob-1 && l >= 0 ; l+=ldf )
1597                         {
1598
1599
1600                                 for( k=0; k<2; k++ )
1601                                 {
1602                                         if( l == locnjob-2 ) k = 1;
1603 #else
1604
1605                         for( jobpos=0; jobpos<nbranch; jobpos++)
1606                         {
1607                                 {
1608                                         myjob = branchtable[jobpos];
1609                                         l = myjob / 2;
1610                                         if( l == locnjob-2 ) k = 1;
1611                                         else k = myjob - l * 2;
1612 #endif
1613         #if 1
1614                                         fprintf( stderr, "STEP %03d-%03d-%d ", iterate+1, l+1, k );
1615                                         fflush( stderr );
1616         #else
1617                                         fprintf( stderr, "STEP %03d-%03d-%d %s", iterate+1, l+1, k, use_fft?"\n":"\n" );
1618         #endif
1619                                         for( i=0; i<locnjob; i++ ) for( j=0; j<locnjob; j++ ) pair[i][j] = 0;
1620
1621                                         OneClusterAndTheOther( locnjob, pair, &s1, &s2, topol, l, k );
1622         #if 0
1623                                         fprintf( stderr, "STEP%d-%d\n", l, k );
1624                                         for( i=0; i<locnjob; i++ )
1625                                         {
1626                                                 for( j=0; j<locnjob; j++ )
1627                                                 {
1628                                                         fprintf( stderr, "%#3d", pair[i][j] );
1629                                                 }
1630                                                 fprintf( stderr, "\n" );
1631                                         }
1632         #endif
1633                                         if( !weight )
1634                                         {
1635                                                 for( i=0; i<locnjob; i++ ) effarr[i] = 1.0;
1636                                                 if( nkozo )
1637                                                 {
1638                                                         for( i=0; i<locnjob; i++ )
1639                                                         {
1640                                                                 if( kozoarivec[i] )
1641                                                                         effarr_kozo[i] = 1.0;
1642                                                                 else
1643                                                                         effarr_kozo[i] = 0.0;
1644                                                         }
1645                                                 }
1646                                         }
1647                                         else if( weight == 2 )
1648                                         {
1649                                                 nodeFromABranch( locnjob, branchnode, node, topol, len, l, k );
1650                                                 node_eff( locnjob, effarr, branchnode );
1651                                         }
1652                                         else if( weight == 4 )
1653                                         {
1654                                                 weightFromABranch( locnjob, effarr, stopol, topol, l, k );
1655         #if 0
1656                                                 if( nkozo )
1657                                                 {
1658                                                         assignstrweight( locnjob, effarr_kozo, stopol, topol, l, k, kozoarivec, effarr );
1659                                                 }
1660
1661         #else
1662                                                 if( nkozo ) // hitomadu single weight.
1663                                                         for( i=0; i<locnjob; i++ )
1664                                                         {
1665                                                                 if( kozoarivec[i] ) effarr_kozo[i] = effarr[i];
1666                                                                 else effarr_kozo[i] = 0.0;
1667                                                         }
1668         #endif
1669         #if 0
1670                                                 fprintf( stderr, "\n" );
1671                                                 fprintf( stderr, "effarr_kozo = \n" );
1672                                                 for( i=0; i<locnjob; i++ ) fprintf( stderr, "%5.3f ", effarr_kozo[i] );
1673                                                 fprintf( stderr, "\n" );
1674                                                 fprintf( stderr, "effarr = \n" );
1675                                                 for( i=0; i<locnjob; i++ ) fprintf( stderr, "%5.3f ", effarr[i] );
1676                                                 fprintf( stderr, "\n\n" );
1677         #endif
1678                                         }
1679
1680                                         for( i=0; i<locnjob; i++ ) strcpy( aseq[i], bseq[i] );
1681                                         length = strlen( aseq[0] );
1682
1683                                         if( nkozo )
1684                                         {
1685         #if 1
1686                                                 double tmptmptmp;
1687                                                 tmptmptmp = 0.0;
1688                                                 clus1 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s1, aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
1689                                                 for( i=0; i<clus1; i++ ) effarr1_kozo[i] *= 1.0; // 0.5 ga sairyo ?
1690                                                 tmptmptmp = 0.0;
1691                                                 clus2 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s2, aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
1692                                                 for( i=0; i<clus2; i++ ) effarr2_kozo[i] *= 1.0; // 0.5 ga sairyo ?
1693
1694         #if 0
1695                                                 fprintf( stderr, "\ngroup1 = %s\n", indication1 );
1696                                                 for( i=0; i<clus1; i++ ) fprintf( stderr, "effarr1_kozo[%d], effarr1[]  = %f, %f\n", i, effarr1_kozo[i], effarr1[i] );
1697                                                 fprintf( stderr, "\ngroup2 = %s\n", indication2 );
1698                                                 for( i=0; i<clus2; i++ ) fprintf( stderr, "effarr2_kozo[%d], effarr2[]  = %f, %f\n", i, effarr2_kozo[i], effarr2[i] );
1699         #endif
1700
1701
1702
1703
1704
1705         #else
1706                                                 clus1 = conjuctionfortbfast_kozo_BUG( pair, s1, aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
1707                                                 clus2 = conjuctionfortbfast_kozo_BUG( pair, s2, aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
1708         #endif
1709                                         }
1710                                         else
1711                                         {
1712                                                 clus1 = conjuctionfortbfast( pair, s1, aseq, mseq1, effarr1, effarr, indication1 );
1713                                                 clus2 = conjuctionfortbfast( pair, s2, aseq, mseq2, effarr2, effarr, indication2 );
1714                                         }
1715
1716
1717
1718                                 if( rnakozo && rnaprediction == 'm' )
1719                                 {
1720                                                 makegrouprnait( grouprna1, singlerna, pair, s1 );
1721                                                 makegrouprnait( grouprna2, singlerna, pair, s2 );
1722                                 }
1723
1724                                         if( score_check == 2 )
1725                                         {
1726                                                 if( constraint )
1727                                                 {
1728         //                                              msshrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
1729                                                         shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
1730                                                         oimpmatch = 0.0;
1731                                                         if( use_fft )
1732                                                         {
1733                                                                 if( alg == 'Q' )
1734                                                                 {
1735                                                                         part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1736                                                                         if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1737                                                                         for(  i=length-1; i>=0; i-- ) oimpmatch += part_imp_match_out_scQ( i, i );
1738                                                                 }
1739                                                                 else
1740                                                                 {
1741                                                                         part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1742                                                                         if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1743                                                                         for(  i=length-1; i>=0; i-- ) oimpmatch += part_imp_match_out_sc( i, i );
1744                                                                 }
1745                                                         }
1746                                                         else
1747                                                         {
1748                                                                 if( alg == 'Q' )
1749                                                                 {
1750                                                                         imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1751                                                                         if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1752                                                                         for(  i=length-1; i>=0; i-- ) oimpmatch += imp_match_out_scQ( i, i );
1753                                                                 }
1754                                                                 else
1755                                                                 {
1756                                                                         imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1757                                                                         fprintf( stderr, "not supported\n" );
1758                                                                         exit( 1 );
1759                                                                 }
1760                                                         }
1761         //                                              fprintf( stderr, "### oimpmatch = %f\n", oimpmatch );
1762                                                 }
1763                                                 else
1764                                                 {
1765                                                         oimpmatch = 0.0;
1766                                                 }
1767                                                 tmpdouble = 0.0;
1768                                                 iu=0;
1769                                                 for( i=s1; i<locnjob; i++ )
1770                                                 {
1771                                                         if( !pair[s1][i] ) continue;
1772                                                         ju=0;
1773                                                         for( j=s2; j<locnjob; j++ )
1774                                                         {
1775                                                                 if( !pair[s2][j] ) continue;
1776         //                                                      fprintf( stderr, "i = %d, j = %d, effarr1=%f, effarr2=%f\n", i, j, effarr1[iu], effarr2[ju] );
1777                                                                 tmpdouble += effarr1[iu] * effarr2[ju] * imanoten[MIN(i,j)][MAX(i,j)];
1778                                                                 ju++;
1779                                                         }
1780                                                         iu++;
1781                                                 }
1782                                                 mscore = oimpmatch + tmpdouble;
1783                                         }
1784                                         else if( score_check )
1785                                         {
1786         #if 1
1787                                                 if( RNAscoremtx == 'r' )
1788                                                         intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
1789                                                 else
1790                                                         intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
1791         #else
1792                                                 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
1793         #endif
1794
1795                                                 if( constraint )
1796                                                 {
1797                                                         shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
1798                 //                                      weightimportance4( clus1, clus2,  effarr1, effarr2, localhomshrink ); // >>>
1799                                                         oimpmatch = 0.0;
1800                                                         if( use_fft )
1801                                                         {
1802                                                                 if( alg == 'Q' )
1803                                                                 {
1804                                                                         part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1805                                                                         if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1806                                                                         for(  i=length-1; i>=0; i-- )
1807                                                                         {
1808                                                                                 oimpmatch += part_imp_match_out_scQ( i, i );
1809         //                                                                      fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1810                                                                         }
1811                                                                 }
1812                                                                 else
1813                                                                 {
1814                                                                         part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1815                                                                         if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1816                                                                         for(  i=length-1; i>=0; i-- )
1817                                                                         {
1818                                                                                 oimpmatch += part_imp_match_out_sc( i, i );
1819                 //                                                              fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1820                                                                         }
1821                                                                 }
1822                 //                                              fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
1823                                                         }
1824                                                         else
1825                                                         {
1826                                                                 if( alg == 'Q' )
1827                                                                 {
1828                                                                         imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1829                                                                         if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1830
1831                                                                         for(  i=length-1; i>=0; i-- )
1832                                                                         {
1833                                                                                 oimpmatch += imp_match_out_scQ( i, i );
1834         //                                                                      fprintf( stderr, "#### i=%d, initial impmatch = %f\n", i, oimpmatch );
1835                                                                         }
1836                                                                 }
1837                                                                 else
1838                                                                 {
1839                                                                         imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1840
1841                                                                         fprintf( stderr, "not supported\n" );
1842                                                                         exit( 1 );
1843
1844                                                                         for(  i=length-1; i>=0; i-- )
1845                                                                         {
1846                                                                                 oimpmatch += imp_match_out_sc( i, i );
1847                 //                                                              fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1848                                                                         }
1849                                                                 }
1850                 //                                              fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
1851                                                         }
1852                 //                                      fprintf( stderr, "#### initial impmatch = %f\n", oimpmatch );
1853                                                 }
1854                                                 else
1855                                                 {
1856                                                         oimpmatch = 0.0;
1857                                                 }
1858
1859
1860         //                                      fprintf( stderr, "#### tmpdouble = %f\n", tmpdouble );
1861                                                 mscore = (double)oimpmatch + tmpdouble;
1862                                         }
1863                                         else
1864                                         {
1865         //                                      fprintf( stderr, "score_check = %d\n", score_check );
1866         /* atode kousokuka */
1867                                                 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1868                                                 mscore = tmpdouble;
1869         /* atode kousokuka */
1870
1871                                                 if( constraint )
1872                                                 {
1873                                                         oimpmatch = 0.0;
1874                                                         shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
1875                                                         if( use_fft )
1876                                                         {
1877                                                                 if( alg == 'Q' )
1878                                                                 {
1879                                                                         part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1880                                                                         if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1881                                                                 }
1882                                                                 else
1883                                                                 {
1884                                                                         part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1885                                                                         if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1886                                                                 }
1887                                                         }
1888                                                         else
1889                                                         {
1890                                                                 if( alg == 'Q' )
1891                                                                 {
1892                                                                         imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1893                                                                         if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1894                                                                 }
1895                                                                 else
1896                                                                 {
1897                                                                         imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1898                                                                         fprintf( stderr, "Not supported\n" );
1899                                                                         exit( 1 );
1900                                                                 }
1901                                                         }
1902                                                 }
1903                                         }
1904
1905         //                              oimpmatch = 0.0;
1906                                         if( constraint )
1907                                         {
1908         #if 0 // iranai
1909                                                 if( alg == 'Q' )
1910                                                 {
1911                                                         imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1912                                                         for(  i=length-1; i>=0; i-- )
1913                                                         {
1914                                                                 oimpmatch += imp_match_out_scQ( i, i );
1915         //                                                      fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1916                                                         }
1917                                                 }
1918                                                 else
1919                                                 {
1920                                                         imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1921                                                         for(  i=length-1; i>=0; i-- )
1922                                                         {
1923                                                                 oimpmatch += imp_match_out_sc( i, i );
1924         //                                                      fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1925                                                         }
1926                                                 }
1927         #endif
1928                                         }
1929         #if 0
1930                                         if( alg == 'H' )
1931                                                 naivescore0 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + oimpmatch;
1932                                         else if( alg == 'Q' )
1933                                                 naivescore0 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + oimpmatch;
1934                                         else if( alg == 'R' )
1935                                                 naivescore0 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + oimpmatch;
1936         #endif
1937
1938         //                              if( rnakozo ) foldalignedrna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, rnapairboth );
1939
1940         //                              if( !use_fft && !rnakozo )
1941                                         if( !use_fft )
1942                                         {
1943                                                 commongappick_record( clus1, mseq1, gapmap1 );
1944                                                 commongappick_record( clus2, mseq2, gapmap2 );
1945                                         }
1946
1947         #if 0
1948                                         fprintf( stderr, "##### mscore = %f\n", mscore );
1949         #endif
1950
1951         #if DEBUG
1952                                         if( !devide )
1953                                         {
1954                                         fprintf( trap_g, "\nSTEP%d-%d-%d\n", iterate+1, l+1, k );
1955                                 fprintf( trap_g, "group1 = %s\n", indication1 );
1956                                         fprintf( trap_g, "group2 = %s\n", indication2 );
1957                                                 fflush( trap_g );
1958                                         }
1959
1960         #endif
1961         #if 0
1962                                         printf( "STEP %d-%d-%d\n", iterate, l, k );
1963                                         for( i=0; i<clus2; i++ ) printf( "%f ", effarr2[i] );
1964                                         printf( "\n" );
1965         #endif
1966                                         if( constraint == 2 )
1967                                         {
1968                                                 if( use_fft )
1969                                                 {
1970         //                                              if( alg == 'Q' )
1971         //                                                      part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1972         //                                              else
1973         //                                                      part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1974                                                         Falign_localhom( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2, NULL, 0, NULL );
1975         //                                              fprintf( stderr, "##### impmatch = %f\n", impmatch );
1976                                                 }
1977                                                 else
1978                                                 {
1979                                                         if( alg == 'Q' )
1980                                                         {
1981                                                                 float wm;
1982         //                                                      imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap  wo tsukuttakara iranai.
1983         //                                                      if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, gapmap1, gapmap2, rnapairboth );
1984
1985                                                                 wm = Q__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL, gapmap1, gapmap2 );
1986                                                                 fprintf( stderr, "wm = %f\n", wm );
1987         #if 0
1988                                                                 fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
1989                                                                 naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
1990                                                                 fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
1991
1992                                                                 if( naivescore1 > naivescore0 )
1993                                                                         fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
1994                                                                 else if( naivescore1 < naivescore0 )
1995                                                                         fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
1996                                                                 else
1997                                                                         fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
1998         #if 0 // chuui
1999                                                                 if( abs( wm - naivescore1 ) > 100 )
2000                                                                 {
2001         //                                                              fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
2002         //                                                              rewind( stdout );
2003         //                                                              for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2004         //                                                              for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2005         //                                                              exit( 1 );
2006                                                                 }
2007         #endif
2008         #endif
2009                                                         }
2010                                                         else if( alg == 'R' )
2011                                                         {
2012                                                                 float wm;
2013                                                                 imp_match_init_strictR( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
2014                                                                 wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
2015         //                                                      fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
2016                                                                 naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
2017         //                                                      fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
2018
2019                                                                 if( naivescore1 > naivescore0 )
2020                                                                         fprintf( stderr, "%d-%d, ns: %f->%f UP!\n", clus1, clus2, naivescore0, naivescore1 );
2021                                                                 else if( naivescore1 < naivescore0 )
2022                                                                         fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2023                                                                 else
2024                                                                         fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2025         #if 0 // chuui
2026                                                                 if( abs( wm - naivescore1 ) > 100 )
2027                                                                 {
2028         //                                                              fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
2029                                                                         rewind( stdout );
2030                                                                         for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2031                                                                         for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2032                                                                         exit( 1 );
2033                                                                 }
2034         #endif
2035                                                         }
2036                                                         else if( alg == 'H' )
2037                                                         {
2038                                                                 float wm;
2039                                                                 imp_match_init_strictH( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
2040                                                                 wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
2041                                                                 fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
2042                                                                 naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
2043                                                                 fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
2044
2045                                                                 if( naivescore1 > naivescore0 )
2046                                                                         fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2047                                                                 else if( naivescore1 < naivescore0 )
2048                                                                         fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2049                                                                 else
2050                                                                         fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2051         #if 0 // chuui
2052                                                                 if( abs( wm - naivescore1 ) > 100 )
2053                                                                 {
2054         //                                                              fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
2055         //                                                              for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2056         //                                                              for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2057         //                                                              exit( 1 );
2058                                                                 }
2059         #endif
2060                                                         }
2061                                                         else
2062                                                         {
2063         //                                                      imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
2064                                                                 A__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2 );
2065         //                                                      fprintf( stderr, "##### impmatch = %f\n", impmatch );
2066                                                         }
2067                                                 }
2068                                         }
2069                                         else if( use_fft )
2070                                         {
2071                                                 float totalwm;
2072                                                 totalwm = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, NULL, 0, NULL );
2073
2074         //                                      fprintf( stderr, "totalwm = %f\n", totalwm );
2075         #if 0
2076                                                 if( alg == 'Q' )
2077                                                 {
2078                                                         fprintf( stderr, "totalwm = %f\n", totalwm );
2079                                                         naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2080
2081                                                         if( naivescore1 > naivescore0 )
2082                                                                 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2083                                                         else if( naivescore1 < naivescore0 )
2084                                                                 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2085                                                         else
2086                                                                 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2087         #if 1 // chuui
2088                                                         if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
2089                                                         {
2090         //                                                      fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
2091         //                                                      for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2092         //                                                      for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2093         //                                                      exit( 1 );
2094                                                         }
2095         #endif
2096                                                 }
2097         #endif
2098                                                 if( alg == 'R' )
2099                                                 {
2100                                                         fprintf( stderr, "totalwm = %f\n", totalwm );
2101                                                         naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2102
2103                                                         if( naivescore1 > naivescore0 )
2104                                                                 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2105                                                         else if( naivescore1 < naivescore0 )
2106                                                                 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2107                                                         else
2108                                                                 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2109         #if 1 // chuui
2110                                                         if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
2111                                                         {
2112         //                                                      fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
2113         //                                                      for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2114         //                                                      for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2115         //                                                      exit( 1 );
2116                                                         }
2117                                                 }
2118         #endif
2119                                         }
2120                                         else
2121                                         {
2122                                                 if( alg == 'M' )
2123                                                 {
2124                                                         MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
2125                                                 }
2126                                                 else if( alg == 'A' )
2127                                                 {
2128                                                         A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL, NULL, 0, NULL, 1, 1 ); // outgap==1
2129                                                 }
2130                                                 else if( alg == 'Q' )
2131                                                 {
2132                                                         float wm;
2133                                                         wm = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
2134                                                         fprintf( stderr, "wm = %f\n", wm );
2135                                                         fprintf( stderr, "impmatch = %f\n", impmatch );
2136                                                         naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2137
2138                                                         if( naivescore1 > naivescore0 )
2139                                                                 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2140                                                         else if( naivescore1 < naivescore0 )
2141                                                                 fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
2142                                                         else
2143                                                                 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2144         #if 1 // chuui
2145                                                         if( abs( wm - naivescore1 ) > 100 )
2146                                                         {
2147         //                                                      fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
2148         //                                                      rewind( stderr );
2149         //                                                      rewind( stdout );
2150         //                                                      for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2151         //                                                      for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2152         //                                                      exit( 1 );
2153                                                         }
2154         #endif
2155                                                 }
2156                                                 else if( alg == 'R' )
2157                                                 {
2158                                                         float wm;
2159                                                         wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
2160                                                         naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2161
2162                                                         if( naivescore1 > naivescore0 )
2163                                                                 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2164                                                         else if( naivescore1 < naivescore0 )
2165                                                                 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2166                                                         else
2167                                                                 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2168         #if 1 // chuui
2169                                                         if( abs( wm - naivescore1 ) > 100 )
2170                                                         {
2171         //                                                      fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
2172         //                                                      rewind( stderr );
2173         //                                                      rewind( stdout );
2174         //                                                      for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2175         //                                                      for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2176         //                                                      exit( 1 );
2177                                                         }
2178         #endif
2179                                                 }
2180                                                 else if( alg == 'H' )
2181                                                 {
2182                                                         float wm;
2183                                                         wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
2184                                                         naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2185
2186                                                         if( naivescore1 > naivescore0 )
2187                                                                 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2188                                                         else if( naivescore1 < naivescore0 )
2189                                                         {
2190                                                                 fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
2191                                                         }
2192                                                         else
2193                                                                 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2194
2195         #if 0 // chuui
2196                                                         if( abs( wm - naivescore1 ) > 100 )
2197                                                         {
2198                                                                 rewind( stdout );
2199                                                                 for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2200                                                                 for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2201                                                                 exit( 1 );
2202                                                         }
2203         #endif
2204                                                 }
2205                                                 else if( alg == 'a' )
2206                                                 {
2207                                                         Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
2208                                                 }
2209                                                 else ErrorExit( "Sorry!" );
2210                                         }
2211         //                              fprintf( stderr, "## impmatch = %f\n", impmatch );
2212
2213                                                 if( checkC )
2214                                                 {
2215                                                         extern double DSPscore();
2216                                                         extern double SSPscore();
2217                                                         static double cur;
2218                                                         static double pre;
2219
2220         /*
2221                                                         pre = SSPscore( locnjob, bseq );
2222                                                         cur = SSPscore( locnjob, aseq );
2223         */
2224                                                         pre = DSPscore( locnjob, bseq );
2225                                                         cur = DSPscore( locnjob, aseq );
2226
2227                                                         fprintf( stderr, "Previous Sscore = %f\n", pre );
2228                                                         fprintf( stderr, "Currnet  Sscore = %f\n\n", cur );
2229                                                 }
2230
2231         //                              fprintf( stderr, "## impmatch = %f\n", impmatch );
2232                                         identity = !strcmp( aseq[s1], bseq[s1] );
2233                                         identity *= !strcmp( aseq[s2], bseq[s2] );
2234
2235
2236         /* Bug?  : idnetitcal but score change when scoreing mtx != JTT  */
2237
2238                                         length = strlen( mseq1[0] );
2239
2240                                         if( identity )
2241                                         {
2242                                                 tscore = mscore;
2243                                                 if( !devide ) fprintf( trap_g, "tscore =  %f   identical.\n", tscore );
2244                                                 fprintf( stderr, " identical." );
2245                                                 converged++;
2246                                         }
2247                                         else
2248                                         {
2249                                                 if( score_check )
2250                                                 {
2251                                                         if( constraint == 2 )
2252                                                         {
2253         #if 1
2254                                                                 if( RNAscoremtx == 'r' )
2255                                                                         intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
2256                                                                 else
2257                                                                         intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
2258         #else
2259                                                                 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
2260         #endif
2261
2262                                                                 tscore = impmatch + tmpdouble;
2263
2264         //                                                      fprintf( stderr, "tmpdouble=%f, impmatch = %f -> %f, tscore = %f\n", tmpdouble, oimpmatch, impmatch, tscore );
2265                                                         }
2266                                                         else
2267                                                         {
2268                                                                 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
2269                                                                 tscore = tmpdouble;
2270                                                         }
2271         //                                              fprintf( stderr, "#######ii=%d, iterate=%d score = %f -> %f \n", ii, iterate , mscore, tscore );
2272                 #if 0
2273                                                         for( i=0; i<1; i++ )
2274                                                                 fprintf( stderr, "%s\n", mseq1[i] );
2275                                                         fprintf( stderr, "+++++++\n" );
2276                                                         for( i=0; i<1; i++ )
2277                                                                 fprintf( stderr, "%s\n", mseq2[i] );
2278                 #endif
2279
2280                                                 }
2281                                                 else
2282                                                 {
2283                                                         tscore = mscore + 1.0;
2284         //                                              tscore = 0.0;
2285         //                                              fprintf( stderr, "in line 705, tscore=%f\n", tscore );
2286         //                                              for( i=0; i<length; i++ )
2287         //                                                      tscore = tscore + (double)mseq1[0][i];
2288         //                                              mscore = tscore - 1.0;
2289                                                 }
2290
2291                                                 if( isnan( mscore ) )
2292                                                 {
2293                                                         fprintf( stderr, "\n\nmscore became NaN\n" );
2294                                                         exit( 1 );
2295                                                 }
2296                                                 if( isnan( tscore ) )
2297                                                 {
2298                                                         fprintf( stderr, "\n\ntscore became NaN\n" );
2299                                                         exit( 1 );
2300                                                 }
2301
2302
2303
2304         //                                      fprintf( stderr, "@@@@@ mscore,tscore = %f,%f\n", mscore, tscore );
2305
2306                                                 if( tscore > mscore - cut/100.0*mscore )
2307                                                 {
2308                                                         writePre( locnjob, name, nlen, aseq, 0 );
2309                                                         for( i=0; i<locnjob; i++ ) strcpy( bseq[i], aseq[i] );
2310                                                         if( score_check == 2 )
2311                                                         {
2312                                                                 effarr1[0] = 1.0;
2313                                                                 effarr2[0] = 1.0;
2314                                                                 for( i=0; i<locnjob-1; i++ )
2315                                                                         for( j=i+1; j<locnjob; j++ )
2316                                                                                 intergroup_score( bseq+i, bseq+j, effarr1, effarr2, 1, 1, length, imanoten[i]+j );
2317                                                         }
2318
2319         #if 0
2320                                                         fprintf( stderr, "tscore =  %f   mscore = %f  accepted.\n", tscore, mscore );
2321         #endif
2322                                                         fprintf( stderr, " accepted." );
2323                                                         converged = 0;
2324
2325                                                 }
2326                                                 else
2327                                                 {
2328         #if 0
2329                                                         fprintf( stderr, "tscore =  %f   mscore = %f  rejected.\n", tscore, mscore );
2330         #endif
2331                                                         fprintf( stderr, " rejected." );
2332                                                         tscore = mscore;
2333                                                         converged++;
2334                                                 }
2335                                         }
2336                                         fprintf( stderr, "\r" );
2337
2338
2339                                         history[iterate][l][k] = (float)tscore;
2340
2341         //                              fprintf( stderr, "tscore = %f\n", tscore );
2342
2343                                         if( converged >= locnjob * 2 )
2344                                         {
2345                                                 fprintf( trap_g, "Converged.\n\n" );
2346                                                 fprintf( stderr, "\nConverged.\n\n" );
2347                                                 if( scoreout )
2348                                                 {
2349                                                         unweightedspscore = plainscore( njob, bseq );
2350                                                         fprintf( stderr, "\nSCORE %d = %.0f, ", iterate * ( (njob-1)*2-1 ), unweightedspscore );
2351                                                         fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( njob * strlen( bseq[0] ) ) );
2352                                                         if( weight || constraint ) fprintf( stderr, " (differs from the objective score)" );
2353                                                         fprintf( stderr, "\n\n" );
2354                                                 }
2355                                                 return( 0 );
2356                                         }
2357                                         if( iterate >= 1 )
2358                                         {
2359                 /*   oscillation check    */
2360                                                 oscillating = 0;
2361                                                 for( ii=iterate-2; ii>=0; ii-=2 )
2362                                                 {
2363                                                         if( (float)tscore == history[ii][l][k] )
2364                                                         {
2365                                                                 oscillating = 1;
2366                                                                 break;
2367                                                         }
2368                                                 }
2369                                                 if( ( oscillating && !cooling ) || ( oscillating && cut < 0.001 && cooling ) )
2370                                                 {
2371                                                         fprintf( trap_g, "Oscillating.\n" );
2372                                                         fprintf( stderr, "\nOscillating.\n\n" );
2373                                                         if( scoreout )
2374                                                         {
2375                                                                 unweightedspscore = plainscore( njob, bseq );
2376                                                                 fprintf( stderr, "\nSCORE %d = %.0f, ", iterate * ( (njob-1)*2-1 ), unweightedspscore );
2377                                                                 fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( njob * strlen( bseq[0] ) ) );
2378                                                                 if( weight || constraint ) fprintf( stderr, " (differs from the objective score)" );
2379                                                                 fprintf( stderr, "\n\n" );
2380                                                         }
2381         #if 1 /* hujuubun */
2382                                                         return( -1 );
2383         #endif
2384                                                 }
2385                                         }      /* if( iterate ) */
2386                                 }          /* for( k ) */
2387                         }              /* for( l ) */
2388                         if( scoreout )
2389                         {
2390                                 unweightedspscore = plainscore( njob, bseq );
2391                                 fprintf( stderr, "\nSCORE %d = %.0f, ", iterate * ( (njob-1)*2-1 ), unweightedspscore );
2392                                 fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( njob * strlen( bseq[0] ) ) );
2393                                 if( weight || constraint ) fprintf( stderr, " (differs from the objective score)" );
2394                                 fprintf( stderr, "\n\n" );
2395                         }
2396                 }                  /* for( iterate ) */
2397         }
2398         return( 2 );
2399 }                          /* int Tree... */