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 //This solves for the value of 'v' using a binary search.
  15 double Boneh::getV(double f1, double n, double rs) {
  16
  17         //cout << "f1 = " << f1 << "\nn = " << n << "\nrs = " << rs << "\n\n";
  18
  19         if(rs == 0)
  20                 return 0;
  21
  22         double accuracy = .0001;
  23         double v = 100000.0;
  24         double step = v/2;
  25         double ls = v * (1 - pow((1 - f1/(n*v)), n));
  26
  27         //cout << "ls = " << ls << "\n";
  28
  29         while(abs(ls - rs) > accuracy) {
  30                 if(ls > rs)
  31                         v -= step;
  32                 else
  33                         v += step;
  34
  35                 ls = v* (1 - pow((1 - f1/(n*v)), n));
  36                 step /= 2;
  37
  38                 //cout << "ls = " << ls << "\n";
  39         }
  40
  41         return v;
  42 }
  43
  44 /***********************************************************************/
  45 EstOutput Boneh::getValues(SAbundVector* rank){
  46
  47         try {
  48                 data.resize(1,0);
  49
  50                 bool valid = false;
  51                 double sum = 0;
  52                 double n = (double)rank->getNumSeqs();
  53                 double f1 = (double)rank->get(1);
  54
  55                 for(int i = 1; i < rank->size(); i++)
  56                         sum += (double)rank->get(i) * exp(-i);
  57
  58                 if(rank->get(1) > sum)
  59                         valid = true;
  60
  61                 sum = 0;
  62                 if(valid) {
  63                         for(int j = 1; j < rank->size(); j++)
  64                                 sum += rank->get(j) * pow((1 - (double)j / n), n);
  65
  66                         double v = getV(f1, n, sum);
  67
  68                         //cout << "v = " << v << "\n";
  69
  70
  71                         sum = 0;
  72                         for(int j = 1; j < rank->size(); j++) {
  73                                 double Xi = 0; //I didn't know what this was, simply replace the 0
  74                                                            //with the appropriate expression for the boneh calculator
  75                                                            //to work.
  76                                 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));
  77                         }
  78                 }
  79
  80                 data[0] = sum;
  81
  82                 return data;
  83         }
  84         catch(exception& e) {
  85                 cout << "Standard Error: " << e.what() << " has occurred in the Coverage class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  86                 exit(1);
  87         }
  88         catch(...) {
  89                 cout << "An unknown error has occurred in the Coverage class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  90                 exit(1);
  91         }
  92 };
  93
  94
  95 /***********************************************************************/