X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=boneh.cpp;h=b4ce2baf81a4371a2f59959f0658f11e50a13fe3;hp=7b52683d36b21b74653e9e074429f07060176e2d;hb=372fb21ea66ced432b109225851a1b80ef0491a3;hpb=37a2c0b7dc68bfa53c8f38a713259929c4c46a9f

diff --git a/boneh.cpp b/boneh.cpp
index 7b52683..b4ce2ba 100644
--- a/boneh.cpp
+++ b/boneh.cpp
@@ -3,7 +3,7 @@
  *  Mothur
  *
  *  Created by Thomas Ryabin on 5/13/09.
- *  Copyright 2009 __MyCompanyName__. All rights reserved.
+ *  Copyright 2009Schloss Lab UMASS Amherst. All rights reserved.
  *
  */
 
@@ -11,62 +11,65 @@
 #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;
 	
 	double accuracy = .0001;
 	double v = 100000.0;
 	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->size() - 1;
-		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++)
-				sum += pow(1 - (double)rank->get(j) / n, n) * (1 - pow(1 - (double)rank->get(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;
@@ -74,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);
-	}	
-};
+}
 
 
 /***********************************************************************/