weighted.h

   1 #ifndef WEIGHTED_H
   2 #define WEIGHTED_H
   3
   4
   5 /*
   6  *  weighted.h
   7  *  Mothur
   8  *
   9  *  Created by Sarah Westcott on 2/9/09.
  10  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
  11  *
  12  */
  13
  14 #include "treecalculator.h"
  15 #include "counttable.h"
  16
  17 /***********************************************************************/
  18
  19 class Weighted : public TreeCalculator  {
  20
  21         public:
  22         Weighted( bool r) : includeRoot(r) {};
  23                 ~Weighted() {};
  24
  25                 EstOutput getValues(Tree*, string, string);
  26                 EstOutput getValues(Tree*, int, string);
  27
  28         private:
  29                 struct linePair {
  30                         int start;
  31                         int num;
  32                         linePair(int i, int j) : start(i), num(j) {}
  33                 };
  34                 vector<linePair> lines;
  35
  36                 EstOutput data;
  37                 map<string, int>::iterator it;
  38                 map<string, double> WScore; //a score for each group combination i.e. AB, AC, BC.
  39                 int processors;
  40                 string outputDir;
  41                 map< vector<string>, set<int> > rootForGrouping;  //maps a grouping combo to the root for that combo
  42                 bool includeRoot;
  43
  44                 EstOutput driver(Tree*, vector< vector<string> >, int, int, CountTable*);
  45                 EstOutput createProcesses(Tree*, vector< vector<string> >, CountTable*);
  46                 double getLengthToRoot(Tree*, int, string, string);
  47 };
  48
  49 /***********************************************************************/
  50 struct weightedData {
  51     int start;
  52         int num;
  53         MothurOut* m;
  54     EstOutput results;
  55     vector< vector<string> > namesOfGroupCombos;
  56     Tree* t;
  57     CountTable* ct;
  58     bool includeRoot;
  59
  60
  61         weightedData(){}
  62         weightedData(MothurOut* mout, int st, int en, vector< vector<string> > ngc, Tree* tree, CountTable* count, bool ir) {
  63         m = mout;
  64                 start = st;
  65                 num = en;
  66         namesOfGroupCombos = ngc;
  67         t = tree;
  68         ct = count;
  69         includeRoot = ir;
  70         }
  71 };
  72
  73 /**************************************************************************************************/
  74 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
  75 #else
  76 static DWORD WINAPI MyWeightedThreadFunction(LPVOID lpParam){
  77         weightedData* pDataArray;
  78         pDataArray = (weightedData*)lpParam;
  79         try {
  80         map<string, int>::iterator it;
  81                 vector<double> D;
  82                 int count = 0;
  83         map< vector<string>, set<int> > rootForGrouping;
  84         map<string, double> WScore;
  85
  86                 for (int h = pDataArray->start; h < (pDataArray->start+pDataArray->num); h++) {
  87
  88                         if (pDataArray->m->control_pressed) { return 0; }
  89
  90                         //initialize weighted score
  91                         string groupA = pDataArray->namesOfGroupCombos[h][0];
  92                         string groupB = pDataArray->namesOfGroupCombos[h][1];
  93
  94                         set<int> validBranches;
  95                         WScore[groupA+groupB] = 0.0;
  96                         D.push_back(0.0000); //initialize a spot in D for each combination
  97
  98                         //adding the wieghted sums from group i
  99                         for (int j = 0; j < pDataArray->t->groupNodeInfo[groupA].size(); j++) { //the leaf nodes that have seqs from group i
 100                                 map<string, int>::iterator it = pDataArray->t->tree[pDataArray->t->groupNodeInfo[groupA][j]].pcount.find(groupA);
 101                                 int numSeqsInGroupI = it->second;
 102
 103                                 //double sum = getLengthToRoot(pDataArray->t, pDataArray->t->groupNodeInfo[groupA][j], groupA, groupB);
 104                 /*************************************************************************************/
 105                 double sum = 0.0;
 106                 int index = pDataArray->t->groupNodeInfo[groupA][j];
 107
 108                 //you are a leaf
 109                 if(pDataArray->t->tree[index].getBranchLength() != -1){ sum += abs(pDataArray->t->tree[index].getBranchLength());       }
 110                 double tempTotal = 0.0;
 111                 index = pDataArray->t->tree[index].getParent();
 112
 113                 vector<string> grouping; grouping.push_back(groupA); grouping.push_back(groupB);
 114
 115                 rootForGrouping[grouping].insert(index);
 116
 117                 //while you aren't at root
 118                 while(pDataArray->t->tree[index].getParent() != -1){
 119
 120                     if (pDataArray->m->control_pressed) {  return 0; }
 121
 122                     int parent = pDataArray->t->tree[index].getParent();
 123
 124                     if (pDataArray->includeRoot) { //add everyone
 125                         if(pDataArray->t->tree[index].getBranchLength() != -1){ sum += abs(pDataArray->t->tree[index].getBranchLength());       }
 126                     }else {
 127
 128                         //am I the root for this grouping? if so I want to stop "early"
 129                         //does my sibling have descendants from the users groups?
 130                         int lc = pDataArray->t->tree[parent].getLChild();
 131                         int rc = pDataArray->t->tree[parent].getRChild();
 132
 133                         int sib = lc;
 134                         if (lc == index) { sib = rc; }
 135
 136                         map<string, int>::iterator itGroup;
 137                         int pcountSize = 0;
 138                         itGroup = pDataArray->t->tree[sib].pcount.find(groupA);
 139                         if (itGroup != pDataArray->t->tree[sib].pcount.end()) { pcountSize++;  }
 140                         itGroup = pDataArray->t->tree[sib].pcount.find(groupB);
 141                         if (itGroup != pDataArray->t->tree[sib].pcount.end()) { pcountSize++;  }
 142
 143                         //if yes, I am not the root so add me
 144                         if (pcountSize != 0) {
 145                             if (pDataArray->t->tree[index].getBranchLength() != -1) {
 146                                 sum += abs(pDataArray->t->tree[index].getBranchLength()) + tempTotal;
 147                                 tempTotal = 0.0;
 148                             }else {
 149                                 sum += tempTotal;
 150                                 tempTotal = 0.0;
 151                             }
 152                             rootForGrouping[grouping].clear();
 153                             rootForGrouping[grouping].insert(parent);
 154                         }else { //if no, I may be the root so add my br to tempTotal until I am proven innocent
 155                             if (pDataArray->t->tree[index].getBranchLength() != -1) {
 156                                 tempTotal += abs(pDataArray->t->tree[index].getBranchLength());
 157                             }
 158                         }
 159                     }
 160
 161                     index = parent;
 162                 }
 163
 164                 //get all nodes above the root to add so we don't add their u values above
 165                 index = *(rootForGrouping[grouping].begin());
 166
 167                 while(pDataArray->t->tree[index].getParent() != -1){
 168                     int parent = pDataArray->t->tree[index].getParent();
 169                     rootForGrouping[grouping].insert(parent);
 170                     index = parent;
 171                 }
 172
 173                 /*************************************************************************************/
 174                                 double weightedSum = ((numSeqsInGroupI * sum) / (double)pDataArray->ct->getGroupCount(groupA));
 175
 176                                 D[count] += weightedSum;
 177                         }
 178
 179                         //adding the wieghted sums from group l
 180                         for (int j = 0; j < pDataArray->t->groupNodeInfo[groupB].size(); j++) { //the leaf nodes that have seqs from group l
 181                                 map<string, int>::iterator it = pDataArray->t->tree[pDataArray->t->groupNodeInfo[groupB][j]].pcount.find(groupB);
 182                                 int numSeqsInGroupL = it->second;
 183
 184                                 //double sum = getLengthToRoot(pDataArray->t, pDataArray->t->groupNodeInfo[groupB][j], groupA, groupB);
 185                 /*************************************************************************************/
 186                 double sum = 0.0;
 187                 int index = pDataArray->t->groupNodeInfo[groupB][j];
 188
 189                 //you are a leaf
 190                 if(pDataArray->t->tree[index].getBranchLength() != -1){ sum += abs(pDataArray->t->tree[index].getBranchLength());       }
 191                 double tempTotal = 0.0;
 192                 index = pDataArray->t->tree[index].getParent();
 193
 194                 vector<string> grouping; grouping.push_back(groupA); grouping.push_back(groupB);
 195
 196                 rootForGrouping[grouping].insert(index);
 197
 198                 //while you aren't at root
 199                 while(pDataArray->t->tree[index].getParent() != -1){
 200
 201                     if (pDataArray->m->control_pressed) {  return 0; }
 202
 203                     int parent = pDataArray->t->tree[index].getParent();
 204
 205                     if (pDataArray->includeRoot) { //add everyone
 206                         if(pDataArray->t->tree[index].getBranchLength() != -1){ sum += abs(pDataArray->t->tree[index].getBranchLength());       }
 207                     }else {
 208
 209                         //am I the root for this grouping? if so I want to stop "early"
 210                         //does my sibling have descendants from the users groups?
 211                         int lc = pDataArray->t->tree[parent].getLChild();
 212                         int rc = pDataArray->t->tree[parent].getRChild();
 213
 214                         int sib = lc;
 215                         if (lc == index) { sib = rc; }
 216
 217                         map<string, int>::iterator itGroup;
 218                         int pcountSize = 0;
 219                         itGroup = pDataArray->t->tree[sib].pcount.find(groupA);
 220                         if (itGroup != pDataArray->t->tree[sib].pcount.end()) { pcountSize++;  }
 221                         itGroup = pDataArray->t->tree[sib].pcount.find(groupB);
 222                         if (itGroup != pDataArray->t->tree[sib].pcount.end()) { pcountSize++;  }
 223
 224                         //if yes, I am not the root so add me
 225                         if (pcountSize != 0) {
 226                             if (pDataArray->t->tree[index].getBranchLength() != -1) {
 227                                 sum += abs(pDataArray->t->tree[index].getBranchLength()) + tempTotal;
 228                                 tempTotal = 0.0;
 229                             }else {
 230                                 sum += tempTotal;
 231                                 tempTotal = 0.0;
 232                             }
 233                             rootForGrouping[grouping].clear();
 234                             rootForGrouping[grouping].insert(parent);
 235                         }else { //if no, I may be the root so add my br to tempTotal until I am proven innocent
 236                             if (pDataArray->t->tree[index].getBranchLength() != -1) {
 237                                 tempTotal += abs(pDataArray->t->tree[index].getBranchLength());
 238                             }
 239                         }
 240                     }
 241
 242                     index = parent;
 243                 }
 244
 245                 //get all nodes above the root to add so we don't add their u values above
 246                 index = *(rootForGrouping[grouping].begin());
 247
 248                 while(pDataArray->t->tree[index].getParent() != -1){
 249                     int parent = pDataArray->t->tree[index].getParent();
 250                     rootForGrouping[grouping].insert(parent);
 251                     index = parent;
 252                 }
 253
 254                 /*************************************************************************************/
 255
 256                                 double weightedSum = ((numSeqsInGroupL * sum) / (double)pDataArray->ct->getGroupCount(groupB));
 257
 258                                 D[count] += weightedSum;
 259                         }
 260                         count++;
 261                 }
 262
 263                 //calculate u for the group comb
 264                 for (int h = pDataArray->start; h < (pDataArray->start+pDataArray->num); h++) {
 265                         //report progress
 266                         //pDataArray->m->mothurOut("Processing combo: " + toString(h)); pDataArray->m->mothurOutEndLine();
 267
 268                         string groupA = pDataArray->namesOfGroupCombos[h][0];
 269                         string groupB = pDataArray->namesOfGroupCombos[h][1];
 270
 271                         //calculate u for the group comb
 272                         for(int i=0;i<pDataArray->t->getNumNodes();i++){
 273
 274                                 if (pDataArray->m->control_pressed) { return 0; }
 275
 276                                 double u;
 277                                 //int pcountSize = 0;
 278                                 //does this node have descendants from groupA
 279                                 it = pDataArray->t->tree[i].pcount.find(groupA);
 280                                 //if it does u = # of its descendants with a certain group / total number in tree with a certain group
 281                                 if (it != pDataArray->t->tree[i].pcount.end()) {
 282                                         u = (double) pDataArray->t->tree[i].pcount[groupA] / (double) pDataArray->ct->getGroupCount(groupA);
 283                                 }else { u = 0.00; }
 284
 285
 286                                 //does this node have descendants from group l
 287                                 it = pDataArray->t->tree[i].pcount.find(groupB);
 288
 289                                 //if it does subtract their percentage from u
 290                                 if (it != pDataArray->t->tree[i].pcount.end()) {
 291                                         u -= (double) pDataArray->t->tree[i].pcount[groupB] / (double) pDataArray->ct->getGroupCount(groupB);
 292                                 }
 293
 294                                 if (pDataArray->includeRoot) {
 295                                         if (pDataArray->t->tree[i].getBranchLength() != -1) {
 296                                                 u = abs(u * pDataArray->t->tree[i].getBranchLength());
 297                                                 WScore[(groupA+groupB)] += u;
 298                                         }
 299                                 }else {
 300                                         //if this is not the root then add it
 301                                         if (rootForGrouping[pDataArray->namesOfGroupCombos[h]].count(i) == 0) {
 302                                                 if (pDataArray->t->tree[i].getBranchLength() != -1) {
 303                                                         u = abs(u * pDataArray->t->tree[i].getBranchLength());
 304                                                         WScore[(groupA+groupB)] += u;
 305                                                 }
 306                                         }
 307                                 }
 308                         }
 309
 310                 }
 311
 312                 /********************************************************/
 313                 //calculate weighted score for the group combination
 314                 double UN;
 315                 count = 0;
 316                 for (int h = pDataArray->start; h < (pDataArray->start+pDataArray->num); h++) {
 317                         UN = (WScore[pDataArray->namesOfGroupCombos[h][0]+pDataArray->namesOfGroupCombos[h][1]] / D[count]);
 318                         if (isnan(UN) || isinf(UN)) { UN = 0; }
 319                         pDataArray->results.push_back(UN);
 320                         count++;
 321                 }
 322
 323                 return 0;
 324
 325     }
 326         catch(exception& e) {
 327                 pDataArray->m->errorOut(e, "Weighted", "MyWeightedThreadFunction");
 328                 exit(1);
 329         }
 330 }
 331 #endif
 332
 333
 334 #endif