binaries/src/clustalw/src/tree/NJTree.cpp

   1 /**
   2  * Author: Mark Larkin
   3  *
   4  * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
   5  */
   6 #ifdef HAVE_CONFIG_H
   7     #include "config.h"
   8 #endif
   9 #include <math.h>
  10 #include "NJTree.h"
  11
  12 namespace clustalw
  13 {
  14
  15     /****************************************************************************
  16      * [ Improvement ideas in fast_nj_tree() ] by DDBJ & FUJITSU Limited.
  17      *                        written by Tadashi Koike
  18      *                        (takoike@genes.nig.ac.jp)
  19      *******************
  20      * <IMPROVEMENT 1> : Store the value of sum of the score to temporary array,
  21      *                   and use again and again.
  22      *
  23      *    In the main cycle, these are calculated again and again :
  24      *        diq = sum of distMat[n][ii]   (n:1 to lastSeq-firstSeq+1),
  25      *        djq = sum of distMat[n][jj]   (n:1 to lastSeq-firstSeq+1),
  26      *        dio = sum of distMat[n][mini] (n:1 to lastSeq-firstSeq+1),
  27      *        djq = sum of distMat[n][minj] (n:1 to lastSeq-firstSeq+1)
  28      *        // 'lastSeq' and 'firstSeq' are both constant values //
  29      *    and the result of above calculations is always same until
  30      *    a best pair of neighbour nodes is joined.
  31      *
  32      *    So, we change the logic to calculate the sum[i] (=sum of distMat[n][i]
  33      *    (n:1 to lastSeq-firstSeq+1)) and store it to array, before
  34      *    beginning to find a best pair of neighbour nodes, and after that
  35      *    we use them again and again.
  36      *
  37      *        tmat[i][j]
  38      *                  1   2   3   4   5
  39      *                +---+---+---+---+---+
  40      *              1 |   |   |   |   |   |
  41      *                +---+---+---+---+---+
  42      *              2 |   |   |   |   |   |  1) calculate sum of distMat[n][i]
  43      *                +---+---+---+---+---+        (n: 1 to lastSeq-firstSeq+1)
  44      *              3 |   |   |   |   |   |  2) store that sum value to sum[i]
  45      *                +---+---+---+---+---+
  46      *              4 |   |   |   |   |   |  3) use sum[i] during finding a best
  47      *                +---+---+---+---+---+     pair of neibour nodes.
  48      *              5 |   |   |   |   |   |
  49      *                +---+---+---+---+---+
  50      *                  |   |   |   |   |
  51      *                  V   V   V   V   V  Calculate sum , and store it to sum[i]
  52      *                +---+---+---+---+---+
  53      *         sum[i] |   |   |   |   |   |
  54      *                +---+---+---+---+---+
  55      *
  56      *    At this time, we thought that we use upper triangle of the matrix
  57      *    because distMat[i][j] is equal to distMat[j][i] and distMat[i][i] is equal
  58      *    to zero. Therefore, we prepared sum_rows[i] and sum_cols[i] instead
  59      *    of sum[i] for storing the sum value.
  60      *
  61      *        distMat[i][j]
  62      *                  1   2   3   4   5     sum_cols[i]
  63      *                +---+---+---+---+---+     +---+
  64      *              1     | # | # | # | # | --> |   | ... sum of distMat[1][2..5]
  65      *                + - +---+---+---+---+     +---+
  66      *              2         | # | # | # | --> |   | ... sum of distMat[2][3..5]
  67      *                + - + - +---+---+---+     +---+
  68      *              3             | # | # | --> |   | ... sum of distMat[3][4..5]
  69      *                + - + - + - +---+---+     +---+
  70      *              4                 | # | --> |   | ... sum of distMat[4][5]
  71      *                + - + - + - + - +---+     +---+
  72      *              5                     | --> |   | ... zero
  73      *                + - + - + - + - + - +     +---+
  74      *                  |   |   |   |   |
  75      *                  V   V   V   V   V  Calculate sum , sotre to sum[i]
  76      *                +---+---+---+---+---+
  77      *    sum_rows[i] |   |   |   |   |   |
  78      *                +---+---+---+---+---+
  79      *                  |   |   |   |   |
  80      *                  |   |   |   |   +----- sum of distMat[1..4][5]
  81      *                  |   |   |   +--------- sum of distMat[1..3][4]
  82      *                  |   |   +------------- sum of distMat[1..2][3]
  83      *                  |   +----------------- sum of distMat[1][2]
  84      *                  +--------------------- zero
  85      *
  86      *    And we use (sum_rows[i] + sum_cols[i]) instead of sum[i].
  87      *
  88      *******************
  89      * <IMPROVEMENT 2> : We manage valid nodes with chain list, instead of
  90      *                   tkill[i] flag array.
  91      *
  92      *    In original logic, invalid(killed?) nodes after nodes-joining
  93      *    are managed with tkill[i] flag array (set to 1 when killed).
  94      *    By this method, it is conspicuous to try next node but skip it
  95      *    at the latter of finding a best pair of neighbor nodes.
  96      *
  97      *    So, we thought that we managed valid nodes by using a chain list
  98      *    as below:
  99      *
 100      *    1) declare the list structure.
 101      *        struct {
 102      *            sint n;        // entry number of node.
 103      *            void *prev;        // pointer to previous entry.
 104      *            void *next;        // pointer to next entry.
 105      *        }
 106      *    2) construct a valid node list.
 107      *
 108      *       +-----+    +-----+    +-----+    +-----+        +-----+
 109      * NULL<-|prev |<---|prev |<---|prev |<---|prev |<- - - -|prev |
 110      *       |  0  |    |  1  |    |  2  |    |  3  |        |  n  |
 111      *       | next|--->| next|--->| next|--->| next|- - - ->| next|->NULL
 112      *       +-----+    +-----+    +-----+    +-----+        +-----+
 113      *
 114      *    3) when finding a best pair of neighbor nodes, we use
 115      *       this chain list as loop counter.
 116      *
 117      *    4) If an entry was killed by node-joining, this chain list is
 118      *       modified to remove that entry.
 119      *
 120      *       EX) remove the entry No 2.
 121      *       +-----+    +-----+               +-----+        +-----+
 122      * NULL<-|prev |<---|prev |<--------------|prev |<- - - -|prev |
 123      *       |  0  |    |  1  |               |  3  |        |  n  |
 124      *       | next|--->| next|-------------->| next|- - - ->| next|->NULL
 125      *       +-----+    +-----+               +-----+        +-----+
 126      *                             +-----+
 127      *                       NULL<-|prev |
 128      *                             |  2  |
 129      *                             | next|->NULL
 130      *                             +-----+
 131      *
 132      *    By this method, speed is up at the latter of finding a best pair of
 133      *    neighbor nodes.
 134      *
 135      *******************
 136      * <IMPROVEMENT 3> : Cut the frequency of division.
 137      *
 138      * At comparison between 'total' and 'tmin' in the main cycle, total is
 139      * divided by (2.0*fnseqs2) before comparison.  If N nodes are available,
 140      * that division happen (N*(N-1))/2 order.
 141      *
 142      * We thought that the comparison relation between tmin and total/(2.0*fnseqs2)
 143      * is equal to the comparison relation between (tmin*2.0*fnseqs2) and total.
 144      * Calculation of (tmin*2.0*fnseqs2) is only one time. so we stop dividing
 145      * a total value and multiply tmin and (tmin*2.0*fnseqs2) instead.
 146      *
 147      *******************
 148      * <IMPROVEMENT 4> : some transformation of the equation (to cut operations).
 149      *
 150      * We transform an equation of calculating 'total' in the main cycle.
 151      *
 152      */
 153
 154 void NJTree::generateTree(clustalw::PhyloTree* phyTree,
 155                           clustalw::DistMatrix* distMat,
 156                           clustalw::SeqInfo* seqInfo,
 157                           ofstream* log)
 158 {
 159     if (log == NULL)
 160     {
 161         verbose = false;
 162     }
 163
 164     register int i;
 165     int l[4], nude, k;
 166     int nc, mini, minj, j, ii, jj;
 167     double fnseqs, fnseqs2 = 0, sumd;
 168     double diq, djq, dij, dio, djo, da;
 169     double tmin, total, dmin;
 170     double bi, bj, b1, b2, b3, branch[4];
 171     int typei, typej; /* 0 = node; 1 = OTU */
 172
 173     int firstSeq = seqInfo->firstSeq;
 174     int lastSeq = seqInfo->lastSeq;
 175     int numSeqs = seqInfo->numSeqs;
 176
 177     /* IMPROVEMENT 1, STEP 0 : declare  variables */
 178     double *sumCols,  *sumRows,  *join;
 179
 180     sumCols = new double[numSeqs + 1];
 181     sumRows = new double[numSeqs + 1];
 182     join = new double[numSeqs + 1];
 183
 184     /* IMPROVEMENT 2, STEP 0 : declare  variables */
 185     int loop_limit;
 186     typedef struct _ValidNodeID
 187     {
 188         int n;
 189         struct _ValidNodeID *prev;
 190         struct _ValidNodeID *next;
 191     } ValidNodeID;
 192     ValidNodeID *tvalid,  *lpi,  *lpj,  *lpii,  *lpjj,  *lp_prev,  *lp_next;
 193
 194     /*
 195      * correspondence of the loop counter variables.
 196      *   i .. lpi->n,    ii .. lpii->n
 197      *   j .. lpj->n,    jj .. lpjj->n
 198      */
 199
 200     fnseqs = (double)lastSeq - firstSeq + 1;
 201
 202     /*********************** First initialisation ***************************/
 203
 204     if (verbose)
 205     {
 206         (*log) << "\n\n\t\t\tNeighbor-joining Method\n"
 207                 << "\n Saitou, N. and Nei, M. (1987)" << " The Neighbor-joining Method:"
 208                 << "\n A New Method for Reconstructing Phylogenetic Trees."
 209                 << "\n Mol. Biol. Evol., 4(4), 406-425\n" << "\n\n This is an UNROOTED tree\n"
 210                 << "\n Numbers in parentheses are branch lengths\n\n";
 211     }
 212
 213     if (fnseqs == 2)
 214     {
 215         if (verbose)
 216         {
 217             (*log) << "Cycle   1     =  SEQ:   1 (" << setw(9) << setprecision(5)
 218                     << (*distMat)(firstSeq, firstSeq + 1)
 219                     << ") joins  SEQ:   2 ("
 220                     << setw(9) << setprecision(5)
 221                     << (*distMat)(firstSeq, firstSeq + 1) << ")";
 222         }
 223         return ;
 224     }
 225
 226     mini = minj = 0;
 227
 228     /* IMPROVEMENT 1, STEP 1 : Allocate memory */
 229     /* IMPROVEMENT 1, STEP 2 : Initialize arrays to 0 */
 230     phyTree->leftBranch.resize(numSeqs + 2, 0.0);
 231     phyTree->rightBranch.resize(numSeqs + 2, 0.0);
 232     tkill.resize(numSeqs + 1, 0);
 233     av.resize(numSeqs + 1, 0.0);
 234
 235     /* IMPROVEMENT 2, STEP 1 : Allocate memory */
 236
 237     tvalid = new ValidNodeID[numSeqs + 1];
 238
 239     /* tvalid[0] is special entry in array. it points a header of valid entry list */
 240     tvalid[0].n = 0;
 241     tvalid[0].prev = NULL;
 242     tvalid[0].next = &tvalid[1];
 243
 244     /* IMPROVEMENT 2, STEP 2 : Construct and initialize the entry chain list */
 245     for (i = 1, loop_limit = lastSeq - firstSeq + 1, lpi = &tvalid[1],
 246         lp_prev = &tvalid[0], lp_next = &tvalid[2]; i <= loop_limit; ++i,
 247         ++lpi, ++lp_prev, ++lp_next)
 248     {
 249         (*distMat)(i, i) = av[i] = 0.0;
 250         tkill[i] = 0;
 251         lpi->n = i;
 252         lpi->prev = lp_prev;
 253         lpi->next = lp_next;
 254     }
 255     tvalid[loop_limit].next = NULL;
 256
 257     /*
 258      * IMPROVEMENT 1, STEP 3 : Calculate the sum of score value that
 259      * is sequence[i] to others.
 260      */
 261     double matValue;
 262     sumd = 0.0;
 263     for (lpj = tvalid[0].next; lpj != NULL; lpj = lpj->next)
 264     {
 265         double tmp_sum = 0.0;
 266         j = lpj->n;
 267         /* calculate sumRows[j] */
 268         for (lpi = tvalid[0].next; lpi->n < j; lpi = lpi->next)
 269         {
 270             i = lpi->n;
 271             matValue = (*distMat)(i, j);
 272             tmp_sum = tmp_sum + matValue;
 273         }
 274         sumRows[j] = tmp_sum;
 275
 276         tmp_sum = 0.0;
 277         /* Set lpi to that lpi->n is greater than j */
 278         if ((lpi != NULL) && (lpi->n == j))
 279         {
 280             lpi = lpi->next;
 281         }
 282         /* calculate sumCols[j] */
 283         for (; lpi != NULL; lpi = lpi->next)
 284         {
 285             i = lpi->n;
 286             tmp_sum += (*distMat)(j, i);
 287         }
 288         sumCols[j] = tmp_sum;
 289     }
 290
 291     /*********************** Enter The Main Cycle ***************************/
 292
 293     for (nc = 1, loop_limit = (lastSeq - firstSeq + 1-3); nc <= loop_limit; ++nc)
 294     {
 295         sumd = 0.0;
 296         /* IMPROVEMENT 1, STEP 4 : use sum value */
 297         for (lpj = tvalid[0].next; lpj != NULL; lpj = lpj->next)
 298         {
 299             sumd += sumCols[lpj->n];
 300         }
 301
 302         /* IMPROVEMENT 3, STEP 0 : multiply tmin and 2*fnseqs2 */
 303         fnseqs2 = fnseqs - 2.0; /* Set fnseqs2 at this point. */
 304         tmin = 99999.0 * 2.0 * fnseqs2;
 305
 306         /*.................compute SMATij values and find the smallest one ........*/
 307
 308         mini = minj = 0;
 309
 310         /* jj must starts at least 2 */
 311         if ((tvalid[0].next != NULL) && (tvalid[0].next->n == 1))
 312         {
 313             lpjj = tvalid[0].next->next;
 314         }
 315         else
 316         {
 317             lpjj = tvalid[0].next;
 318         }
 319
 320         for (; lpjj != NULL; lpjj = lpjj->next)
 321         {
 322             jj = lpjj->n;
 323             for (lpii = tvalid[0].next; lpii->n < jj; lpii = lpii->next)
 324             {
 325                 ii = lpii->n;
 326                 diq = djq = 0.0;
 327
 328                 /* IMPROVEMENT 1, STEP 4 : use sum value */
 329                 diq = sumCols[ii] + sumRows[ii];
 330                 djq = sumCols[jj] + sumRows[jj];
 331                 /*
 332                  * always ii < jj in this point. Use upper
 333                  * triangle of score matrix.
 334                  */
 335                 dij = (*distMat)(ii, jj);
 336                 /*
 337                  * IMPROVEMENT 3, STEP 1 : fnseqs2 is
 338                  * already calculated.
 339                  */
 340                 /* fnseqs2 = fnseqs - 2.0 */
 341
 342                 /* IMPROVEMENT 4 : transform the equation */
 343                 /*-------------------------------------------------------------------*
 344                  * OPTIMIZE of expression 'total = d2r + fnseqs2*dij + dr*2.0'       *
 345                  * total = d2r + fnseq2*dij + 2.0*dr                                 *
 346                  *       = d2r + fnseq2*dij + 2(sumd - dij - d2r)                    *
 347                  *       = d2r + fnseq2*dij + 2*sumd - 2*dij - 2*d2r                 *
 348                  *       =       fnseq2*dij + 2*sumd - 2*dij - 2*d2r + d2r           *
 349                  *       = fnseq2*dij + 2*sumd - 2*dij - d2r                         *
 350                  *       = fnseq2*dij + 2*sumd - 2*dij - (diq + djq - 2*dij)         *
 351                  *       = fnseq2*dij + 2*sumd - 2*dij - diq - djq + 2*dij           *
 352                  *       = fnseq2*dij + 2*sumd - 2*dij + 2*dij - diq - djq           *
 353                  *       = fnseq2*dij + 2*sumd  - diq - djq                          *
 354                  *-------------------------------------------------------------------*/
 355                 total = fnseqs2 * dij + 2.0 * sumd - diq - djq;
 356                 /*
 357                  * IMPROVEMENT 3, STEP 2 : abbrevlate
 358                  * the division on comparison between
 359                  * total and tmin.
 360                  */
 361                 /* total = total / (2.0*fnseqs2); */
 362
 363                 if (total < tmin)
 364                 {
 365                     tmin = total;
 366                     mini = ii;
 367                     minj = jj;
 368                 }
 369             }
 370         }
 371
 372         /* MEMO: always ii < jj in avobe loop, so mini < minj */
 373
 374         /*.................compute branch lengths and print the results ........*/
 375
 376
 377         dio = djo = 0.0;
 378
 379         /* IMPROVEMENT 1, STEP 4 : use sum value */
 380         dio = sumCols[mini] + sumRows[mini];
 381         djo = sumCols[minj] + sumRows[minj];
 382
 383         dmin = (*distMat)(mini, minj);
 384         dio = (dio - dmin) / fnseqs2;
 385         djo = (djo - dmin) / fnseqs2;
 386         bi = (dmin + dio - djo) *0.5;
 387         bj = dmin - bi;
 388         bi = bi - av[mini];
 389         bj = bj - av[minj];
 390
 391 #if 0
 392         (*log) << endl << "***  cycle " << nc << endl;
 393         (*log) << "dmin     = " << setw(9) << right << setprecision(9) << dmin << endl;
 394         (*log) << "dio      = " << setw(9) << right << setprecision(9) <<  dio << endl;
 395         (*log) << "djo      = " << setw(9) << right << setprecision(9) <<  djo << endl;
 396         (*log) << "bi       = " << setw(9) << right << setprecision(9) << bi << endl;
 397         (*log) << "bj       = " << setw(9) << right << setprecision(9) <<  bj << endl;
 398         (*log) << "mini     = " << setw(9) << mini << endl;
 399         (*log) << "minj     = " << setw(9) << minj << endl;
 400         (*log) << "av[minj] = " << setw(9) << right << setprecision(9) <<  av[minj] << endl;
 401         (*log) << "av[minj] = " << setw(9) << right << setprecision(9) << av[minj] << endl;
 402 #endif
 403         if (av[mini] > 0.0)
 404         {
 405             typei = 0;
 406         }
 407         else
 408         {
 409             typei = 1;
 410         }
 411         if (av[minj] > 0.0)
 412         {
 413             typej = 0;
 414         }
 415         else
 416         {
 417             typej = 1;
 418         }
 419
 420         if (verbose)
 421         {
 422             (*log) <<  "\n Cycle" << setw(4) << nc << "     = ";
 423         }
 424
 425         /*
 426          set (tiny? (AW&FS)) negative branch lengths to zero.  Also set any tiny positive
 427          branch lengths to zero.
 428          */
 429         if (fabs(bi) < 0.0001)
 430         {
 431             bi = 0.0;
 432         }
 433         if (fabs(bj) < 0.0001)
 434         {
 435             bj = 0.0;
 436         }
 437
 438         if (verbose)
 439         {
 440             if (typei == 0)
 441             {
 442                 (*log) <<  "Node:" << setw(4) << mini << " (" << setw(9) << setprecision(5)
 443                         << bi << ") joins ";
 444             }
 445             else
 446             {
 447                 (*log) << " SEQ:" << setw(4) << mini << " (" << setw(9) << setprecision(5)
 448                         << bi << ") joins ";
 449             }
 450
 451             if (typej == 0)
 452             {
 453                 (*log) << "Node:" << setw(4) << minj << " (" << setw(9) << setprecision(5)
 454                         << bj << ")";
 455             }
 456             else
 457             {
 458                 (*log) << " SEQ:" << setw(4) << minj << " (" << setw(9) << setprecision(5)
 459                         << bj << ")";
 460             }
 461
 462             (*log) << "\n";
 463         }
 464
 465         phyTree->leftBranch[nc] = bi;
 466         phyTree->rightBranch[nc] = bj;
 467
 468         for (i = 1; i <= lastSeq - firstSeq + 1; i++)
 469         {
 470             phyTree->treeDesc[nc][i] = 0;
 471         }
 472
 473         if (typei == 0)
 474         {
 475             for (i = nc - 1; i >= 1; i--)
 476                 if (phyTree->treeDesc[i][mini] == 1)
 477                 {
 478                     for (j = 1; j <= lastSeq - firstSeq + 1; j++)
 479                         if (phyTree->treeDesc[i][j] == 1)
 480                         {
 481                             phyTree->treeDesc[nc][j] = 1;
 482                         }
 483                     break;
 484                 }
 485         }
 486         else
 487         {
 488             phyTree->treeDesc[nc][mini] = 1;
 489         }
 490
 491         if (typej == 0)
 492         {
 493             for (i = nc - 1; i >= 1; i--)
 494                 if (phyTree->treeDesc[i][minj] == 1)
 495                 {
 496                     for (j = 1; j <= lastSeq - firstSeq + 1; j++)
 497                         if (phyTree->treeDesc[i][j] == 1)
 498                         {
 499                             phyTree->treeDesc[nc][j] = 1;
 500                         }
 501                     break;
 502                 }
 503         }
 504         else
 505         {
 506             phyTree->treeDesc[nc][minj] = 1;
 507         }
 508
 509
 510         /*
 511         Here is where the -0.00005 branch lengths come from for 3 or more
 512         identical seqs.
 513          */
 514         /*        if(dmin <= 0.0) dmin = 0.0001; */
 515         if (dmin <= 0.0)
 516         {
 517             dmin = 0.000001;
 518         }
 519         av[mini] = dmin * 0.5;
 520
 521         /*........................Re-initialisation................................*/
 522
 523         fnseqs = fnseqs - 1.0;
 524         tkill[minj] = 1;
 525
 526         /* IMPROVEMENT 2, STEP 3 : Remove tvalid[minj] from chain list. */
 527         /* [ Before ]
 528          *  +---------+        +---------+        +---------+
 529          *  |prev     |<-------|prev     |<-------|prev     |<---
 530          *  |    n    |        | n(=minj)|        |    n    |
 531          *  |     next|------->|     next|------->|     next|----
 532          *  +---------+        +---------+        +---------+
 533          *
 534          * [ After ]
 535          *  +---------+                           +---------+
 536          *  |prev     |<--------------------------|prev     |<---
 537          *  |    n    |                           |    n    |
 538          *  |     next|-------------------------->|     next|----
 539          *  +---------+                           +---------+
 540          *                     +---------+
 541          *              NULL---|prev     |
 542          *                     | n(=minj)|
 543          *                     |     next|---NULL
 544          *                     +---------+
 545          */
 546         (tvalid[minj].prev)->next = tvalid[minj].next;
 547         if (tvalid[minj].next != NULL)
 548         {
 549             (tvalid[minj].next)->prev = tvalid[minj].prev;
 550         }
 551         tvalid[minj].prev = tvalid[minj].next = NULL;
 552
 553         /* IMPROVEMENT 1, STEP 5 : re-calculate sum values. */
 554         for (lpj = tvalid[0].next; lpj != NULL; lpj = lpj->next)
 555         {
 556             double tmp_di = 0.0;
 557             double tmp_dj = 0.0;
 558             j = lpj->n;
 559
 560             /*
 561              * subtrace a score value related with 'minj' from
 562              * sum arrays .
 563              */
 564             if (j < minj)
 565             {
 566                 tmp_dj = (*distMat)(j, minj);
 567                 sumCols[j] -= tmp_dj;
 568             }
 569             else if (j > minj)
 570             {
 571                 tmp_dj = (*distMat)(minj, j);
 572                 sumRows[j] -= tmp_dj;
 573             } /* nothing to do when j is equal to minj. */
 574
 575
 576             /*
 577              * subtrace a score value related with 'mini' from
 578              * sum arrays .
 579              */
 580             if (j < mini)
 581             {
 582                 tmp_di = (*distMat)(j, mini);
 583                 sumCols[j] -= tmp_di;
 584             }
 585             else if (j > mini)
 586             {
 587                 tmp_di = (*distMat)(mini, j);
 588                 sumRows[j] -= tmp_di;
 589             } /* nothing to do when j is equal to mini. */
 590
 591             /*
 592              * calculate a score value of the new inner node.
 593              * then, store it temporary to join[] array.
 594              */
 595             join[j] = (tmp_dj + tmp_di) *0.5;
 596         }
 597
 598         /*
 599          * 1)
 600          * Set the score values (stored in join[]) into the matrix,
 601          * row/column position is 'mini'.
 602          * 2)
 603          * Add a score value of the new inner node to sum arrays.
 604          */
 605         for (lpj = tvalid[0].next; lpj != NULL; lpj = lpj->next)
 606         {
 607             j = lpj->n;
 608             if (j < mini)
 609             {
 610                 distMat->SetAt(j, mini, join[j]);
 611                 sumCols[j] += join[j];
 612             }
 613             else if (j > mini)
 614             {
 615                 distMat->SetAt(mini, j, join[j]);
 616                 sumRows[j] += join[j];
 617             } /* nothing to do when j is equal to mini. */
 618         }
 619
 620         /* Re-calculate sumRows[mini],sumCols[mini]. */
 621         sumCols[mini] = sumRows[mini] = 0.0;
 622
 623         /* calculate sumRows[mini] */
 624         da = 0.0;
 625         for (lpj = tvalid[0].next; lpj->n < mini; lpj = lpj->next)
 626         {
 627             da = da + join[lpj->n];
 628         }
 629         sumRows[mini] = da;
 630
 631         /* skip if 'lpj->n' is equal to 'mini' */
 632         if ((lpj != NULL) && (lpj->n == mini))
 633         {
 634             lpj = lpj->next;
 635         }
 636
 637         /* calculate sumCols[mini] */
 638         da = 0.0;
 639         for (; lpj != NULL; lpj = lpj->next)
 640         {
 641             da = da + join[lpj->n];
 642         }
 643         sumCols[mini] = da;
 644
 645         /*
 646          * Clean up sumRows[minj], sumCols[minj] and score matrix
 647          * related with 'minj'.
 648          */
 649         sumCols[minj] = sumRows[minj] = 0.0;
 650         for (j = 1; j <= lastSeq - firstSeq + 1; ++j)
 651         {
 652             distMat->SetAt(minj, j, 0.0);
 653             distMat->SetAt(j, minj, 0.0);
 654             join[j] = 0.0;
 655         }
 656
 657     } /** end main cycle **/
 658
 659     /******************************Last Cycle (3 Seqs. left)********************/
 660
 661     nude = 1;
 662
 663     for (lpi = tvalid[0].next; lpi != NULL; lpi = lpi->next)
 664     {
 665         l[nude] = lpi->n;
 666         ++nude;
 667     }
 668
 669     b1 = ((*distMat)(l[1], l[2]) + (*distMat)(l[1], l[3]) - (*distMat)(l[2], l[3])) *0.5;
 670     b2 = (*distMat)(l[1], l[2]) - b1;
 671     b3 = (*distMat)(l[1], l[3]) - b1;
 672
 673     branch[1] = b1 - av[l[1]];
 674     branch[2] = b2 - av[l[2]];
 675     branch[3] = b3 - av[l[3]];
 676
 677     /* Reset tiny negative and positive branch lengths to zero */
 678     if (fabs(branch[1]) < 0.0001)
 679     {
 680         branch[1] = 0.0;
 681     }
 682     if (fabs(branch[2]) < 0.0001)
 683     {
 684         branch[2] = 0.0;
 685     }
 686     if (fabs(branch[3]) < 0.0001)
 687     {
 688         branch[3] = 0.0;
 689     }
 690
 691     phyTree->leftBranch[lastSeq - firstSeq + 1-2] = branch[1];
 692     phyTree->leftBranch[lastSeq - firstSeq + 1-1] = branch[2];
 693     phyTree->leftBranch[lastSeq - firstSeq + 1] = branch[3];
 694
 695     for (i = 1; i <= lastSeq - firstSeq + 1; i++)
 696     {
 697         phyTree->treeDesc[lastSeq - firstSeq + 1-2][i] = 0;
 698     }
 699
 700     if (verbose)
 701     {
 702         (*log) << "\n Cycle" << setw(4) << nc << " (Last cycle, trichotomy):\n";
 703     }
 704
 705     for (i = 1; i <= 3; ++i)
 706     {
 707         if (av[l[i]] > 0.0)
 708         {
 709             if (verbose)
 710             {
 711                 (*log) << "\n\t\t Node:" << setw(4) << l[i] <<" (" << setw(9)
 712                         << setprecision(5) << branch[i] << ") ";
 713             }
 714             for (k = lastSeq - firstSeq + 1-3; k >= 1; k--)
 715                 if (phyTree->treeDesc[k][l[i]] == 1)
 716                 {
 717                     for (j = 1; j <= lastSeq - firstSeq + 1; j++)
 718                         if (phyTree->treeDesc[k][j] == 1)
 719                         {
 720                             phyTree->treeDesc[lastSeq - firstSeq + 1-2][j] = i;
 721                         }
 722                     break;
 723                 }
 724         }
 725         else
 726         {
 727             if (verbose)
 728             {
 729                 (*log) << "\n\t\t  SEQ:" << setw(4) << l[i] << " (" << setw(9)
 730                         << setprecision(5) << branch[i] << ") ";
 731             }
 732             phyTree->treeDesc[lastSeq - firstSeq + 1-2][l[i]] = i;
 733         }
 734         if (i < 3)
 735         {
 736             if (verbose)
 737             {
 738                 (*log) << "joins";
 739             }
 740         }
 741     }
 742
 743     if (verbose)
 744     {
 745         (*log) << "\n";
 746     }
 747
 748     /* IMPROVEMENT 2, STEP 4 : release memory area */
 749
 750     delete [] tvalid;
 751     delete [] sumCols;
 752     delete [] sumRows;
 753     delete [] join;
 754
 755
 756 }
 757 }