+++ /dev/null
-/**
- * Author: Mark Larkin
- *
- * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
- */
-#ifdef HAVE_CONFIG_H
- #include "config.h"
-#endif
-#include <math.h>
-#include "NJTree.h"
-
-namespace clustalw
-{
-
- /****************************************************************************
- * [ Improvement ideas in fast_nj_tree() ] by DDBJ & FUJITSU Limited.
- * written by Tadashi Koike
- * (takoike@genes.nig.ac.jp)
- *******************
- * <IMPROVEMENT 1> : Store the value of sum of the score to temporary array,
- * and use again and again.
- *
- * In the main cycle, these are calculated again and again :
- * diq = sum of distMat[n][ii] (n:1 to lastSeq-firstSeq+1),
- * djq = sum of distMat[n][jj] (n:1 to lastSeq-firstSeq+1),
- * dio = sum of distMat[n][mini] (n:1 to lastSeq-firstSeq+1),
- * djq = sum of distMat[n][minj] (n:1 to lastSeq-firstSeq+1)
- * // 'lastSeq' and 'firstSeq' are both constant values //
- * and the result of above calculations is always same until
- * a best pair of neighbour nodes is joined.
- *
- * So, we change the logic to calculate the sum[i] (=sum of distMat[n][i]
- * (n:1 to lastSeq-firstSeq+1)) and store it to array, before
- * beginning to find a best pair of neighbour nodes, and after that
- * we use them again and again.
- *
- * tmat[i][j]
- * 1 2 3 4 5
- * +---+---+---+---+---+
- * 1 | | | | | |
- * +---+---+---+---+---+
- * 2 | | | | | | 1) calculate sum of distMat[n][i]
- * +---+---+---+---+---+ (n: 1 to lastSeq-firstSeq+1)
- * 3 | | | | | | 2) store that sum value to sum[i]
- * +---+---+---+---+---+
- * 4 | | | | | | 3) use sum[i] during finding a best
- * +---+---+---+---+---+ pair of neibour nodes.
- * 5 | | | | | |
- * +---+---+---+---+---+
- * | | | | |
- * V V V V V Calculate sum , and store it to sum[i]
- * +---+---+---+---+---+
- * sum[i] | | | | | |
- * +---+---+---+---+---+
- *
- * At this time, we thought that we use upper triangle of the matrix
- * because distMat[i][j] is equal to distMat[j][i] and distMat[i][i] is equal
- * to zero. Therefore, we prepared sum_rows[i] and sum_cols[i] instead
- * of sum[i] for storing the sum value.
- *
- * distMat[i][j]
- * 1 2 3 4 5 sum_cols[i]
- * +---+---+---+---+---+ +---+
- * 1 | # | # | # | # | --> | | ... sum of distMat[1][2..5]
- * + - +---+---+---+---+ +---+
- * 2 | # | # | # | --> | | ... sum of distMat[2][3..5]
- * + - + - +---+---+---+ +---+
- * 3 | # | # | --> | | ... sum of distMat[3][4..5]
- * + - + - + - +---+---+ +---+
- * 4 | # | --> | | ... sum of distMat[4][5]
- * + - + - + - + - +---+ +---+
- * 5 | --> | | ... zero
- * + - + - + - + - + - + +---+
- * | | | | |
- * V V V V V Calculate sum , sotre to sum[i]
- * +---+---+---+---+---+
- * sum_rows[i] | | | | | |
- * +---+---+---+---+---+
- * | | | | |
- * | | | | +----- sum of distMat[1..4][5]
- * | | | +--------- sum of distMat[1..3][4]
- * | | +------------- sum of distMat[1..2][3]
- * | +----------------- sum of distMat[1][2]
- * +--------------------- zero
- *
- * And we use (sum_rows[i] + sum_cols[i]) instead of sum[i].
- *
- *******************
- * <IMPROVEMENT 2> : We manage valid nodes with chain list, instead of
- * tkill[i] flag array.
- *
- * In original logic, invalid(killed?) nodes after nodes-joining
- * are managed with tkill[i] flag array (set to 1 when killed).
- * By this method, it is conspicuous to try next node but skip it
- * at the latter of finding a best pair of neighbor nodes.
- *
- * So, we thought that we managed valid nodes by using a chain list
- * as below:
- *
- * 1) declare the list structure.
- * struct {
- * sint n; // entry number of node.
- * void *prev; // pointer to previous entry.
- * void *next; // pointer to next entry.
- * }
- * 2) construct a valid node list.
- *
- * +-----+ +-----+ +-----+ +-----+ +-----+
- * NULL<-|prev |<---|prev |<---|prev |<---|prev |<- - - -|prev |
- * | 0 | | 1 | | 2 | | 3 | | n |
- * | next|--->| next|--->| next|--->| next|- - - ->| next|->NULL
- * +-----+ +-----+ +-----+ +-----+ +-----+
- *
- * 3) when finding a best pair of neighbor nodes, we use
- * this chain list as loop counter.
- *
- * 4) If an entry was killed by node-joining, this chain list is
- * modified to remove that entry.
- *
- * EX) remove the entry No 2.
- * +-----+ +-----+ +-----+ +-----+
- * NULL<-|prev |<---|prev |<--------------|prev |<- - - -|prev |
- * | 0 | | 1 | | 3 | | n |
- * | next|--->| next|-------------->| next|- - - ->| next|->NULL
- * +-----+ +-----+ +-----+ +-----+
- * +-----+
- * NULL<-|prev |
- * | 2 |
- * | next|->NULL
- * +-----+
- *
- * By this method, speed is up at the latter of finding a best pair of
- * neighbor nodes.
- *
- *******************
- * <IMPROVEMENT 3> : Cut the frequency of division.
- *
- * At comparison between 'total' and 'tmin' in the main cycle, total is
- * divided by (2.0*fnseqs2) before comparison. If N nodes are available,
- * that division happen (N*(N-1))/2 order.
- *
- * We thought that the comparison relation between tmin and total/(2.0*fnseqs2)
- * is equal to the comparison relation between (tmin*2.0*fnseqs2) and total.
- * Calculation of (tmin*2.0*fnseqs2) is only one time. so we stop dividing
- * a total value and multiply tmin and (tmin*2.0*fnseqs2) instead.
- *
- *******************
- * <IMPROVEMENT 4> : some transformation of the equation (to cut operations).
- *
- * We transform an equation of calculating 'total' in the main cycle.
- *
- */
-
-void NJTree::generateTree(clustalw::PhyloTree* phyTree,
- clustalw::DistMatrix* distMat,
- clustalw::SeqInfo* seqInfo,
- ofstream* log)
-{
- if (log == NULL)
- {
- verbose = false;
- }
-
- register int i;
- int l[4], nude, k;
- int nc, mini, minj, j, ii, jj;
- double fnseqs, fnseqs2 = 0, sumd;
- double diq, djq, dij, dio, djo, da;
- double tmin, total, dmin;
- double bi, bj, b1, b2, b3, branch[4];
- int typei, typej; /* 0 = node; 1 = OTU */
-
- int firstSeq = seqInfo->firstSeq;
- int lastSeq = seqInfo->lastSeq;
- int numSeqs = seqInfo->numSeqs;
-
- /* IMPROVEMENT 1, STEP 0 : declare variables */
- double *sumCols, *sumRows, *join;
-
- sumCols = new double[numSeqs + 1];
- sumRows = new double[numSeqs + 1];
- join = new double[numSeqs + 1];
-
- /* IMPROVEMENT 2, STEP 0 : declare variables */
- int loop_limit;
- typedef struct _ValidNodeID
- {
- int n;
- struct _ValidNodeID *prev;
- struct _ValidNodeID *next;
- } ValidNodeID;
- ValidNodeID *tvalid, *lpi, *lpj, *lpii, *lpjj, *lp_prev, *lp_next;
-
- /*
- * correspondence of the loop counter variables.
- * i .. lpi->n, ii .. lpii->n
- * j .. lpj->n, jj .. lpjj->n
- */
-
- fnseqs = (double)lastSeq - firstSeq + 1;
-
- /*********************** First initialisation ***************************/
-
- if (verbose)
- {
- (*log) << "\n\n\t\t\tNeighbor-joining Method\n"
- << "\n Saitou, N. and Nei, M. (1987)" << " The Neighbor-joining Method:"
- << "\n A New Method for Reconstructing Phylogenetic Trees."
- << "\n Mol. Biol. Evol., 4(4), 406-425\n" << "\n\n This is an UNROOTED tree\n"
- << "\n Numbers in parentheses are branch lengths\n\n";
- }
-
- if (fnseqs == 2)
- {
- if (verbose)
- {
- (*log) << "Cycle 1 = SEQ: 1 (" << setw(9) << setprecision(5)
- << (*distMat)(firstSeq, firstSeq + 1)
- << ") joins SEQ: 2 ("
- << setw(9) << setprecision(5)
- << (*distMat)(firstSeq, firstSeq + 1) << ")";
- }
- return ;
- }
-
- mini = minj = 0;
-
- /* IMPROVEMENT 1, STEP 1 : Allocate memory */
- /* IMPROVEMENT 1, STEP 2 : Initialize arrays to 0 */
- phyTree->leftBranch.resize(numSeqs + 2, 0.0);
- phyTree->rightBranch.resize(numSeqs + 2, 0.0);
- tkill.resize(numSeqs + 1, 0);
- av.resize(numSeqs + 1, 0.0);
-
- /* IMPROVEMENT 2, STEP 1 : Allocate memory */
-
- tvalid = new ValidNodeID[numSeqs + 1];
-
- /* tvalid[0] is special entry in array. it points a header of valid entry list */
- tvalid[0].n = 0;
- tvalid[0].prev = NULL;
- tvalid[0].next = &tvalid[1];
-
- /* IMPROVEMENT 2, STEP 2 : Construct and initialize the entry chain list */
- for (i = 1, loop_limit = lastSeq - firstSeq + 1, lpi = &tvalid[1],
- lp_prev = &tvalid[0], lp_next = &tvalid[2]; i <= loop_limit; ++i,
- ++lpi, ++lp_prev, ++lp_next)
- {
- (*distMat)(i, i) = av[i] = 0.0;
- tkill[i] = 0;
- lpi->n = i;
- lpi->prev = lp_prev;
- lpi->next = lp_next;
- }
- tvalid[loop_limit].next = NULL;
-
- /*
- * IMPROVEMENT 1, STEP 3 : Calculate the sum of score value that
- * is sequence[i] to others.
- */
- double matValue;
- sumd = 0.0;
- for (lpj = tvalid[0].next; lpj != NULL; lpj = lpj->next)
- {
- double tmp_sum = 0.0;
- j = lpj->n;
- /* calculate sumRows[j] */
- for (lpi = tvalid[0].next; lpi->n < j; lpi = lpi->next)
- {
- i = lpi->n;
- matValue = (*distMat)(i, j);
- tmp_sum = tmp_sum + matValue;
- }
- sumRows[j] = tmp_sum;
-
- tmp_sum = 0.0;
- /* Set lpi to that lpi->n is greater than j */
- if ((lpi != NULL) && (lpi->n == j))
- {
- lpi = lpi->next;
- }
- /* calculate sumCols[j] */
- for (; lpi != NULL; lpi = lpi->next)
- {
- i = lpi->n;
- tmp_sum += (*distMat)(j, i);
- }
- sumCols[j] = tmp_sum;
- }
-
- /*********************** Enter The Main Cycle ***************************/
-
- for (nc = 1, loop_limit = (lastSeq - firstSeq + 1-3); nc <= loop_limit; ++nc)
- {
- sumd = 0.0;
- /* IMPROVEMENT 1, STEP 4 : use sum value */
- for (lpj = tvalid[0].next; lpj != NULL; lpj = lpj->next)
- {
- sumd += sumCols[lpj->n];
- }
-
- /* IMPROVEMENT 3, STEP 0 : multiply tmin and 2*fnseqs2 */
- fnseqs2 = fnseqs - 2.0; /* Set fnseqs2 at this point. */
- tmin = 99999.0 * 2.0 * fnseqs2;
-
- /*.................compute SMATij values and find the smallest one ........*/
-
- mini = minj = 0;
-
- /* jj must starts at least 2 */
- if ((tvalid[0].next != NULL) && (tvalid[0].next->n == 1))
- {
- lpjj = tvalid[0].next->next;
- }
- else
- {
- lpjj = tvalid[0].next;
- }
-
- for (; lpjj != NULL; lpjj = lpjj->next)
- {
- jj = lpjj->n;
- for (lpii = tvalid[0].next; lpii->n < jj; lpii = lpii->next)
- {
- ii = lpii->n;
- diq = djq = 0.0;
-
- /* IMPROVEMENT 1, STEP 4 : use sum value */
- diq = sumCols[ii] + sumRows[ii];
- djq = sumCols[jj] + sumRows[jj];
- /*
- * always ii < jj in this point. Use upper
- * triangle of score matrix.
- */
- dij = (*distMat)(ii, jj);
- /*
- * IMPROVEMENT 3, STEP 1 : fnseqs2 is
- * already calculated.
- */
- /* fnseqs2 = fnseqs - 2.0 */
-
- /* IMPROVEMENT 4 : transform the equation */
- /*-------------------------------------------------------------------*
- * OPTIMIZE of expression 'total = d2r + fnseqs2*dij + dr*2.0' *
- * total = d2r + fnseq2*dij + 2.0*dr *
- * = d2r + fnseq2*dij + 2(sumd - dij - d2r) *
- * = d2r + fnseq2*dij + 2*sumd - 2*dij - 2*d2r *
- * = fnseq2*dij + 2*sumd - 2*dij - 2*d2r + d2r *
- * = fnseq2*dij + 2*sumd - 2*dij - d2r *
- * = fnseq2*dij + 2*sumd - 2*dij - (diq + djq - 2*dij) *
- * = fnseq2*dij + 2*sumd - 2*dij - diq - djq + 2*dij *
- * = fnseq2*dij + 2*sumd - 2*dij + 2*dij - diq - djq *
- * = fnseq2*dij + 2*sumd - diq - djq *
- *-------------------------------------------------------------------*/
- total = fnseqs2 * dij + 2.0 * sumd - diq - djq;
- /*
- * IMPROVEMENT 3, STEP 2 : abbrevlate
- * the division on comparison between
- * total and tmin.
- */
- /* total = total / (2.0*fnseqs2); */
-
- if (total < tmin)
- {
- tmin = total;
- mini = ii;
- minj = jj;
- }
- }
- }
-
- /* MEMO: always ii < jj in avobe loop, so mini < minj */
-
- /*.................compute branch lengths and print the results ........*/
-
-
- dio = djo = 0.0;
-
- /* IMPROVEMENT 1, STEP 4 : use sum value */
- dio = sumCols[mini] + sumRows[mini];
- djo = sumCols[minj] + sumRows[minj];
-
- dmin = (*distMat)(mini, minj);
- dio = (dio - dmin) / fnseqs2;
- djo = (djo - dmin) / fnseqs2;
- bi = (dmin + dio - djo) *0.5;
- bj = dmin - bi;
- bi = bi - av[mini];
- bj = bj - av[minj];
-
-#if 0
- (*log) << endl << "*** cycle " << nc << endl;
- (*log) << "dmin = " << setw(9) << right << setprecision(9) << dmin << endl;
- (*log) << "dio = " << setw(9) << right << setprecision(9) << dio << endl;
- (*log) << "djo = " << setw(9) << right << setprecision(9) << djo << endl;
- (*log) << "bi = " << setw(9) << right << setprecision(9) << bi << endl;
- (*log) << "bj = " << setw(9) << right << setprecision(9) << bj << endl;
- (*log) << "mini = " << setw(9) << mini << endl;
- (*log) << "minj = " << setw(9) << minj << endl;
- (*log) << "av[minj] = " << setw(9) << right << setprecision(9) << av[minj] << endl;
- (*log) << "av[minj] = " << setw(9) << right << setprecision(9) << av[minj] << endl;
-#endif
- if (av[mini] > 0.0)
- {
- typei = 0;
- }
- else
- {
- typei = 1;
- }
- if (av[minj] > 0.0)
- {
- typej = 0;
- }
- else
- {
- typej = 1;
- }
-
- if (verbose)
- {
- (*log) << "\n Cycle" << setw(4) << nc << " = ";
- }
-
- /*
- set (tiny? (AW&FS)) negative branch lengths to zero. Also set any tiny positive
- branch lengths to zero.
- */
- if (fabs(bi) < 0.0001)
- {
- bi = 0.0;
- }
- if (fabs(bj) < 0.0001)
- {
- bj = 0.0;
- }
-
- if (verbose)
- {
- if (typei == 0)
- {
- (*log) << "Node:" << setw(4) << mini << " (" << setw(9) << setprecision(5)
- << bi << ") joins ";
- }
- else
- {
- (*log) << " SEQ:" << setw(4) << mini << " (" << setw(9) << setprecision(5)
- << bi << ") joins ";
- }
-
- if (typej == 0)
- {
- (*log) << "Node:" << setw(4) << minj << " (" << setw(9) << setprecision(5)
- << bj << ")";
- }
- else
- {
- (*log) << " SEQ:" << setw(4) << minj << " (" << setw(9) << setprecision(5)
- << bj << ")";
- }
-
- (*log) << "\n";
- }
-
- phyTree->leftBranch[nc] = bi;
- phyTree->rightBranch[nc] = bj;
-
- for (i = 1; i <= lastSeq - firstSeq + 1; i++)
- {
- phyTree->treeDesc[nc][i] = 0;
- }
-
- if (typei == 0)
- {
- for (i = nc - 1; i >= 1; i--)
- if (phyTree->treeDesc[i][mini] == 1)
- {
- for (j = 1; j <= lastSeq - firstSeq + 1; j++)
- if (phyTree->treeDesc[i][j] == 1)
- {
- phyTree->treeDesc[nc][j] = 1;
- }
- break;
- }
- }
- else
- {
- phyTree->treeDesc[nc][mini] = 1;
- }
-
- if (typej == 0)
- {
- for (i = nc - 1; i >= 1; i--)
- if (phyTree->treeDesc[i][minj] == 1)
- {
- for (j = 1; j <= lastSeq - firstSeq + 1; j++)
- if (phyTree->treeDesc[i][j] == 1)
- {
- phyTree->treeDesc[nc][j] = 1;
- }
- break;
- }
- }
- else
- {
- phyTree->treeDesc[nc][minj] = 1;
- }
-
-
- /*
- Here is where the -0.00005 branch lengths come from for 3 or more
- identical seqs.
- */
- /* if(dmin <= 0.0) dmin = 0.0001; */
- if (dmin <= 0.0)
- {
- dmin = 0.000001;
- }
- av[mini] = dmin * 0.5;
-
- /*........................Re-initialisation................................*/
-
- fnseqs = fnseqs - 1.0;
- tkill[minj] = 1;
-
- /* IMPROVEMENT 2, STEP 3 : Remove tvalid[minj] from chain list. */
- /* [ Before ]
- * +---------+ +---------+ +---------+
- * |prev |<-------|prev |<-------|prev |<---
- * | n | | n(=minj)| | n |
- * | next|------->| next|------->| next|----
- * +---------+ +---------+ +---------+
- *
- * [ After ]
- * +---------+ +---------+
- * |prev |<--------------------------|prev |<---
- * | n | | n |
- * | next|-------------------------->| next|----
- * +---------+ +---------+
- * +---------+
- * NULL---|prev |
- * | n(=minj)|
- * | next|---NULL
- * +---------+
- */
- (tvalid[minj].prev)->next = tvalid[minj].next;
- if (tvalid[minj].next != NULL)
- {
- (tvalid[minj].next)->prev = tvalid[minj].prev;
- }
- tvalid[minj].prev = tvalid[minj].next = NULL;
-
- /* IMPROVEMENT 1, STEP 5 : re-calculate sum values. */
- for (lpj = tvalid[0].next; lpj != NULL; lpj = lpj->next)
- {
- double tmp_di = 0.0;
- double tmp_dj = 0.0;
- j = lpj->n;
-
- /*
- * subtrace a score value related with 'minj' from
- * sum arrays .
- */
- if (j < minj)
- {
- tmp_dj = (*distMat)(j, minj);
- sumCols[j] -= tmp_dj;
- }
- else if (j > minj)
- {
- tmp_dj = (*distMat)(minj, j);
- sumRows[j] -= tmp_dj;
- } /* nothing to do when j is equal to minj. */
-
-
- /*
- * subtrace a score value related with 'mini' from
- * sum arrays .
- */
- if (j < mini)
- {
- tmp_di = (*distMat)(j, mini);
- sumCols[j] -= tmp_di;
- }
- else if (j > mini)
- {
- tmp_di = (*distMat)(mini, j);
- sumRows[j] -= tmp_di;
- } /* nothing to do when j is equal to mini. */
-
- /*
- * calculate a score value of the new inner node.
- * then, store it temporary to join[] array.
- */
- join[j] = (tmp_dj + tmp_di) *0.5;
- }
-
- /*
- * 1)
- * Set the score values (stored in join[]) into the matrix,
- * row/column position is 'mini'.
- * 2)
- * Add a score value of the new inner node to sum arrays.
- */
- for (lpj = tvalid[0].next; lpj != NULL; lpj = lpj->next)
- {
- j = lpj->n;
- if (j < mini)
- {
- distMat->SetAt(j, mini, join[j]);
- sumCols[j] += join[j];
- }
- else if (j > mini)
- {
- distMat->SetAt(mini, j, join[j]);
- sumRows[j] += join[j];
- } /* nothing to do when j is equal to mini. */
- }
-
- /* Re-calculate sumRows[mini],sumCols[mini]. */
- sumCols[mini] = sumRows[mini] = 0.0;
-
- /* calculate sumRows[mini] */
- da = 0.0;
- for (lpj = tvalid[0].next; lpj->n < mini; lpj = lpj->next)
- {
- da = da + join[lpj->n];
- }
- sumRows[mini] = da;
-
- /* skip if 'lpj->n' is equal to 'mini' */
- if ((lpj != NULL) && (lpj->n == mini))
- {
- lpj = lpj->next;
- }
-
- /* calculate sumCols[mini] */
- da = 0.0;
- for (; lpj != NULL; lpj = lpj->next)
- {
- da = da + join[lpj->n];
- }
- sumCols[mini] = da;
-
- /*
- * Clean up sumRows[minj], sumCols[minj] and score matrix
- * related with 'minj'.
- */
- sumCols[minj] = sumRows[minj] = 0.0;
- for (j = 1; j <= lastSeq - firstSeq + 1; ++j)
- {
- distMat->SetAt(minj, j, 0.0);
- distMat->SetAt(j, minj, 0.0);
- join[j] = 0.0;
- }
-
- } /** end main cycle **/
-
- /******************************Last Cycle (3 Seqs. left)********************/
-
- nude = 1;
-
- for (lpi = tvalid[0].next; lpi != NULL; lpi = lpi->next)
- {
- l[nude] = lpi->n;
- ++nude;
- }
-
- b1 = ((*distMat)(l[1], l[2]) + (*distMat)(l[1], l[3]) - (*distMat)(l[2], l[3])) *0.5;
- b2 = (*distMat)(l[1], l[2]) - b1;
- b3 = (*distMat)(l[1], l[3]) - b1;
-
- branch[1] = b1 - av[l[1]];
- branch[2] = b2 - av[l[2]];
- branch[3] = b3 - av[l[3]];
-
- /* Reset tiny negative and positive branch lengths to zero */
- if (fabs(branch[1]) < 0.0001)
- {
- branch[1] = 0.0;
- }
- if (fabs(branch[2]) < 0.0001)
- {
- branch[2] = 0.0;
- }
- if (fabs(branch[3]) < 0.0001)
- {
- branch[3] = 0.0;
- }
-
- phyTree->leftBranch[lastSeq - firstSeq + 1-2] = branch[1];
- phyTree->leftBranch[lastSeq - firstSeq + 1-1] = branch[2];
- phyTree->leftBranch[lastSeq - firstSeq + 1] = branch[3];
-
- for (i = 1; i <= lastSeq - firstSeq + 1; i++)
- {
- phyTree->treeDesc[lastSeq - firstSeq + 1-2][i] = 0;
- }
-
- if (verbose)
- {
- (*log) << "\n Cycle" << setw(4) << nc << " (Last cycle, trichotomy):\n";
- }
-
- for (i = 1; i <= 3; ++i)
- {
- if (av[l[i]] > 0.0)
- {
- if (verbose)
- {
- (*log) << "\n\t\t Node:" << setw(4) << l[i] <<" (" << setw(9)
- << setprecision(5) << branch[i] << ") ";
- }
- for (k = lastSeq - firstSeq + 1-3; k >= 1; k--)
- if (phyTree->treeDesc[k][l[i]] == 1)
- {
- for (j = 1; j <= lastSeq - firstSeq + 1; j++)
- if (phyTree->treeDesc[k][j] == 1)
- {
- phyTree->treeDesc[lastSeq - firstSeq + 1-2][j] = i;
- }
- break;
- }
- }
- else
- {
- if (verbose)
- {
- (*log) << "\n\t\t SEQ:" << setw(4) << l[i] << " (" << setw(9)
- << setprecision(5) << branch[i] << ") ";
- }
- phyTree->treeDesc[lastSeq - firstSeq + 1-2][l[i]] = i;
- }
- if (i < 3)
- {
- if (verbose)
- {
- (*log) << "joins";
- }
- }
- }
-
- if (verbose)
- {
- (*log) << "\n";
- }
-
- /* IMPROVEMENT 2, STEP 4 : release memory area */
-
- delete [] tvalid;
- delete [] sumCols;
- delete [] sumRows;
- delete [] join;
-
-
-}
-}