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 += 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 += 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                 cout << "Standard Error: " << e.what() << " has occurred in the Weighted class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 113                 exit(1);
 114         }
 115         catch(...) {
 116                 cout << "An unknown error has occurred in the Weighted class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 117                 exit(1);
 118         }
 119
 120 }
 121
 122 /**************************************************************************************************/
 123 EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) {
 124  try {
 125                 globaldata = GlobalData::getInstance();
 126
 127                 data.clear(); //clear out old values
 128
 129                 //initialize weighted score
 130                 WScore[(groupA+groupB)] = 0.0;
 131                 float D = 0.0;
 132
 133
 134                 /********************************************************/
 135                 //calculate a D value for the group combo
 136                 for(int v=0;v<t->getNumLeaves();v++){
 137                         int index = v;
 138                         double sum = 0.0000;
 139
 140                         //while you aren't at root
 141                         while(t->tree[index].getParent() != -1){
 142
 143                                 //if you have a BL
 144                                 if(t->tree[index].getBranchLength() != -1){
 145                                         sum += t->tree[index].getBranchLength();
 146                                 }
 147                                 index = t->tree[index].getParent();
 148                         }
 149
 150                         //get last breanch length added
 151                         if(t->tree[index].getBranchLength() != -1){
 152                                 sum += t->tree[index].getBranchLength();
 153                         }
 154
 155                         if ((t->tree[v].getGroup() == groupA) || (t->tree[v].getGroup() == groupB)) {
 156                                 sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
 157                                 D += sum;
 158                         }
 159                 }
 160                 /********************************************************/
 161
 162                 //calculate u for the group comb
 163                 for(int i=0;i<t->getNumNodes();i++){
 164                         double u;
 165                         //does this node have descendants from groupA
 166                         it = t->tree[i].pcount.find(groupA);
 167                         //if it does u = # of its descendants with a certain group / total number in tree with a certain group
 168                         if (it != t->tree[i].pcount.end()) {
 169                                 u = (double) t->tree[i].pcount[groupA] / (double) tmap->seqsPerGroup[groupA];
 170                         }else { u = 0.00; }
 171
 172
 173                         //does this node have descendants from group l
 174                         it = t->tree[i].pcount.find(groupB);
 175                         //if it does subtract their percentage from u
 176                         if (it != t->tree[i].pcount.end()) {
 177                                 u -= (double) t->tree[i].pcount[groupB] / (double) tmap->seqsPerGroup[groupB];
 178                         }
 179
 180                         u = abs(u) * t->tree[i].getBranchLength();
 181
 182                         //save groupcombs u value
 183                         WScore[(groupA+groupB)] += u;
 184                 }
 185
 186                 /********************************************************/
 187
 188                 //calculate weighted score for the group combination
 189                 double UN;
 190                 UN = (WScore[(groupA+groupB)] / D);
 191
 192                 if (isnan(UN) || isinf(UN)) { UN = 0; }
 193                 data.push_back(UN);
 194
 195                 return data;
 196         }
 197         catch(exception& e) {
 198                 cout << "Standard Error: " << e.what() << " has occurred in the Weighted class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 199                 exit(1);
 200         }
 201         catch(...) {
 202                 cout << "An unknown error has occurred in the Weighted class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 203                 exit(1);
 204         }
 205
 206 }
 207
 208
 209
 210
 211
 212
 213
 214
 215