binaries/src/clustalw/src/tree/UPGMA/RootedGuideTree.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 <iostream>
  10 #include "RootedGuideTree.h"
  11 #include "../../general/userparams.h"
  12 using namespace std;
  13
  14 namespace clustalw
  15 {
  16
  17 RootedGuideTree::RootedGuideTree()
  18  : root(0){ }
  19
  20 RootedGuideTree::RootedGuideTree(Node* r)
  21  : root(r) { }
  22
  23 RootedGuideTree::~RootedGuideTree()
  24 {
  25     root->makeEmpty();
  26 }
  27
  28 void RootedGuideTree::setRoot(Node* r)
  29 {
  30     makeEmpty();
  31
  32     root = r;
  33 }
  34
  35 void RootedGuideTree::makeEmpty()
  36 {
  37     root->makeEmpty();
  38 }
  39
  40 void RootedGuideTree::orderNodes()
  41 {
  42     calcOrderNode(root);
  43 }
  44
  45 int RootedGuideTree::calcOrderNode(Node* node)
  46 {
  47     if(node != 0)
  48     {
  49         if(node->getLeft() == 0 && node->getRight() == 0) // Leaf Node
  50         {
  51             node->setOrder(1);
  52             return 1;
  53         }
  54         else // Internal node
  55         {
  56             node->setOrder(calcOrderNode(node->getLeft()) + calcOrderNode(node->getRight()));
  57             return node->getOrder();
  58         }
  59     }
  60     return 0;
  61 }
  62
  63 void RootedGuideTree::calcSeqWeights(int firstSeq, int lastSeq, vector<int>* seqWeights)
  64 {
  65     if((int)seqWeights->size() < lastSeq - 1)
  66     {
  67         seqWeights->resize(lastSeq - 1);
  68     }
  69
  70     int i = 0, _nSeqs = 0;
  71     int temp = 0, sum = 0;
  72     //
  73     // If there are more than three sequences....
  74     //
  75     _nSeqs = lastSeq - firstSeq;
  76     if ((_nSeqs >= 2) && (userParameters->getDistanceTree() == true) &&
  77         (userParameters->getNoWeights() == false))
  78     {
  79         //
  80         // Calculate sequence weights based on Phylip tree.
  81         //
  82         orderNodes();
  83         calcWeights(seqWeights);
  84
  85         //
  86         // Normalise the weights, such that the sum of the weights = INT_SCALE_FACTOR
  87         //
  88
  89         sum = 0;
  90         for (i = firstSeq; i < lastSeq; i++)
  91         {
  92             sum += (*seqWeights)[i];
  93         }
  94
  95         if (sum == 0)
  96         {
  97             for (i = firstSeq; i < lastSeq; i++)
  98             {
  99                 (*seqWeights)[i] = 1;
 100             }
 101             sum = i;
 102         }
 103
 104         for (i = firstSeq; i < lastSeq; i++)
 105         {
 106             (*seqWeights)[i] = ((*seqWeights)[i] * INT_SCALE_FACTOR) / sum;
 107             if ((*seqWeights)[i] < 1)
 108             {
 109                 (*seqWeights)[i] = 1;
 110             }
 111         }
 112     }
 113     else
 114     {
 115         //
 116         // Otherwise, use identity weights.
 117         //
 118         temp = INT_SCALE_FACTOR / _nSeqs;
 119         // AW 2009-07-09: goes wrong if we have more than
 120         // INT_SCALE_FACTOR seqs. if so, set to 1, just as above
 121         // same as in Tree.cpp
 122         if (temp < 1)
 123             temp = 1;
 124
 125
 126         for (i = firstSeq; i < lastSeq; i++)
 127         {
 128             (*seqWeights)[i] = temp;
 129         }
 130     }
 131 }
 132
 133 void RootedGuideTree::calcWeights(vector<int>* seqWeights)
 134 {
 135     vector<float> weights;
 136     int sizeSeqWeights = seqWeights->size();
 137     weights.resize(sizeSeqWeights, 0.0);
 138
 139     doWeightCalc(0.0, &weights, root);
 140
 141     for(int i = 0; i < sizeSeqWeights; i++)
 142     {
 143         (*seqWeights)[i] = static_cast<int>(weights[i] * 100);
 144     }
 145 }
 146
 147 void RootedGuideTree::doWeightCalc(float weightSoFar, vector<float>* weights, Node* t)
 148 {
 149     if(t != 0)
 150     {
 151         if(t->getLeft() == 0 && t->getRight() == 0) // Leaf Node
 152         {
 153             (*weights)[t->getSeqNum() - 1] = weightSoFar;
 154         }
 155         else // Internal node
 156         {
 157             float w = weightSoFar + (t->getHeight() / t->getOrder());
 158             doWeightCalc(w, weights, t->getLeft());
 159             doWeightCalc(w, weights, t->getRight());
 160         }
 161     }
 162 }
 163
 164 }