added shen calculator

diff --git a/boneh.cpp b/boneh.cpp
index 7b52683d36b21b74653e9e074429f07060176e2d..5174f7519477b9c5ad8301ed0fa2d3471e2d1700 100644 (file)
--- a/boneh.cpp
+++ b/boneh.cpp
@@ -11,10 +11,14 @@
 #include <math.h>
 
 /***********************************************************************/
+//This solves for the value of 'v' using a binary search.
 double Boneh::getV(double f1, double n, double rs) {
 
        //cout << "f1 = " << f1 << "\nn = " << n << "\nrs = " << rs << "\n\n";
        
+       if(rs == 0)
+               return 0;
+       
        double accuracy = .0001;
        double v = 100000.0;
        double step = v/2;
@@ -31,7 +35,7 @@ double Boneh::getV(double f1, double n, double rs) {
                ls = v* (1 - pow((1 - f1/(n*v)), n));
                step /= 2;
                
-       //      cout << "ls = " << ls << "\n";
+               //cout << "ls = " << ls << "\n";
        }
        
        return v;
@@ -45,7 +49,7 @@ EstOutput Boneh::getValues(SAbundVector* rank){
                
                bool valid = false;
                double sum = 0;
-               double n = (double)rank->size() - 1;
+               double n = (double)rank->getNumSeqs();
                double f1 = (double)rank->get(1);
                
                for(int i = 1; i < rank->size(); i++)
@@ -61,12 +65,16 @@ EstOutput Boneh::getValues(SAbundVector* rank){
                        
                        double v = getV(f1, n, sum);
                        
-               //      cout << "v = " << v << "\n";
+                       //cout << "v = " << v << "\n";
+                       
                        
                        sum = 0;
-                       for(int j = 1; j < rank->size(); j++)
-                               sum += pow(1 - (double)rank->get(j) / n, n) * (1 - pow(1 - (double)rank->get(j) / n, m)) + v * pow(1 - f1/(n*v), n) * (1 - pow(1 - f1/(n*v), m));
-                               
+                       for(int j = 1; j < rank->size(); j++) {
+                               double Xi = 0; //I didn't know what this was, simply replace the 0
+                                                          //with the appropriate expression for the boneh calculator
+                                                          //to work.
+                               sum += pow(1 - Xi / n, n) * (1 - pow(1 - Xi / n, m)) + v * pow(1 - f1/(n*v), n) * (1 - pow(1 - f1/(n*v), m));
+                       }
                }
 
                data[0] = sum;

diff --git a/collectcommand.cpp b/collectcommand.cpp
index e429dffb19721c6f8c124734f78f7ff2bf937e3c..d317a06aff4ebf68557a081475cd15a4df58c86a 100644 (file)
--- a/collectcommand.cpp
+++ b/collectcommand.cpp
@@ -26,7 +26,7 @@
 #include "efron.h"
 #include "boneh.h"
 #include "solow.h"
-
+#include "shen.h"
 #include "coverage.h"
 
 
@@ -82,6 +82,9 @@ CollectCommand::CollectCommand(){
                                }else if (globaldata->Estimators[i] == "solow") {
                                        convert(globaldata->getSize(), size); 
                                        cDisplays.push_back(new CollectDisplay(new Solow(size), new OneColumnFile(fileNameRoot+"solow")));
+                               }else if (globaldata->Estimators[i] == "shen") {
+                                       convert(globaldata->getSize(), size); 
+                                       cDisplays.push_back(new CollectDisplay(new Shen(size), new OneColumnFile(fileNameRoot+"shen")));
                                }
                        }
                }

