1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
|
/*
* boneh.cpp
* Mothur
*
* Created by Thomas Ryabin on 5/13/09.
* Copyright 2009Schloss Lab UMASS Amherst. All rights reserved.
*
*/
#include "boneh.h"
#include <math.h>
/***********************************************************************/
//This solves for the value of 'v' using a binary search.
double Boneh::getV(double f1, double n, double rs) {
if(util.isEqual(rs, 0.0)) { return 0; }
double accuracy = .0001;
double v = 100000.0;
double step = v/2;
double ls = v * (1 - pow((1 - f1/(n*v)), n));
while(abs(ls - rs) > accuracy) {
if(ls > rs) { v -= step; }
else { v += step; }
ls = v * (1 - pow((1 - f1/(n * v)), n));
step /= 2;
}
return v;
}
/***********************************************************************/
EstOutput Boneh::getValues(SAbundVector* sabund){
try {
data.resize(1,0);
bool valid = false;
double sum = 0;
double n = (double)sabund->getNumSeqs();
if(util.isEqual(f,0.0)){ f=n; }
double f1 = (double)sabund->get(1);
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 < sabund->size(); j++){
sum += sabund->get(j) * pow((1 - (double)j / n), n);
}
double v = getV(f1, n, sum);
sum = 0;
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) {
m->errorOut(e, "Boneh", "getValues");
exit(1);
}
}
/***********************************************************************/
|
