X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=geom.cpp;h=5bd150ef7f8201254a392a52a7583b07b83e0b00;hp=31580db15b7bbb08b648cd10ff3c0bcf4cd98a10;hb=615301e57c25e241356a9c2380648d117709458d;hpb=e4827e0945cbda536064e5a345996b2a7dfcbb81 diff --git a/geom.cpp b/geom.cpp index 31580db..5bd150e 100644 --- a/geom.cpp +++ b/geom.cpp @@ -21,8 +21,7 @@ RAbundVector Geom::getRAbundVector(SAbundVector* rank){ int nb = 0; int ns = 0; - for(int i = rank->size()-1; i > 0; i--) - { + for(int i = rank->size()-1; i > 0; i--) { int cur = rank->get(i); if(mr == 1 && cur > 0) mr = i; @@ -45,14 +44,13 @@ EstOutput Geom::getValues(SAbundVector* rank){ data.resize(3,0); rdata = getRAbundVector(rank); - int numInd = rdata.getNumSeqs(); - int numSpec = rdata.getNumBins(); - int min = rdata.get(rdata.size()-1); + double numInd = rdata.getNumSeqs(); + double numSpec = rdata.getNumBins(); + double min = rdata.get(rdata.size()-1); double k = .5; double step = .49999; - while(fabs(min - numInd*kEq(k, (double)numSpec)) > .0001) //This uses a binary search to find the value of k. - { + while(fabs(min - numInd*kEq(k, (double)numSpec)) > .0001) { //This uses a binary search to find the value of k. if(numInd*kEq(k, numSpec) > min) k += step; else @@ -75,25 +73,11 @@ EstOutput Geom::getValues(SAbundVector* rank){ } - double DStatistic = maxDiff/numInd; - double critVal = 0; - /*cout << "Geom:\n"; - cout << "D-Statistic = " << DStatistic << "\n"; - cout << "Critical value for 95% confidence interval = ";*/ - if(rdata.size() > 20) - { - data[1] = 1.031/sqrt(rdata.size()); - data[2] = 0.886/sqrt(rdata.size()); - } - else - { - KOSTable table; - critVal = table.getConfLimit(numSpec); - } - /*cout << critVal << "\n"; - cout << "If D-Statistic is less than the critical value then the data fits the Geometric Series model w/ 95% confidence.\n\n";*/ - - data[0] = DStatistic; + + data[0] = maxDiff/numInd; + data[1] = 0.886/sqrt(numSpec); + data[2] = 1.031/sqrt(numSpec); + if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; } @@ -102,13 +86,9 @@ EstOutput Geom::getValues(SAbundVector* rank){ return data; } catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + m->errorOut(e, "Geom", "getValues"); exit(1); } - catch(...) { - cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } } /***********************************************************************/