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                 //if the user has not entered specific groups to analyze then do them all
  21                 if (globaldata->Groups.size() == 0) {
  22                         numGroups = tmap->getNumGroups();
  23                 }else {
  24                         numGroups = globaldata->Groups.size();
  25                 }
  26
  27                 //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
  28                 int n = 1;
  29                 int count = 0;
  30                 for (int i=1; i<numGroups; i++) {
  31                         for (int l = n; l < numGroups; l++) {
  32                                 //initialize weighted scores
  33                                 if (globaldata->Groups.size() == 0) {
  34                                         WScore[tmap->namesOfGroups[i-1]+tmap->namesOfGroups[l]] = 0.0;
  35                                 }else {
  36                                         WScore[globaldata->Groups[i-1]+globaldata->Groups[l]] = 0.0;
  37                                 }
  38
  39                                 D.push_back(0.0000); //initialize a spot in D for each combination
  40
  41                                 /********************************************************/
  42                                 //calculate a D value for each group combo
  43                                 for(int v=0;v<t->getNumLeaves();v++){
  44                                         int index = v;
  45                                         double sum = 0.0000;
  46
  47                                         //while you aren't at root
  48                                         while(t->tree[index].getParent() != -1){
  49
  50                                                 //if you have a BL
  51                                                 if(t->tree[index].getBranchLength() != -1){
  52                                                         sum += t->tree[index].getBranchLength();
  53                                                 }
  54                                                 index = t->tree[index].getParent();
  55                                         }
  56
  57                                         //get last breanch length added
  58                                         if(t->tree[index].getBranchLength() != -1){
  59                                                 sum += t->tree[index].getBranchLength();
  60                                         }
  61
  62                                         if (globaldata->Groups.size() == 0) {
  63                                                 //is this sum from a sequence which is in one of the users groups
  64                                                 if (inUsersGroups(t->tree[v].getGroup(), tmap->namesOfGroups) == true) {
  65                                                         //is this sum from a sequence which is in this groupCombo
  66                                                         if ((t->tree[v].getGroup() == tmap->namesOfGroups[i-1]) || (t->tree[v].getGroup() == tmap->namesOfGroups[l])) {
  67                                                                 sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
  68                                                                 D[count] += sum;
  69                                                         }
  70                                                 }
  71                                         }else {
  72                                                 //is this sum from a sequence which is in one of the users groups
  73                                                 if (inUsersGroups(t->tree[v].getGroup(), globaldata->Groups) == true) {
  74                                                         //is this sum from a sequence which is in this groupCombo
  75                                                         if ((t->tree[v].getGroup() == globaldata->Groups[i-1]) || (t->tree[v].getGroup() == globaldata->Groups[l])) {
  76                                                                 sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
  77                                                                 D[count] += sum;
  78                                                         }
  79                                                 }
  80                                         }
  81                                 }
  82                                 /*********************************************************/
  83                                 count++;
  84                         }
  85                         n++;
  86                 }
  87
  88                 data.clear(); //clear out old values
  89
  90                 for(int i=0;i<t->getNumNodes();i++){
  91                         //calculate weighted score for each of the group comb i.e. with groups A,B,C = AB, AC, BC.
  92                         n = 1;
  93                         for (int b=1; b<numGroups; b++) {
  94                                 for (int l = n; l < numGroups; l++) {
  95                                         //calculate a u value for each combo
  96                                         double u;
  97                                         //the user has not entered specific groups
  98                                         if (globaldata->Groups.size() == 0) {
  99                                                 //does this node have descendants from group b-1
 100                                                 it = t->tree[i].pcount.find(tmap->namesOfGroups[b-1]);
 101                                                 //if it does u = # of its descendants with a certain group / total number in tree with a certain group
 102                                                 if (it != t->tree[i].pcount.end()) {
 103                                                         u = (double) t->tree[i].pcount[tmap->namesOfGroups[b-1]] / (double) tmap->seqsPerGroup[tmap->namesOfGroups[b-1]];
 104                                                 }else { u = 0.00; }
 105
 106                                                 //does this node have descendants from group l
 107                                                 it = t->tree[i].pcount.find(tmap->namesOfGroups[l]);
 108                                                 //if it does subtract their percentage from u
 109                                                 if (it != t->tree[i].pcount.end()) {
 110                                                         u -= (double) t->tree[i].pcount[tmap->namesOfGroups[l]] / (double) tmap->seqsPerGroup[tmap->namesOfGroups[l]];
 111                                                 }
 112
 113                                                 u = abs(u) * t->tree[i].getBranchLength();
 114
 115                                                 //save groupcombs u value
 116                                                 WScore[tmap->namesOfGroups[b-1]+tmap->namesOfGroups[l]] += u;
 117
 118                                         //the user has entered specific groups
 119                                         }else {
 120                                                 //does this node have descendants from group b-1
 121                                                 it = t->tree[i].pcount.find(globaldata->Groups[b-1]);
 122                                                 //if it does u = # of its descendants with a certain group / total number in tree with a certain group
 123                                                 if (it != t->tree[i].pcount.end()) {
 124                                                         u = (double) t->tree[i].pcount[globaldata->Groups[b-1]] / (double) tmap->seqsPerGroup[globaldata->Groups[b-1]];
 125                                                 }else { u = 0.00; }
 126
 127                                                 //does this node have descendants from group l
 128                                                 it = t->tree[i].pcount.find(globaldata->Groups[l]);
 129                                                 //if it does subtract their percentage from u
 130                                                 if (it != t->tree[i].pcount.end()) {
 131                                                         u -= (double) t->tree[i].pcount[globaldata->Groups[l]] / (double) tmap->seqsPerGroup[globaldata->Groups[l]];
 132                                                 }
 133
 134                                                 u = abs(u) * t->tree[i].getBranchLength();
 135
 136                                                 //save groupcombs u value
 137                                                 WScore[globaldata->Groups[b-1]+globaldata->Groups[l]] += u;
 138                                         }
 139                                 /*********************************************************/
 140                                 }
 141                                 n++;
 142                         }
 143                 }
 144
 145                 //calculate weighted score for each group combination
 146                 double UN;
 147                 n = 1;
 148                 count = 0;
 149                 for (int i=1; i<numGroups; i++) {
 150                         for (int l = n; l < numGroups; l++) {
 151                                 //the user has not entered specific groups
 152                                 if (globaldata->Groups.size() == 0) {
 153                                         UN = (WScore[tmap->namesOfGroups[i-1]+tmap->namesOfGroups[l]] / D[count]);
 154                                 }else {//they have entered specific groups
 155                                         UN = (WScore[globaldata->Groups[i-1]+globaldata->Groups[l]] / D[count]);
 156                                 }
 157                                 if (isnan(UN) || isinf(UN)) { UN = 0; }
 158                                 data.push_back(UN);
 159                                 count++;
 160                         }
 161                         n++;
 162                 }
 163                 return data;
 164         }
 165         catch(exception& e) {
 166                 cout << "Standard Error: " << e.what() << " has occurred in the Weighted class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 167                 exit(1);
 168         }
 169         catch(...) {
 170                 cout << "An unknown error has occurred in the Weighted class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 171                 exit(1);
 172         }
 173
 174 }
 175
 176 /**************************************************************************************************/