shannonrange.cpp

   1 //
   2 //  shannonrange.cpp
   3 //  Mothur
   4 //
   5 //  Created by SarahsWork on 1/3/14.
   6 //  Copyright (c) 2014 Schloss Lab. All rights reserved.
   7 //
   8
   9 #include "shannonrange.h"
  10
  11 /***********************************************************************/
  12
  13 EstOutput RangeShannon::getValues(SAbundVector* rank){
  14         try {
  15         data.resize(3,0);
  16
  17         double commSize = 1e20;
  18         double sampleSize = rank->getNumSeqs();
  19
  20         vector<int> freqx;
  21         vector<int> freqy;
  22         for (int i = 1; i <=rank->getMaxRank(); i++) {
  23             int abund = rank->get(i);
  24             if (abund != 0) {
  25                 freqx.push_back(i);
  26                 freqy.push_back(abund);
  27             }
  28         }
  29
  30         double aux = ceil(pow((sampleSize+1), (1/(double)3)));
  31         double est0 = max(freqy[0]+1, aux);
  32
  33         vector<double> ests;
  34         double numr = 0.0;
  35         double denr = 0.0;
  36         for (int i = 0; i < freqx.size()-1; i++) {
  37
  38             if (m->control_pressed) { break; }
  39
  40             if (freqx[i+1] == freqx[i]+1)   { numr = max(freqy[i+1]+1, aux);    }
  41             else                            { numr = aux;                       }
  42
  43             denr = max(freqy[i], aux);
  44             ests.push_back((freqx[i]+1)*numr/(double)denr);
  45         }
  46         numr = aux;
  47         denr = max(freqy[freqy.size()-1], aux);
  48         ests.push_back((freqx[freqx.size()-1]+1)*numr/(double)denr);
  49
  50         double sum = 0.0;
  51         for (int i = 0; i < freqy.size(); i++) {  sum += (ests[i]*freqy[i]); }
  52         double nfac = est0 + sum;
  53         est0 /= nfac;
  54
  55         for (int i = 0; i < ests.size(); i++) {  ests[i] /= nfac;   }
  56
  57         double abunup = 1 / commSize;
  58         double nbrup = est0 / abunup;
  59         double abunlow = ests[0];
  60         double nbrlow = est0 / abunlow;
  61
  62         if (alpha == 1) {
  63             double sum = 0.0;
  64             for (int i = 0; i < freqy.size(); i++) {
  65                 if (m->control_pressed) { break; }
  66                 sum += (freqy[i] * ests[i] * log(ests[i]));
  67             }
  68             data[0] = -sum;
  69             data[1] = exp(data[0]+nbrlow*(-abunlow*log(abunlow)));
  70             data[2] = exp(data[0]+nbrup*(-abunup*log(abunup)));
  71         }else {
  72             for (int i = 0; i < freqy.size(); i++) {
  73                 if (m->control_pressed) { break; }
  74                 data[0] += (freqy[i] * (pow(ests[i],alpha)));
  75             }
  76             data[1] = pow(data[0]+nbrup*pow(abunup,alpha), (1/(1-alpha)));
  77             data[2] = pow(data[0]+nbrlow*pow(abunlow,alpha), (1/(1-alpha)));
  78         }
  79
  80         //this calc has no data[0], just a lower and upper estimate. set data[0] to lower estimate.
  81         data[0] = data[1];
  82         if (data[1] > data[2]) { data[1] = data[2]; data[2] = data[0]; }
  83         data[0] = data[1];
  84
  85         if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
  86
  87                 return data;
  88         }
  89         catch(exception& e) {
  90                 m->errorOut(e, "RangeShannon", "getValues");
  91                 exit(1);
  92         }
  93 }
  94 /***********************************************************************/