5 * Created by Sarah Westcott on 11/18/08.
6 * Copyright 2008 Schloss Lab UMASS Amherst. All rights reserved.
11 //#include "ordervector.hpp"
13 /***********************************************************************/
15 void Rarefact::getCurve(int increment = 1, int nIters = 1000){
17 RarefactionCurveData* rcd = new RarefactionCurveData();
18 for(int i=0;i<displays.size();i++){
19 rcd->registerDisplay(displays[i]);
22 for(int iter=0;iter<nIters;iter++){
24 for(int i=0;i<displays.size();i++){
25 displays[i]->init(label);
28 RAbundVector* lookup = new RAbundVector(order->getNumBins());
29 SAbundVector* rank = new SAbundVector(order->getMaxRank()+1);
30 random_shuffle(order->begin(), order->end());
32 for(int i=0;i<numSeqs;i++){
34 int binNumber = order->get(i);
35 int abundance = lookup->get(binNumber);
37 rank->set(abundance, rank->get(abundance)-1);
40 lookup->set(binNumber, abundance);
41 rank->set(abundance, rank->get(abundance)+1);
43 if((i == 0) || (i+1) % increment == 0){
44 rcd->updateRankData(rank);
48 if(numSeqs % increment != 0){
49 rcd->updateRankData(rank);
52 for(int i=0;i<displays.size();i++){
60 for(int i=0;i<displays.size();i++){
66 errorOut(e, "Rarefact", "getCurve");
71 /***********************************************************************/
73 void Rarefact::getSharedCurve(int increment = 1, int nIters = 1000){
75 SharedRarefactionCurveData* rcd = new SharedRarefactionCurveData();
77 label = lookup[0]->getLabel();
79 //register the displays
80 for(int i=0;i<displays.size();i++){
81 rcd->registerDisplay(displays[i]);
84 for(int iter=0;iter<nIters;iter++){
86 for(int i=0;i<displays.size();i++){
87 displays[i]->init(label);
90 //randomize the groups
91 random_shuffle(lookup.begin(), lookup.end());
93 //make merge the size of lookup[0]
94 SharedRAbundVector* merge = new SharedRAbundVector(lookup[0]->size());
96 //make copy of lookup zero
97 for(int i = 0; i<lookup[0]->size(); i++) {
98 merge->set(i, lookup[0]->getAbundance(i), "merge");
101 vector<SharedRAbundVector*> subset;
102 //send each group one at a time
103 for (int k = 0; k < lookup.size(); k++) {
104 subset.clear(); //clears out old pair of sharedrabunds
105 //add in new pair of sharedrabunds
106 subset.push_back(merge); subset.push_back(lookup[k]);
108 rcd->updateSharedData(subset, k+1, numGroupComb);
109 mergeVectors(merge, lookup[k]);
112 //resets output files
113 for(int i=0;i<displays.size();i++){
114 displays[i]->reset();
120 for(int i=0;i<displays.size();i++){
121 displays[i]->close();
126 catch(exception& e) {
127 errorOut(e, "Rarefact", "getSharedCurve");
132 /**************************************************************************************/
133 void Rarefact::mergeVectors(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
135 for (int k = 0; k < shared1->size(); k++) {
136 //merge new species into shared1
137 shared1->set(k, (shared1->getAbundance(k) + shared2->getAbundance(k)), "combo"); //set to 'combo' since this vector now contains multiple groups
140 catch(exception& e) {
141 errorOut(e, "Rarefact", "mergeVectors");