]> git.donarmstrong.com Git - mothur.git/blob - bstick.cpp
added square format for phylipfile. made command work with list files containing...
[mothur.git] / bstick.cpp
1 /*
2  *  bstick.cpp
3  *  Mothur
4  *
5  *  Created by Thomas Ryabin on 3/6/09.
6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
7  *
8  */
9
10 #include "bstick.h"
11
12
13 /***********************************************************************/
14 double BStick::invSum(int index, double numSpec)
15 {
16         double sum = 0;
17         for(int i = index; i <= numSpec; i++)
18                 sum += 1/(double)i;
19         return sum;
20 }
21 /***********************************************************************/
22 RAbundVector BStick::getRAbundVector(SAbundVector* rank){
23                 vector <int> rData;
24                 int mr = 1;
25                 int nb = 0;
26                 int ns = 0;
27                 
28                 for(int i = rank->size()-1; i > 0; i--) {
29                         double cur = rank->get(i);
30                         if(mr == 1 && cur > 0)
31                                 mr = i;
32                         nb += cur;
33                         ns += i*cur;
34                         for(int j = 0; j < cur; j++)
35                                 rData.push_back(i);
36                 }
37                 
38                 RAbundVector rav = RAbundVector(rData, mr, nb, ns);
39                 rav.setLabel(rank->getLabel());
40                 return rav;
41 }
42
43 /***************************************************************************/
44
45 /***************************************************************************/
46 EstOutput BStick::getValues(SAbundVector* rank){
47         try {
48                 data.resize(3,0);
49                 rdata = getRAbundVector(rank);
50                 double numInd = (double)rdata.getNumSeqs();
51                 double numSpec = (double)rdata.getNumBins();
52                 
53                 double sumExp = 0;
54                 double sumObs = 0;
55                 double maxDiff = 0;
56
57                 for(int i = 0; i < rdata.size(); i++) {
58                         sumObs += rdata.get(i);
59                         sumExp += numInd/numSpec*invSum(i+1,numSpec);
60                         double diff = fabs(sumObs-sumExp);
61                         if(diff > maxDiff)
62                                 maxDiff = diff;
63                 }
64                 
65
66                 data[0] = maxDiff/numInd;
67                 data[1] = 0.886/sqrt(rdata.size());
68                 data[2] = 1.031/sqrt(rdata.size());
69
70                 /*m->mothurOut(critVal); m->mothurOutEndLine();
71                 m->mothurOut("If D-Statistic is less than the critical value then the data fits the Broken Stick model w/ 95% confidence.\n");*/
72                 
73
74                 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
75                 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
76                 if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
77                 
78                 return data;
79         }
80         catch(exception& e) {
81                 m->errorOut(e, "BStick", "getValues");
82                 exit(1);
83         }
84 }
85
86 /***********************************************************************/
87
88