changed how rarefaction commands find the lci and hci. they now save values from...

diff --git a/mothurout.cpp b/mothurout.cpp
index 2900c7ec4fc8e0df952163554f6fefb7e622cd64..8a3b6313072a158f611e1a69a84d0a9713add18c 100644 (file)
--- a/mothurout.cpp
+++ b/mothurout.cpp
@@ -3272,6 +3272,26 @@ vector<double> MothurOut::getAverages(vector< vector<double> >& dists) {
                exit(1);
        }
 }
+/**************************************************************************************************/
+double MothurOut::getAverage(vector<double> 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<double> MothurOut::getStandardDeviation(vector< vector<double> >& dists) {
        try{

diff --git a/mothurout.h b/mothurout.h
index 845e6dd79fc5a7143b238d48bbc9deb1ae401df6..f99eac1917cbe2ce5a1ce278c66d27d1df71737e 100644 (file)
--- a/mothurout.h
+++ b/mothurout.h
@@ -163,6 +163,7 @@ class MothurOut {
         vector<double> getStandardDeviation(vector< vector<double> >&);
         vector<double> getStandardDeviation(vector< vector<double> >&, vector<double>&);
         vector<double> getAverages(vector< vector<double> >&);
+        double getAverage(vector<double>);
         vector< vector<seqDist> > getStandardDeviation(vector< vector< vector<seqDist> > >&);
         vector< vector<seqDist> > getStandardDeviation(vector< vector< vector<seqDist> > >&, vector< vector<seqDist> >&);
         vector< vector<seqDist> > getAverages(vector< vector< vector<seqDist> > >&, string);

diff --git a/raredisplay.cpp b/raredisplay.cpp
index b82127c51e6a64b50480a0f9b359cda3cd8e75f0..da22e74e90731d5294c6250d4c19f1a84045c06d 100644 (file)
--- a/raredisplay.cpp
+++ b/raredisplay.cpp
@@ -28,25 +28,14 @@ void RareDisplay::update(SAbundVector* rank){
                int newNSeqs = rank->getNumSeqs();
                vector<double> 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<int, vector<double> >::iterator it = results.find(newNSeqs);
+        if (it == results.end()) { //first iter for this count
+            vector<double> 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<SharedRAbundVector*> shared, int numSeqs, int numGroupComb) {
        try {
                vector<double> 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<int, vector<double> >::iterator it = results.find(numSeqs);
+        if (it == results.end()) { //first iter for this count
+            vector<double> 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<SharedRAbundVector*> 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<int, vector<double> >::iterator it = results.begin(); it != results.end(); it++) {
                
                        vector<double> 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<int, vector<double> >::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<double> 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<int, vector<double> >::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 a87596225b818c525fdd97156f5baf64b1dceb71..6d07efc55ab9ab677bb9f3604f46862d02984e16 100644 (file)
--- 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<int> seqs;  
-       vector<double> results;
-       vector<double> var;
-       int index, nIters;
+       map<int, vector<double> > results; //maps seqCount to results for that number of sequences
+       int nIters;
 };
 
 #endif