weighted.cpp

   1 /*
   2  *  weighted.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 2/9/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "weighted.h"
  11
  12 /**************************************************************************************************/
  13
  14 EstOutput Weighted::getValues(Tree* t) {
  15     try {
  16                 globaldata = GlobalData::getInstance();
  17                 int numGroups;
  18                 vector<double> D;
  19
  20                 numGroups = globaldata->Groups.size();
  21
  22                 //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
  23                 int count = 0;
  24                 for (int i=0; i<numGroups; i++) {
  25                         for (int l = i+1; l < numGroups; l++) {
  26                                 //initialize weighted scores
  27                                 WScore[globaldata->Groups[i]+globaldata->Groups[l]] = 0.0;
  28
  29                                 D.push_back(0.0000); //initialize a spot in D for each combination
  30
  31                                 /********************************************************/
  32                                 //calculate a D value for each group combo
  33                                 for(int v=0;v<t->getNumLeaves();v++){
  34                                         int index = v;
  35                                         double sum = 0.0000;
  36
  37                                         //while you aren't at root
  38                                         while(t->tree[index].getParent() != -1){
  39
  40                                                 //if you have a BL
  41                                                 if(t->tree[index].getBranchLength() != -1){
  42                                                         sum += abs(t->tree[index].getBranchLength());
  43                                                 }
  44                                                 index = t->tree[index].getParent();
  45                                         }
  46
  47                                         //get last breanch length added
  48                                         if(t->tree[index].getBranchLength() != -1){
  49                                                 sum += abs(t->tree[index].getBranchLength());
  50                                         }
  51
  52                                         //is this sum from a sequence which is in one of the users groups
  53                                         if (inUsersGroups(t->tree[v].getGroup(), globaldata->Groups) == true) {
  54                                                 //is this sum from a sequence which is in this groupCombo
  55                                                 if ((t->tree[v].getGroup() == globaldata->Groups[i]) || (t->tree[v].getGroup() == globaldata->Groups[l])) {
  56                                                         sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
  57                                                         D[count] += sum;
  58                                                 }
  59                                         }
  60                                 }
  61                                 /*********************************************************/
  62                                 count++;
  63                         }
  64                 }
  65
  66                 data.clear(); //clear out old values
  67
  68                 for(int i=0;i<t->getNumNodes();i++){
  69                         //calculate weighted score for each of the group comb i.e. with groups A,B,C = AB, AC, BC.
  70                         for (int b=0; b<numGroups; b++) {
  71                                 for (int l = b+1; l < numGroups; l++) {
  72                                         //calculate a u value for each combo
  73                                         double u;
  74                                         //does this node have descendants from group b-1
  75                                         it = t->tree[i].pcount.find(globaldata->Groups[b]);
  76                                         //if it does u = # of its descendants with a certain group / total number in tree with a certain group
  77                                         if (it != t->tree[i].pcount.end()) {
  78                                                 u = (double) t->tree[i].pcount[globaldata->Groups[b]] / (double) tmap->seqsPerGroup[globaldata->Groups[b]];
  79                                         }else { u = 0.00; }
  80
  81                                         //does this node have descendants from group l
  82                                         it = t->tree[i].pcount.find(globaldata->Groups[l]);
  83                                         //if it does subtract their percentage from u
  84                                         if (it != t->tree[i].pcount.end()) {
  85                                                 u -= (double) t->tree[i].pcount[globaldata->Groups[l]] / (double) tmap->seqsPerGroup[globaldata->Groups[l]];
  86                                         }
  87
  88                                         u = abs(u * t->tree[i].getBranchLength());
  89
  90                                         //save groupcombs u value
  91                                         WScore[globaldata->Groups[b]+globaldata->Groups[l]] += u;
  92                                 /*********************************************************/
  93                                 }
  94                         }
  95                 }
  96
  97                 //calculate weighted score for each group combination
  98                 double UN;
  99                 count = 0;
 100                 for (int i=0; i<numGroups; i++) {
 101                         for (int l = i+1; l < numGroups; l++) {
 102                                 UN = (WScore[globaldata->Groups[i]+globaldata->Groups[l]] / D[count]);
 103
 104                                 if (isnan(UN) || isinf(UN)) { UN = 0; }
 105                                 data.push_back(UN);
 106                                 count++;
 107                         }
 108                 }
 109                 return data;
 110         }
 111         catch(exception& e) {
 112                 errorOut(e, "Weighted", "getValues");
 113                 exit(1);
 114         }
 115 }
 116
 117 /**************************************************************************************************/
 118 EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) {
 119  try {
 120                 globaldata = GlobalData::getInstance();
 121
 122                 data.clear(); //clear out old values
 123
 124                 //initialize weighted score
 125                 WScore[(groupA+groupB)] = 0.0;
 126                 float D = 0.0;
 127
 128
 129                 /********************************************************/
 130                 //calculate a D value for the group combo
 131                 for(int v=0;v<t->getNumLeaves();v++){
 132                         int index = v;
 133                         double sum = 0.0000;
 134
 135                         //while you aren't at root
 136                         while(t->tree[index].getParent() != -1){
 137
 138                                 //if you have a BL
 139                                 if(t->tree[index].getBranchLength() != -1){
 140                                         sum += abs(t->tree[index].getBranchLength());
 141                                 }
 142                                 index = t->tree[index].getParent();
 143                         }
 144
 145                         //get last breanch length added
 146                         if(t->tree[index].getBranchLength() != -1){
 147                                 sum += abs(t->tree[index].getBranchLength());
 148                         }
 149
 150                         if ((t->tree[v].getGroup() == groupA) || (t->tree[v].getGroup() == groupB)) {
 151                                 sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
 152                                 D += sum;
 153                         }
 154                 }
 155                 /********************************************************/
 156
 157                 //calculate u for the group comb
 158                 for(int i=0;i<t->getNumNodes();i++){
 159                         double u;
 160                         //does this node have descendants from groupA
 161                         it = t->tree[i].pcount.find(groupA);
 162                         //if it does u = # of its descendants with a certain group / total number in tree with a certain group
 163                         if (it != t->tree[i].pcount.end()) {
 164                                 u = (double) t->tree[i].pcount[groupA] / (double) tmap->seqsPerGroup[groupA];
 165                         }else { u = 0.00; }
 166
 167
 168                         //does this node have descendants from group l
 169                         it = t->tree[i].pcount.find(groupB);
 170                         //if it does subtract their percentage from u
 171                         if (it != t->tree[i].pcount.end()) {
 172                                 u -= (double) t->tree[i].pcount[groupB] / (double) tmap->seqsPerGroup[groupB];
 173                         }
 174
 175                         u = abs(u * t->tree[i].getBranchLength());
 176
 177                         //save groupcombs u value
 178                         WScore[(groupA+groupB)] += u;
 179                 }
 180
 181                 /********************************************************/
 182
 183                 //calculate weighted score for the group combination
 184                 double UN;
 185                 UN = (WScore[(groupA+groupB)] / D);
 186
 187                 if (isnan(UN) || isinf(UN)) { UN = 0; }
 188                 data.push_back(UN);
 189
 190                 return data;
 191         }
 192         catch(exception& e) {
 193                 errorOut(e, "Weighted", "getValues");
 194                 exit(1);
 195         }
 196 }
 197
 198
 199
 200
 201
 202
 203
 204
 205