--- /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;
+
+
+}
+}
git://source.jalview.org / jabaws.git / blobdiff