diff --git a/efron.cpp b/efron.cpp
index e84364160e2a48dd6e87927ade301c942565088d..0441f265ad45c2b929e0fa493b789b33559b41fe 100644 (file)
--- a/efron.cpp
+++ b/efron.cpp
@@ -16,15 +16,16 @@ EstOutput Efron::getValues(SAbundVector* rank){
        try {
                data.resize(1,0);
                
-               if(m > rank->size()-1) {
-                       cout << "Error in the 'efron' calculator. 'size' must be less than the length of the smalled sabund vector.\n";
+               double n = (double)rank->getNumSeqs();
+               if(m > n) {
+                       cout << "Error in the 'efron' calculator. 'size' must be less than the length of the smallest sabund vector.\n";
                        data[0] = 0;
                        return data;
                }
                
                double sum = 0;
                for(int i = 1; i < rank->size(); i++)
-                       sum += pow(-1, (double)(i+1)) * pow(((double)m / (double)(rank->size()-1)), i) * (double)(rank->get(i));
+                       sum += pow(-1, (double)(i+1)) * pow(((double)m / n), i) * (double)(rank->get(i));
 
                data[0] = sum;
                

diff --git a/globaldata.cpp b/globaldata.cpp
index 88c5a13aa364239131f0d86ecb2b99dfe25d5930..8b8f9e39b083d9d4329ff500f335d9f4f04cb786 100644 (file)
--- a/globaldata.cpp
+++ b/globaldata.cpp
@@ -203,7 +203,7 @@ void GlobalData::parseGlobalData(string commandString, string optionText){
                //input defaults for calculators
                if (commandName == "collect.single") {
 
-                       if ((calc == "default") || (calc == "")) { calc = "sobs-chao-ace-jack-shannon-npshannon-simpson-efron-boneh-solow"; }
+                       if ((calc == "default") || (calc == "")) { calc = "sobs-chao-ace-jack-shannon-npshannon-simpson-efron-boneh-solow-shen"; }
                        Estimators.clear();
                        splitAtDash(calc, Estimators); 
                }
@@ -219,7 +219,7 @@ void GlobalData::parseGlobalData(string commandString, string optionText){
                        splitAtDash(calc, Estimators); 
                }
                if (commandName == "summary.single") {
-                       if ((calc == "default") || (calc == "")) { calc = "sobs-chao-ace-jack-shannon-npshannon-simpson-efron-boneh-solow"; }
+                       if ((calc == "default") || (calc == "")) { calc = "sobs-chao-ace-jack-shannon-npshannon-simpson-efron-boneh-solow-shen"; }
                        Estimators.clear();
                        splitAtDash(calc, Estimators); 
                }

diff --git a/shen.cpp b/shen.cpp
new file mode 100644 (file)
index 0000000..6bdb9b6
--- /dev/null
+++ b/shen.cpp
@@ -0,0 +1,44 @@
+/*
+ *  shen.cpp
+ *  Mothur
+ *
+ *  Created by Thomas Ryabin on 5/18/09.
+ *  Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "shen.h"
+#include <math.h>
+
+       
+/***********************************************************************/      
+EstOutput Shen::getValues(SAbundVector* rank){
+
+       try {
+               data.resize(1,0);
+               
+               double n = (double)rank->getNumSeqs();
+               double f1 = (double)rank->get(1);
+               
+               double D_rare = 0; //I didn't know what this was. Simply replace the '0' with the appropriate expression.
+               double C_rare = 1; //I didn't know what this was. Simply replace the '1' with the appropriate expression.
+               double Y_rare = 1; //I didn't know what this was. Simply replace the '1' with the appropriate expression.
+               
+               double f0 = D_rare/C_rare + f1/C_rare * Y_rare*Y_rare - D_rare;
+               
+               data[0] = f0 * (1 - pow(1 - f1/n/f0, m));
+
+               return data;
+       }
+       catch(exception& e) {
+               cout << "Standard Error: " << e.what() << " has occurred in the Coverage class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }
+       catch(...) {
+               cout << "An unknown error has occurred in the Coverage class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }       
+};
+
+
+/***********************************************************************/

diff --git a/shen.h b/shen.h
new file mode 100644 (file)
index 0000000..4a97805
--- /dev/null
+++ b/shen.h
@@ -0,0 +1,33 @@
+#ifndef SHEN_H
+#define SHEN_H
+
+/*
+ *  shen.h
+ *  Mothur
+ *
+ *  Created by Thomas Ryabin on 5/18/09.
+ *  Copyright 2009 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "calculator.h"
+
+/* This class implements the shen calculator on single group. 
+ It is a child of the calculator class. */
+
+/***********************************************************************/
+
+class Shen : public Calculator  {
+       
+public: 
+       Shen(int size) : m(size), Calculator("shen", 1, false) {};
+       EstOutput getValues(SAbundVector*);     
+       EstOutput getValues(vector<SharedRAbundVector*>) {return data;};
+private:
+       int m;
+};
+
+
+/***********************************************************************/
+
+#endif

diff --git a/solow.cpp b/solow.cpp
index 3a15fdfd68dea1fe0f91d5ca3fd1e87df10cb745..8c3d155add993c93a33719686baf3911641acde3 100644 (file)
--- a/solow.cpp
+++ b/solow.cpp
@@ -17,7 +17,7 @@ EstOutput Solow::getValues(SAbundVector* rank){
        try {
                data.resize(1,0);
                
-               double n = (double)rank->size() - 1;
+               double n = (double)rank->getNumSeqs();
                double f1 = (double)rank->get(1);
                double f2 = (double)rank->get(2);
 

diff --git a/summarycommand.cpp b/summarycommand.cpp
index d0124d8516033800804e16cf5a092af8517e7fe7..b649582758636370f812a51cf182ce7813bc3e02 100644 (file)
--- a/summarycommand.cpp
+++ b/summarycommand.cpp
@@ -27,6 +27,7 @@
 #include "efron.h"
 #include "boneh.h"
 #include "solow.h"
+#include "shen.h"
 
 //**********************************************************************************************************************
 
@@ -82,6 +83,9 @@ SummaryCommand::SummaryCommand(){
                                }else if (globaldata->Estimators[i] == "solow") { 
                                        convert(globaldata->getSize(), size);
                                        sumCalculators.push_back(new Solow(size));
+                               }else if (globaldata->Estimators[i] == "shen") { 
+                                       convert(globaldata->getSize(), size);
+                                       sumCalculators.push_back(new Shen(size));
                                }
                        }
                }

diff --git a/validcalculator.cpp b/validcalculator.cpp
index 7a07fd1af9e7eab2b7e7551cb170a852b62aa031..edcf8c7a33269ece0758903f4d1fc425da558e97 100644 (file)
--- a/validcalculator.cpp
+++ b/validcalculator.cpp
@@ -204,6 +204,7 @@ void ValidCalculators::initialSingle() {
                single["efron"]         = "efron";
                single["boneh"]         = "boneh";
                single["solow"]         = "solow";
+               single["shen"]          = "shen";
                single["default"]           = "default";
        }
        catch(exception& e) {
@@ -300,6 +301,7 @@ void ValidCalculators::initialSummary() {
                summary["efron"]        = "efron";
                summary["boneh"]        = "boneh";
                summary["solow"]        = "solow";
+               summary["shen"]         = "shen";
                summary["default"]          = "default";
        }
        catch(exception& e) {