]> git.donarmstrong.com Git - mothur.git/blob - sharedjackknife.cpp
added new class forest
[mothur.git] / sharedjackknife.cpp
1 /*
2  *  sharedjackknife.cpp
3  *  Mothur
4  *
5  *  Created by Thomas Ryabin on 3/30/09.
6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
7  *
8  */
9
10 #include "sharedjackknife.h"
11
12 /***************************************************************************************
13
14 ***************************************************************************************/
15 double SharedJackknife::simpson(vector<int> abunds, double numInd, int numBins){
16         double denom = numInd*(numInd-1);
17         double sum = 0;
18         for(int i = 0; i < numBins; i++)
19                 sum += (double)abunds[i]*((double)abunds[i]-1)/denom;
20         
21         return sum;
22 }
23
24 /*****************************************************************************************/
25
26 double* SharedJackknife::jackknife(){           
27         int numBins = groups.at(0)->getNumBins()-1;
28         vector<int> cArray(numBins);
29         for(int i = 0; i < numBins; i++)
30                 cArray[i] = 0;
31
32         double numInd = 0;
33         for(int i = 0; i < numGroups; i++)
34                 for(int j = 0; j < numBins; j++) {
35                         int curAbund = groups.at(i)->get(j+1).abundance;
36                         cArray[j] += curAbund;
37                         numInd += (double)curAbund;
38                 }
39
40         double baseD = 1/simpson(cArray, numInd, numBins);
41         
42         vector<double> pseudoVals(numBins);
43         double jackknifeEstimate = 0;
44         for(int i = 0; i < numGroups; i++) {
45                 for(int j = 0; j < numBins-1; j++) {
46                         int abundDiff = -groups.at(i)->get(j+1).abundance;
47                         if(i > 0)
48                                 abundDiff += groups.at(i-1)->get(j+1).abundance;
49
50                         cArray[j] += abundDiff;
51                         numInd += abundDiff;    
52                 }
53                 
54                 double curD = 1/simpson(cArray, numInd, numBins);
55                 pseudoVals[i] = (double)numGroups*(baseD - curD) + curD;
56                 jackknifeEstimate += pseudoVals[i];
57         }
58         jackknifeEstimate /= (double)numGroups;
59                 
60         double variance = 0;
61         for(int i = 0; i < numGroups; i++)
62                 variance += pow(pseudoVals[i]-jackknifeEstimate, 2);
63                 
64         variance /= (double)numGroups*((double)numGroups-1);
65         double stErr = sqrt(variance);
66         TDTable table;
67         double confLimit = 0;
68         if(numGroups <= 30)
69                 confLimit = table.getConfLimit(numGroups-1, 1);
70         else
71                 confLimit = 1.645;
72         
73         confLimit *= stErr;
74         
75         double* rdata = new double[3];
76         rdata[0] = baseD;
77         rdata[1] = jackknifeEstimate - confLimit;
78         rdata[2] = jackknifeEstimate + confLimit;
79         
80         return rdata;
81 }
82
83 /************************************************************************************************
84 ************************************************************************************************/
85
86 EstOutput SharedJackknife::getValues(vector<SharedRAbundVector*> vectorShared){ //Fix this for collect, mistake was that it was made with summary in mind.
87         try {
88                 SharedRAbundVector* shared1 = vectorShared[0];
89                 SharedRAbundVector* shared2 = vectorShared[1];
90                 if(numGroups == -1) {
91                         numGroups = m->getNumGroups();
92                 }
93
94                 if(callCount == numGroups*(numGroups-1)/2) {
95                         currentCallDone = true;
96                         callCount = 0;
97                 }
98                 callCount++;
99
100                 if(currentCallDone) {
101                         groups.clear(); 
102                         currentCallDone = false;
103                 }
104                 
105                 if(groups.size() != numGroups) {        
106                         if(groups.size() == 0)
107                                 groups.push_back(shared1);
108                         groups.push_back(shared2);
109                 }
110                 
111                 if(groups.size() == numGroups && callCount < numGroups) {
112                         data.resize(3,0);
113
114                         double* rdata = jackknife();
115                         data[0] = rdata[0];
116                         data[1] = rdata[1];
117                         data[2] = rdata[2];
118                 
119                         if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
120                         if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
121                         if (isnan(data[2]) || isinf(data[0])) { data[2] = 0; }
122                         
123                         return data;
124                 }
125                 
126                 data.resize(3,0);
127                 data[0] = 0;
128                 data[1] = 0;
129                 data[2] = 0;
130                 
131                 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
132                 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
133                 if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
134                 return data;    
135         }
136                 
137         catch(exception& e) {
138                 m->errorOut(e, "SharedJackknife", "getValues");
139                 exit(1);
140         }
141 }
142
143 /***********************************************************************/
144