boneh.cpp

   1 /*
   2  *  boneh.cpp
   3  *  Mothur
   4  *
   5  *  Created by Thomas Ryabin on 5/13/09.
   6  *  Copyright 2009Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "boneh.h"
  11 #include <math.h>
  12
  13 /***********************************************************************/
  14
  15 //This solves for the value of 'v' using a binary search.
  16 double Boneh::getV(double f1, double n, double rs) {
  17
  18         if(rs == 0)
  19                 return 0;
  20
  21         double accuracy = .0001;
  22         double v = 100000.0;
  23         double step = v/2;
  24         double ls = v * (1 - pow((1 - f1/(n*v)), n));
  25
  26         while(abs(ls - rs) > accuracy) {
  27                 if(ls > rs)     {       v -= step;      }
  28                 else            {       v += step;      }
  29
  30                 ls = v * (1 - pow((1 - f1/(n * v)), n));
  31                 step /= 2;
  32         }
  33
  34         return v;
  35 }
  36
  37 /***********************************************************************/
  38 EstOutput Boneh::getValues(SAbundVector* sabund){
  39
  40         try {
  41                 data.resize(1,0);
  42
  43
  44                 bool valid = false;
  45                 double sum = 0;
  46                 double n = (double)sabund->getNumSeqs();
  47                 if(f==0){       f=n;    }
  48
  49                 double f1 = (double)sabund->get(1);
  50
  51                 for(int i = 1; i < sabund->size(); i++){
  52                         sum += (double)sabund->get(i) * exp(-i);
  53                 }
  54
  55                 if(sabund->get(1) > sum)
  56                         valid = true;
  57
  58                 sum = 0;
  59                 if(valid) {
  60                         for(int j = 1; j < sabund->size(); j++){
  61                                 sum += sabund->get(j) * pow((1 - (double)j / n), n);
  62                         }
  63
  64                         double v = getV(f1, n, sum);
  65
  66                         sum = 0;
  67                         for(int j = 1; j < sabund->size(); j++) {
  68                                 for (int i = 0; i < sabund->get(j); i++) {
  69                                         sum += pow(1 - j / n, n) * (1 - pow(1 - j / n, f));
  70                                 }
  71                         }
  72                         sum +=  v * pow(1 - f1/(n*v), n) * (1 - pow(1 - f1/(n*v), f));
  73                 }
  74
  75                 data[0] = sum;
  76
  77                 return data;
  78         }
  79         catch(exception& e) {
  80                 m->errorOut(e, "Boneh", "getValues");
  81                 exit(1);
  82         }
  83 }
  84
  85
  86 /***********************************************************************/