changes while testing

[mothur.git] / boneh.cpp

diff --git a/boneh.cpp b/boneh.cpp
index 2025ae70d10311c6184ee723097be4e67b22b318..b4ce2baf81a4371a2f59959f0658f11e50a13fe3 100644 (file)
--- a/boneh.cpp
+++ b/boneh.cpp
@@ -11,11 +11,10 @@
 #include <math.h>
 
 /***********************************************************************/
+
 //This solves for the value of 'v' using a binary search.
 double Boneh::getV(double f1, double n, double rs) {
 
-       //cout << "f1 = " << f1 << "\nn = " << n << "\nrs = " << rs << "\n\n";
-       
        if(rs == 0)
                return 0;
        
@@ -24,56 +23,53 @@ double Boneh::getV(double f1, double n, double rs) {
        double step = v/2;
        double ls = v * (1 - pow((1 - f1/(n*v)), n));
        
-       //cout << "ls = " << ls << "\n";
-       
        while(abs(ls - rs) > accuracy) {
-               if(ls > rs)
-                       v -= step;
-               else
-                       v += step;
+               if(ls > rs)     {       v -= step;      }
+               else            {       v += step;      }
                
-               ls = v* (1 - pow((1 - f1/(n*v)), n));
+               ls = v * (1 - pow((1 - f1/(n * v)), n));
                step /= 2;
-               
-               //cout << "ls = " << ls << "\n";
        }
-       
+
        return v;
 }
        
 /***********************************************************************/      
-EstOutput Boneh::getValues(SAbundVector* rank){
+EstOutput Boneh::getValues(SAbundVector* sabund){
 
        try {
                data.resize(1,0);
                
+
                bool valid = false;
                double sum = 0;
-               double n = (double)rank->getNumSeqs();
-               double f1 = (double)rank->get(1);
+               double n = (double)sabund->getNumSeqs();
+               if(f==0){       f=n;    }
                
-               for(int i = 1; i < rank->size(); i++)
-                       sum += (double)rank->get(i) * exp(-i);
+               double f1 = (double)sabund->get(1);
                
-               if(rank->get(1) > sum)
+               for(int i = 1; i < sabund->size(); i++){
+                       sum += (double)sabund->get(i) * exp(-i);
+               }
+
+               if(sabund->get(1) > sum)
                        valid = true;
                
                sum = 0;
                if(valid) {
-                       for(int j = 1; j < rank->size(); j++)
-                               sum += rank->get(j) * pow((1 - (double)j / n), n);
+                       for(int j = 1; j < sabund->size(); j++){
+                               sum += sabund->get(j) * pow((1 - (double)j / n), n);
+                       }
                        
                        double v = getV(f1, n, sum);
                        
-                       //cout << "v = " << v << "\n";
-                       
-                       
                        sum = 0;
-                       for(int j = 1; j < rank->size(); j++) {
-                               for (int i = 0; i < rank->get(j); i++) {
-                                       sum += pow(1 - j / n, n) * (1 - pow(1 - j / n, m)) + v * pow(1 - f1/(n*v), n) * (1 - pow(1 - f1/(n*v), m));
+                       for(int j = 1; j < sabund->size(); j++) {
+                               for (int i = 0; i < sabund->get(j); i++) {
+                                       sum += pow(1 - j / n, n) * (1 - pow(1 - j / n, f));
                                }
                        }
+                       sum +=  v * pow(1 - f1/(n*v), n) * (1 - pow(1 - f1/(n*v), f));
                }
 
                data[0] = sum;
@@ -81,14 +77,10 @@ EstOutput Boneh::getValues(SAbundVector* rank){
                return data;
        }
        catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the Coverage class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               m->errorOut(e, "Boneh", "getValues");
                exit(1);
        }
-       catch(...) {
-               cout << "An unknown error has occurred in the Coverage class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }       
-};
+}
 
 
 /***********************************************************************/