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