website/archive/binaries/mac/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 #if 0
 306 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 )
 307 #else
 308 int conjuctionforgaln( int s0, int s1, char **seq, char **aseq, double *peff, double *eff, char **name, char **aname, char *d )
 309 #endif
 310 {
 311         int m, k;
 312         char b[B];
 313         double total;
 314
 315 #if 0
 316         fprintf( stderr, "s0 = %d, s1 = %d\n", s0, s1 );
 317 #endif
 318
 319         total = 0.0;
 320         d[0] = 0;
 321         for( m=s0, k=0; m<s1; m++ )
 322         {
 323                 {
 324                         sprintf( b, " %d", m+1 );
 325 #if 1
 326                         if( strlen( d ) < 100 ) strcat( d, b );
 327 #else
 328                         strcat( d, b );
 329 #endif
 330                         aseq[k] = seq[m];
 331                         peff[k] = eff[m];
 332                         total += peff[k];
 333 #if 0
 334                         strcpy( aname[k], name[m] );
 335 #endif
 336                         k++;
 337                 }
 338         }
 339 #if 1
 340         for( m=0; m<k; m++ )
 341         {
 342                 peff[m] /= total;
 343 //              fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
 344         }
 345 #endif
 346         return( k );
 347 }
 348
 349 void makegrouprna( RNApair ***group, RNApair ***all, int *memlist )
 350 {
 351         int k, m;
 352         for( k=0; (m=*memlist)!=-1; memlist++, k++ )
 353         {
 354                 group[k] = all[m];
 355         }
 356 }
 357
 358 void makegrouprnait( RNApair ***group, RNApair ***all, char **pair, int s )
 359 {
 360         int k, m;
 361         for( m=s, k=0; m<njob; m++ )
 362         {
 363                 if( pair[s][m] != 0 )
 364                 {
 365                         group[k++] = all[m];
 366                 }
 367         }
 368 }
 369
 370 int fastconjuction_noweight( int *memlist, char **seq, char **aseq, double *peff, char *d )
 371 {
 372         int m, k, dln;
 373         char b[B];
 374         double total;
 375
 376 #if DEBUG
 377         fprintf( stderr, "s = %d\n", s );
 378 #endif
 379
 380         total = 0.0;
 381         d[0] = 0;
 382         dln = 0;
 383         for( k=0; *memlist!=-1; memlist++, k++ )
 384         {
 385                 m = *memlist;
 386                 dln += sprintf( b, " %d", m+1 );
 387 #if 1
 388                 if( dln < 100 ) strcat( d, b );
 389 #else
 390                 strcat( d, b );
 391 #endif
 392                 aseq[k] = seq[m];
 393                 peff[k] = 1.0;
 394                 total += peff[k];
 395         }
 396 #if 1
 397         for( m=0; m<k; m++ )
 398                 peff[m] /= total;
 399 #endif
 400         return( k );
 401 }
 402
 403 int fastconjuction_noname_kozo( int *memlist, char **seq, char **aseq, double *peff, double *eff, double *peff_kozo, double *eff_kozo, char *d )
 404 {
 405         int m, k, dln;
 406         char b[B];
 407         double total;
 408         double total_kozo;
 409
 410 #if DEBUG
 411         fprintf( stderr, "s = %d\n", s );
 412 #endif
 413
 414         total = 0.0;
 415         total_kozo = 0.0;
 416         d[0] = 0;
 417         dln = 0;
 418         for( k=0; *memlist!=-1; memlist++, k++ )
 419         {
 420                 m = *memlist;
 421                 dln += sprintf( b, " %d", m+1 );
 422 #if 1
 423                 if( dln < 100 ) strcat( d, b );
 424 #else
 425                 strcat( d, b );
 426 #endif
 427                 aseq[k] = seq[m];
 428                 peff[k] = eff[m];
 429                 peff_kozo[k] = eff_kozo[m];
 430                 total += peff[k];
 431                 total_kozo += peff_kozo[k];
 432         }
 433 #if 1
 434         for( m=0; m<k; m++ )
 435         {
 436 //              fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
 437                 peff[m] /= total;
 438         }
 439         if( total_kozo )
 440         {
 441                 for( m=0; m<k; m++ )
 442                 {
 443 //                      fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
 444                         peff_kozo[m] /= total_kozo;
 445                         if( peff_kozo[m] > 0.0 ) peff_kozo[m] += peff[m];
 446                 }
 447         }
 448         else  //iranai
 449         {
 450                 for( m=0; m<k; m++ )
 451                 {
 452                         peff_kozo[m] = 0.0;
 453                 }
 454         }
 455 #endif
 456
 457 //      fprintf( stderr, "\n\ntbfast_total_kozo = %f\n\n", total_kozo );
 458         return( k );
 459 }
 460
 461
 462 int fastconjuction_noname( int *memlist, char **seq, char **aseq, double *peff, double *eff, char *d )
 463 {
 464         int m, k, dln;
 465         char b[B];
 466         double total;
 467
 468 #if DEBUG
 469         fprintf( stderr, "s = %d\n", s );
 470 #endif
 471
 472         total = 0.0;
 473         d[0] = 0;
 474         dln = 0;
 475         for( k=0; *memlist!=-1; memlist++, k++ )
 476         {
 477                 m = *memlist;
 478                 dln += sprintf( b, " %d", m+1 );
 479 #if 1
 480                 if( dln < 100 ) strcat( d, b );
 481 #else
 482                 strcat( d, b );
 483 #endif
 484                 aseq[k] = seq[m];
 485                 peff[k] = eff[m];
 486                 total += peff[k];
 487         }
 488 #if 1
 489         for( m=0; m<k; m++ )
 490         {
 491 //              fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
 492                 peff[m] /= total;
 493         }
 494 #endif
 495         return( k );
 496 }
 497
 498 int fastconjuction( int *memlist, char **seq, char **aseq, double *peff, double *eff, char name[M][B], char aname[M][B], char *d )
 499 {
 500         int m, k, dln;
 501         char b[B];
 502         double total;
 503
 504 #if DEBUG
 505         fprintf( stderr, "s = %d\n", s );
 506 #endif
 507
 508         total = 0.0;
 509         d[0] = 0;
 510         dln = 0;
 511         for( k=0; *memlist!=-1; memlist++, k++ )
 512         {
 513                 m = *memlist;
 514                 dln += sprintf( b, " %d", m+1 );
 515 #if 1
 516                 if( dln < 100 ) strcat( d, b );
 517 #else
 518                 strcat( d, b );
 519 #endif
 520                 aseq[k] = seq[m];
 521                 peff[k] = eff[m];
 522                 total += peff[k];
 523 #if 0
 524                         strcpy( aname[k], name[m] );
 525 #endif
 526         }
 527 #if 1
 528         for( m=0; m<k; m++ )
 529                 peff[m] /= total;
 530 #endif
 531         return( k );
 532 }
 533
 534
 535 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 )
 536 {
 537         int m, k;
 538         char b[B];
 539         double total;
 540         double total_kozo;
 541
 542 #if DEBUG
 543         fprintf( stderr, "s = %d\n", s );
 544 #endif
 545
 546         total = 0.0;
 547 //      total_kozo = 0.0;
 548         total_kozo = *tmptmptmp; // masaka
 549         d[0] = 0;
 550         for( m=s, k=0; m<njob; m++ )
 551         {
 552                 if( pair[s][m] != 0 )
 553                 {
 554                         sprintf( b, " %d", m+1 );
 555 #if 1
 556                         if( strlen( d ) < 100 ) strcat( d, b );
 557 #else
 558                         strcat( d, b );
 559 #endif
 560                         aseq[k] = seq[m];
 561                         peff[k] = eff[m];
 562                         peff_kozo[k] = eff_kozo[m];
 563                         total += peff[k];
 564                         total_kozo += peff_kozo[k];
 565 #if 0
 566                         strcpy( aname[k], name[m] );
 567 #endif
 568                         k++;
 569                 }
 570         }
 571 #if 1
 572         for( m=0; m<k; m++ )
 573                 peff[m] /= total;
 574         if( total_kozo > 0.0 )
 575         {
 576                 for( m=0; m<k; m++ )
 577                 {
 578                         peff_kozo[m] /= total_kozo;
 579                         if( peff_kozo[m] > 0.0 ) peff_kozo[m] += peff[m];
 580                 }
 581         }
 582         else //iranai
 583         {
 584                 for( m=0; m<k; m++ ) peff_kozo[m] = 0.0;
 585         }
 586 #endif
 587 //      fprintf( stderr, "\n\ndvtditr_total_kozo = %f\n\n", total_kozo );
 588         *tmptmptmp = total_kozo;
 589         return( k );
 590 }
 591
 592 int conjuctionfortbfast( char **pair, int s, char **seq, char **aseq, double *peff, double *eff, char *d )
 593 {
 594         int m, k;
 595         char *b;
 596         double total;
 597
 598         b = calloc( B, sizeof( char ) );
 599 #if DEBUG
 600         fprintf( stderr, "s = %d\n", s );
 601 #endif
 602
 603         total = 0.0;
 604         d[0] = 0;
 605         for( m=s, k=0; m<njob; m++ )
 606         {
 607                 if( pair[s][m] != 0 )
 608                 {
 609                         sprintf( b, " %d", m+1 );
 610 #if 1
 611                         if( strlen( d ) < 100 ) strcat( d, b );
 612 #else
 613                         strcat( d, b );
 614 #endif
 615                         aseq[k] = seq[m];
 616                         peff[k] = eff[m];
 617                         total += peff[k];
 618 #if 0
 619                         strcpy( aname[k], name[m] );
 620 #endif
 621                         k++;
 622                 }
 623         }
 624 #if 1
 625         for( m=0; m<k; m++ )
 626                 peff[m] /= total;
 627 #endif
 628         free( b );
 629         return( k );
 630 }
 631 int conjuction( char **pair, int s, char **seq, char **aseq, double *peff, double *eff, char **name, char **aname, char *d )
 632 {
 633         int m, k;
 634         char b[B];
 635         double total;
 636
 637 #if DEBUG
 638         fprintf( stderr, "s = %d\n", s );
 639 #endif
 640
 641         total = 0.0;
 642         d[0] = 0;
 643         for( m=s, k=0; m<njob; m++ )
 644         {
 645                 if( pair[s][m] != 0 )
 646                 {
 647                         sprintf( b, " %d", m+1 );
 648 #if 1
 649                         if( strlen( d ) < 100 ) strcat( d, b );
 650 #else
 651                         strcat( d, b );
 652 #endif
 653                         aseq[k] = seq[m];
 654                         peff[k] = eff[m];
 655                         total += peff[k];
 656 #if 0
 657                         strcpy( aname[k], name[m] );
 658 #endif
 659                         k++;
 660                 }
 661         }
 662 #if 0
 663         for( m=0; m<k; m++ )
 664                 peff[m] /= total;
 665 #endif
 666         return( k );
 667 }
 668
 669 void floatdelete( float **cpmx, int d, int len )
 670 {
 671         int i, j;
 672
 673         for( i=d; i<len-1; i++ )
 674         {
 675                 for( j=0; j<26; j++ )
 676                 {
 677                         cpmx[j][i] = cpmx[j][i+1];
 678                 }
 679         }
 680 }
 681
 682 void chardelete( char *seq, int d )
 683 {
 684         char b[N];
 685
 686                 strcpy( b, seq+d+1 );
 687                 strcpy( seq+d, b );
 688 }
 689
 690 int RootBranchNode( int nseq, int ***topol, int step, int branch )
 691 {
 692         int i, j, s1, s2, value;
 693
 694         value = 1;
 695         for( i=step+1; i<nseq-2; i++ )
 696         {
 697                 for( j=0; (s1=topol[i][0][j])>-1; j++ )
 698                         if( s1 == topol[step][branch][0] ) value++;
 699                 for( j=0; (s2=topol[i][1][j])>-1; j++ )
 700                         if( s2 == topol[step][branch][0] ) value++;
 701         }
 702         return( value );
 703 }
 704
 705 void BranchLeafNode( int nseq, int ***topol, int *node, int step, int branch )
 706 {
 707         int i, j, k, s;
 708
 709         for( i=0; i<nseq; i++ ) node[i] = 0;
 710         for( i=0; i<step-1; i++ )
 711                 for( k=0; k<2; k++ )
 712                         for( j=0; (s=topol[i][k][j])>-1; j++ )
 713                                 node[s]++;
 714         for( k=0; k<branch+1; k++ )
 715                 for( j=0; (s=topol[step][k][j])>-1; j++ )
 716                         node[s]++;
 717 }
 718
 719 void RootLeafNode( int nseq, int ***topol, int *node )
 720 {
 721         int i, j, k, s;
 722         for( i=0; i<nseq; i++ ) node[i] = 0;
 723         for( i=0; i<nseq-2; i++ )
 724                 for( k=0; k<2; k++ )
 725                         for( j=0; (s=topol[i][k][j])>-1; j++ )
 726                                 node[s]++;
 727 }
 728
 729 void nodeFromABranch( int nseq, int *result, int **pairwisenode, int ***topol, double **len, int step, int num )
 730 {
 731         int i, s, count;
 732         int *innergroup;
 733         int *outergroup1;
 734 #if 0
 735         int outergroup2[nseq];
 736         int table[nseq];
 737 #else
 738         static int *outergroup2 = NULL;
 739         static int *table = NULL;
 740         if( outergroup2 == NULL )
 741         {
 742                 outergroup2 = AllocateIntVec( nseq );
 743                 table = AllocateIntVec( nseq );
 744         }
 745 #endif
 746         innergroup = topol[step][num];
 747         outergroup1 = topol[step][!num];
 748         for( i=0; i<nseq; i++ ) table[i] = 1;
 749         for( i=0; (s=innergroup[i])>-1; i++ ) table[s] = 0;
 750         for( i=0; (s=outergroup1[i])>-1; i++ ) table[s] = 0;
 751         for( i=0, count=0; i<nseq; i++ )
 752         {
 753                 if( table[i] )
 754                 {
 755                         outergroup2[count] = i;
 756                         count++;
 757                 }
 758         }
 759         outergroup2[count] = -1;
 760
 761         for( i=0; (s=innergroup[i])>-1; i++ )
 762         {
 763                 result[s] = pairwisenode[s][outergroup1[0]]
 764                                   + pairwisenode[s][outergroup2[0]]
 765                                   - pairwisenode[outergroup1[0]][outergroup2[0]] - 1;
 766                 result[s] /= 2;
 767         }
 768         for( i=0; (s=outergroup1[i])>-1; i++ )
 769         {
 770                 result[s] = pairwisenode[s][outergroup2[0]]
 771                                   + pairwisenode[s][innergroup[0]]
 772                                   - pairwisenode[innergroup[0]][outergroup2[0]] + 1;
 773                 result[s] /= 2;
 774         }
 775         for( i=0; (s=outergroup2[i])>-1; i++ )
 776         {
 777                 result[s] = pairwisenode[s][outergroup1[0]]
 778                                   + pairwisenode[s][innergroup[0]]
 779                                   - pairwisenode[innergroup[0]][outergroup1[0]] + 1;
 780                 result[s] /= 2;
 781         }
 782
 783 #if 0
 784         for( i=0; i<nseq; i++ )
 785                 fprintf( stderr, "result[%d] = %d\n", i+1, result[i] );
 786 #endif
 787 }
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803 #if 0
 804 void OneClusterAndTheOther( int locnjob, char pair[njob][njob], int *s1, int *s2, int ***topol, int step, int branch )
 805 #else
 806 void OneClusterAndTheOther( int locnjob, char **pair, int *s1, int *s2, int ***topol, int step, int branch )
 807 #endif
 808 {
 809         int i;
 810         int r1;
 811
 812     *s1 = topol[step][branch][0];
 813     for( i=0; (r1=topol[step][branch][i])>-1; i++ )
 814         pair[*s1][r1] = 1;
 815     for( i=0; i<locnjob; i++ )
 816     {
 817         if( !pair[*s1][i] )
 818         {
 819             *s2 = i;
 820             break;
 821         }
 822     }
 823     for( i=*s2; i<locnjob; i++ )
 824     {
 825         if( !pair[*s1][i] )
 826             pair[*s2][i] = 1;
 827     }
 828 }
 829
 830 void makeEffMtx( int nseq, double **mtx, double *vec )
 831 {
 832         int i, j;
 833         for( i=0; i<nseq; i++ ) for( j=0; j<nseq; j++ )
 834                 mtx[i][j] = vec[i] * vec[j];
 835 }
 836
 837 void node_eff( int nseq, double *eff, int *node )
 838 {
 839         int i;
 840         extern double ipower( double, int );
 841         for( i=0; i<nseq; i++ )
 842                 eff[i] = ipower( 0.5, node[i] ) + geta2;
 843         /*
 844                 eff[i] = ipower( 0.5, node[i] ) + 0;
 845         */
 846 #if DEBUG
 847         for( i=0; i<nseq; i++ )
 848 #endif
 849 }
 850
 851
 852 int shrinklocalhom( char **pair, int s1, int s2, LocalHom **localhom, LocalHom ***localhomshrink )
 853 {
 854         int m1, k1, m2, k2;
 855
 856         for( m1=s1, k1=0; m1<njob; m1++ )
 857         {
 858                 if( pair[s1][m1] != 0 )
 859                 {
 860                         for( m2=s2, k2=0; m2<njob; m2++ )
 861                         {
 862                                 if( pair[s2][m2] != 0 )
 863                                 {
 864                                         if( localhom[m1][m2].opt == -1 )
 865                                                 localhomshrink[k1][k2] = NULL;
 866                                         else
 867                                                 localhomshrink[k1][k2] = localhom[m1]+m2;
 868                                         k2++;
 869                                 }
 870                         }
 871                         k1++;
 872                 }
 873         }
 874         return( 0 );
 875 }
 876 int msshrinklocalhom( char **pair, int s1, int s2, LocalHom **localhom, LocalHom ***localhomshrink )
 877 {
 878         int m1, k1, n1, m2, k2, n2;
 879
 880         for( m1=s1, k1=0; m1<njob; m1++ )
 881         {
 882                 if( pair[s1][m1] != 0 )
 883                 {
 884                         for( m2=s2, k2=0; m2<njob; m2++ )
 885                         {
 886                                 if( pair[s2][m2] != 0 )
 887                                 {
 888                                         n1 = MIN(m1,m2); n2=MAX(m1,m2);
 889                                         if( localhom[m1][m2].opt == -1 )
 890                                                 localhomshrink[k1][k2] = NULL;
 891                                         else
 892                                                 localhomshrink[k1][k2] = localhom[n1]+n2;
 893                                         k2++;
 894                                 }
 895                         }
 896                         k1++;
 897                 }
 898         }
 899         return( 0 );
 900 }
 901
 902 int fastshrinklocalhom( int *mem1, int *mem2, LocalHom **localhom, LocalHom ***localhomshrink )
 903 {
 904         int k1, k2;
 905         int *intpt1, *intpt2;
 906
 907
 908         for( intpt1=mem1, k1=0; *intpt1!=-1; intpt1++, k1++ )
 909         {
 910                 for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
 911                 {
 912                         if( localhom[*intpt1][*intpt2].opt == -1 )
 913                                 localhomshrink[k1][k2] = NULL;
 914                         else
 915                                 localhomshrink[k1][k2] = localhom[*intpt1]+*intpt2;
 916                 }
 917         }
 918         return( 0 );
 919 }
 920
 921 int msfastshrinklocalhom( int *mem1, int *mem2, LocalHom **localhom, LocalHom ***localhomshrink )
 922 {
 923         int k1, k2;
 924         int *intpt1, *intpt2;
 925         int m1, m2;
 926
 927         for( intpt1=mem1, k1=0; *intpt1!=-1; intpt1++, k1++ )
 928         {
 929                 for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
 930                 {
 931                         m1 = MIN(*intpt1,*intpt2); m2 = MAX(*intpt1,*intpt2);
 932                         if( localhom[m1][m2].opt == -1 )
 933                                 localhomshrink[k1][k2] = NULL;
 934                         else
 935                                 localhomshrink[k1][k2] = localhom[m1]+m2;
 936                 }
 937         }
 938         return( 0 );
 939 }
 940