+++ /dev/null
-#ifdef HAVE_CONFIG_H
- #include "config.h"
-#endif
-#include "UPGMAAlgorithm.h"
-#include <cmath>
-#include <ctime>
-#include <iomanip>
-#include <sstream>
-
-#include "../../general/SymMatrix.h"
-#include "../../general/debuglogObject.h"
-#include "../../general/clustalw.h"
-#include "upgmadata.h"
-
-namespace clustalw
-{
-
-UPGMAAlgorithm::UPGMAAlgorithm()
-: overwriteMatrix(false),
- numSeqs(0),
- verbose(false),
- orderNode1(0),
- orderNode2(0),
- orderNewNode(0)
-{}
-
-auto_ptr<AlignmentSteps> UPGMAAlgorithm::generateTree(RootedGuideTree* phyTree,
- DistMatrix* distMat,
- SeqInfo* seqInfo, bool overwrite,
- ofstream* tree)
-{
- if (tree == 0 || !tree->is_open())
- {
- verbose = false;
- }
-
- if (verbose)
- {
- (*tree) << "\n\n\t\t\tUPGMA Method\n"
- << "\n\n This is a ROOTED tree\n"
- << "\n Numbers in parentheses are branch lengths\n\n";
- }
-
- progSteps.reset(new AlignmentSteps);
-
- Node** clusters;
- Node* root;
- numSeqs = seqInfo->numSeqs;
- const int sizeDistMat = ((numSeqs + 1) * (numSeqs + 2)) / 2;
-
- double* elements = overwrite ?
- distMat->getDistMatrix(seqInfo->firstSeq, seqInfo->numSeqs) :
- (double *)memcpy(new double[sizeDistMat],
- distMat->getDistMatrix(seqInfo->firstSeq, seqInfo->numSeqs),
- sizeDistMat * sizeof(double));
-
- clusters = initialiseNodes(elements, seqInfo->firstSeq);
- root = doUPGMA(clusters, tree);
-
- phyTree->setRoot(root);
- delete [] clusters;
-
- if(!overwrite)
- {
- delete [] elements;
- }
- distMat->clearSubArray();
-
- return progSteps;
-}
-
-
-Node** UPGMAAlgorithm::initialiseNodes(double* distanceMatrix, int fSeq)
-{
- int firstSeq = fSeq;
-
-
- Node** nodes = new Node*[numSeqs];
- Node** nodeIter = nodes;
-
- *nodes = new Node(firstSeq, 0, 0);
-
- distanceMatrix++;
-
- // Move to first position in distanceMatrix.
- for(int elementIndex = 1, e = numSeqs; elementIndex < e;
- distanceMatrix += ++elementIndex)
- {
- Node* newcluster = new Node(elementIndex + firstSeq,
- distanceMatrix, elementIndex);
- (*nodeIter++)->next = newcluster;
- *nodeIter = newcluster;
- }
-
- return nodes;
-}
-
-void UPGMAAlgorithm::printAllNodes(Node** nodes)
-{
- int numNodes = 0;
- for(Node* nodeIter = *nodes; nodeIter; nodeIter = nodeIter->next)
- {
- numNodes++;
- cout << "Node " << numNodes << "\n";
- nodeIter->printElements();
- cout << "\n\n";
- }
- cout << "There are " << numNodes << " nodes\n";
-}
-
-Node* UPGMAAlgorithm::doUPGMA(Node** nodes, ofstream* tree)
-{
- if (tree == 0 || !tree->is_open())
- {
- verbose = false;
- }
-
- string type1, type2;
- int step = 0;
-
- while((*nodes)->next) // While there is more than 1 node.
- {
- step++;
- if (verbose)
- {
- (*tree) << "\n Cycle" << setw(4) << step << " = ";
- }
- Node** ptrNodeWithMin;
-
- /**
- * STEP 1 find node with min distance
- */
- ptrNodeWithMin = getNodeWithMinDist(nodes);
- Node* nodeToJoin2 = *ptrNodeWithMin;
- const int indexToMinDist = nodeToJoin2->indexToMinDist;
- Node* const nodeToJoin1 = nodes[indexToMinDist];
-
- orderNode1 = nodeToJoin1->size;
- orderNode2 = nodeToJoin2->size;
- orderNewNode = orderNode1 + orderNode2;
-
- /**
- * STEP 2 Recompute the dist Matrix row for nodeToJoin1
- * example: Join nodes 2 and 6
- *
- * 1 0
- * 2 X 0
- * 3 0 0 0
- * 4 0 0 0 0
- * 5 0 0 0 0 0
- * 6 0 0 0 0 0 0
- * 7 0 0 0 0 0 0 0
- */
- double* nodeToJoin2DistIter = nodeToJoin2->ptrToDistMatRow;
-
- if (indexToMinDist != 0)
- {
- recomputeNodeToJoin1DistMatRow(nodeToJoin1, &nodeToJoin2DistIter);
- }
-
- /**
- * STEP 3 Recompute all distances in column
- * example: Join nodes 2 and 6
- *
- * 1 0
- * 2 C 0
- * 3 0 X 0
- * 4 0 X 0 0
- * 5 0 X 0 0 0
- * 6 0 0 0 0 0 0
- * 7 0 X 0 0 0 0 0
- */
-
- computeAllOtherDistsToNewNode(nodeToJoin1, nodeToJoin2, &nodeToJoin2DistIter);
-
- /**
- * STEP 4 Add the step to the progSteps.
- * This creates the multiple alignment steps.
- */
- addAlignmentStep(&nodeToJoin1->allElements, &nodeToJoin2->allElements);
-
-
- double minDistance = (*ptrNodeWithMin)->minDist;
-
- double height = 0.0;
- height = minDistance / 2.0;
-
- if(verbose)
- {
- if(nodeToJoin1->allElements.size() > 1)
- {
- type1 = "NODE: ";
- }
- else
- {
- type1 = "SEQ: ";
- }
- if(nodeToJoin2->allElements.size() > 1)
- {
- type2 = "NODE: ";
- }
- else
- {
- type2 = "SEQ: ";
- }
- (*tree) << type1 << nodeToJoin1->allElements[0] << " (" << setw(9)
- << setprecision(5) << height << ") joins " << type2
- << setw(4) << nodeToJoin2->allElements[0] << " ("
- << setw(9) << setprecision(5) << height << ")";
- }
- /**
- * STEP 5 merge 2 nodes
- */
- nodeToJoin1->merge(ptrNodeWithMin, height);
-
- }
-
- return *nodes;
-}
-
-/**
- * This function returns a Node object that has the minimum distance of
- * all the nodes left.
- */
-Node** UPGMAAlgorithm::getNodeWithMinDist(Node** nodes)
-{
- Node** ptrMinNode = NULL;
- double minDistance = numeric_limits<double>::max();
- //Start at node 1, check see if it points to something, then move on the next one.
- Node** nodeIter;
- for(nodeIter = &((*nodes)->next); *nodeIter;
- nodeIter = &(*nodeIter)->next)
- {
- if ((*nodeIter)->getMinDist() < minDistance)
- {
- minDistance = (*nodeIter)->getMinDist();
- ptrMinNode = nodeIter; // ptrMinNode will hold a ptr to node with sm'st dist
- }
- }
- return ptrMinNode;
-}
-
-/**
- * This function is used to recompute the nodeToJoin1 dist mat row. Each row in the distance
- * matrix has the distances to all nodes before it, not after. For example the dist mat row
- * for Node 3 would have the distances to nodes 1 and 2.
- * 1 0
- * 2 0 0
- * 3 X X 0 Dists to node 1 and 2.
- * 4 0 0 0 0
- * 5 0 0 0 0 0
- * 6 0 0 0 0 0 0
- * 7 0 0 0 0 0 0 0
- */
-void UPGMAAlgorithm::recomputeNodeToJoin1DistMatRow(Node* nodeToJoin1,
- double** nodeToJoin2DistIter)
-{
- double* nodeToJoin1DistIter = nodeToJoin1->ptrToDistMatRow;
- // calculate new distance
- *nodeToJoin1DistIter = calcNewDist(*nodeToJoin1DistIter, **nodeToJoin2DistIter);
-
- const double* minIndex1 = nodeToJoin1DistIter;
- nodeToJoin1DistIter++;
- (*nodeToJoin2DistIter)++;
- int numDistToUpdate = nodeToJoin1->numDists - 1;
-
- // For each of the distances in nodeToJoin1
- while(numDistToUpdate > 0)
- {
- if (*nodeToJoin1DistIter >= 0)
- {
- // Calculate the average
- *nodeToJoin1DistIter = calcNewDist(*nodeToJoin1DistIter, **nodeToJoin2DistIter);
-
- if (*nodeToJoin1DistIter < *minIndex1)
- {
- minIndex1 = nodeToJoin1DistIter;
- }
- }
- nodeToJoin1DistIter++;
- (*nodeToJoin2DistIter)++;
- numDistToUpdate--;
- }
-
- // We have found the minimum distance
- nodeToJoin1->minDist = *minIndex1;
- nodeToJoin1->indexToMinDist = minIndex1 - nodeToJoin1->ptrToDistMatRow;
-}
-
-/**
- * This function is used to recompute all the other distances. It does this by calling
- * two other functions. We first compute all the distances until we get to node 2.
- * Then we call another function to do the rest. This is because nodes after node 2 may
- * have had EITHER node1 of node2 as their min distance.
- */
-void UPGMAAlgorithm::computeAllOtherDistsToNewNode(Node* nodeToJoin1, Node* nodeToJoin2,
- double** nodeToJoin2DistIter)
-{
- computeDistsUpToNodeToJoin2(nodeToJoin1, nodeToJoin2, nodeToJoin2DistIter);
- computeDistsForNodesAfterNode2(nodeToJoin2);
-}
-
-/**
- * STEP 3A:
- * This function recomputes the distances in column until we get to nodeToJoin2
- * example: Join nodes 2 and 6
- *
- * 1 0
- * 2 C 0
- * 3 0 X 0
- * 4 0 X 0 0
- * 5 0 X 0 0 0
- * 6 0 0 0 0 0 0
- * 7 0 0 0 0 0 0 0
- *
- * For each node until we get to nodeToJoin2
- * If newdistance is less than the old min distance, set it to this.
- * else if its greater than old minDist and indexTominDist is the same, recompute min
- * else leave the min the same as it hasnt changed.
- */
-void UPGMAAlgorithm::computeDistsUpToNodeToJoin2(Node* nodeToJoin1, Node* nodeToJoin2, double** nodeToJoin2DistIter)
-{
- const int indexToMinDist = nodeToJoin2->indexToMinDist;
-
- movePtrPastUnusedDistances(nodeToJoin2DistIter);
-
- Node* nodeIter;
- // For each node until we get to the second node we are joining
- for(nodeIter = nodeToJoin1->next; nodeIter != nodeToJoin2; nodeIter = nodeIter->next)
- {
- (*nodeToJoin2DistIter)++; // Skip the distance to the node we are joining with
- movePtrPastUnusedDistances(nodeToJoin2DistIter);
-
- double distToNode = nodeIter->ptrToDistMatRow[indexToMinDist];
-
- double newDistToNode = calcNewDist(distToNode, **nodeToJoin2DistIter);
- nodeIter->ptrToDistMatRow[indexToMinDist] = newDistToNode;
-
-
- if (newDistToNode < nodeIter->minDist)
- { /** new value is smaller than current min. */
- nodeIter->minDist = newDistToNode;
- nodeIter->indexToMinDist = indexToMinDist;
- }
- else if ((newDistToNode > nodeIter->minDist) &&
- (nodeIter->indexToMinDist == indexToMinDist))
- { /** indexToMinDist was the min dist, but it is now a larger num, recompute */
- nodeIter->findMinDist();
- }
-
- }
-
-}
-
-/**
- * STEP 3B Recompute distance for nodes after nodeToJoin2
- * example: Join nodes 2 and 6
- *
- * 1 0
- * 2 C 0
- * 3 0 C 0
- * 4 0 C 0 0
- * 5 0 C 0 0 0
- * 6 0 0 0 0 0 0
- * 7 0 X 0 0 0 0 0
- *
- * For each node until we get to the end.
- * if dist is less than minDist, set mindist to new distance, set index to index
- * else if dist is greater than mindist and the index was either nodetojoin1 or
- * nodetojoin2, recompute distance, set entry for nodetojoin2 to NOTUSED
- * else set the distance to unused.
- */
-
-void UPGMAAlgorithm::computeDistsForNodesAfterNode2(Node* nodeToJoin2)
-{
- int indexToNode2 = nodeToJoin2->numDists;
- const int indexToMinDist = nodeToJoin2->indexToMinDist;
-
- Node* nodeIter;
-
- for(nodeIter = nodeToJoin2->next; nodeIter; nodeIter = nodeIter->next)
- {
- double &distUpdate = nodeIter->ptrToDistMatRow[indexToMinDist];
-
- distUpdate =
- ((distUpdate * orderNode1) +
- (nodeIter->ptrToDistMatRow[indexToNode2] * orderNode2))
- / orderNewNode;
-
- /* The comparison (distUpdate > nodeIter->minDist) is unsafe.
- * Specifically, we get different results for 32/64bit machines,
- * which leads to different branching in the if/else statement
- * and nasty behaviour down the line.
- * Using all of Pfam as a benchmark distUpdate can 'wobble' by 40 ULPs
- * (Unit of Least Precision) which is difficult to translate into
- * a maximum relative error, so we pick COMPARE_REL_EPSILON
- * phenomenologically: approx 2E-15 was the biggest we saw in Pfam,
- * but suggest 1E-14 (for good measure).
- * Using ULPs, eg, *(long long int*)&(distUpdate),
- * would be better (more elegant?) but gives problems
- * with aliasing and/or performance reduction.
- * FS, 2009-04-30
- */
-#define COMPARE_REL_EPSILON 1E-14
- if ( (distUpdate < nodeIter->minDist) &&
- ((nodeIter->minDist-distUpdate) /
- nodeIter->minDist > COMPARE_REL_EPSILON) )
- { /** new distance is smaller */
- nodeIter->minDist = distUpdate;
- nodeIter->indexToMinDist = indexToMinDist;
- }
- else if (((distUpdate > nodeIter->minDist) &&
- ((distUpdate-nodeIter->minDist) /
- distUpdate > COMPARE_REL_EPSILON) &&
- (nodeIter->indexToMinDist == indexToMinDist)) ||
- (nodeIter->indexToMinDist == indexToNode2))
- { /** if old min dist was to either nodeToJoin1 or nodeToJoin2 */
- nodeIter->ptrToDistMatRow[indexToNode2] = BLANKDIST;
- nodeIter->findMinDist();
- }
- else
- {
- nodeIter->ptrToDistMatRow[indexToNode2] = BLANKDIST;
- }
- }
-
-}
-
-void UPGMAAlgorithm::addAlignmentStep(vector<int>* group1, vector<int>* group2)
-{
- int sizeGroup1 = group1->size();
- int sizeGroup2 = group2->size();
-
- vector<int> groups;
- groups.resize(numSeqs + 1, 0);
- int sizeGroup = groups.size();
-
- for(int i = 0; i < sizeGroup1 && (*group1)[i] < sizeGroup; i++)
- {
- groups[(*group1)[i]] = 1; // Set to be apart of group1
- }
-
- for(int i = 0; i < sizeGroup2 && (*group2)[i] < sizeGroup; i++)
- {
- groups[(*group2)[i]] = 2; // Set to be apart of group1
- }
-
- progSteps->saveSet(&groups);
-}
-
-/**
- * At the moment we are only using average distance, UPGMA, but we can also use this
- * function to have different distance measures, min, max etc.
- */
-double UPGMAAlgorithm::calcNewDist(double dist1, double dist2)
-{
- return ((dist1 * orderNode1) + (dist2 * orderNode2)) / orderNewNode;
-}
-
-}
git://source.jalview.org / jabaws.git / blobdiff