pcr.seqs bug fix. working on shannon range calc

[mothur.git] / shannonrange.cpp

diff --git a/shannonrange.cpp b/shannonrange.cpp
index f26f08a1f541177ae607d60be6aec0964fb19241..2b8468bd2d096f71c6c33289d371829d0a69dcba 100644 (file)
--- a/shannonrange.cpp
+++ b/shannonrange.cpp
@@ -10,18 +10,40 @@
 
 /***********************************************************************/
 
-EstOutput RangeShannon::getValues(vector<SharedRAbundVector*> shared) {
+EstOutput RangeShannon::getValues(SAbundVector* rank){
        try {
         data.resize(3,0);
         
         double commSize = 1e20;
+        double sampleSize = rank->getNumSeqs();
         
-        SAbundVector sabund1 = shared[0]->getSAbundVector();
-        SAbundVector sabund2 = shared[1]->getSAbundVector();
+        double aux = ceil(pow(sampleSize+1, (1/(double)3)));
+        double est0 = rank->get(1)+1;
+        if (aux > est0) { est0 = aux; } //est0 = max(rank->get(1)+1, aux)
+        
+        vector<double> ests;
+        double numr = 0.0;
+        for (int i = 1; i < rank->getNumBins()-1; i++) {
+            
+            if (m->control_pressed) { break; }
+            
+            int abund = rank->get(i);
+            
+            if (abund != 0) {
+            
+                int abundNext = rank->get(i+1);
+                if (abundNext == 0) {  numr = aux;  }
+                else {
+                    if (abundNext+1 > aux) { numr = abundNext+1; } //numr = max(abundNext+1, aux)
+                    else { numr = aux; }
+                }
+                double denr = aux;
+                if (abund > aux) { denr = abund; } //denr = max(abund, aux)
+                ests.push_back((abund+1)*numr/denr);
+             }
+        }
+        numr = aux;
         
-        double sampleSize = 0;
-        for (int i = 0; i < sabund1.getNumBins(); i++) {  sampleSize += (sabund1.get(i) * sabund2.get(i));  }
-        int aux = ceil(pow((sampleSize+1), 0.33333));
                
                if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }