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