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