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 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
29 driver(rcd, increment, nIters);
31 vector<int> procIters;
33 int numItersPerProcessor = nIters / processors;
35 //divide iters between processes
36 for (int i = 0; i < processors; i++) {
37 if(i == processors - 1){
38 numItersPerProcessor = nIters - i * numItersPerProcessor;
40 procIters.push_back(numItersPerProcessor);
43 createProcesses(procIters, rcd, increment);
47 driver(rcd, increment, nIters);
50 for(int i=0;i<displays.size();i++){
58 m->errorOut(e, "Rarefact", "getCurve");
62 /***********************************************************************/
63 int Rarefact::driver(RarefactionCurveData* rcd, int increment, int nIters = 1000){
66 for(int iter=0;iter<nIters;iter++){
68 for(int i=0;i<displays.size();i++){
69 displays[i]->init(label);
72 RAbundVector* lookup = new RAbundVector(order->getNumBins());
73 SAbundVector* rank = new SAbundVector(order->getMaxRank()+1);
74 random_shuffle(order->begin(), order->end());
76 for(int i=0;i<numSeqs;i++){
78 if (m->control_pressed) { delete lookup; delete rank; delete rcd; return 0; }
80 int binNumber = order->get(i);
81 int abundance = lookup->get(binNumber);
83 rank->set(abundance, rank->get(abundance)-1);
86 lookup->set(binNumber, abundance);
87 rank->set(abundance, rank->get(abundance)+1);
89 if((i == 0) || ((i+1) % increment == 0) || (ends.count(i+1) != 0)){
90 rcd->updateRankData(rank);
94 if((numSeqs % increment != 0) || (ends.count(numSeqs) != 0)){
95 rcd->updateRankData(rank);
98 for(int i=0;i<displays.size();i++){
108 catch(exception& e) {
109 m->errorOut(e, "Rarefact", "driver");
113 /**************************************************************************************************/
115 int Rarefact::createProcesses(vector<int>& procIters, RarefactionCurveData* rcd, int increment) {
117 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
120 vector<int> processIDS;
124 //loop through and create all the processes you want
125 while (process != processors) {
129 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
132 driver(rcd, increment, procIters[process]);
134 //pass numSeqs to parent
135 for(int i=0;i<displays.size();i++){
136 string tempFile = toString(getpid()) + toString(i) + ".rarefact.temp";
137 displays[i]->outputTempFiles(tempFile);
141 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
142 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
147 driver(rcd, increment, procIters[0]);
149 //force parent to wait until all the processes are done
150 for (int i=0;i<(processors-1);i++) {
151 int temp = processIDS[i];
155 //get data created by processes
156 for (int i=0;i<(processors-1);i++) {
157 for(int j=0;j<displays.size();j++){
158 string s = toString(processIDS[i]) + toString(j) + ".rarefact.temp";
159 displays[j]->inputTempFiles(s);
167 catch(exception& e) {
168 m->errorOut(e, "Rarefact", "createProcesses");
172 /***********************************************************************/
174 int Rarefact::getSharedCurve(float percentFreq = 0.01, int nIters = 1000){
176 SharedRarefactionCurveData* rcd = new SharedRarefactionCurveData();
178 label = lookup[0]->getLabel();
180 //register the displays
181 for(int i=0;i<displays.size();i++){
182 rcd->registerDisplay(displays[i]);
185 //if jumble is false all iters will be the same
186 if (m->jumble == false) { nIters = 1; }
188 //convert freq percentage to number
190 if (percentFreq < 1.0) { increment = numSeqs * percentFreq; }
191 else { increment = percentFreq; }
193 for(int iter=0;iter<nIters;iter++){
195 for(int i=0;i<displays.size();i++){
196 displays[i]->init(label);
199 if (m->jumble == true) {
200 //randomize the groups
201 random_shuffle(lookup.begin(), lookup.end());
204 //make merge the size of lookup[0]
205 SharedRAbundVector* merge = new SharedRAbundVector(lookup[0]->size());
207 //make copy of lookup zero
208 for(int i = 0; i<lookup[0]->size(); i++) {
209 merge->set(i, lookup[0]->getAbundance(i), "merge");
212 vector<SharedRAbundVector*> subset;
213 //send each group one at a time
214 for (int k = 0; k < lookup.size(); k++) {
215 if (m->control_pressed) { delete merge; delete rcd; return 0; }
217 subset.clear(); //clears out old pair of sharedrabunds
218 //add in new pair of sharedrabunds
219 subset.push_back(merge); subset.push_back(lookup[k]);
221 rcd->updateSharedData(subset, k+1, numGroupComb);
222 mergeVectors(merge, lookup[k]);
225 //resets output files
226 for(int i=0;i<displays.size();i++){
227 displays[i]->reset();
233 for(int i=0;i<displays.size();i++){
234 displays[i]->close();
240 catch(exception& e) {
241 m->errorOut(e, "Rarefact", "getSharedCurve");
246 /**************************************************************************************/
247 void Rarefact::mergeVectors(SharedRAbundVector* shared1, SharedRAbundVector* shared2) {
249 for (int k = 0; k < shared1->size(); k++) {
250 //merge new species into shared1
251 shared1->set(k, (shared1->getAbundance(k) + shared2->getAbundance(k)), "combo"); //set to 'combo' since this vector now contains multiple groups
254 catch(exception& e) {
255 m->errorOut(e, "Rarefact", "mergeVectors");