rarefact.cpp

   1 /*
   2  *  rarefact.cpp
   3  *  Dotur
   4  *
   5  *  Created by Sarah Westcott on 11/18/08.
   6  *  Copyright 2008 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "rarefact.h"
  11 //#include "ordervector.hpp"
  12
  13 /***********************************************************************/
  14
  15 int Rarefact::getCurve(float percentFreq = 0.01, int nIters = 1000){
  16         try {
  17                 RarefactionCurveData* rcd = new RarefactionCurveData();
  18                 for(int i=0;i<displays.size();i++){
  19                         rcd->registerDisplay(displays[i]);
  20                 }
  21
  22                 //convert freq percentage to number
  23                 int increment = 1;
  24                 if (percentFreq < 1.0) {  increment = numSeqs * percentFreq;  }
  25                 else { increment = percentFreq;  }
  26
  27                 for(int iter=0;iter<nIters;iter++){
  28
  29                         for(int i=0;i<displays.size();i++){
  30                                 displays[i]->init(label);
  31                         }
  32
  33                         RAbundVector* lookup    = new RAbundVector(order->getNumBins());
  34                         SAbundVector* rank      = new SAbundVector(order->getMaxRank()+1);
  35                         random_shuffle(order->begin(), order->end());
  36
  37                         for(int i=0;i<numSeqs;i++){
  38
  39                                 if (m->control_pressed) { delete lookup; delete rank; delete rcd; return 0;  }
  40
  41                                 int binNumber = order->get(i);
  42                                 int abundance = lookup->get(binNumber);
  43
  44                                 rank->set(abundance, rank->get(abundance)-1);
  45                                 abundance++;
  46
  47                                 lookup->set(binNumber, abundance);
  48                                 rank->set(abundance, rank->get(abundance)+1);
  49
  50                                 if((i == 0) || (i+1) % increment == 0){
  51                                         rcd->updateRankData(rank);
  52                                 }
  53                         }
  54
  55                         if(numSeqs % increment != 0){
  56                                 rcd->updateRankData(rank);
  57                         }
  58
  59                         for(int i=0;i<displays.size();i++){
  60                                 displays[i]->reset();
  61                         }
  62
  63                         delete lookup;
  64                         delete rank;
  65                 }
  66
  67                 for(int i=0;i<displays.size();i++){
  68                         displays[i]->close();
  69                 }
  70                 delete rcd;
  71                 return 0;
  72         }
  73         catch(exception& e) {
  74                 m->errorOut(e, "Rarefact", "getCurve");
  75                 exit(1);
  76         }
  77 }
  78
  79 /***********************************************************************/
  80
  81 int Rarefact::getSharedCurve(float percentFreq = 0.01, int nIters = 1000){
  82 try {
  83                 SharedRarefactionCurveData* rcd = new SharedRarefactionCurveData();
  84
  85                 label = lookup[0]->getLabel();
  86
  87                 //register the displays
  88                 for(int i=0;i<displays.size();i++){
  89                         rcd->registerDisplay(displays[i]);
  90                 }
  91
  92                 //if jumble is false all iters will be the same
  93                 if (globaldata->jumble == false)  {  nIters = 1;  }
  94
  95                 //convert freq percentage to number
  96                 int increment = 1;
  97                 if (percentFreq < 1.0) {  increment = numSeqs * percentFreq;  }
  98                 else { increment = percentFreq;  }
  99
 100                 for(int iter=0;iter<nIters;iter++){
 101
 102                         for(int i=0;i<displays.size();i++){
 103                                 displays[i]->init(label);
 104                         }
 105
 106                         if (globaldata->jumble == true)  {
 107                                 //randomize the groups
 108                                 random_shuffle(lookup.begin(), lookup.end());
 109                         }
 110
 111                         //make merge the size of lookup[0]
 112                         SharedRAbundVector* merge = new SharedRAbundVector(lookup[0]->size());
 113
 114                         //make copy of lookup zero
 115                         for(int i = 0; i<lookup[0]->size(); i++) {
 116                                 merge->set(i, lookup[0]->getAbundance(i), "merge");
 117                         }
 118
 119                         vector<SharedRAbundVector*> subset;
 120                         //send each group one at a time
 121                         for (int k = 0; k < lookup.size(); k++) {
 122                                 if (m->control_pressed) {  delete merge; delete rcd; return 0;  }
 123
 124                                 subset.clear(); //clears out old pair of sharedrabunds
 125                                 //add in new pair of sharedrabunds
 126                                 subset.push_back(merge); subset.push_back(lookup[k]);
 127
 128                                 rcd->updateSharedData(subset, k+1, numGroupComb);
 129                                 mergeVectors(merge, lookup[k]);
 130                         }
 131
 132                         //resets output files
 133                         for(int i=0;i<displays.size();i++){
 134                                 displays[i]->reset();
 135                         }
 136
 137                         delete merge;
 138                 }
 139
 140                 for(int i=0;i<displays.size();i++){
 141                         displays[i]->close();
 142                 }
 143
 144                 delete rcd;
 145                 return 0;
 146         }
 147         catch(exception& e) {
 148                 m->errorOut(e, "Rarefact", "getSharedCurve");
 149                 exit(1);
 150         }
 151 }
 152
 153 /**************************************************************************************/
 154 void Rarefact::mergeVectors(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
 155         try{
 156                 for (int k = 0; k < shared1->size(); k++) {
 157                         //merge new species into shared1
 158                         shared1->set(k, (shared1->getAbundance(k) + shared2->getAbundance(k)), "combo");  //set to 'combo' since this vector now contains multiple groups
 159                 }
 160         }
 161         catch(exception& e) {
 162                 m->errorOut(e, "Rarefact", "mergeVectors");
 163                 exit(1);
 164         }
 165 }
 166