From de3d202264ab8a55fc91e0c2776aa5bed92bbf55 Mon Sep 17 00:00:00 2001 From: Sarah Westcott Date: Fri, 17 May 2013 14:44:40 -0400 Subject: [PATCH] changed how rarefaction commands find the lci and hci. they now save values from each iter at each frequency, report values at the .025 and .975. --- mothurout.cpp | 20 ++++++++ mothurout.h | 1 + raredisplay.cpp | 128 ++++++++++++++++++++++-------------------------- raredisplay.h | 8 ++- 4 files changed, 83 insertions(+), 74 deletions(-) diff --git a/mothurout.cpp b/mothurout.cpp index 2900c7e..8a3b631 100644 --- a/mothurout.cpp +++ b/mothurout.cpp @@ -3272,6 +3272,26 @@ vector MothurOut::getAverages(vector< vector >& dists) { exit(1); } } +/**************************************************************************************************/ +double MothurOut::getAverage(vector dists) { + try{ + double average = 0; + + for (int i = 0; i < dists.size(); i++) { + average += dists[i]; + } + + //finds average. + average /= (double) dists.size(); + + return average; + } + catch(exception& e) { + errorOut(e, "MothurOut", "getAverage"); + exit(1); + } +} + /**************************************************************************************************/ vector MothurOut::getStandardDeviation(vector< vector >& dists) { try{ diff --git a/mothurout.h b/mothurout.h index 845e6dd..f99eac1 100644 --- a/mothurout.h +++ b/mothurout.h @@ -163,6 +163,7 @@ class MothurOut { vector getStandardDeviation(vector< vector >&); vector getStandardDeviation(vector< vector >&, vector&); vector getAverages(vector< vector >&); + double getAverage(vector); vector< vector > getStandardDeviation(vector< vector< vector > >&); vector< vector > getStandardDeviation(vector< vector< vector > >&, vector< vector >&); vector< vector > getAverages(vector< vector< vector > >&, string); diff --git a/raredisplay.cpp b/raredisplay.cpp index b82127c..da22e74 100644 --- a/raredisplay.cpp +++ b/raredisplay.cpp @@ -28,25 +28,14 @@ void RareDisplay::update(SAbundVector* rank){ int newNSeqs = rank->getNumSeqs(); vector data = estimate->getValues(rank); - if(nIters != 1){ - - double oldS = var[index] * ( nIters - 2 ); - double delta = data[0] - results[index]; - double newMean = results[index] + delta / nIters; - double newS = oldS + delta * ( data[0] - newMean ); - double newVar = newS / ( nIters - 1 ); - - seqs[index] = newNSeqs; - results[index] = newMean; - var[index] = newVar; - - index++; - } - else{ - seqs.push_back(newNSeqs); - results.push_back(data[0]); - var.push_back(0.0); - } + map >::iterator it = results.find(newNSeqs); + if (it == results.end()) { //first iter for this count + vector temp; + temp.push_back(data[0]); + results[newNSeqs] = temp; + }else { + it->second.push_back(data[0]); + } } catch(exception& e) { m->errorOut(e, "RareDisplay", "update"); @@ -58,28 +47,15 @@ void RareDisplay::update(SAbundVector* rank){ void RareDisplay::update(vector shared, int numSeqs, int numGroupComb) { try { vector data = estimate->getValues(shared); - double newNSeqs = data[0]; - if(nIters != 1){ - - double oldS = var[index] * ( nIters - 2 ); - double delta = data[0] - results[index]; - double newMean = results[index] + delta / nIters; - double newS = oldS + delta * ( data[0] - newMean ); - double newVar = newS / ( nIters - 1 ); - seqs[index] = newNSeqs; - results[index] = newMean; - var[index] = newVar; - - index++; - } - else{ - - seqs.push_back(newNSeqs); - results.push_back(data[0]); - var.push_back(0.0); - - } + map >::iterator it = results.find(numSeqs); + if (it == results.end()) { //first iter for this count + vector temp; + temp.push_back(data[0]); + results[numSeqs] = temp; + }else { + it->second.push_back(data[0]); + } } catch(exception& e) { m->errorOut(e, "RareDisplay", "update"); @@ -92,7 +68,6 @@ void RareDisplay::update(vector shared, int numSeqs, int nu void RareDisplay::reset(){ try { nIters++; - index = 0; } catch(exception& e) { m->errorOut(e, "RareDisplay", "reset"); @@ -104,29 +79,28 @@ void RareDisplay::reset(){ void RareDisplay::close(){ try { - output->initFile(label); - for (int i = 0; i < seqs.size(); i++) { + for (map >::iterator it = results.begin(); it != results.end(); it++) { vector data(3,0); - double variance = var[i]; - - data[0] = results[i]; - - double ci = 1.96 * pow(variance, 0.5); - data[1] = data[0] - ci; - data[2] = data[0] + ci; + + sort((it->second).begin(), (it->second).end()); + + //lci=results[int(0.025*iter)] and hci=results[int(0.975*iter)] + data[0] = (it->second)[(int)(0.50*(nIters-1))]; + //data[0] = m->getAverage(it->second); + data[1] = (it->second)[(int)(0.025*(nIters-1))]; + data[2] = (it->second)[(int)(0.975*(nIters-1))]; + //cout << nIters << '\t' << (int)(0.025*(nIters-1)) << '\t' << (int)(0.975*(nIters-1)) << endl; + + //cout << it->first << '\t' << data[1] << '\t' << data[2] << endl; - output->output(seqs[i], data); + output->output(it->first, data); } nIters = 1; - index = 0; - - seqs.clear(); - results.clear(); - var.clear(); + results.clear(); output->resetFile(); } @@ -142,17 +116,29 @@ void RareDisplay::inputTempFiles(string filename){ ifstream in; m->openInputFile(filename, in); - int thisIters; - in >> thisIters; m->gobble(in); + int thisIters, size; + in >> thisIters >> size; m->gobble(in); + nIters += thisIters; - for (int i = 0; i < seqs.size(); i++) { - double tempresult, tempvar; - in >> tempresult >> tempvar; m->gobble(in); - - //find weighted result - results[i] = ((nIters * results[i]) + (thisIters * tempresult)) / (float)(nIters + thisIters); - - var[i] = ((nIters * var[i]) + (thisIters * tempvar)) / (float)(nIters + thisIters); + for (int i = 0; i < size; i++) { + int tempCount; + in >> tempCount; m->gobble(in); + map >::iterator it = results.find(tempCount); + if (it != results.end()) { + for (int j = 0; j < thisIters; j++) { + double value; + in >> value; m->gobble(in); + (it->second).push_back(value); + } + }else { + vector tempValues; + for (int j = 0; j < thisIters; j++) { + double value; + in >> value; m->gobble(in); + tempValues.push_back(value); + } + results[tempCount] = tempValues; + } } in.close(); @@ -170,10 +156,14 @@ void RareDisplay::outputTempFiles(string filename){ ofstream out; m->openOutputFile(filename, out); - out << nIters << endl; + out << nIters-1 << '\t' << results.size() << endl; - for (int i = 0; i < seqs.size(); i++) { - out << results[i] << '\t' << var[i] << endl; + for (map >::iterator it = results.begin(); it != results.end(); it++) { + out << it->first << '\t'; + for(int i = 0; i < (it->second).size(); i++) { + out << (it->second)[i] << '\t'; + } + out << endl; } out.close(); diff --git a/raredisplay.h b/raredisplay.h index a875962..6d07efc 100644 --- a/raredisplay.h +++ b/raredisplay.h @@ -11,7 +11,7 @@ class RareDisplay : public Display { public: - RareDisplay(Calculator* calc, FileOutput* file) : estimate(calc), output(file), nIters(1), index(0) {}; + RareDisplay(Calculator* calc, FileOutput* file) : estimate(calc), output(file), nIters(1) {}; ~RareDisplay() { delete estimate; delete output; }; void init(string); void reset(); @@ -27,10 +27,8 @@ private: Calculator* estimate; FileOutput* output; string label; - vector seqs; - vector results; - vector var; - int index, nIters; + map > results; //maps seqCount to results for that number of sequences + int nIters; }; #endif -- 2.39.2