rarefact.cpp

   1 /*
   2  *  rarefact.cpp
   3  *  Dotur
   4  *
   5  *  Created by Sarah Westcott on 11/18/08.
   6  *  Copyright 2008 __MyCompanyName__. 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
  57                 for(int i=0;i<displays.size();i++){
  58                         displays[i]->close();
  59                 }
  60         }
  61         catch(exception& e) {
  62                 cout << "Standard Error: " << e.what() << " has occurred in the Rarefact class Function getCurve. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  63                 exit(1);
  64         }
  65         catch(...) {
  66                 cout << "An unknown error has occurred in the Rarefact class function getCurve. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  67                 exit(1);
  68         }
  69 }
  70
  71 /***********************************************************************/
  72
  73 void Rarefact::getSharedCurve(int increment = 1, int nIters = 1000){
  74 try {
  75                 globaldata = GlobalData::getInstance();
  76                 SharedRarefactionCurveData* rcd = new SharedRarefactionCurveData();
  77
  78                 //register the displays
  79                 for(int i=0;i<displays.size();i++){
  80                         rcd->registerDisplay(displays[i]);
  81                 }
  82
  83                 for(int iter=0;iter<nIters;iter++){
  84                         //clear out old values for new iter
  85                         lookup.clear();
  86
  87                         //create and initialize vector of sharedvectors, one for each group
  88                         for (int i = 0; i < globaldata->gGroupmap->getNumGroups(); i++) {
  89                                 SharedRAbundVector* temp = new SharedRAbundVector(sharedorder->getNumBins());
  90                                 temp->setLabel(sharedorder->getLabel());
  91                                 temp->setGroup(globaldata->gGroupmap->namesOfGroups[i]);
  92                                 lookup.push_back(temp);
  93                         }
  94
  95                         for(int i=0;i<displays.size();i++){
  96                                 displays[i]->init(label);
  97                         }
  98
  99                         //sample all the members
 100                         for(int i=0;i<numSeqs;i++){
 101                                 //get first sample
 102                                 individual chosen = sharedorder->get(i);
 103                                 int abundance;
 104
 105                                 //set info for sharedvector in chosens group
 106                                 for (int j = 0; j < lookup.size(); j++) {
 107                                         if (chosen.group == lookup[j]->getGroup()) {
 108                                                 abundance = lookup[j]->getAbundance(chosen.bin);
 109                                                 lookup[j]->set(chosen.bin, (abundance + 1), chosen.group);
 110                                                 break;
 111                                         }
 112                                 }
 113                         }
 114
 115                         //randomize the groups
 116                         random_shuffle(lookup.begin(), lookup.end());
 117
 118                         //send the first group
 119                         rcd->updateSharedData(lookup[0], lookup[0], 1, numGroupComb);
 120
 121                         //send each additional group one at a time
 122                         int n = 1;
 123                         for (int k = 0; k < (lookup.size() - 1); k++) {
 124                                 for (int l = n; l < lookup.size(); l++) {
 125                                         rcd->updateSharedData(lookup[k], lookup[l], l+1, numGroupComb);
 126                                         mergeVectors(lookup[0], lookup[l]);
 127                                 }
 128                                 n++;
 129                         }
 130
 131                         //resets output files
 132                         for(int i=0;i<displays.size();i++){
 133                                 displays[i]->reset();
 134                         }
 135                 }
 136
 137                 for(int i=0;i<displays.size();i++){
 138                         displays[i]->close();
 139                 }
 140
 141         }
 142         catch(exception& e) {
 143                 cout << "Standard Error: " << e.what() << " has occurred in the Rarefact class Function getSharedCurve. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 144                 exit(1);
 145         }
 146         catch(...) {
 147                 cout << "An unknown error has occurred in the Rarefact class function getSharedCurve. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 148                 exit(1);
 149         }
 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                         if ((shared1->getAbundance(k) == 0) && (shared2->getAbundance(k) != 0)) {
 159                                 shared1->set(k, shared2->getAbundance(k), "combo");  //set to 'combo' since this vector now contains multiple groups
 160                         }
 161                 }
 162         }
 163         catch(exception& e) {
 164                 cout << "Standard Error: " << e.what() << " has occurred in the Rarefact class Function mergeVectors. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 165                 exit(1);
 166         }
 167         catch(...) {
 168                 cout << "An unknown error has occurred in the Rarefact class function mergeVectors. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 169                 exit(1);
 170         }
 171 }
 172