X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=geom.cpp;fp=geom.cpp;h=0000000000000000000000000000000000000000;hb=4a877efa127e56e81a21f53cfdbbfd3bfbe8c4ff;hp=5bd150ef7f8201254a392a52a7583b07b83e0b00;hpb=a6cf29fa4dac0909c7582cb1094151d34093ee76;p=mothur.git diff --git a/geom.cpp b/geom.cpp deleted file mode 100644 index 5bd150e..0000000 --- a/geom.cpp +++ /dev/null @@ -1,101 +0,0 @@ -/* - * geom.cpp - * Mothur - * - * Created by Thomas Ryabin on 2/23/09. - * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved. - * - */ - -#include "geom.h" - -/***********************************************************************/ - -double Geom::kEq(double k, double spec){ - return k/(1-k)*pow(1-k, spec)/(1-pow(1-k, spec)); -} - -RAbundVector Geom::getRAbundVector(SAbundVector* rank){ - vector rData; - int mr = 1; - int nb = 0; - int ns = 0; - - for(int i = rank->size()-1; i > 0; i--) { - int cur = rank->get(i); - if(mr == 1 && cur > 0) - mr = i; - nb += cur; - ns += i*cur; - for(int j = 0; j < cur; j++) - rData.push_back(i); - } - - RAbundVector rav = RAbundVector(rData, mr, nb, ns); - rav.setLabel(rank->getLabel()); - return rav; -} - -/***********************************************************************************/ - -/***********************************************************************************/ -EstOutput Geom::getValues(SAbundVector* rank){ - try { - data.resize(3,0); - - rdata = getRAbundVector(rank); - 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. - if(numInd*kEq(k, numSpec) > min) - k += step; - else - k -= step; - step /= 2; - } - - double cK = 1/(1-pow(1-k, numSpec)); - double sumExp = 0; - double sumObs = 0; - double maxDiff = 0; - - for(int i = 0; i < numSpec; i++) - { - sumObs += rdata.get(i); - sumExp += numInd*cK*k*pow(1-k, i); - - double diff = fabs(sumObs-sumExp); - if(diff > maxDiff) { maxDiff = diff; } - - } - - - 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; } - if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; } - - return data; - } - catch(exception& e) { - m->errorOut(e, "Geom", "getValues"); - exit(1); - } -} - -/***********************************************************************/ - - - - - - -