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 void Rarefact::getCurve(int increment = 1, 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                 for(int iter=0;iter<nIters;iter++){
  23
  24                         for(int i=0;i<displays.size();i++){
  25                                 displays[i]->init(label);
  26                         }
  27
  28                         RAbundVector* lookup    = new RAbundVector(order->getNumBins());
  29                         SAbundVector* rank      = new SAbundVector(order->getMaxRank()+1);
  30                         random_shuffle(order->begin(), order->end());
  31
  32                         for(int i=0;i<numSeqs;i++){
  33
  34                                 int binNumber = order->get(i);
  35                                 int abundance = lookup->get(binNumber);
  36
  37                                 rank->set(abundance, rank->get(abundance)-1);
  38                                 abundance++;
  39
  40                                 lookup->set(binNumber, abundance);
  41                                 rank->set(abundance, rank->get(abundance)+1);
  42
  43                                 if((i == 0) || (i+1) % increment == 0){
  44                                         rcd->updateRankData(rank);
  45                                 }
  46                         }
  47
  48                         if(numSeqs % increment != 0){
  49                                 rcd->updateRankData(rank);
  50                         }
  51
  52                         for(int i=0;i<displays.size();i++){
  53                                 displays[i]->reset();
  54                         }
  55
  56                         delete lookup;
  57                         delete rank;
  58                 }
  59
  60                 for(int i=0;i<displays.size();i++){
  61                         displays[i]->close();
  62                 }
  63                 delete rcd;
  64         }
  65         catch(exception& e) {
  66                 cout << "Standard Error: " << e.what() << " has occurred in the Rarefact class Function getCurve. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  67                 exit(1);
  68         }
  69         catch(...) {
  70                 cout << "An unknown error has occurred in the Rarefact class function getCurve. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  71                 exit(1);
  72         }
  73 }
  74
  75 /***********************************************************************/
  76
  77 void Rarefact::getSharedCurve(int increment = 1, int nIters = 1000){
  78 try {
  79                 SharedRarefactionCurveData* rcd = new SharedRarefactionCurveData();
  80
  81                 label = lookup[0]->getLabel();
  82
  83                 //register the displays
  84                 for(int i=0;i<displays.size();i++){
  85                         rcd->registerDisplay(displays[i]);
  86                 }
  87
  88                 for(int iter=0;iter<nIters;iter++){
  89
  90                         for(int i=0;i<displays.size();i++){
  91                                 displays[i]->init(label);
  92                         }
  93
  94                         //randomize the groups
  95                         random_shuffle(lookup.begin(), lookup.end());
  96
  97                         //make merge the size of lookup[0]
  98                         SharedRAbundVector* merge = new SharedRAbundVector(lookup[0]->size());
  99
 100                         //make copy of lookup zero
 101                         for(int i = 0; i<lookup[0]->size(); i++) {
 102                                 merge->set(i, lookup[0]->getAbundance(i), "merge");
 103                         }
 104
 105                         vector<SharedRAbundVector*> subset;
 106                         //send each group one at a time
 107                         for (int k = 0; k < lookup.size(); k++) {
 108                                 subset.clear(); //clears out old pair of sharedrabunds
 109                                 //add in new pair of sharedrabunds
 110                                 subset.push_back(merge); subset.push_back(lookup[k]);
 111
 112                                 rcd->updateSharedData(subset, k+1, numGroupComb);
 113                                 mergeVectors(merge, lookup[k]);
 114                         }
 115
 116                         //resets output files
 117                         for(int i=0;i<displays.size();i++){
 118                                 displays[i]->reset();
 119                         }
 120
 121                         delete merge;
 122                 }
 123
 124                 for(int i=0;i<displays.size();i++){
 125                         displays[i]->close();
 126                 }
 127
 128                 delete rcd;
 129         }
 130         catch(exception& e) {
 131                 cout << "Standard Error: " << e.what() << " has occurred in the Rarefact class Function getSharedCurve. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 132                 exit(1);
 133         }
 134         catch(...) {
 135                 cout << "An unknown error has occurred in the Rarefact class function getSharedCurve. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 136                 exit(1);
 137         }
 138
 139 }
 140
 141 /**************************************************************************************/
 142 void Rarefact::mergeVectors(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
 143         try{
 144                 for (int k = 0; k < shared1->size(); k++) {
 145                         //merge new species into shared1
 146                         shared1->set(k, (shared1->getAbundance(k) + shared2->getAbundance(k)), "combo");  //set to 'combo' since this vector now contains multiple groups
 147                 }
 148         }
 149         catch(exception& e) {
 150                 cout << "Standard Error: " << e.what() << " has occurred in the Rarefact class Function mergeVectors. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 151                 exit(1);
 152         }
 153         catch(...) {
 154                 cout << "An unknown error has occurred in the Rarefact class function mergeVectors. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 155                 exit(1);
 156         }
 157 }
 158