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