From 2bb51be6ba6a3fa741f2131d328cf21c395b506a Mon Sep 17 00:00:00 2001 From: westcott Date: Mon, 18 May 2009 19:56:29 +0000 Subject: [PATCH] added shen calculator --- boneh.cpp | 20 ++++++++++++++------ collectcommand.cpp | 5 ++++- efron.cpp | 7 ++++--- globaldata.cpp | 4 ++-- shen.cpp | 44 ++++++++++++++++++++++++++++++++++++++++++++ shen.h | 33 +++++++++++++++++++++++++++++++++ solow.cpp | 2 +- summarycommand.cpp | 4 ++++ validcalculator.cpp | 2 ++ 9 files changed, 108 insertions(+), 13 deletions(-) create mode 100644 shen.cpp create mode 100644 shen.h diff --git a/boneh.cpp b/boneh.cpp index 7b52683..5174f75 100644 --- a/boneh.cpp +++ b/boneh.cpp @@ -11,10 +11,14 @@ #include /***********************************************************************/ +//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 e429dff..d317a06 100644 --- 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 e843641..0441f26 100644 --- 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 88c5a13..8b8f9e3 100644 --- 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 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 + + +/***********************************************************************/ +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 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) {return data;}; +private: + int m; +}; + + +/***********************************************************************/ + +#endif diff --git a/solow.cpp b/solow.cpp index 3a15fdf..8c3d155 100644 --- 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 d0124d8..b649582 100644 --- 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 7a07fd1..edcf8c7 100644 --- 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) { -- 2.39.2