#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 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++) {
- double Xi = 0; //I didn't know what this was, simply replace the 0
- //with the appropriate expression for the boneh calculator
- //to work.
- sum += pow(1 - Xi / n, n) * (1 - pow(1 - Xi / 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;
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);
- }
-};
+}
/***********************************************************************/