binaries/src/mafft/core/tddis.c

   1 #include "mltaln.h"
   2
   3 #define DEBUG 0
   4
   5 #if 0
   6 void mdfymtx( char pair[njob][njob], int s1, double **partialmtx, double **mtx )
   7 #else
   8 void mdfymtx( char **pair, int s1, double **partialmtx, double **mtx )
   9 #endif
  10 {
  11         int i, j;
  12         int icount, jcount;
  13 #if DEBUG
  14         FILE *fp;
  15         static char name[M][B];
  16
  17         for( i=0; i<M; i++ ) name[i][0] = 0;
  18         fprintf( stdout, "s1 = %d\n", s1 );
  19         for( i=0; i<njob; i++ )
  20         {
  21                 for( j=0; j<njob; j++ )
  22                 {
  23                         printf( "%#2d", pair[i][j] );
  24                 }
  25                 printf( "\n" );
  26         }
  27 #endif
  28
  29         for( i=0, icount=0; i<njob-1; i++ )
  30         {
  31                 if( !pair[s1][i] ) continue;
  32                 for( j=i+1, jcount=icount+1; j<njob; j++ )
  33                 {
  34                         if( !pair[s1][j] ) continue;
  35                         partialmtx[icount][jcount] = mtx[i][j];
  36                         jcount++;
  37                 }
  38                 icount++;
  39         }
  40 #if DEBUG
  41         fp = fopen( "hat2.org", "w" );
  42         WriteHat2( fp, njob, name, mtx );
  43         fclose( fp );
  44         fp = fopen( "hat2.mdf", "w" );
  45         WriteHat2( fp, icount, name, partialmtx );
  46         fclose( fp );
  47 #endif
  48
  49 }
  50
  51
  52 float score_calc( char **seq, int s )  /* method 3  */
  53 {
  54     int i, j, k, c;
  55     int len = strlen( seq[0] );
  56     float score;
  57     int tmpscore;
  58     char *mseq1, *mseq2;
  59
  60     score = 0.0;
  61     for( i=0; i<s-1; i++ )
  62     {
  63         for( j=i+1; j<s; j++ )
  64         {
  65             mseq1 = seq[i];
  66             mseq2 = seq[j];
  67             tmpscore = 0;
  68             c = 0;
  69             for( k=0; k<len; k++ )
  70             {
  71                 if( mseq1[k] == '-' && mseq2[k] == '-' ) continue;
  72                 c++;
  73                 tmpscore += amino_dis[(int)mseq1[k]][(int)mseq2[k]];
  74                 if( mseq1[k] == '-' )
  75                 {
  76                     tmpscore += penalty;
  77                     while( mseq1[++k] == '-' )
  78                         ;
  79                     k--;
  80                     if( k > len-2 ) break;
  81                     continue;
  82                 }
  83                 if( mseq2[k] == '-' )
  84                 {
  85                     tmpscore += penalty;
  86                     while( mseq2[++k] == '-' )
  87                         ;
  88                     k--;
  89                     if( k > len-2 ) break;
  90                     continue;
  91                 }
  92             }
  93             /*
  94             if( mseq1[0] == '-' || mseq2[0] == '-' )
  95             {
  96                 for( k=0; k<len; k++ )
  97                 {
  98                     if( mseq1[k] == '-' && mseq2[k] == '-' ) continue;
  99                     if( !( mseq1[k] != '-' && mseq2[k] != '-' ) )
 100                     {
 101                         c--;
 102                         tmpscore -= penalty;
 103                         break;
 104                     }
 105                     else break;
 106                 }
 107             }
 108             if( mseq1[len-1] == '-' || mseq2[len-1] == '-' )
 109             {
 110                 for( k=0; k<len; k++ )
 111                 {
 112                     if( mseq1[k] == '-' && mseq2[k] == '-' ) continue;
 113                     if( !( mseq1[k] != '-' && mseq2[k] != '-' ) )
 114                     {
 115                         c--;
 116                         tmpscore -= penalty;
 117                         break;
 118                     }
 119                     else break;
 120                 }
 121             }
 122             */
 123             score += (double)tmpscore / (double)c;
 124         }
 125     }
 126     score = (float)score / ( ( (double)s * ((double)s-1.0) ) / 2.0 );
 127         fprintf( stderr, "score in score_calc = %f\n", score );
 128     return( score );
 129 }
 130
 131 void cpmx_calc( char **seq, float **cpmx, double *eff, int lgth, int clus )
 132 {
 133         int  i, j, k;
 134         double totaleff = 0.0;
 135
 136         for( i=0; i<clus; i++ ) totaleff += eff[i];
 137         for( i=0; i<26; i++ ) for( j=0; j<lgth; j++ ) cpmx[i][j] = 0.0;
 138         for( j=0; j<lgth; j++ ) for( k=0; k<clus; k++ )
 139                         cpmx[(int)amino_n[(int)seq[k][j]]][j] += (float)eff[k] / totaleff;
 140 }
 141
 142
 143 void cpmx_calc_new_bk( char **seq, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
 144 {
 145     int  i, j, k;
 146     float feff;
 147
 148     for( i=0; i<26; i++ ) for( j=0; j<lgth; j++ ) cpmx[i][j] = 0.0;
 149     for( k=0; k<clus; k++ )
 150     {
 151         feff = (float)eff[k];
 152         for( j=0; j<lgth; j++ )
 153         {
 154             cpmx[(int)amino_n[(int)seq[k][j]]][j] += feff;
 155         }
 156     }
 157 }
 158
 159 void cpmx_calc_new( char **seq, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
 160 {
 161         int  i, j, k;
 162         float feff;
 163         float *cpmxpt, **cpmxptpt;
 164         char *seqpt;
 165
 166         j = 26;
 167         cpmxptpt = cpmx;
 168         while( j-- )
 169         {
 170                 cpmxpt = *cpmxptpt++;
 171                 i = lgth;
 172                 while( i-- )
 173                         *cpmxpt++ = 0.0;
 174         }
 175         for( k=0; k<clus; k++ )
 176         {
 177                 feff = (float)eff[k];
 178                 seqpt = seq[k];
 179 //              fprintf( stderr, "seqpt = %s, lgth=%d\n", seqpt, lgth );
 180                 for( j=0; j<lgth; j++ )
 181                 {
 182                         cpmx[(int)amino_n[(int)*seqpt++]][j] += feff;
 183                 }
 184         }
 185 }
 186 void MScpmx_calc_new( char **seq, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
 187 {
 188         int  i, j, k;
 189         float feff;
 190         float **cpmxptpt, *cpmxpt;
 191         char *seqpt;
 192
 193         j = lgth;
 194         cpmxptpt = cpmx;
 195         while( j-- )
 196         {
 197                 cpmxpt = *cpmxptpt++;
 198                 i = 26;
 199                 while( i-- )
 200                         *cpmxpt++ = 0.0;
 201         }
 202         for( k=0; k<clus; k++ )
 203         {
 204                 feff = (float)eff[k];
 205                 seqpt = seq[k];
 206                 cpmxptpt = cpmx;
 207                 j = lgth;
 208                 while( j-- )
 209                         (*cpmxptpt++)[(int)amino_n[(int)*seqpt++]] += feff;
 210         }
 211 #if 0
 212         for( j=0; j<lgth; j++ ) for( i=0; i<26; i++ ) cpmx[j][i] = 0.0;
 213         for( k=0; k<clus; k++ )
 214         {
 215                 feff = (float)eff[k];
 216                 for( j=0; j<lgth; j++ )
 217                         cpmx[j][(int)amino_n[(int)seq[k][j]]] += feff;
 218         }
 219 #endif
 220 }
 221
 222 void cpmx_ribosum( char **seq, char **seqr, char *dir, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
 223 {
 224         int  i, j, k;
 225         float feff;
 226         float **cpmxptpt, *cpmxpt;
 227         char *seqpt, *seqrpt, *dirpt;
 228         int code, code1, code2;
 229
 230         j = lgth;
 231         cpmxptpt = cpmx;
 232         while( j-- )
 233         {
 234                 cpmxpt = *cpmxptpt++;
 235                 i = 37;
 236                 while( i-- )
 237                         *cpmxpt++ = 0.0;
 238         }
 239         for( k=0; k<clus; k++ )
 240         {
 241                 feff = (float)eff[k];
 242                 seqpt = seq[k];
 243                 seqrpt = seqr[k];
 244                 dirpt = dir;
 245                 cpmxptpt = cpmx;
 246                 j = lgth;
 247                 while( j-- )
 248                 {
 249 #if 0
 250                         code1 = amino_n[(int)*seqpt];
 251                         if( code1 > 3 ) code = 36;
 252                         else
 253                                 code = code1;
 254 #else
 255                         code1 = amino_n[(int)*seqpt];
 256                         code2 = amino_n[(int)*seqrpt];
 257                         if( code1 > 3 )
 258                         {
 259                                 code = 36;
 260                         }
 261                         else if( code2 > 3 )
 262                         {
 263                                 code = code1;
 264                         }
 265                         else if( *dirpt == '5' )
 266                         {
 267                                 code = 4 + code2 * 4  + code1;
 268                         }
 269                         else if( *dirpt == '3' )
 270                         {
 271                                 code = 20 + code2 * 4  + code1;
 272                         }
 273                         else // if( *dirpt == 'o' ) // nai
 274                         {
 275                                 code = code1;
 276                         }
 277 #endif
 278
 279 //                      fprintf( stderr, "%c -> code=%d toa=%d, tog=%d, toc=%d, tot=%d, ton=%d, efee=%f\n", *seqpt, code%4, ribosumdis[code][4+0], ribosumdis[code][4+1], ribosumdis[code][4+2], ribosumdis[code][20+3], ribosumdis[code][36], feff );
 280
 281                         seqpt++;
 282                         seqrpt++;
 283                         dirpt++;
 284
 285                         (*cpmxptpt++)[code] += feff;
 286                 }
 287         }
 288 }
 289
 290 void mseqcat( char **seq1, char **seq2, double **eff, double *effarr1, double *effarr2, char name1[M][B], char name2[M][B], int clus1, int clus2 )
 291 {
 292         int i, j;
 293     for( i=0; i<clus2; i++ )
 294         seq1[clus1+i] = seq2[i];
 295     for( i=0; i<clus2; i++ ) strcpy( name1[clus1+i], name2[i] );
 296
 297         for( i=0; i<clus1; i++ ) for( j=0; j<clus1; j++ ) eff[i][j] = effarr1[i]* effarr1[j];
 298         for( i=0; i<clus1; i++ ) for( j=clus1; j<clus1+clus2; j++ ) eff[i][j] = effarr1[i]* effarr2[j-clus1];
 299         for( i=clus1; i<clus1+clus2; i++ ) for( j=0; j<clus1; j++ ) eff[i][j] = effarr2[i-clus1] * effarr1[j];
 300         for( i=clus1; i<clus1+clus2; i++ ) for( j=clus1; j<clus1+clus2; j++ ) eff[i][j] = effarr2[i-clus1] * effarr2[j-clus1];
 301 }
 302
 303
 304
 305 void strnbcat( char *s1, char *s2, int m )   /* s1 + s2 ---> s2  */
 306 {
 307         static char b[N];
 308         strncpy( b, s1, m ); b[m] = 0;
 309         strcat( b, s2 );
 310         strcpy( s2, b );
 311 }
 312
 313 #if 0
 314 int conjuction( char pair[njob][njob], int s, char **seq, char **aseq, double *peff, double *eff, char name[M][B], char aname[M][B], char *d )
 315 #else
 316 int conjuctionforgaln( int s0, int s1, char **seq, char **aseq, double *peff, double *eff, char **name, char **aname, char *d )
 317 #endif
 318 {
 319         int m, k;
 320         char b[B];
 321         double total;
 322
 323 #if 0
 324         fprintf( stderr, "s0 = %d, s1 = %d\n", s0, s1 );
 325 #endif
 326
 327         total = 0.0;
 328         d[0] = 0;
 329         for( m=s0, k=0; m<s1; m++ )
 330         {
 331                 {
 332                         sprintf( b, " %d", m+1 );
 333 #if 1
 334                         if( strlen( d ) < 100 ) strcat( d, b );
 335 #else
 336                         strcat( d, b );
 337 #endif
 338                         aseq[k] = seq[m];
 339                         peff[k] = eff[m];
 340                         total += peff[k];
 341 #if 0
 342                         strcpy( aname[k], name[m] );
 343 #endif
 344                         k++;
 345                 }
 346         }
 347 #if 1
 348         for( m=0; m<k; m++ )
 349         {
 350                 peff[m] /= total;
 351 //              fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
 352         }
 353 #endif
 354         return( k );
 355 }
 356
 357 void makegrouprna( RNApair ***group, RNApair ***all, int *memlist )
 358 {
 359         int k, m;
 360         for( k=0; (m=*memlist)!=-1; memlist++, k++ )
 361         {
 362                 group[k] = all[m];
 363         }
 364 }
 365
 366 void makegrouprnait( RNApair ***group, RNApair ***all, char **pair, int s )
 367 {
 368         int k, m;
 369         for( m=s, k=0; m<njob; m++ )
 370         {
 371                 if( pair[s][m] != 0 )
 372                 {
 373                         group[k++] = all[m];
 374                 }
 375         }
 376 }
 377
 378 int fastconjuction_noweight( int *memlist, char **seq, char **aseq, double *peff, char *d )
 379 {
 380         int m, k, dln;
 381         char b[B];
 382         double total;
 383
 384 #if DEBUG
 385         fprintf( stderr, "s = %d\n", s );
 386 #endif
 387
 388         total = 0.0;
 389         d[0] = 0;
 390         dln = 0;
 391         for( k=0; *memlist!=-1; memlist++, k++ )
 392         {
 393                 m = *memlist;
 394                 dln += sprintf( b, " %d", m+1 );
 395 #if 1
 396                 if( dln < 100 ) strcat( d, b );
 397 #else
 398                 strcat( d, b );
 399 #endif
 400                 aseq[k] = seq[m];
 401                 peff[k] = 1.0;
 402                 total += peff[k];
 403         }
 404 #if 1
 405         for( m=0; m<k; m++ )
 406                 peff[m] /= total;
 407 #endif
 408         return( k );
 409 }
 410
 411 int fastconjuction_noname_kozo( int *memlist, char **seq, char **aseq, double *peff, double *eff, double *peff_kozo, double *eff_kozo, char *d )
 412 {
 413         int m, k, dln;
 414         char b[B];
 415         double total;
 416         double total_kozo;
 417
 418 #if DEBUG
 419         fprintf( stderr, "s = %d\n", s );
 420 #endif
 421
 422         total = 0.0;
 423         total_kozo = 0.0;
 424         d[0] = 0;
 425         dln = 0;
 426         for( k=0; *memlist!=-1; memlist++, k++ )
 427         {
 428                 m = *memlist;
 429                 dln += sprintf( b, " %d", m+1 );
 430 #if 1
 431                 if( dln < 100 ) strcat( d, b );
 432 #else
 433                 strcat( d, b );
 434 #endif
 435                 aseq[k] = seq[m];
 436                 peff[k] = eff[m];
 437                 peff_kozo[k] = eff_kozo[m];
 438                 total += peff[k];
 439                 total_kozo += peff_kozo[k];
 440         }
 441 #if 1
 442         for( m=0; m<k; m++ )
 443         {
 444 //              fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
 445                 peff[m] /= total;
 446         }
 447         if( total_kozo )
 448         {
 449                 for( m=0; m<k; m++ )
 450                 {
 451 //                      fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
 452                         peff_kozo[m] /= total_kozo;
 453                         if( peff_kozo[m] > 0.0 ) peff_kozo[m] += peff[m];
 454                 }
 455         }
 456         else  //iranai
 457         {
 458                 for( m=0; m<k; m++ )
 459                 {
 460                         peff_kozo[m] = 0.0;
 461                 }
 462         }
 463 #endif
 464
 465 //      fprintf( stderr, "\n\ntbfast_total_kozo = %f\n\n", total_kozo );
 466         return( k );
 467 }
 468
 469
 470 int fastconjuction_noname( int *memlist, char **seq, char **aseq, double *peff, double *eff, char *d )
 471 {
 472         int m, k, dln;
 473         char b[B];
 474         double total;
 475
 476 #if DEBUG
 477         fprintf( stderr, "s = %d\n", s );
 478 #endif
 479
 480         total = 0.0;
 481         d[0] = 0;
 482         dln = 0;
 483         for( k=0; *memlist!=-1; memlist++, k++ )
 484         {
 485                 m = *memlist;
 486                 dln += sprintf( b, " %d", m+1 );
 487 #if 1
 488                 if( dln < 100 ) strcat( d, b );
 489 #else
 490                 strcat( d, b );
 491 #endif
 492                 aseq[k] = seq[m];
 493                 peff[k] = eff[m];
 494                 total += peff[k];
 495         }
 496 #if 1
 497         for( m=0; m<k; m++ )
 498         {
 499 //              fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
 500                 peff[m] /= total;
 501         }
 502 #endif
 503         return( k );
 504 }
 505
 506 int fastconjuction( int *memlist, char **seq, char **aseq, double *peff, double *eff, char name[M][B], char aname[M][B], char *d )
 507 {
 508         int m, k, dln;
 509         char b[B];
 510         double total;
 511
 512 #if DEBUG
 513         fprintf( stderr, "s = %d\n", s );
 514 #endif
 515
 516         total = 0.0;
 517         d[0] = 0;
 518         dln = 0;
 519         for( k=0; *memlist!=-1; memlist++, k++ )
 520         {
 521                 m = *memlist;
 522                 dln += sprintf( b, " %d", m+1 );
 523 #if 1
 524                 if( dln < 100 ) strcat( d, b );
 525 #else
 526                 strcat( d, b );
 527 #endif
 528                 aseq[k] = seq[m];
 529                 peff[k] = eff[m];
 530                 total += peff[k];
 531 #if 0
 532                         strcpy( aname[k], name[m] );
 533 #endif
 534         }
 535 #if 1
 536         for( m=0; m<k; m++ )
 537                 peff[m] /= total;
 538 #endif
 539         return( k );
 540 }
 541
 542
 543 int conjuctionfortbfast_kozo( double *tmptmptmp, char **pair, int s, char **seq, char **aseq, double *peff, double *eff, double *peff_kozo, double *eff_kozo, char *d )
 544 {
 545         int m, k;
 546         char b[B];
 547         double total;
 548         double total_kozo;
 549
 550 #if DEBUG
 551         fprintf( stderr, "s = %d\n", s );
 552 #endif
 553
 554         total = 0.0;
 555 //      total_kozo = 0.0;
 556         total_kozo = *tmptmptmp; // masaka
 557         d[0] = 0;
 558         for( m=s, k=0; m<njob; m++ )
 559         {
 560                 if( pair[s][m] != 0 )
 561                 {
 562                         sprintf( b, " %d", m+1 );
 563 #if 1
 564                         if( strlen( d ) < 100 ) strcat( d, b );
 565 #else
 566                         strcat( d, b );
 567 #endif
 568                         aseq[k] = seq[m];
 569                         peff[k] = eff[m];
 570                         peff_kozo[k] = eff_kozo[m];
 571                         total += peff[k];
 572                         total_kozo += peff_kozo[k];
 573 #if 0
 574                         strcpy( aname[k], name[m] );
 575 #endif
 576                         k++;
 577                 }
 578         }
 579 #if 1
 580         for( m=0; m<k; m++ )
 581                 peff[m] /= total;
 582         if( total_kozo > 0.0 )
 583         {
 584                 for( m=0; m<k; m++ )
 585                 {
 586                         peff_kozo[m] /= total_kozo;
 587                         if( peff_kozo[m] > 0.0 ) peff_kozo[m] += peff[m];
 588                 }
 589         }
 590         else //iranai
 591         {
 592                 for( m=0; m<k; m++ ) peff_kozo[m] = 0.0;
 593         }
 594 #endif
 595 //      fprintf( stderr, "\n\ndvtditr_total_kozo = %f\n\n", total_kozo );
 596         *tmptmptmp = total_kozo;
 597         return( k );
 598 }
 599
 600 int conjuctionfortbfast( char **pair, int s, char **seq, char **aseq, double *peff, double *eff, char *d )
 601 {
 602         int m, k;
 603         char b[B];
 604         double total;
 605
 606 #if DEBUG
 607         fprintf( stderr, "s = %d\n", s );
 608 #endif
 609
 610         total = 0.0;
 611         d[0] = 0;
 612         for( m=s, k=0; m<njob; m++ )
 613         {
 614                 if( pair[s][m] != 0 )
 615                 {
 616                         sprintf( b, " %d", m+1 );
 617 #if 1
 618                         if( strlen( d ) < 100 ) strcat( d, b );
 619 #else
 620                         strcat( d, b );
 621 #endif
 622                         aseq[k] = seq[m];
 623                         peff[k] = eff[m];
 624                         total += peff[k];
 625 #if 0
 626                         strcpy( aname[k], name[m] );
 627 #endif
 628                         k++;
 629                 }
 630         }
 631 #if 1
 632         for( m=0; m<k; m++ )
 633                 peff[m] /= total;
 634 #endif
 635         return( k );
 636 }
 637 int conjuction( char **pair, int s, char **seq, char **aseq, double *peff, double *eff, char name[M][B], char aname[M][B], char *d )
 638 {
 639         int m, k;
 640         char b[B];
 641         double total;
 642
 643 #if DEBUG
 644         fprintf( stderr, "s = %d\n", s );
 645 #endif
 646
 647         total = 0.0;
 648         d[0] = 0;
 649         for( m=s, k=0; m<njob; m++ )
 650         {
 651                 if( pair[s][m] != 0 )
 652                 {
 653                         sprintf( b, " %d", m+1 );
 654 #if 1
 655                         if( strlen( d ) < 100 ) strcat( d, b );
 656 #else
 657                         strcat( d, b );
 658 #endif
 659                         aseq[k] = seq[m];
 660                         peff[k] = eff[m];
 661                         total += peff[k];
 662 #if 0
 663                         strcpy( aname[k], name[m] );
 664 #endif
 665                         k++;
 666                 }
 667         }
 668 #if 0
 669         for( m=0; m<k; m++ )
 670                 peff[m] /= total;
 671 #endif
 672         return( k );
 673 }
 674
 675 void floatdelete( float **cpmx, int d, int len )
 676 {
 677         int i, j;
 678
 679         for( i=d; i<len-1; i++ )
 680         {
 681                 for( j=0; j<26; j++ )
 682                 {
 683                         cpmx[j][i] = cpmx[j][i+1];
 684                 }
 685         }
 686 }
 687
 688 void chardelete( char *seq, int d )
 689 {
 690         char b[N];
 691
 692                 strcpy( b, seq+d+1 );
 693                 strcpy( seq+d, b );
 694 }
 695
 696 int RootBranchNode( int nseq, int ***topol, int step, int branch )
 697 {
 698         int i, j, s1, s2, value;
 699
 700         value = 1;
 701         for( i=step+1; i<nseq-2; i++ )
 702         {
 703                 for( j=0; (s1=topol[i][0][j])>-1; j++ )
 704                         if( s1 == topol[step][branch][0] ) value++;
 705                 for( j=0; (s2=topol[i][1][j])>-1; j++ )
 706                         if( s2 == topol[step][branch][0] ) value++;
 707         }
 708         return( value );
 709 }
 710
 711 void BranchLeafNode( int nseq, int ***topol, int *node, int step, int branch )
 712 {
 713         int i, j, k, s;
 714
 715         for( i=0; i<nseq; i++ ) node[i] = 0;
 716         for( i=0; i<step-1; i++ )
 717                 for( k=0; k<2; k++ )
 718                         for( j=0; (s=topol[i][k][j])>-1; j++ )
 719                                 node[s]++;
 720         for( k=0; k<branch+1; k++ )
 721                 for( j=0; (s=topol[step][k][j])>-1; j++ )
 722                         node[s]++;
 723 }
 724
 725 void RootLeafNode( int nseq, int ***topol, int *node )
 726 {
 727         int i, j, k, s;
 728         for( i=0; i<nseq; i++ ) node[i] = 0;
 729         for( i=0; i<nseq-2; i++ )
 730                 for( k=0; k<2; k++ )
 731                         for( j=0; (s=topol[i][k][j])>-1; j++ )
 732                                 node[s]++;
 733 }
 734
 735 void nodeFromABranch( int nseq, int *result, int **pairwisenode, int ***topol, double **len, int step, int num )
 736 {
 737         int i, s, count;
 738         int *innergroup;
 739         int *outergroup1;
 740 #if 0
 741         int outergroup2[nseq];
 742         int table[nseq];
 743 #else
 744         static int *outergroup2 = NULL;
 745         static int *table = NULL;
 746         if( outergroup2 == NULL )
 747         {
 748                 outergroup2 = AllocateIntVec( nseq );
 749                 table = AllocateIntVec( nseq );
 750         }
 751 #endif
 752         innergroup = topol[step][num];
 753         outergroup1 = topol[step][!num];
 754         for( i=0; i<nseq; i++ ) table[i] = 1;
 755         for( i=0; (s=innergroup[i])>-1; i++ ) table[s] = 0;
 756         for( i=0; (s=outergroup1[i])>-1; i++ ) table[s] = 0;
 757         for( i=0, count=0; i<nseq; i++ )
 758         {
 759                 if( table[i] )
 760                 {
 761                         outergroup2[count] = i;
 762                         count++;
 763                 }
 764         }
 765         outergroup2[count] = -1;
 766
 767         for( i=0; (s=innergroup[i])>-1; i++ )
 768         {
 769                 result[s] = pairwisenode[s][outergroup1[0]]
 770                                   + pairwisenode[s][outergroup2[0]]
 771                                   - pairwisenode[outergroup1[0]][outergroup2[0]] - 1;
 772                 result[s] /= 2;
 773         }
 774         for( i=0; (s=outergroup1[i])>-1; i++ )
 775         {
 776                 result[s] = pairwisenode[s][outergroup2[0]]
 777                                   + pairwisenode[s][innergroup[0]]
 778                                   - pairwisenode[innergroup[0]][outergroup2[0]] + 1;
 779                 result[s] /= 2;
 780         }
 781         for( i=0; (s=outergroup2[i])>-1; i++ )
 782         {
 783                 result[s] = pairwisenode[s][outergroup1[0]]
 784                                   + pairwisenode[s][innergroup[0]]
 785                                   - pairwisenode[innergroup[0]][outergroup1[0]] + 1;
 786                 result[s] /= 2;
 787         }
 788
 789 #if 0
 790         for( i=0; i<nseq; i++ )
 791                 fprintf( stderr, "result[%d] = %d\n", i+1, result[i] );
 792 #endif
 793 }
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809 #if 0
 810 void OneClusterAndTheOther( int locnjob, char pair[njob][njob], int *s1, int *s2, int ***topol, int step, int branch )
 811 #else
 812 void OneClusterAndTheOther( int locnjob, char **pair, int *s1, int *s2, int ***topol, int step, int branch )
 813 #endif
 814 {
 815         int i;
 816         int r1;
 817
 818     *s1 = topol[step][branch][0];
 819     for( i=0; (r1=topol[step][branch][i])>-1; i++ )
 820         pair[*s1][r1] = 1;
 821     for( i=0; i<locnjob; i++ )
 822     {
 823         if( !pair[*s1][i] )
 824         {
 825             *s2 = i;
 826             break;
 827         }
 828     }
 829     for( i=*s2; i<locnjob; i++ )
 830     {
 831         if( !pair[*s1][i] )
 832             pair[*s2][i] = 1;
 833     }
 834 }
 835
 836 void makeEffMtx( int nseq, double **mtx, double *vec )
 837 {
 838         int i, j;
 839         for( i=0; i<nseq; i++ ) for( j=0; j<nseq; j++ )
 840                 mtx[i][j] = vec[i] * vec[j];
 841 }
 842
 843 void node_eff( int nseq, double *eff, int *node )
 844 {
 845         int i;
 846         extern double ipower( double, int );
 847         for( i=0; i<nseq; i++ )
 848                 eff[i] = ipower( 0.5, node[i] ) + geta2;
 849         /*
 850                 eff[i] = ipower( 0.5, node[i] ) + 0;
 851         */
 852 #if DEBUG
 853         for( i=0; i<nseq; i++ )
 854 #endif
 855 }
 856
 857
 858 int shrinklocalhom( char **pair, int s1, int s2, LocalHom **localhom, LocalHom ***localhomshrink )
 859 {
 860         int m1, k1, m2, k2;
 861
 862         for( m1=s1, k1=0; m1<njob; m1++ )
 863         {
 864                 if( pair[s1][m1] != 0 )
 865                 {
 866                         for( m2=s2, k2=0; m2<njob; m2++ )
 867                         {
 868                                 if( pair[s2][m2] != 0 )
 869                                 {
 870                                         if( localhom[m1][m2].opt == -1 )
 871                                                 localhomshrink[k1][k2] = NULL;
 872                                         else
 873                                                 localhomshrink[k1][k2] = localhom[m1]+m2;
 874                                         k2++;
 875                                 }
 876                         }
 877                         k1++;
 878                 }
 879         }
 880         return( 0 );
 881 }
 882 int msshrinklocalhom( char **pair, int s1, int s2, LocalHom **localhom, LocalHom ***localhomshrink )
 883 {
 884         int m1, k1, n1, m2, k2, n2;
 885
 886         for( m1=s1, k1=0; m1<njob; m1++ )
 887         {
 888                 if( pair[s1][m1] != 0 )
 889                 {
 890                         for( m2=s2, k2=0; m2<njob; m2++ )
 891                         {
 892                                 if( pair[s2][m2] != 0 )
 893                                 {
 894                                         n1 = MIN(m1,m2); n2=MAX(m1,m2);
 895                                         if( localhom[m1][m2].opt == -1 )
 896                                                 localhomshrink[k1][k2] = NULL;
 897                                         else
 898                                                 localhomshrink[k1][k2] = localhom[n1]+n2;
 899                                         k2++;
 900                                 }
 901                         }
 902                         k1++;
 903                 }
 904         }
 905         return( 0 );
 906 }
 907
 908 int fastshrinklocalhom( int *mem1, int *mem2, LocalHom **localhom, LocalHom ***localhomshrink )
 909 {
 910         int k1, k2;
 911         int *intpt1, *intpt2;
 912
 913
 914         for( intpt1=mem1, k1=0; *intpt1!=-1; intpt1++, k1++ )
 915         {
 916                 for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
 917                 {
 918                         if( localhom[*intpt1][*intpt2].opt == -1 )
 919                                 localhomshrink[k1][k2] = NULL;
 920                         else
 921                                 localhomshrink[k1][k2] = localhom[*intpt1]+*intpt2;
 922                 }
 923         }
 924         return( 0 );
 925 }
 926
 927 int msfastshrinklocalhom( int *mem1, int *mem2, LocalHom **localhom, LocalHom ***localhomshrink )
 928 {
 929         int k1, k2;
 930         int *intpt1, *intpt2;
 931         int m1, m2;
 932
 933         for( intpt1=mem1, k1=0; *intpt1!=-1; intpt1++, k1++ )
 934         {
 935                 for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
 936                 {
 937                         m1 = MIN(*intpt1,*intpt2); m2 = MAX(*intpt1,*intpt2);
 938                         if( localhom[m1][m2].opt == -1 )
 939                                 localhomshrink[k1][k2] = NULL;
 940                         else
 941                                 localhomshrink[k1][k2] = localhom[m1]+m2;
 942                 }
 943         }
 944         return( 0 );
 945 }
 946