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 int Rarefact::getCurve(float percentFreq = 0.01, int nIters = 1000){
17 RarefactionCurveData* rcd = new RarefactionCurveData();
18 for(int i=0;i<displays.size();i++){
19 rcd->registerDisplay(displays[i]);
22 //convert freq percentage to number
24 if (percentFreq < 1.0) { increment = numSeqs * percentFreq; }
25 else { increment = percentFreq; }
27 for(int iter=0;iter<nIters;iter++){
29 for(int i=0;i<displays.size();i++){
30 displays[i]->init(label);
33 RAbundVector* lookup = new RAbundVector(order->getNumBins());
34 SAbundVector* rank = new SAbundVector(order->getMaxRank()+1);
35 random_shuffle(order->begin(), order->end());
37 for(int i=0;i<numSeqs;i++){
39 if (m->control_pressed) { delete lookup; delete rank; delete rcd; return 0; }
41 int binNumber = order->get(i);
42 int abundance = lookup->get(binNumber);
44 rank->set(abundance, rank->get(abundance)-1);
47 lookup->set(binNumber, abundance);
48 rank->set(abundance, rank->get(abundance)+1);
50 if((i == 0) || (i+1) % increment == 0){
51 rcd->updateRankData(rank);
55 if(numSeqs % increment != 0){
56 rcd->updateRankData(rank);
59 for(int i=0;i<displays.size();i++){
67 for(int i=0;i<displays.size();i++){
74 m->errorOut(e, "Rarefact", "getCurve");
79 /***********************************************************************/
81 int Rarefact::getSharedCurve(float percentFreq = 0.01, int nIters = 1000){
83 SharedRarefactionCurveData* rcd = new SharedRarefactionCurveData();
85 label = lookup[0]->getLabel();
87 //register the displays
88 for(int i=0;i<displays.size();i++){
89 rcd->registerDisplay(displays[i]);
92 //if jumble is false all iters will be the same
93 if (globaldata->jumble == false) { nIters = 1; }
95 //convert freq percentage to number
97 if (percentFreq < 1.0) { increment = numSeqs * percentFreq; }
98 else { increment = percentFreq; }
100 for(int iter=0;iter<nIters;iter++){
102 for(int i=0;i<displays.size();i++){
103 displays[i]->init(label);
106 if (globaldata->jumble == true) {
107 //randomize the groups
108 random_shuffle(lookup.begin(), lookup.end());
111 //make merge the size of lookup[0]
112 SharedRAbundVector* merge = new SharedRAbundVector(lookup[0]->size());
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");
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; }
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]);
128 rcd->updateSharedData(subset, k+1, numGroupComb);
129 mergeVectors(merge, lookup[k]);
132 //resets output files
133 for(int i=0;i<displays.size();i++){
134 displays[i]->reset();
140 for(int i=0;i<displays.size();i++){
141 displays[i]->close();
147 catch(exception& e) {
148 m->errorOut(e, "Rarefact", "getSharedCurve");
153 /**************************************************************************************/
154 void Rarefact::mergeVectors(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
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
161 catch(exception& e) {
162 m->errorOut(e, "Rarefact", "mergeVectors");