geom.cpp

   1 /*
   2  *  geom.cpp
   3  *  Mothur
   4  *
   5  *  Created by Thomas Ryabin on 2/23/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "geom.h"
  11
  12 /***********************************************************************/
  13
  14 double Geom::kEq(double k, double spec){
  15         return k/(1-k)*pow(1-k, spec)/(1-pow(1-k, spec));
  16 }
  17
  18 RAbundVector Geom::getRAbundVector(SAbundVector* rank){
  19                 vector <int> rData;
  20                 int mr = 1;
  21                 int nb = 0;
  22                 int ns = 0;
  23
  24                 for(int i = rank->size()-1; i > 0; i--) {
  25                         int cur = rank->get(i);
  26                         if(mr == 1 && cur > 0)
  27                                 mr = i;
  28                         nb += cur;
  29                         ns += i*cur;
  30                         for(int j = 0; j < cur; j++)
  31                                 rData.push_back(i);
  32                 }
  33
  34                 RAbundVector rav = RAbundVector(rData, mr, nb, ns);
  35                 rav.setLabel(rank->getLabel());
  36                 return rav;
  37 }
  38
  39 /***********************************************************************************/
  40
  41 /***********************************************************************************/
  42 EstOutput Geom::getValues(SAbundVector* rank){
  43         try {
  44                 data.resize(3,0);
  45
  46                 rdata = getRAbundVector(rank);
  47                 double numInd = rdata.getNumSeqs();
  48                 double numSpec = rdata.getNumBins();
  49                 double min = rdata.get(rdata.size()-1);
  50                 double k = .5;
  51                 double step = .49999;
  52
  53                 while(fabs(min - numInd*kEq(k, (double)numSpec)) > .0001) { //This uses a binary search to find the value of k.
  54                         if(numInd*kEq(k, numSpec) > min)
  55                                 k += step;
  56                         else
  57                                 k -= step;
  58                         step /= 2;
  59                 }
  60
  61                 double cK = 1/(1-pow(1-k, numSpec));
  62                 double sumExp = 0;
  63                 double sumObs = 0;
  64                 double maxDiff = 0;
  65
  66                 for(int i = 0; i < numSpec; i++)
  67                 {
  68                         sumObs += rdata.get(i);
  69                         sumExp += numInd*cK*k*pow(1-k, i);
  70
  71                         double diff = fabs(sumObs-sumExp);
  72                         if(diff > maxDiff)      {       maxDiff = diff;         }
  73
  74                 }
  75
  76
  77                 data[0] = maxDiff/numInd;
  78                 data[1] = 0.886/sqrt(numSpec);
  79                 data[2] = 1.031/sqrt(numSpec);
  80
  81
  82                 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
  83                 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
  84                 if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
  85
  86                 return data;
  87         }
  88         catch(exception& e) {
  89                 m->errorOut(e, "Geom", "getValues");
  90                 exit(1);
  91         }
  92 }
  93
  94 /***********************************************************************/
  95
  96
  97
  98
  99
 100
 101