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
|
//
// lnabundace.cpp
// Mothur
//
// Created by Sarah Westcott on 5/8/19.
// Copyright © 2019 Schloss Lab. All rights reserved.
//
#include "lnabundance.hpp"
/***********************************************************************/
LNAbundance::LNAbundance() : DiversityCalculator(true) {}
/***********************************************************************/
vector<double> LNAbundance::getValues(int maxRank, vector<mcmcSample>& sampling) { //rank->getMaxRank(); //nMax
try {
maxRank = floor(pow(2.0,ceil(log((double) maxRank)/log(2.0)) + 2.0) + 1.0e-7);
results.resize(maxRank, 0.0);
int nSamples = sampling.size();
if (nSamples == 0) { return results; }
#ifdef USE_GSL
gsl_set_error_handler_off();
DiversityUtils dutils("lnabund");
for(int i = 0; i < sampling.size(); i++) {
if (m->getControl_pressed()) { break; }
for (int j = 1; j <= maxRank; j++) {
int nA = j;
double dLog = 0.0, dP = 0.0;
if(nA < 100){ //MAX_QUAD
dLog = dutils.logLikelihoodQuad(nA, sampling[i].alpha, sampling[i].beta); //nA, dMDash, dV
}
else{
dLog = dutils.logLikelihoodRampal(nA, sampling[i].alpha, sampling[i].beta); //nA, dMDash, dV
}
dP = exp(dLog);
results[j - 1] += dP*sampling[i].ns;
}
}
for (int i = 1; i<=maxRank; i++) {
results[i-1] /= (double)nSamples;
if (isnan(results[i-1]) || isinf(results[i-1])) { results[i-1] = 0.0; }
}
#endif
return results;
}
catch(exception& e) {
m->errorOut(e, "LNAbundance", "getValues");
exit(1);
}
}
/***********************************************************************/
|
