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 #define DEBUG 0
  12 #define RECORD 0
  13
  14 extern char **seq_g;
  15 extern char **res_g;
  16
  17 static void Writeoption2( FILE *fp, int cycle, double cut )
  18 {
  19         fprintf( fp, "%dth cycle\n", cycle );
  20     fprintf( fp, "marginal score to search : current score * (100-%d) / 100\n", (int)cut );
  21 }
  22
  23 static void Writeoptions( FILE *fp )
  24 {
  25         fprintf( fp, "Tree-dependent-iteration\n" );
  26     if( scoremtx == 1 )
  27         fprintf( fp, "Blosum %d\n", nblosum );
  28     else if( scoremtx == -1 )
  29         fprintf( fp, "DNA\n" );
  30     else if( scoremtx == 2 )
  31         fprintf( fp, "Miyata-Yasunaga\n" );
  32         else
  33                 fprintf( fp, "JTT %dPAM\n", pamN );
  34
  35         if( scoremtx == 0 || scoremtx == 1 )
  36         fprintf( fp, "Gap Penalty = %+5.3f, %5.2f, %+5.3f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
  37         else
  38                 fprintf( fp, "Gap Penalty = %+5.3f\n", (double)penalty/1000 );
  39
  40
  41     if( scmtd == 3 )
  42         fprintf( fp, "score of rnd or sco\n" );
  43     else if( scmtd == 4 )
  44         fprintf( fp, "score = sigma( score for a pair of homologous amino acids ) / ( number of amino acids pairs )\n" );
  45     else if( scmtd == 5 )
  46         fprintf( fp, "score : SP\n" );
  47     if( mix )
  48         fprintf( fp, "?\n" );
  49     else
  50     {
  51         if( weight == 2 )
  52             fprintf( fp, "weighted rationale-1,  geta2 = %f\n", geta2 );
  53         else if( weight == 3 )
  54             fprintf( fp, "weighted like ClustalW," );
  55         else if( weight == 4 )
  56             fprintf( fp, "weighted rationale-2,  geta2 = %f\n", geta2 );
  57         else
  58             fprintf( fp, "unweighted\n" );
  59     }
  60     if( weight && utree )
  61         fprintf( fp, "using tree defined by the file hat2.\n" );
  62     if( weight && !utree )
  63         fprintf( fp, "using temporary tree.\n" );
  64
  65         if( treemethod == 'n' )
  66                 fprintf( fp, "Tree is calculated with Neighbor-Joining Method.\n" );
  67         else if( treemethod == 'q' )
  68                 fprintf( fp, "Tree is calculated with Minimum linkage.\n" );
  69         else if( treemethod == 'X' )
  70                 fprintf( fp, "Tree is calculated with simplified UPG Method and UPG Method.\n" );
  71         else if( treemethod == 'E' )
  72                 fprintf( fp, "Tree is calculated with UPG Method.\n" );
  73         else
  74                 fprintf( fp, "Tree is calculated with unknown Method.\n" );
  75
  76         if( alg == 'C' )
  77                 fprintf( fp, "Algorithm A+ / C\n" );
  78         else if( alg == 'S' )
  79                 fprintf( fp, "Algorithm S \n" );
  80         else if( alg == 'A' )
  81                 fprintf( fp, "Algorithm A+ \n" );
  82         else if( alg == 'a' )
  83                 fprintf( fp, "Algorithm A \n" );
  84         else
  85                 fprintf( fp, "Algorithm ? \n" );
  86
  87         if( use_fft )
  88         {
  89                 if( scoremtx == -1 )
  90                 {
  91                         fprintf( fp, "Basis : 4 nucleotides\n" );
  92                 }
  93                 else
  94                 {
  95                         if( fftscore )
  96                                 fprintf( fp, "Basis : Polarity and Volume\n" );
  97                         else
  98                                 fprintf( fp, "Basis : 20 amino acids\n" );
  99                 }
 100                 fprintf( fp, "Threshold   of anchors = %d%%\n", fftThreshold );
 101                 fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
 102         }
 103 }
 104
 105
 106 int TreeDependentIteration( int locnjob, char name[M][B], int nlen[M],
 107                              char **aseq, char **bseq, int ***topol, double **len,
 108                              int alloclen, LocalHom **localhomtable,
 109                                                          RNApair ***singlerna,
 110                                                          int nkozo, char *kozoarivec )
 111 {
 112         int i, j, k, l, iterate, ii, iu, ju;
 113         int lin, ldf, length;
 114         int clus1, clus2;
 115         int s1, s2;
 116         static double **imanoten;
 117         static Node *stopol;
 118         static double *effarrforlocalhom = NULL;
 119         static double *effarr = NULL;
 120         static double *effarr1 = NULL;
 121         static double *effarr2 = NULL;
 122         static double *effarr_kozo = NULL;
 123         static double *effarr1_kozo = NULL;
 124         static double *effarr2_kozo = NULL;
 125         static double **mtx = NULL;
 126         static int **node = NULL;
 127         static int *branchnode = NULL;
 128         static double **branchWeight = NULL;
 129         static char **mseq1, **mseq2;
 130         static float ***history;
 131         FILE *trap;
 132         double tscore, mscore;
 133         int identity;
 134         int converged;
 135         int oscillating;
 136         float naivescore0 = 0.0; // by D.Mathog, a guess
 137         float naivescore1;
 138 #if 0
 139         char pair[njob][njob];
 140 #else
 141         static char **pair;
 142 #endif
 143 #if DEBUG + RECORD
 144         double score_for_check0, score_for_check1;
 145         static double **effmtx = NULL;
 146         extern double score_calc0();
 147 #endif
 148         static char *indication1, *indication2;
 149         static LocalHom ***localhomshrink = NULL;
 150         float impmatch, oimpmatch;
 151         static int *gapmap1;
 152         static int *gapmap2;
 153         double tmpdouble;
 154         int intdum;
 155         static RNApair *rnapairboth;
 156         RNApair ***grouprna1, ***grouprna2;
 157
 158         if( rnakozo && rnaprediction == 'm' )
 159         {
 160                 grouprna1 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
 161                 grouprna2 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
 162         }
 163         else
 164         {
 165                 grouprna1 = grouprna2 = NULL;
 166         }
 167
 168         Writeoptions( trap_g );
 169         fflush( trap_g );
 170
 171         if( effarr == NULL ) /* locnjob == njob ni kagiru */
 172         {
 173                 indication1 = AllocateCharVec( njob*3+50 );
 174                 indication2 = AllocateCharVec( njob*3+50 );
 175                 effarr = AllocateDoubleVec( locnjob );
 176                 effarrforlocalhom = AllocateDoubleVec( locnjob );
 177                 effarr1 = AllocateDoubleVec( locnjob );
 178                 effarr2 = AllocateDoubleVec( locnjob );
 179                 mseq1 = AllocateCharMtx( locnjob, 0 );
 180                 mseq2 = AllocateCharMtx( locnjob, 0 );
 181                 mtx = AllocateDoubleMtx( locnjob, locnjob );
 182                 node = AllocateIntMtx( locnjob, locnjob );
 183                 branchnode = AllocateIntVec( locnjob );
 184                 branchWeight = AllocateDoubleMtx( locnjob, 2 );
 185                 history = AllocateFloatCub( niter, locnjob, 2 );
 186                 stopol = (Node *)calloc( locnjob * 2, sizeof( Node ) );
 187                 gapmap1 = AllocateIntVec( alloclen );
 188                 gapmap2 = AllocateIntVec( alloclen );
 189                 if( score_check == 2 ) imanoten = AllocateDoubleMtx( njob, njob );
 190
 191                 effarr1_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
 192                 effarr2_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
 193                 effarr_kozo = AllocateDoubleVec( locnjob );
 194                 for( i=0; i<locnjob; i++ )
 195                         effarr_kozo[i] = effarr1_kozo[i] = effarr2_kozo[i] = 0.0;
 196
 197 #if 0
 198 #else
 199                 pair = AllocateCharMtx( locnjob, locnjob );
 200                 if( rnakozo ) rnapairboth = (RNApair *)calloc( alloclen, sizeof( RNApair ) );
 201
 202                 if( constraint )
 203                 {
 204                         localhomshrink = (LocalHom ***)calloc( njob, sizeof( LocalHom ** ) );
 205                         for( i=0; i<njob; i++)
 206                         {
 207                                 localhomshrink[i] = (LocalHom **)calloc( njob, sizeof( LocalHom * ) );
 208                         }
 209                 }
 210 #endif
 211         }
 212 #if DEBUG + RECORD
 213         if( !effmtx ) effmtx = AllocateDoubleMtx( locnjob, locnjob );
 214         for( i=0; i<locnjob; i++ ) for( j=0; j<locnjob; j++ ) effmtx[i][j] = 1.0;
 215 #endif
 216
 217         for( i=0; i<locnjob; i++ ) strcpy( bseq[i], aseq[i] );
 218
 219         writePre( locnjob, name, nlen, aseq, 0 );
 220
 221         if( utree )
 222         {
 223                 if( constraint )
 224                 {
 225                         counteff_simple( locnjob, topol, len, effarrforlocalhom );
 226                         calcimportance( locnjob, effarrforlocalhom, aseq, localhomtable );
 227                 }
 228
 229                 if( weight == 2 )
 230                 {
 231                         countnode_int( locnjob, topol, node );
 232                         if( nkozo )
 233                         {
 234                                 fprintf( stderr, "Not supported, weight=%d nkozo=%d.\n", weight, nkozo );
 235                         }
 236                 }
 237                 else if( weight == 4 )
 238                 {
 239                         treeCnv( stopol, locnjob, topol, len, branchWeight );
 240                         calcBranchWeight( branchWeight, locnjob, stopol, topol, len );
 241                 }
 242         }
 243
 244         converged = 0;
 245         if( cooling ) cut *= 2.0;
 246         for( iterate = 0; iterate<niter; iterate++ )
 247         {
 248                 if( cooling ) cut *= 0.5; /* ... */
 249
 250                 fprintf( trap_g, "\n" );
 251                 Writeoption2( trap_g, iterate, cut );
 252                 fprintf( trap_g, "\n" );
 253
 254                 if( utree == 0 )
 255                 {
 256                         if( nkozo )
 257                         {
 258                                 fprintf( stderr, "Dynamic tree and kozo is not supported.\n" );
 259                                 exit( 1 );
 260                         }
 261                         if( devide )
 262                         {
 263                                 static char *buff1 = NULL;
 264                                 static char *buff2 = NULL;
 265                                 if( !buff1 )
 266                                 {
 267                                         buff1 = AllocateCharVec( alloclen );
 268                                         buff2 = AllocateCharVec( alloclen );
 269                                 }
 270
 271                                 for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ )
 272                                 {
 273                                         buff1[0] = buff2[0] = 0;
 274                                         strcat( buff1, res_g[i] ); strcat( buff2, res_g[j] );
 275                                         strcat( buff1,  bseq[i] ); strcat( buff2,  bseq[j] );
 276                                         strcat( buff1, seq_g[i] ); strcat( buff2, seq_g[j] );
 277
 278                                         mtx[i][j] = (double)substitution_hosei( buff1, buff2 );
 279                                 }
 280                         }
 281                         else
 282                         {
 283                                 for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ )
 284                                         mtx[i][j] = (double)substitution_hosei( bseq[i], bseq[j] );
 285                         }
 286
 287                         if     ( treemethod == 'n' )
 288                                 nj( locnjob, mtx, topol, len );
 289                         else if( treemethod == 's' )
 290                                 spg( locnjob, mtx, topol, len );
 291                         else if( treemethod == 'X' )
 292                                 supg( locnjob, mtx, topol, len );
 293                         else if( treemethod == 'p' )
 294                                 upg2( locnjob, mtx, topol, len );
 295                         /* veryfastsupg\e$B$O!":#$N$H$3$m;H$($^$;$s!#\e(B*/
 296                         /* \e$B=gHV$NLdBj$,$"$k$N$G\e(B                  */
 297
 298                         if( weight == 2 )
 299                                 countnode_int( locnjob, topol, node );
 300                         else if( weight == 4 )
 301                         {
 302                                 treeCnv( stopol, locnjob, topol, len, branchWeight );
 303                                 calcBranchWeight( branchWeight, locnjob, stopol, topol, len );
 304                         }
 305                         trap = fopen( "hat2", "w" );
 306                         if( !trap ) ErrorExit( "Cannot open hat2." );
 307                         WriteHat2( trap, locnjob, name, mtx );
 308                         fclose( trap );
 309                         if( constraint )
 310                         {
 311                                 counteff_simple( locnjob, topol, len, effarrforlocalhom );
 312                                 calcimportance( locnjob, effarrforlocalhom, aseq, localhomtable );
 313                         }
 314                 }
 315
 316                 if( iterate % 2 == 0 )
 317                 {
 318                         lin = 0; ldf = +1;
 319                 }
 320                 else
 321                 {
 322                         lin = locnjob - 2; ldf = -1;
 323                 }
 324
 325                 if( score_check == 2 )
 326                 {
 327                         effarr1[0] = 1.0;
 328                         effarr2[0] = 1.0;
 329                         length = strlen( bseq[0] );
 330                         for( i=0; i<locnjob-1; i++ )
 331                                 for( j=i+1; j<locnjob; j++ )
 332                                         intergroup_score( bseq+i, bseq+j, effarr1, effarr2, 1, 1, length, imanoten[i]+j );
 333                 }
 334
 335                 for( l=lin; l < locnjob-1 && l >= 0 ; l+=ldf )
 336                 {
 337                         for( k=0; k<2; k++ )
 338                         {
 339                                 if( l == locnjob-2 ) k = 1;
 340 #if 1
 341                                 fprintf( stderr, "STEP %03d-%03d-%d ", iterate+1, l+1, k );
 342 #else
 343                                 fprintf( stderr, "STEP %03d-%03d-%d %s", iterate+1, l+1, k, use_fft?"\n":"\n" );
 344 #endif
 345                                 for( i=0; i<locnjob; i++ ) for( j=0; j<locnjob; j++ ) pair[i][j] = 0;
 346
 347                                 OneClusterAndTheOther( locnjob, pair, &s1, &s2, topol, l, k );
 348 #if 0
 349                                 fprintf( stderr, "STEP%d-%d\n", l, k );
 350                                 for( i=0; i<locnjob; i++ )
 351                                 {
 352                                         for( j=0; j<locnjob; j++ )
 353                                         {
 354                                                 fprintf( stderr, "%#3d", pair[i][j] );
 355                                         }
 356                                         fprintf( stderr, "\n" );
 357                                 }
 358 #endif
 359                                 if( !weight )
 360                                 {
 361                                         for( i=0; i<locnjob; i++ ) effarr[i] = 1.0;
 362                                         if( nkozo )
 363                                         {
 364                                                 for( i=0; i<locnjob; i++ )
 365                                                 {
 366                                                         if( kozoarivec[i] )
 367                                                                 effarr_kozo[i] = 1.0;
 368                                                         else
 369                                                                 effarr_kozo[i] = 0.0;
 370                                                 }
 371                                         }
 372                                 }
 373                                 else if( weight == 2 )
 374                                 {
 375                                         nodeFromABranch( locnjob, branchnode, node, topol, len, l, k );
 376                                         node_eff( locnjob, effarr, branchnode );
 377                                 }
 378                                 else if( weight == 4 )
 379                                 {
 380                                         weightFromABranch( locnjob, effarr, stopol, topol, l, k );
 381 #if 0
 382                                         if( nkozo )
 383                                         {
 384                                                 assignstrweight( locnjob, effarr_kozo, stopol, topol, l, k, kozoarivec, effarr );
 385                                         }
 386
 387 #else
 388                                         if( nkozo ) // hitomadu single weight.
 389                                                 for( i=0; i<locnjob; i++ )
 390                                                 {
 391                                                         if( kozoarivec[i] ) effarr_kozo[i] = effarr[i];
 392                                                         else effarr_kozo[i] = 0.0;
 393                                                 }
 394 #endif
 395 #if 0
 396                                         fprintf( stderr, "\n" );
 397                                         fprintf( stderr, "effarr_kozo = \n" );
 398                                         for( i=0; i<locnjob; i++ ) fprintf( stderr, "%5.3f ", effarr_kozo[i] );
 399                                         fprintf( stderr, "\n" );
 400                                         fprintf( stderr, "effarr = \n" );
 401                                         for( i=0; i<locnjob; i++ ) fprintf( stderr, "%5.3f ", effarr[i] );
 402                                         fprintf( stderr, "\n\n" );
 403 #endif
 404                                 }
 405
 406                                 for( i=0; i<locnjob; i++ ) strcpy( aseq[i], bseq[i] );
 407                                 length = strlen( aseq[0] );
 408
 409                                 if( nkozo )
 410                                 {
 411 #if 1
 412                                         double tmptmptmp;
 413                                         tmptmptmp = 0.0;
 414                                         clus1 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s1, aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
 415                                         for( i=0; i<clus1; i++ ) effarr1_kozo[i] *= 1.0; // 0.5 ga sairyo ?
 416                                         tmptmptmp = 0.0;
 417                                         clus2 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s2, aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
 418                                         for( i=0; i<clus2; i++ ) effarr2_kozo[i] *= 1.0; // 0.5 ga sairyo ?
 419
 420 #if 0
 421                                         fprintf( stderr, "\ngroup1 = %s\n", indication1 );
 422                                         for( i=0; i<clus1; i++ ) fprintf( stderr, "effarr1_kozo[%d], effarr1[]  = %f, %f\n", i, effarr1_kozo[i], effarr1[i] );
 423                                         fprintf( stderr, "\ngroup2 = %s\n", indication2 );
 424                                         for( i=0; i<clus2; i++ ) fprintf( stderr, "effarr2_kozo[%d], effarr2[]  = %f, %f\n", i, effarr2_kozo[i], effarr2[i] );
 425 #endif
 426
 427
 428
 429
 430
 431 #else
 432                                         clus1 = conjuctionfortbfast_kozo_BUG( pair, s1, aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
 433                                         clus2 = conjuctionfortbfast_kozo_BUG( pair, s2, aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
 434 #endif
 435                                 }
 436                                 else
 437                                 {
 438                                         clus1 = conjuctionfortbfast( pair, s1, aseq, mseq1, effarr1, effarr, indication1 );
 439                                         clus2 = conjuctionfortbfast( pair, s2, aseq, mseq2, effarr2, effarr, indication2 );
 440                                 }
 441
 442
 443
 444                         if( rnakozo && rnaprediction == 'm' )
 445                         {
 446                                         makegrouprnait( grouprna1, singlerna, pair, s1 );
 447                                         makegrouprnait( grouprna2, singlerna, pair, s2 );
 448                         }
 449
 450                                 if( score_check == 2 )
 451                                 {
 452                                         if( constraint )
 453                                         {
 454 //                                              msshrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
 455                                                 shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
 456                                                 oimpmatch = 0.0;
 457                                                 if( use_fft )
 458                                                 {
 459                                                         if( alg == 'Q' )
 460                                                         {
 461                                                                 part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 462                                                                 if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
 463                                                                 for(  i=length-1; i>=0; i-- ) oimpmatch += part_imp_match_out_scQ( i, i );
 464                                                         }
 465                                                         else
 466                                                         {
 467                                                                 part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 468                                                                 if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
 469                                                                 for(  i=length-1; i>=0; i-- ) oimpmatch += part_imp_match_out_sc( i, i );
 470                                                         }
 471                                                 }
 472                                                 else
 473                                                 {
 474                                                         if( alg == 'Q' )
 475                                                         {
 476                                                                 imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 477                                                                 if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
 478                                                                 for(  i=length-1; i>=0; i-- ) oimpmatch += imp_match_out_scQ( i, i );
 479                                                         }
 480                                                         else
 481                                                         {
 482                                                                 imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 483                                                                 fprintf( stderr, "not supported\n" );
 484                                                                 exit( 1 );
 485                                                         }
 486                                                 }
 487 //                                              fprintf( stderr, "### oimpmatch = %f\n", oimpmatch );
 488                                         }
 489                                         else
 490                                         {
 491                                                 oimpmatch = 0.0;
 492                                         }
 493                                         tmpdouble = 0.0;
 494                                         iu=0;
 495                                         for( i=s1; i<locnjob; i++ )
 496                                         {
 497                                                 if( !pair[s1][i] ) continue;
 498                                                 ju=0;
 499                                                 for( j=s2; j<locnjob; j++ )
 500                                                 {
 501                                                         if( !pair[s2][j] ) continue;
 502 //                                                      fprintf( stderr, "i = %d, j = %d, effarr1=%f, effarr2=%f\n", i, j, effarr1[iu], effarr2[ju] );
 503                                                         tmpdouble += effarr1[iu] * effarr2[ju] * imanoten[MIN(i,j)][MAX(i,j)];
 504                                                         ju++;
 505                                                 }
 506                                                 iu++;
 507                                         }
 508                                         mscore = oimpmatch + tmpdouble;
 509                                 }
 510                                 else if( score_check )
 511                                 {
 512 #if 1
 513                                         if( RNAscoremtx == 'r' )
 514                                                 intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
 515                                         else
 516                                                 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
 517 #else
 518                                         intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
 519 #endif
 520
 521                                         if( constraint )
 522                                         {
 523                                                 shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
 524         //                                      weightimportance4( clus1, clus2,  effarr1, effarr2, localhomshrink ); // >>>
 525                                                 oimpmatch = 0.0;
 526                                                 if( use_fft )
 527                                                 {
 528                                                         if( alg == 'Q' )
 529                                                         {
 530                                                                 part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 531                                                                 if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
 532                                                                 for(  i=length-1; i>=0; i-- )
 533                                                                 {
 534                                                                         oimpmatch += part_imp_match_out_scQ( i, i );
 535 //                                                                      fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
 536                                                                 }
 537                                                         }
 538                                                         else
 539                                                         {
 540                                                                 part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 541                                                                 if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
 542                                                                 for(  i=length-1; i>=0; i-- )
 543                                                                 {
 544                                                                         oimpmatch += part_imp_match_out_sc( i, i );
 545         //                                                              fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
 546                                                                 }
 547                                                         }
 548         //                                              fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
 549                                                 }
 550                                                 else
 551                                                 {
 552                                                         if( alg == 'Q' )
 553                                                         {
 554                                                                 imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 555                                                                 if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
 556
 557                                                                 for(  i=length-1; i>=0; i-- )
 558                                                                 {
 559                                                                         oimpmatch += imp_match_out_scQ( i, i );
 560 //                                                                      fprintf( stderr, "#### i=%d, initial impmatch = %f\n", i, oimpmatch );
 561                                                                 }
 562                                                         }
 563                                                         else
 564                                                         {
 565                                                                 imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 566
 567                                                                 fprintf( stderr, "not supported\n" );
 568                                                                 exit( 1 );
 569
 570                                                                 for(  i=length-1; i>=0; i-- )
 571                                                                 {
 572                                                                         oimpmatch += imp_match_out_sc( i, i );
 573         //                                                              fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
 574                                                                 }
 575                                                         }
 576         //                                              fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
 577                                                 }
 578         //                                      fprintf( stderr, "#### initial impmatch = %f\n", oimpmatch );
 579                                         }
 580                                         else
 581                                         {
 582                                                 oimpmatch = 0.0;
 583                                         }
 584
 585
 586 //                                      fprintf( stderr, "#### tmpdouble = %f\n", tmpdouble );
 587                                         mscore = (double)oimpmatch + tmpdouble;
 588                                 }
 589                                 else
 590                                 {
 591 //                                      fprintf( stderr, "score_check = %d\n", score_check );
 592 /* atode kousokuka */
 593                                         intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
 594                                         mscore = tmpdouble;
 595 /* atode kousokuka */
 596
 597                                         if( constraint )
 598                                         {
 599                                                 oimpmatch = 0.0;
 600                                                 shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
 601                                                 if( use_fft )
 602                                                 {
 603                                                         if( alg == 'Q' )
 604                                                         {
 605                                                                 part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 606                                                                 if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
 607                                                         }
 608                                                         else
 609                                                         {
 610                                                                 part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 611                                                                 if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
 612                                                         }
 613                                                 }
 614                                                 else
 615                                                 {
 616                                                         if( alg == 'Q' )
 617                                                         {
 618                                                                 imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 619                                                                 if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
 620                                                         }
 621                                                         else
 622                                                         {
 623                                                                 imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 624                                                                 fprintf( stderr, "Not supported\n" );
 625                                                                 exit( 1 );
 626                                                         }
 627                                                 }
 628                                         }
 629                                 }
 630
 631 //                              oimpmatch = 0.0;
 632                                 if( constraint )
 633                                 {
 634 #if 0 // iranai
 635                                         if( alg == 'Q' )
 636                                         {
 637                                                 imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 638                                                 for(  i=length-1; i>=0; i-- )
 639                                                 {
 640                                                         oimpmatch += imp_match_out_scQ( i, i );
 641 //                                                      fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
 642                                                 }
 643                                         }
 644                                         else
 645                                         {
 646                                                 imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 647                                                 for(  i=length-1; i>=0; i-- )
 648                                                 {
 649                                                         oimpmatch += imp_match_out_sc( i, i );
 650 //                                                      fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
 651                                                 }
 652                                         }
 653 #endif
 654                                 }
 655 #if 0
 656                                 if( alg == 'H' )
 657                                         naivescore0 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + oimpmatch;
 658                                 else if( alg == 'Q' )
 659                                         naivescore0 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + oimpmatch;
 660                                 else if( alg == 'R' )
 661                                         naivescore0 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + oimpmatch;
 662 #endif
 663
 664 //                              if( rnakozo ) foldalignedrna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, rnapairboth );
 665
 666 //                              if( !use_fft && !rnakozo )
 667                                 if( !use_fft )
 668                                 {
 669                                         commongappick_record( clus1, mseq1, gapmap1 );
 670                                         commongappick_record( clus2, mseq2, gapmap2 );
 671                                 }
 672
 673 #if 0
 674                                 fprintf( stderr, "##### mscore = %f\n", mscore );
 675 #endif
 676
 677 #if DEBUG
 678                                 if( !devide )
 679                                 {
 680                                 fprintf( trap_g, "\nSTEP%d-%d-%d\n", iterate+1, l+1, k );
 681                         fprintf( trap_g, "group1 = %s\n", indication1 );
 682                                 fprintf( trap_g, "group2 = %s\n", indication2 );
 683                                         fflush( trap_g );
 684                                 }
 685
 686 #endif
 687 #if 0
 688                                 printf( "STEP %d-%d-%d\n", iterate, l, k );
 689                                 for( i=0; i<clus2; i++ ) printf( "%f ", effarr2[i] );
 690                                 printf( "\n" );
 691 #endif
 692                                 if( constraint == 2 )
 693                                 {
 694                                         if( use_fft )
 695                                         {
 696 //                                              if( alg == 'Q' )
 697 //                                                      part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
 698 //                                              else
 699 //                                                      part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 700                                                 Falign_localhom( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2 );
 701 //                                              fprintf( stderr, "##### impmatch = %f\n", impmatch );
 702                                         }
 703                                         else
 704                                         {
 705                                                 if( alg == 'Q' )
 706                                                 {
 707                                                         float wm;
 708 //                                                      imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap  wo tsukuttakara iranai.
 709 //                                                      if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, gapmap1, gapmap2, rnapairboth );
 710
 711                                                         wm = Q__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL, gapmap1, gapmap2 );
 712                                                         fprintf( stderr, "wm = %f\n", wm );
 713 #if 0
 714                                                         fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
 715                                                         naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
 716                                                         fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
 717
 718                                                         if( naivescore1 > naivescore0 )
 719                                                                 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
 720                                                         else if( naivescore1 < naivescore0 )
 721                                                                 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
 722                                                         else
 723                                                                 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 724 #if 0 // chuui
 725                                                         if( abs( wm - naivescore1 ) > 100 )
 726                                                         {
 727 //                                                              fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
 728 //                                                              rewind( stdout );
 729 //                                                              for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 730 //                                                              for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 731 //                                                              exit( 1 );
 732                                                         }
 733 #endif
 734 #endif
 735                                                 }
 736                                                 else if( alg == 'R' )
 737                                                 {
 738                                                         float wm;
 739                                                         imp_match_init_strictR( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
 740                                                         wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
 741 //                                                      fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
 742                                                         naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
 743 //                                                      fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
 744
 745                                                         if( naivescore1 > naivescore0 )
 746                                                                 fprintf( stderr, "%d-%d, ns: %f->%f UP!\n", clus1, clus2, naivescore0, naivescore1 );
 747                                                         else if( naivescore1 < naivescore0 )
 748                                                                 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
 749                                                         else
 750                                                                 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 751 #if 0 // chuui
 752                                                         if( abs( wm - naivescore1 ) > 100 )
 753                                                         {
 754 //                                                              fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
 755                                                                 rewind( stdout );
 756                                                                 for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 757                                                                 for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 758                                                                 exit( 1 );
 759                                                         }
 760 #endif
 761                                                 }
 762                                                 else if( alg == 'H' )
 763                                                 {
 764                                                         float wm;
 765                                                         imp_match_init_strictH( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
 766                                                         wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
 767                                                         fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
 768                                                         naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
 769                                                         fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
 770
 771                                                         if( naivescore1 > naivescore0 )
 772                                                                 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
 773                                                         else if( naivescore1 < naivescore0 )
 774                                                                 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
 775                                                         else
 776                                                                 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 777 #if 0 // chuui
 778                                                         if( abs( wm - naivescore1 ) > 100 )
 779                                                         {
 780 //                                                              fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
 781 //                                                              for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 782 //                                                              for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 783 //                                                              exit( 1 );
 784                                                         }
 785 #endif
 786                                                 }
 787                                                 else
 788                                                 {
 789 //                                                      imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
 790                                                         A__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2 );
 791 //                                                      fprintf( stderr, "##### impmatch = %f\n", impmatch );
 792                                                 }
 793                                         }
 794                                 }
 795                                 else if( use_fft )
 796                                 {
 797                                         float totalwm;
 798                                         totalwm = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum );
 799
 800 //                                      fprintf( stderr, "totalwm = %f\n", totalwm );
 801 #if 0
 802                                         if( alg == 'Q' )
 803                                         {
 804                                                 fprintf( stderr, "totalwm = %f\n", totalwm );
 805                                                 naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
 806
 807                                                 if( naivescore1 > naivescore0 )
 808                                                         fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
 809                                                 else if( naivescore1 < naivescore0 )
 810                                                         fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
 811                                                 else
 812                                                         fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 813 #if 1 // chuui
 814                                                 if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
 815                                                 {
 816 //                                                      fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
 817 //                                                      for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 818 //                                                      for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 819 //                                                      exit( 1 );
 820                                                 }
 821 #endif
 822                                         }
 823 #endif
 824                                         if( alg == 'R' )
 825                                         {
 826                                                 fprintf( stderr, "totalwm = %f\n", totalwm );
 827                                                 naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
 828
 829                                                 if( naivescore1 > naivescore0 )
 830                                                         fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
 831                                                 else if( naivescore1 < naivescore0 )
 832                                                         fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
 833                                                 else
 834                                                         fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 835 #if 1 // chuui
 836                                                 if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
 837                                                 {
 838 //                                                      fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
 839 //                                                      for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 840 //                                                      for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 841 //                                                      exit( 1 );
 842                                                 }
 843                                         }
 844 #endif
 845                                 }
 846                                 else
 847                                 {
 848                                         if( alg == 'M' )
 849                                         {
 850                                                 MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL );
 851                                         }
 852                                         else if( alg == 'A' )
 853                                         {
 854                                                 A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
 855                                         }
 856                                         else if( alg == 'Q' )
 857                                         {
 858                                                 float wm;
 859                                                 wm = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
 860                                                 fprintf( stderr, "wm = %f\n", wm );
 861                                                 fprintf( stderr, "impmatch = %f\n", impmatch );
 862                                                 naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
 863
 864                                                 if( naivescore1 > naivescore0 )
 865                                                         fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
 866                                                 else if( naivescore1 < naivescore0 )
 867                                                         fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
 868                                                 else
 869                                                         fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 870 #if 1 // chuui
 871                                                 if( abs( wm - naivescore1 ) > 100 )
 872                                                 {
 873 //                                                      fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
 874 //                                                      rewind( stderr );
 875 //                                                      rewind( stdout );
 876 //                                                      for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 877 //                                                      for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 878 //                                                      exit( 1 );
 879                                                 }
 880 #endif
 881                                         }
 882                                         else if( alg == 'R' )
 883                                         {
 884                                                 float wm;
 885                                                 wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
 886                                                 naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
 887
 888                                                 if( naivescore1 > naivescore0 )
 889                                                         fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
 890                                                 else if( naivescore1 < naivescore0 )
 891                                                         fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
 892                                                 else
 893                                                         fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 894 #if 1 // chuui
 895                                                 if( abs( wm - naivescore1 ) > 100 )
 896                                                 {
 897 //                                                      fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
 898 //                                                      rewind( stderr );
 899 //                                                      rewind( stdout );
 900 //                                                      for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 901 //                                                      for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 902 //                                                      exit( 1 );
 903                                                 }
 904 #endif
 905                                         }
 906                                         else if( alg == 'H' )
 907                                         {
 908                                                 float wm;
 909                                                 wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
 910                                                 naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
 911
 912                                                 if( naivescore1 > naivescore0 )
 913                                                         fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
 914                                                 else if( naivescore1 < naivescore0 )
 915                                                 {
 916                                                         fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
 917                                                 }
 918                                                 else
 919                                                         fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
 920
 921 #if 0 // chuui
 922                                                 if( abs( wm - naivescore1 ) > 100 )
 923                                                 {
 924                                                         rewind( stdout );
 925                                                         for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
 926                                                         for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
 927                                                         exit( 1 );
 928                                                 }
 929 #endif
 930                                         }
 931                                         else if( alg == 'a' )
 932                                         {
 933                                                 Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
 934                                         }
 935                                         else ErrorExit( "Sorry!" );
 936                                 }
 937 //                              fprintf( stderr, "## impmatch = %f\n", impmatch );
 938
 939                                         if( checkC )
 940                                         {
 941                                                 extern double DSPscore();
 942                                                 extern double SSPscore();
 943                                                 static double cur;
 944                                                 static double pre;
 945
 946 /*
 947                                                 pre = SSPscore( locnjob, bseq );
 948                                                 cur = SSPscore( locnjob, aseq );
 949 */
 950                                                 pre = DSPscore( locnjob, bseq );
 951                                                 cur = DSPscore( locnjob, aseq );
 952
 953                                                 fprintf( stderr, "Previous Sscore = %f\n", pre );
 954                                                 fprintf( stderr, "Currnet  Sscore = %f\n\n", cur );
 955                                         }
 956
 957 //                              fprintf( stderr, "## impmatch = %f\n", impmatch );
 958                                 identity = !strcmp( aseq[s1], bseq[s1] );
 959                                 identity *= !strcmp( aseq[s2], bseq[s2] );
 960
 961
 962 /* Bug?  : idnetitcal but score change when scoreing mtx != JTT  */
 963
 964                                 length = strlen( mseq1[0] );
 965
 966                                 if( identity )
 967                                 {
 968                                         tscore = mscore;
 969                                         if( !devide ) fprintf( trap_g, "tscore =  %f   identical.\n", tscore );
 970                                         fprintf( stderr, " identical." );
 971                                         converged++;
 972                                 }
 973                                 else
 974                                 {
 975                                         if( score_check )
 976                                         {
 977                                                 if( constraint == 2 )
 978                                                 {
 979 #if 1
 980                                                         if( RNAscoremtx == 'r' )
 981                                                                 intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
 982                                                         else
 983                                                                 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
 984 #else
 985                                                         intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
 986 #endif
 987
 988                                                         tscore = impmatch + tmpdouble;
 989
 990 //                                                      fprintf( stderr, "tmpdouble=%f, impmatch = %f -> %f, tscore = %f\n", tmpdouble, oimpmatch, impmatch, tscore );
 991                                                 }
 992                                                 else
 993                                                 {
 994                                                         intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
 995                                                         tscore = tmpdouble;
 996                                                 }
 997 //                                              fprintf( stderr, "#######ii=%d, iterate=%d score = %f -> %f \n", ii, iterate , mscore, tscore );
 998         #if 0
 999                                                 for( i=0; i<1; i++ )
1000                                                         fprintf( stderr, "%s\n", mseq1[i] );
1001                                                 fprintf( stderr, "+++++++\n" );
1002                                                 for( i=0; i<1; i++ )
1003                                                         fprintf( stderr, "%s\n", mseq2[i] );
1004         #endif
1005
1006                                         }
1007                                         else
1008                                         {
1009                                                 tscore = mscore + 1.0;
1010 //                                              tscore = 0.0;
1011 //                                              fprintf( stderr, "in line 705, tscore=%f\n", tscore );
1012 //                                              for( i=0; i<length; i++ )
1013 //                                                      tscore = tscore + (double)mseq1[0][i];
1014 //                                              mscore = tscore - 1.0;
1015                                         }
1016
1017                                         if( isnan( mscore ) )
1018                                         {
1019                                                 fprintf( stderr, "\n\nmscore became NaN\n" );
1020                                                 exit( 1 );
1021                                         }
1022                                         if( isnan( tscore ) )
1023                                         {
1024                                                 fprintf( stderr, "\n\ntscore became NaN\n" );
1025                                                 exit( 1 );
1026                                         }
1027
1028
1029
1030 //                                      fprintf( stderr, "@@@@@ mscore,tscore = %f,%f\n", mscore, tscore );
1031
1032                                         if( tscore > mscore - cut/100.0*mscore )
1033                                         {
1034                                                 writePre( locnjob, name, nlen, aseq, 0 );
1035                                                 for( i=0; i<locnjob; i++ ) strcpy( bseq[i], aseq[i] );
1036                                                 if( score_check == 2 )
1037                                                 {
1038                                                         effarr1[0] = 1.0;
1039                                                         effarr2[0] = 1.0;
1040                                                         for( i=0; i<locnjob-1; i++ )
1041                                                                 for( j=i+1; j<locnjob; j++ )
1042                                                                         intergroup_score( bseq+i, bseq+j, effarr1, effarr2, 1, 1, length, imanoten[i]+j );
1043                                                 }
1044
1045 #if 0
1046                                                 fprintf( stderr, "tscore =  %f   mscore = %f  accepted.\n", tscore, mscore );
1047 #endif
1048                                                 fprintf( stderr, " accepted." );
1049                                                 converged = 0;
1050
1051                                         }
1052                                         else
1053                                         {
1054 #if 0
1055                                                 fprintf( stderr, "tscore =  %f   mscore = %f  rejected.\n", tscore, mscore );
1056 #endif
1057                                                 fprintf( stderr, " rejected." );
1058                                                 tscore = mscore;
1059                                                 converged++;
1060                                         }
1061                                 }
1062                                 fprintf( stderr, "\r" );
1063
1064
1065                                 history[iterate][l][k] = (float)tscore;
1066
1067 //                              fprintf( stderr, "tscore = %f\n", tscore );
1068
1069                                 if( converged >= locnjob * 2 )
1070                                 {
1071                                         fprintf( trap_g, "Converged.\n\n" );
1072                                         fprintf( stderr, "\nConverged.\n\n" );
1073                                         return( 0 );
1074                                 }
1075                                 if( iterate >= 1 )
1076                                 {
1077         /*   oscillation check    */
1078                                         oscillating = 0;
1079                                         for( ii=iterate-2; ii>=0; ii-=2 )
1080                                         {
1081                                                 if( (float)tscore == history[ii][l][k] )
1082                                                 {
1083                                                         oscillating = 1;
1084                                                         break;
1085                                                 }
1086                                         }
1087                                         if( ( oscillating && !cooling ) || ( oscillating && cut < 0.001 && cooling ) )
1088                                         {
1089                                                 fprintf( trap_g, "Oscillating.\n" );
1090                                                 fprintf( stderr, "\nOscillating.\n\n" );
1091 #if 1 /* hujuubun */
1092                                                 return( -1 );
1093 #endif
1094                                         }
1095                                 }      /* if( iterate ) */
1096                         }          /* for( k ) */
1097                 }              /* for( l ) */
1098         }                  /* for( iterate ) */
1099         return( 2 );
1100 }                          /* int Tree... */