5 // Created by SarahsWork on 1/3/14.
6 // Copyright (c) 2014 Schloss Lab. All rights reserved.
9 #include "shannonrange.h"
11 /***********************************************************************/
13 EstOutput RangeShannon::getValues(SAbundVector* rank){
17 double commSize = 1e20;
18 double sampleSize = rank->getNumSeqs();
22 for (int i = 1; i <=rank->getMaxRank(); i++) {
23 int abund = rank->get(i);
26 freqy.push_back(abund);
30 double aux = ceil(pow((sampleSize+1), (1/(double)3)));
31 double est0 = max(freqy[0]+1, aux);
36 for (int i = 0; i < freqx.size()-1; i++) {
38 if (m->control_pressed) { break; }
40 if (freqx[i+1] == freqx[i]+1) { numr = max(freqy[i+1]+1, aux); }
43 denr = max(freqy[i], aux);
44 ests.push_back((freqx[i]+1)*numr/(double)denr);
47 denr = max(freqy[freqy.size()-1], aux);
48 ests.push_back((freqx[freqx.size()-1]+1)*numr/(double)denr);
51 for (int i = 0; i < freqy.size(); i++) { sum += (ests[i]*freqy[i]); }
52 double nfac = est0 + sum;
55 for (int i = 0; i < ests.size(); i++) { ests[i] /= nfac; }
57 double abunup = 1 / commSize;
58 double nbrup = est0 / abunup;
59 double abunlow = ests[0];
60 double nbrlow = est0 / abunlow;
64 for (int i = 0; i < freqy.size(); i++) {
65 if (m->control_pressed) { break; }
66 sum += (freqy[i] * ests[i] * log(ests[i]));
69 data[1] = exp(data[0]+nbrlow*(-abunlow*log(abunlow)));
70 data[2] = exp(data[0]+nbrup*(-abunup*log(abunup)));
72 for (int i = 0; i < freqy.size(); i++) {
73 if (m->control_pressed) { break; }
74 data[0] += (freqy[i] * (pow(ests[i],alpha)));
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)));
80 //this calc has no data[0], just a lower and upper estimate. set data[0] to lower estimate.
82 if (data[1] > data[2]) { data[1] = data[2]; data[2] = data[0]; }
85 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
90 m->errorOut(e, "RangeShannon", "getValues");
94 /***********************************************************************/