boneh.cpp

   1 /*
   2  *  boneh.cpp
   3  *  Mothur
   4  *
   5  *  Created by Thomas Ryabin on 5/13/09.
   6  *  Copyright 2009 __MyCompanyName__. All rights reserved.
   7  *
   8  */
   9
  10 #include "boneh.h"
  11 #include <math.h>
  12
  13 /***********************************************************************/
  14 double Boneh::getV(double f1, double n, double rs) {
  15
  16         //cout << "f1 = " << f1 << "\nn = " << n << "\nrs = " << rs << "\n\n";
  17
  18         double accuracy = .0001;
  19         double v = 100000.0;
  20         double step = v/2;
  21         double ls = v * (1 - pow((1 - f1/(n*v)), n));
  22
  23         //cout << "ls = " << ls << "\n";
  24
  25         while(abs(ls - rs) > accuracy) {
  26                 if(ls > rs)
  27                         v -= step;
  28                 else
  29                         v += step;
  30
  31                 ls = v* (1 - pow((1 - f1/(n*v)), n));
  32                 step /= 2;
  33
  34         //      cout << "ls = " << ls << "\n";
  35         }
  36
  37         return v;
  38 }
  39
  40 /***********************************************************************/
  41 EstOutput Boneh::getValues(SAbundVector* rank){
  42
  43         try {
  44                 data.resize(1,0);
  45
  46                 bool valid = false;
  47                 double sum = 0;
  48                 double n = (double)rank->size() - 1;
  49                 double f1 = (double)rank->get(1);
  50
  51                 for(int i = 1; i < rank->size(); i++)
  52                         sum += (double)rank->get(i) * exp(-i);
  53
  54                 if(rank->get(1) > sum)
  55                         valid = true;
  56
  57                 sum = 0;
  58                 if(valid) {
  59                         for(int j = 1; j < rank->size(); j++)
  60                                 sum += rank->get(j) * pow((1 - (double)j / n), n);
  61
  62                         double v = getV(f1, n, sum);
  63
  64                 //      cout << "v = " << v << "\n";
  65
  66                         sum = 0;
  67                         for(int j = 1; j < rank->size(); j++)
  68                                 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));
  69
  70                 }
  71
  72                 data[0] = sum;
  73
  74                 return data;
  75         }
  76         catch(exception& e) {
  77                 cout << "Standard Error: " << e.what() << " has occurred in the Coverage class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  78                 exit(1);
  79         }
  80         catch(...) {
  81                 cout << "An unknown error has occurred in the Coverage class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  82                 exit(1);
  83         }
  84 };
  85
  86
  87 /***********************************************************************/