1 //
2 //  diversityutils.cpp
3 //  Mothur
4 //
5 //  Created by Sarah Westcott on 4/11/19.
6 //  Copyright © 2019 Schloss Lab. All rights reserved.
7 //
8 
9 #include "diversityutils.hpp"
10 /***********************************************************************/
f1_2(double x,void * pvParams)11 double f1_2(double x, void *pvParams)
12 {
13     t_LSParams *ptLSParams = (t_LSParams *) pvParams;
14     double dMDash = ptLSParams->dMDash, dV = ptLSParams->dV, dNu = ptLSParams->dNu;
15     int n = ptLSParams->n;
16     double t = ((x - dMDash)*(x - dMDash))/dV;
17     double dExp  = x*((double) n) - exp(x);
18     double dF  = pow(1.0 + t/dNu, -0.5*(dNu + 1.0));
19 
20     return exp(dExp)*dF;
21 }
22 /***********************************************************************/
f1Log(double x,void * pvParams)23 double f1Log(double x, void *pvParams)
24 {
25     t_LNParams *ptLNParams = (t_LNParams *) pvParams;
26     double dMDash = ptLNParams->dMDash, dV = ptLNParams->dV;
27     int n = ptLNParams->n;
28     double dTemp = (x - dMDash);
29     double dExp  = x*((double) n) - exp(x) - 0.5*((dTemp*dTemp)/dV);
30     double dRet  = exp(dExp);
31 
32     return dRet;
33 }
34 /***********************************************************************/
f1(double x,void * pvParams)35 double f1(double x, void *pvParams)
36 {
37     t_LNParams *ptLNParams = (t_LNParams *) pvParams;
38     double dMDash = ptLNParams->dMDash, dV = ptLNParams->dV, n = ptLNParams->n;
39     double dTemp = (x - dMDash);
40     double dExp  = x*((double) n) - exp(x) - 0.5*((dTemp*dTemp)/dV);
41     double dRet  = exp(dExp);
42 
43     return dRet;
44 }
45 
46 /***********************************************************************/
derivExponent(double x,void * pvParams)47 double derivExponent(double x, void *pvParams)
48 {
49     t_LNParams *ptLNParams = (t_LNParams *) pvParams;
50     double dMDash = ptLNParams->dMDash, dV = ptLNParams->dV, n = ptLNParams->n;
51     double dTemp = (x - dMDash)/dV, dRet = 0.0;
52 
53     dRet = ((double) n) - exp(x) - dTemp;
54 
55     return dRet;
56 }
57 /***********************************************************************/
loadAbundance(t_Data * ptData,SAbundVector * rank)58 void DiversityUtils::loadAbundance(t_Data *ptData, SAbundVector* rank) {
59     try {
60         int nNA = 0; int  nL = 0, nJ = 0;
61         int maxRank = rank->getMaxRank();
62 
63         for(int i = 1; i <= maxRank; i++){ if (rank->get(i) != 0) { nNA++; } }
64 
65         int **aanAbund = NULL;
66         aanAbund = (int **) malloc(nNA*sizeof(int*));
67 
68         int count = 0;
69         for(int i = 1; i <= maxRank; i++){
70 
71             if (m->getControl_pressed()) { break; }
72 
73             if (rank->get(i) != 0) {
74                 aanAbund[count] = (int *) malloc(sizeof(int)*2);
75 
76                 int nA = i;
77                 int nC = rank->get(i);
78 
79                 nL += nC;
80                 nJ += nC*nA;
81 
82                 aanAbund[count][0]  = nA;
83                 aanAbund[count][1]  = nC;
84 
85                 count++;
86             }
87 
88         }
89 
90         ptData->nJ          = nJ;
91         ptData->nL          = nL;
92         ptData->aanAbund    = aanAbund;
93         ptData->nNA         = nNA;
94     }
95     catch(exception& e) {
96         m->errorOut(e, "DiversityUtils", "loadAbundance");
97         exit(1);
98     }
99 }
100 /***********************************************************************/
freeAbundance(t_Data * ptData)101 void DiversityUtils::freeAbundance(t_Data *ptData) {
102     try {
103         for(int i = 0; i < ptData->nNA; i++){ free(ptData->aanAbund[i]); }
104         free(ptData->aanAbund);
105     }
106     catch(exception& e) {
107         m->errorOut(e, "DiversityUtils", "freeAbundance");
108         exit(1);
109     }
110 }
111 /***********************************************************************/
chao(t_Data * ptData)112 double DiversityUtils::chao(t_Data *ptData) {
113     try {
114         double n1 = 0.0, n2 = 0.0;
115         int **aanAbund = ptData->aanAbund;
116 
117         if(aanAbund[0][0] == 1 && aanAbund[1][0] == 2){
118 
119             n1 = (double) aanAbund[0][1];
120             n2 = (double) aanAbund[1][1];
121 
122             return ((double) ptData->nL) + 0.5*((n1*n1)/n2);
123         }
124         else{ return -1.0; }
125     }
126     catch(exception& e) {
127         m->errorOut(e, "DiversityUtils", "chao");
128         exit(1);
129     }
130 }
131 /***********************************************************************/
fX(double x,double dA,double dB,double dNDash)132 double DiversityUtils::fX(double x, double dA, double dB, double dNDash){
133     try {
134         double dTemp1 = (dA*(x - dB)*(x - dB))/x;
135 
136         return log(x) - (1.0/dNDash)*(x + dTemp1);
137     }
138     catch(exception& e) {
139         m->errorOut(e, "DiversityUtils", "fX");
140         exit(1);
141     }
142 }
143 /***********************************************************************/
f2X(double x,double dA,double dB,double dNDash)144 double DiversityUtils::f2X(double x, double dA, double dB, double dNDash) {
145     try {
146         double dRet = 0.0, dTemp = 2.0*dA*dB*dB;
147 
148         dRet = (1.0/(x*x))*(1.0 + (1.0/dNDash)*(dTemp/x));
149 
150         return -dRet;
151     }
152     catch(exception& e) {
153         m->errorOut(e, "DiversityUtils", "f2X");
154         exit(1);
155     }
156 }
157  #ifdef USE_GSL
158 
159 /***********************************************************************/
fMu_sirarefaction(double x,void * pvParams)160 double fMu_sirarefaction(double x, void* pvParams)
161 {
162     DiversityUtils dutils("sirarefaction");
163 
164     t_LSParams* ptSIParams = (t_LSParams *) pvParams;
165     double dAlphaDD = ptSIParams->dMDash*sqrt(x);
166     double dBetaDD  = ptSIParams->dV*x;
167     double dLogP0   = dutils.logLikelihood(0, dAlphaDD, dBetaDD, ptSIParams->dNu);
168 
169     return (1.0 - exp(dLogP0)) - ptSIParams->dC;
170 }
171 /***********************************************************************/
fMu_igrarefaction(double x,void * pvParams)172 double fMu_igrarefaction(double x, void* pvParams)
173 {
174     t_IGParams* ptIGParams = (t_IGParams *) pvParams;
175     // double tx = x / 1667.0;
176 
177     DiversityUtils dutils("igrarefaction");
178 
179     double dAlphaDD = ptIGParams->dAlpha*sqrt(x);
180     double dBetaDD  = ptIGParams->dBeta*x;
181     double dLogP0   = dutils.logLikelihood(0, dAlphaDD, dBetaDD);
182 
183     // printf("dAlpha %f dBeta %f x %f dAlphaDD %f  dBetaDD %f dLofP0 %f", ptIGParams->dAlpha, ptIGParams->dBeta, x, dAlphaDD, dBetaDD, dLogP0);
184 
185     return (1.0 - exp(dLogP0)) - ptIGParams->dC;
186 }
187 /***********************************************************************/
fMu_lsrarefaction(double x,void * pvParams)188 double fMu_lsrarefaction(double x, void* pvParams)
189 {
190 
191     DiversityUtils dutils("lsrarefaction");
192 
193     t_LSParams* ptLSParams = (t_LSParams *) pvParams;
194     double dMDD            = ptLSParams->dMDash + x;
195     double dLogP0          = dutils.logLikelihoodQuad(0, dMDD, ptLSParams->dV, ptLSParams->dNu);
196 
197     return (1.0 - exp(dLogP0)) - ptLSParams->dC;
198 }
199 /***********************************************************************/
fMu_lnrarefaction(double x,void * pvParams)200 double fMu_lnrarefaction(double x, void* pvParams)
201 {
202     t_IGParams* ptIGParams = (t_IGParams *) pvParams;
203 
204     DiversityUtils dutils("lnrarefaction");
205 
206     double dMDD = ptIGParams->dAlpha + x;
207     double dLogP0 = dutils.logLikelihoodQuad(0, dMDD, ptIGParams->dBeta);
208 
209     return (1.0 - exp(dLogP0)) - ptIGParams->dC;
210 }
211 /***********************************************************************/
calcMu(void * pvParams)212 double DiversityUtils::calcMu(void *pvParams){
213     try {
214         double dLogMu = 0.0;
215 
216         if (method == "lnrarefaction") {
217             t_IGParams* ptIGParams = (t_IGParams *) pvParams;
218             solveF(0, 1.0e7, ptIGParams, 1.0e-7, &dLogMu);
219             return exp(dLogMu);
220         }else if ((method == "igrarefaction") || (method == "sirarefaction")) {
221             t_IGParams* ptIGParams = (t_IGParams *) pvParams;
222             solveF(1.0, 1.0e10, ptIGParams, 1.0e-7, &dLogMu);
223             return dLogMu;
224         }else if (method == "lsrarefaction") {
225             t_LSParams *ptLSParams = (t_LSParams *) pvParams;
226             solveF(0, 1.0e7, ptLSParams, 1.0e-7, &dLogMu);
227             return exp(dLogMu);
228         }
229 
230         return dLogMu;
231     }
232     catch(exception& e) {
233         m->errorOut(e, "DiversityUtils", "calcMu");
234         exit(1);
235     }
236 }
237 //***********************************************************************/
logStirlingsGamma(double dZ)238 double DiversityUtils::logStirlingsGamma(double dZ) { return 0.5*log(2.0*M_PI) + (dZ - 0.5)*log(dZ) - dZ; }
239 /***********************************************************************/
logLikelihoodQuad(int n,double dMDash,double dV)240 double DiversityUtils::logLikelihoodQuad(int n, double dMDash, double dV){
241     try {
242         gsl_integration_workspace *ptGSLWS =
243         gsl_integration_workspace_alloc(1000);
244         double dLogFac1   = 0.0, dLogFacN  = 0.0;
245         double dResult = 0.0, dError = 0.0, dPrecision = 0.0;
246         gsl_function tGSLF;
247         double dEst = dMDash + ((double)n)*dV, dA = 0.0, dB = 0.0;
248 
249         t_LNParams tLNParams; tLNParams.n = n; tLNParams.dMDash = dMDash; tLNParams.dV = dV;
250 
251         tGSLF.function = &f1;
252         if (method == "metrols") {  tGSLF.function = &f1Log;  }
253         tGSLF.params   = (void *) &tLNParams;
254 
255         dLogFac1 = log(2.0*M_PI*dV);
256 
257         if(n < 50){
258             dLogFacN = gsl_sf_fact(n);
259             dLogFacN = log(dLogFacN);
260         }
261         else{
262             dLogFacN = gsl_sf_lngamma(((double) n) + 1.0);
263         }
264 
265         if(dEst > dV){
266             double dMax = 0.0;
267             double dUpper = (((double) n) + (dMDash/dV) - 1.0)/(1.0 + 1.0/dV);
268             double dVar   = 0.0;
269 
270             if ((method == "metroln") || (method == "metrols")) { if(fabs(dUpper) > 1.0e-7){  solveF(0.0, dUpper, (void *) &tLNParams, 1.0e-5, &dMax); } }
271             else {  solveF(0.0, dUpper, derivExponent, (void *) &tLNParams, 1.0e-5, &dMax);  } //lnabund, lnshift, lnrarefact
272 
273             dVar = sqrt(1.0/((1.0/dV) + exp(dMax)));
274 
275             dA = dMax - V_MULT*dVar; dB = dMax + V_MULT*dVar;
276         }
277         else{
278             double dMax = 0.0;
279             double dLower = dEst - dV;
280             double dUpper = (((double) n) + (dMDash/dV) - 1.0)/(1.0 + 1.0/dV);
281             double dVar   = 0.0;
282 
283             if ((method == "metroln") || (method == "metrols")) {
284                 if(fabs(dUpper - dLower) > 1.0e-7){ solveF(dLower, dUpper, (void *) &tLNParams, 1.0e-5, &dMax); }
285                 else{ dMax = 0.5*(dLower + dUpper); }
286             }else {  solveF(dLower, dUpper, derivExponent, (void *) &tLNParams, 1.0e-5, &dMax); } //lnabund, lnshift, lnrarefact
287 
288             dVar = sqrt(1.0/((1.0/dV) + exp(dMax)));
289 
290             dA = dMax - V_MULT*dVar; dB = dMax + V_MULT*dVar;
291         }
292 
293         if(n < 10)  {  dPrecision = HI_PRECISION;   }
294         else        {  dPrecision = LO_PRECISION;   }
295 
296         gsl_integration_qag(&tGSLF, dA, dB, dPrecision, 0.0, 1000, GSL_INTEG_GAUSS61, ptGSLWS, &dResult, &dError);
297 
298         gsl_integration_workspace_free(ptGSLWS);
299 
300         return log(dResult) - dLogFacN -0.5*dLogFac1;
301     }
302     catch(exception& e) {
303         m->errorOut(e, "DiversityUtils", "logLikelihoodQuad");
304         exit(1);
305     }
306 }
307 //***********************************************************************/
logLikelihoodQuad(int n,double dMDash,double dV,double dNu)308 double DiversityUtils::logLikelihoodQuad(int n, double dMDash, double dV, double dNu){
309     try {
310         gsl_integration_workspace *ptGSLWS =
311         gsl_integration_workspace_alloc(1000);
312         double dLogFac1   = 0.0, dLogFacN  = 0.0;
313         double dN = (double) n, dResult = 0.0, dError = 0.0, dPrecision = 0.0;
314         gsl_function tGSLF;
315         t_LSParams tLSParams;
316         double dA = 0.0, dB = 0.0;
317 
318         tLSParams.n = n; tLSParams.dMDash = dMDash; tLSParams.dV = dV; tLSParams.dNu = dNu;
319 
320         tGSLF.function = &f1_2;
321         tGSLF.params   = (void *) &tLSParams;
322 
323         if(dNu < 100){ //MAX_MU_GAMMA
324             dLogFac1 = gsl_sf_lngamma(0.5*(dNu + 1.0)) - gsl_sf_lngamma(0.5*dNu) - 0.5*log(M_PI*dNu);
325         }
326         else{
327             dLogFac1 = 0.5*dNu*(log(0.5*(dNu + 1.0)) - log(0.5*dNu)) -0.5*log(2.0*M_PI) - 0.5;
328         }
329 
330         if(n < 50){
331             dLogFacN = gsl_sf_fact(n);
332             dLogFacN = log(dLogFacN);
333         }
334         else if(n < 100){
335             dLogFacN = gsl_sf_lngamma(dN + 1.0);
336         }
337         else{
338             dLogFacN = logStirlingsGamma(dN + 1.0);
339         }
340 
341         dA = -100.0; dB = 100.0;
342 
343         if(n < 10)  { dPrecision = HI_PRECISION; }
344         else        { dPrecision = LO_PRECISION; }
345 
346         gsl_integration_qag(&tGSLF, dA, dB, dPrecision, 0.0, 1000, GSL_INTEG_GAUSS61, ptGSLWS, &dResult, &dError);
347 
348         //printf("%f %f\n", dResult, dError);
349 
350         gsl_integration_workspace_free(ptGSLWS);
351 
352         return log(dResult) - dLogFacN + dLogFac1 - 0.5*log(dV);
353     }
354     catch(exception& e) {
355         m->errorOut(e, "DiversityUtils", "logLikelihoodQuad");
356         exit(1);
357     }
358 }
359 /***********************************************************************/
logLikelihoodRampal(int n,double dMDash,double dV)360 double DiversityUtils::logLikelihoodRampal(int n, double dMDash, double dV){
361     try {
362         double dN = (double) n;
363         double dLogLik = 0.0, dTemp = gsl_pow_int(log(dN) - dMDash,2), dTemp3 = gsl_pow_int(log(dN) - dMDash,3);
364 
365         dLogLik = -0.5*log(2.0*M_PI*dV) - log(dN) - (dTemp/(2.0*dV));
366 
367         dLogLik += log(1.0 + 1.0/(2.0*dN*dV)*(dTemp/dV + log(dN) - dMDash - 1.0)
368                        + 1.0/(6.0*dN*dN*dV*dV*dV)*(3.0*dV*dV - (3.0*dV - 2.0*dV*dV)*(dMDash - log(dN))
369                                                    - 3.0*dV*dTemp + dTemp3));
370 
371         return dLogLik;
372     }
373     catch(exception& e) {
374         m->errorOut(e, "DiversityUtils", "logLikelihoodRampal");
375         exit(1);
376     }
377 }
378 //***********************************************************************/
logLikelihoodRampal(int n,double dMDash,double dV,double dNu)379 double DiversityUtils::logLikelihoodRampal(int n, double dMDash, double dV, double dNu){
380     try {
381         double dGamma = 0.5*(dNu + 1.0), dN = (double) n, dRN = 1.0/dN, dRSV = 1.0/(sqrt(dV)*sqrt(dNu));
382         double dZ = (log(dN) - dMDash)*dRSV;
383         double dDZDX = dRN*dRSV, dDZDX2 = -dRN*dRN*dRSV;
384         double dF = (1.0 + dZ*dZ);
385         double dA = 0.0, dB = 0.0, dTemp = 0.0;
386         double dLogFac1 = 0.0;
387 
388         if(dNu < 100){ //MAX_MU_GAMMA
389             dLogFac1 = gsl_sf_lngamma(0.5*(dNu + 1.0)) - gsl_sf_lngamma(0.5*dNu) - 0.5*log(M_PI*dNu);
390         }
391         else{
392             dLogFac1 = 0.5*dNu*(log(0.5*(dNu + 1.0)) - log(0.5*dNu)) -0.5*log(2.0*M_PI) - 0.5;
393         }
394 
395         dA = 4.0*dZ*dZ*dDZDX*dDZDX*dGamma*(dGamma + 1.0);
396         dA /= dF*dF;
397 
398         dB = -2.0*dGamma*(dDZDX*dDZDX + dZ*dDZDX2);
399         dB /= dF;
400 
401         dTemp = dRN + dA + dB;
402 
403         return -dGamma*log(dF) + log(dTemp) + dLogFac1 - 0.5*log(dV);
404     }
405     catch(exception& e) {
406         m->errorOut(e, "DiversityUtils", "logLikelihoodRampal");
407         exit(1);
408     }
409 }
410 //***********************************************************************/
solveF(double x_lo,double x_hi,double (* f)(double,void *),void * params,double tol,double * xsolve)411 int DiversityUtils::solveF(double x_lo, double x_hi, double (*f)(double, void*), void* params, double tol, double *xsolve){
412     try {
413         int status, iter = 0, max_iter = 100;
414         const gsl_root_fsolver_type *T;
415         gsl_root_fsolver *s;
416         double r = 0;
417         gsl_function F;
418 
419         F.function = f;
420         F.params = params;
421 
422         T = gsl_root_fsolver_brent;
423         s = gsl_root_fsolver_alloc (T);
424         gsl_root_fsolver_set (s, &F, x_lo, x_hi);
425 
426         do{
427             iter++;
428 
429             if (m->getControl_pressed()) { break; }
430 
431             status = gsl_root_fsolver_iterate (s);
432             r = gsl_root_fsolver_root (s);
433             x_lo = gsl_root_fsolver_x_lower (s);
434             x_hi = gsl_root_fsolver_x_upper (s);
435 
436             status = gsl_root_test_interval (x_lo, x_hi, 0, tol);
437         }
438         while (status == GSL_CONTINUE && iter < max_iter);
439 
440         (*xsolve) = gsl_root_fsolver_root (s);
441         gsl_root_fsolver_free (s);
442 
443         return status;
444     }
445     catch(exception& e) {
446         m->errorOut(e, "DiversityUtils", "solveF");
447         exit(1);
448     }
449 }
450 //***********************************************************************/
solveF(double x_lo,double x_hi,void * params,double tol,double * xsolve)451 int DiversityUtils::solveF(double x_lo, double x_hi, void* params, double tol, double *xsolve){
452     try {
453         int status, iter = 0, max_iter = 100;
454         const gsl_root_fsolver_type *T;
455         gsl_root_fsolver *s;
456         double r = 0;
457         gsl_function F;
458 
459         F.function = &derivExponent;
460         if (method == "igrarefaction") {  F.function = &fMu_igrarefaction; }
461         else if (method == "lnrarefaction") {  F.function = &fMu_lnrarefaction; }
462         else if (method == "lsrarefaction") {  F.function = &fMu_lsrarefaction; }
463         else if (method == "sirarefaction") {  F.function = &fMu_sirarefaction; }
464         F.params = params;
465 
466         T = gsl_root_fsolver_brent;
467         s = gsl_root_fsolver_alloc (T);
468         gsl_root_fsolver_set (s, &F, x_lo, x_hi);
469 
470         do{
471             iter++;
472 
473             if (m->getControl_pressed()) { break; }
474 
475             status = gsl_root_fsolver_iterate (s);
476             r = gsl_root_fsolver_root (s);
477             x_lo = gsl_root_fsolver_x_lower (s);
478             x_hi = gsl_root_fsolver_x_upper (s);
479 
480             status = gsl_root_test_interval (x_lo, x_hi, 0, tol);
481         }
482         while (status == GSL_CONTINUE && iter < max_iter);
483 
484         (*xsolve) = gsl_root_fsolver_root (s);
485         gsl_root_fsolver_free (s);
486 
487         return status;
488     }
489     catch(exception& e) {
490         m->errorOut(e, "DiversityUtils", "solveF");
491         exit(1);
492     }
493 }
494 
495 /***********************************************************************/
logLikelihood(int n,double dAlpha,double dBeta)496 double DiversityUtils::logLikelihood(int n, double dAlpha, double dBeta){
497     try {
498         double dLogFacN = 0.0;
499         bool status      = false;
500         double dRet     = 0.0;
501 
502         if(n < 50){
503             dLogFacN = gsl_sf_fact(n);
504             dLogFacN = log(dLogFacN);
505         }
506         else{
507             dLogFacN = gsl_sf_lngamma(((double) n) + 1.0);
508         }
509 
510         status = bessel(&dRet,n, dAlpha,dBeta);
511 
512         if(!status){ dRet = sd(n, dAlpha,dBeta); }
513 
514         return dRet - dLogFacN;
515     }
516     catch(exception& e) {
517         m->errorOut(e, "DiversityUtils", "logLikelihood");
518         exit(1);
519     }
520 }
521 /***********************************************************************/
logLikelihood(int n,double dAlpha,double dBeta,double dGamma)522 double DiversityUtils::logLikelihood(int n, double dAlpha, double dBeta, double dGamma){
523     try {
524         double dLogFacN = 0.0;
525         bool status      = false;
526         double dRet     = 0.0;
527 
528         if(n < 50){
529             dLogFacN = gsl_sf_fact(n);
530             dLogFacN = log(dLogFacN);
531         }
532         else{
533             dLogFacN = gsl_sf_lngamma(((double) n) + 1.0);
534         }
535 
536         status = bessel(&dRet,n, dAlpha,dBeta,dGamma);
537 
538         if(!status){ dRet = sd(n, dAlpha,dBeta,dGamma); }
539 
540         return dRet - dLogFacN;
541     }
542     catch(exception& e) {
543         m->errorOut(e, "DiversityUtils", "logLikelihood");
544         exit(1);
545     }
546 }
547 /***********************************************************************/
sd(int n,double dAlpha,double dBeta)548 double DiversityUtils::sd(int n, double dAlpha, double dBeta){
549     try {
550         double dGamma = -0.5;
551         double dA = 0.5*(-1.0 + sqrt(1.0 + (dAlpha*dAlpha)/(dBeta*dBeta)));
552         double dN = (double) n, dNDash = dN + dGamma - 1.0, dRN = 1.0/dN;
553         double dTemp1 = (0.5*dN)/(1.0 + dA), dTemp2 = 4.0*dRN*dRN*(1.0 + dA)*dA*dBeta*dBeta;
554         double dXStar = dTemp1*(1.0 + sqrt(1.0 + dTemp2));
555         double dFX = fX(dXStar, dA, dBeta, dNDash);
556         double d2FX = -dNDash*f2X(dXStar, dA, dBeta, dNDash);
557         double dLogK = 0.0, dGamma1 = dGamma;
558 
559         if(dGamma1 < 0.0){ dGamma1 *= -1.0; } //invert sign
560 
561         dLogK = gsl_sf_bessel_lnKnu(dGamma1,2.0*dA*dBeta);
562 
563         return -2.0*dA*dBeta -log(2.0) -dLogK -dGamma*log(dBeta) + dNDash*dFX + 0.5*log(2.0*M_PI) - 0.5*log(d2FX);
564     }
565     catch(exception& e) {
566         m->errorOut(e, "DiversityUtils", "sd");
567         exit(1);
568     }
569 }
570 /***********************************************************************/
sd(int n,double dAlpha,double dBeta,double dGamma)571 double DiversityUtils::sd(int n, double dAlpha, double dBeta, double dGamma){
572     try {
573         double dA = 0.5*(-1.0 + sqrt(1.0 + (dAlpha*dAlpha)/(dBeta*dBeta)));
574         double dN = (double) n, dNDash = dN + dGamma - 1.0, dRN = 1.0/dN;
575         double dTemp1 = (0.5*dN)/(1.0 + dA), dTemp2 = 4.0*dRN*dRN*(1.0 + dA)*dA*dBeta*dBeta;
576         double dXStar = dTemp1*(1.0 + sqrt(1.0 + dTemp2));
577         double dFX = fX(dXStar, dA, dBeta, dNDash);
578         double d2FX = -dNDash*f2X(dXStar, dA, dBeta, dNDash);
579         double dLogK = 0.0, dGamma1 = dGamma;
580 
581         if(dGamma1 < 0.0){ dGamma1 *= -1.0; } //invert sign
582 
583         dLogK = gsl_sf_bessel_lnKnu(dGamma1,2.0*dA*dBeta);
584 
585         return -2.0*dA*dBeta -log(2.0) -dLogK -dGamma*log(dBeta) + dNDash*dFX + 0.5*log(2.0*M_PI) - 0.5*log(d2FX);
586     }
587     catch(exception& e) {
588         m->errorOut(e, "DiversityUtils", "sd");
589         exit(1);
590     }
591 }
592 /***********************************************************************/
bessel(double * pdResult,int n,double dAlpha,double dBeta)593 bool DiversityUtils::bessel(double* pdResult, int n, double dAlpha, double dBeta){
594     try {
595         double dGamma  = -0.5;
596         double dResult = 0.0;
597         double dOmega = 0.0, dGamma2 = 0.0;
598         double dLogK1 = 0.0, dLogK2 = 0.0;
599         double dN = (double) n, dNu = dGamma + dN;
600         double dTemp1 = 0.0;
601 
602         if(dNu < 0.0){ dNu = -dNu; }
603 
604         if(dGamma < 0.0)    { dGamma2 = -dGamma;    }
605         else                { dGamma2 = dGamma;     }
606 
607         dOmega = sqrt(dBeta*dBeta + dAlpha*dAlpha) - dBeta;
608 
609         dLogK2 = gsl_sf_bessel_lnKnu(dNu, dAlpha);
610 
611         if(!gsl_finite(dLogK2)){
612             if(dAlpha < 0.1*sqrt(dNu + 1.0)){
613                 dLogK2 = gsl_sf_lngamma(dNu) + (dNu - 1.0)*log(2.0) - dNu*log(dAlpha);
614             }else{
615                 (*pdResult) = dResult;
616                 return false;
617             }
618         }
619 
620         dLogK1 = dGamma*log(dOmega/dAlpha) -gsl_sf_bessel_lnKnu(dGamma2,dOmega);
621 
622         dTemp1 = log((dBeta*dOmega)/dAlpha);
623 
624         dResult = dN*dTemp1 + dLogK2 + dLogK1;
625 
626         (*pdResult) = dResult;
627 
628         return true;
629     }
630     catch(exception& e) {
631         m->errorOut(e, "DiversityUtils", "bessel");
632         exit(1);
633     }
634 }
635 /***********************************************************************/
bessel(double * pdResult,int n,double dAlpha,double dBeta,double dGamma)636 bool DiversityUtils::bessel(double* pdResult, int n, double dAlpha, double dBeta, double dGamma){
637     try {
638         double dResult = 0.0;
639         double dOmega = 0.0, dGamma2 = 0.0;
640         double dLogK1 = 0.0, dLogK2 = 0.0;
641         double dN = (double) n, dNu = dGamma + dN;
642         double dTemp1 = 0.0;
643 
644         if(dNu < 0.0){ dNu = -dNu; }
645 
646         if(dGamma < 0.0){ dGamma2 = -dGamma;    }
647         else            { dGamma2 = dGamma;     }
648 
649         dOmega = sqrt(dBeta*dBeta + dAlpha*dAlpha) - dBeta;
650 
651         dLogK2 = gsl_sf_bessel_lnKnu(dNu, dAlpha);
652 
653         if(!gsl_finite(dLogK2)){
654             if(dAlpha < 0.1*sqrt(dNu + 1.0)){
655                 dLogK2 = gsl_sf_lngamma(dNu) + (dNu - 1.0)*log(2.0) - dNu*log(dAlpha);
656             }else{
657                 (*pdResult) = dResult;
658                 return false;
659             }
660         }
661 
662         dLogK1 = dGamma*log(dOmega/dAlpha) -gsl_sf_bessel_lnKnu(dGamma2,dOmega);
663 
664         dTemp1 = log((dBeta*dOmega)/dAlpha);
665 
666         dResult = dN*dTemp1 + dLogK2 + dLogK1;
667 
668         (*pdResult) = dResult;
669 
670         return true;
671     }
672     catch(exception& e) {
673         m->errorOut(e, "DiversityUtils", "bessel");
674         exit(1);
675     }
676 }
677 
678 /***********************************************************************/
679 
outputResults(gsl_vector * ptX,t_Data * ptData,double (* f)(const gsl_vector *,void * params))680 vector<double> DiversityUtils::outputResults(gsl_vector *ptX, t_Data *ptData, double (*f)(const gsl_vector*, void* params)){
681     try {
682         vector<double> results;
683 
684         double dAlpha = 0.0, dBeta = 0.0, dS = 0.0, dL = 0.0, dGamma = 0.0;
685 
686         dAlpha = gsl_vector_get(ptX, 0);
687 
688         dBeta  = gsl_vector_get(ptX, 1);
689 
690         if ((method == "metrols") || (method == "metrosichel"))  {
691             dGamma = gsl_vector_get(ptX, 2);
692             dS = gsl_vector_get(ptX, 3);
693         }else {
694             dS = gsl_vector_get(ptX, 2);
695         }
696 
697         dL = f(ptX, ptData);
698 
699         results.push_back(dAlpha); results.push_back(dBeta);
700 
701         if (method == "metroig")        {
702             m->mothurOut("\nMetroIG - ML simplex: a = " + toString(dAlpha) +  " b = " + toString(dBeta) +  " S = " + toString(dS) +  " NLL = " + toString(dL) + "\n");
703             results.push_back(dS); results.push_back(dL);
704         }
705         else if (method == "metroln")   {
706             m->mothurOut("\nMetroLogNormal - ML simplex: M = " + toString(dAlpha) +  " V = " + toString(dBeta) +  " S = " + toString(dS) +  " NLL = " + toString(dL) + "\n");
707             results.push_back(dS); results.push_back(dL);
708         }
709         else if (method == "metrols")   {
710             m->mothurOut("\nMetroLogStudent - ML simplex: M = " + toString(dAlpha) +  " V = " + toString(dBeta) +  " Nu = " + toString(dGamma) +  " S = " + toString(dS) +  " NLL = " + toString(dL) + "\n");
711             results.push_back(dGamma); results.push_back(dS); results.push_back(dL);
712         }
713         else if (method == "metrosichel")   {
714             m->mothurOut("\nMetroSichel - ML simplex: a = " + toString(dAlpha) +  " b = " + toString(dBeta) +  " g = " + toString(dGamma) +  " S = " + toString(dS) +  " NLL = " + toString(dL) + "\n");
715             results.push_back(dGamma); results.push_back(dS); results.push_back(dL);
716         }
717 
718         return results;
719 
720     }
721     catch(exception& e) {
722         m->errorOut(e, "DiversityUtils", "outputResults");
723         exit(1);
724     }
725 }
726 /***********************************************************************/
minimiseSimplex(gsl_vector * ptX,size_t nP,void * pvData,double (* f)(const gsl_vector *,void * params),double initSimplexSize,double minSimplexSize,double maxSimplexSize)727 int DiversityUtils::minimiseSimplex(gsl_vector* ptX, size_t nP, void* pvData, double (*f)(const gsl_vector*, void* params), double initSimplexSize, double minSimplexSize, double maxSimplexSize){
728     try {
729         const gsl_multimin_fminimizer_type *T =
730         gsl_multimin_fminimizer_nmsimplex;
731         gsl_multimin_fminimizer *s = NULL;
732         gsl_vector *ss;
733         gsl_multimin_function minex_func;
734         size_t iter = 0;
735         int status;
736         double size;
737 
738         /* Initial vertex size vector */
739         ss = gsl_vector_alloc (nP);
740 
741         if (method == "metroig")        {
742             gsl_vector_set_all(ss, initSimplexSize);
743             gsl_vector_set(ss,nP - 1,0.1*gsl_vector_get(ptX,0));
744         }
745         else if (method == "metroln")   {
746             gsl_vector_set_all(ss, initSimplexSize);
747             gsl_vector_set(ss,2,0.1*gsl_vector_get(ptX,2));
748         }
749         else if (method == "metrols" ) {
750             for(int i = 0; i < nP; i++){
751                 gsl_vector_set(ss, i,initSimplexSize*fabs(gsl_vector_get(ptX,i)));
752             }
753         }else if (method == "metrosichel" ) {
754             for(int i = 0; i < nP; i++){
755                 gsl_vector_set(ss, i,initSimplexSize*gsl_vector_get(ptX,i));
756             }
757         }
758 
759         /* Initialize method and iterate */
760         minex_func.f = f;
761         minex_func.n = nP;
762         minex_func.params = pvData;
763 
764         s = gsl_multimin_fminimizer_alloc (T, nP);
765         gsl_multimin_fminimizer_set(s, &minex_func, ptX, ss);
766 
767         do{
768             iter++;
769 
770             if (m->getControl_pressed()) { break; }
771 
772             status = gsl_multimin_fminimizer_iterate(s);
773 
774             if(status) { break; }
775 
776             size = gsl_multimin_fminimizer_size(s);
777             status = gsl_multimin_test_size(size, minSimplexSize);
778 
779             if(status == GSL_SUCCESS){
780                 for(int i = 0; i < nP; i++){
781                     gsl_vector_set(ptX, i, gsl_vector_get(s->x, i));
782                 }
783             }
784         }
785         while(status == GSL_CONTINUE && iter < maxSimplexSize);
786 
787         if(status == GSL_CONTINUE){
788             for(int i = 0; i < nP; i++){
789                 gsl_vector_set(ptX, i, gsl_vector_get(s->x, i));
790             }
791         }
792 
793 
794         gsl_vector_free(ss);
795         gsl_multimin_fminimizer_free (s);
796 
797         return status;
798     }
799     catch(exception& e) {
800         m->errorOut(e, "DiversityUtils", "minimiseSimplex");
801         exit(1);
802     }
803 }
804 /***********************************************************************/
getProposal(gsl_rng * ptGSLRNG,gsl_vector * ptXDash,gsl_vector * ptX,int * pnSDash,int nS,t_Params * ptParams)805 void DiversityUtils::getProposal(gsl_rng *ptGSLRNG, gsl_vector *ptXDash, gsl_vector *ptX, int* pnSDash, int nS, t_Params *ptParams){
806     try {
807         double dDeltaS =  gsl_ran_gaussian(ptGSLRNG, ptParams->dSigmaS);
808         double dDeltaA =  gsl_ran_gaussian(ptGSLRNG, ptParams->dSigmaX);
809         double dDeltaB =  gsl_ran_gaussian(ptGSLRNG, ptParams->dSigmaY);
810 
811         double dDeltaN =  0;
812         int    nSDash = 0;
813         if ((method == "metrols") || (method == "metrosichel")) { dDeltaN = gsl_ran_gaussian(ptGSLRNG, ptParams->dSigmaN); }
814 
815         gsl_vector_set(ptXDash, 0, gsl_vector_get(ptX,0) + dDeltaA);
816         gsl_vector_set(ptXDash, 1, gsl_vector_get(ptX,1) + dDeltaB);
817         if ((method == "metrols") || (method == "metrosichel")) {  gsl_vector_set(ptXDash, 2, gsl_vector_get(ptX,2) + dDeltaN); }
818 
819         nSDash = nS + (int) floor(dDeltaS);
820 
821         if(nSDash < 1){ nSDash = 1; }
822 
823         (*pnSDash) = nSDash;
824     }
825     catch(exception& e) {
826         m->errorOut(e, "DiversityUtils", "getProposal");
827         exit(1);
828     }
829 }
830 /***********************************************************************/
fitSigma(vector<double> acceptanceRates,vector<double> parameterResults,int fitIters,t_Params * ptParams,t_Data * ptData,gsl_vector * ptX,void * f (void * pvInitMetro))831 int DiversityUtils::fitSigma(vector<double> acceptanceRates, vector<double> parameterResults, int fitIters, t_Params *ptParams, t_Data *ptData, gsl_vector* ptX, void* f (void * pvInitMetro)){
832     try {
833         double sigmaA = 0.1;
834         vector<double> defaults;
835         defaults.push_back(ptParams->dSigmaX); defaults.push_back(ptParams->dSigmaY); defaults.push_back(ptParams->dSigmaN);
836         acceptRatioPos defaultRatio = findBest(acceptanceRates);
837 
838         if (defaultRatio.acceptRatio <= 0.05) { return defaultRatio.pos;  }
839 
840         for (int i = 0; i < parameterResults.size(); i++) { parameterResults[i] /= 10.0; }
841 
842         int numTries = 1;
843         map<vector<double>, acceptRatioPos> sigmaToAccept; //sigma value -> acceptance ratio
844         map<acceptRatioPos, vector<double>> acceptToSigma; //acceptance ratio -> sigma value
845 
846         acceptRatioPos temp; //1.0 and pos 0 be default
847         sigmaToAccept[parameterResults] = temp;     //0.01
848         vector<double> testTries = parameterResults;
849         for (int i = 0; i < testTries.size(); i++) { testTries[i] /= 10.0; }
850         sigmaToAccept[testTries] = temp;        //0.01
851         testTries = parameterResults;
852         for (int i = 0; i < testTries.size(); i++) { testTries[i] /= 100.0; }
853         sigmaToAccept[testTries] = temp;       //0.001
854         testTries = parameterResults;
855         for (int i = 0; i < testTries.size(); i++) { testTries[i] /= 1000.0; }
856         sigmaToAccept[testTries] = temp;       //0.001
857         testTries = parameterResults;
858         for (int i = 0; i < testTries.size(); i++) { testTries[i] /= 10000.0; }
859         sigmaToAccept[testTries] = temp;       //0.001
860 
861         double newSigmaA = sigmaA + (sigmaA/2.0);         //0.15
862         testTries = parameterResults;
863         for (int i = 0; i < testTries.size(); i++) { testTries[i] += newSigmaA; }
864         sigmaToAccept[testTries] = temp;
865         testTries = parameterResults;
866         for (int i = 0; i < testTries.size(); i++) { testTries[i] += sigmaA+sigmaA; }
867         sigmaToAccept[testTries] = temp;
868 
869 
870         //adjust around closest "high" and closest "low" values
871         acceptRatioPos thisBestHigh, thisBestLow;
872         map<acceptRatioPos, vector<double>> acceptToSigmaHigh; //acceptance ratio -> sigma value
873         map<acceptRatioPos, vector<double>> acceptToSigmaLow; //acceptance ratio -> sigma value
874 
875         //set iters to 1000, get close to value then run with nIters
876         int savedIters = ptParams->nIter;
877         ptParams->nIter = 1000;
878 
879         for (map<vector<double>, acceptRatioPos>::iterator it = sigmaToAccept.begin(); it != sigmaToAccept.end(); it++) {
880             if (m->getControl_pressed()) { break; }
881 
882             ptParams->dSigmaX = it->first[0]; ptParams->dSigmaY = it->first[1]; ptParams->dSigmaN = it->first[2];
883 
884             acceptanceRates = mcmc(ptParams, ptData, ptX, f);
885 
886             it->second = findBest(acceptanceRates);
887 
888             if (it->second.high) { //high
889                 if (it->second.acceptRatio < thisBestHigh.acceptRatio) {  thisBestHigh = it->second; }
890                 acceptToSigmaHigh[it->second] = it->first;
891             }else { //low
892                 if (it->second.acceptRatio < thisBestLow.acceptRatio) {  thisBestLow = it->second; }
893                 acceptToSigmaLow[it->second] = it->first;
894             }
895             acceptToSigma[it->second] = it->first;
896 
897             if (it->second.acceptRatio <= 0.05) {
898                 //try with nIters to confirm
899                 ptParams->nIter = savedIters;
900 
901                 acceptanceRates = mcmc(ptParams, ptData, ptX, f);
902 
903                 it->second = findBest(acceptanceRates);
904 
905                 if (it->second.high) { //high
906                     if (it->second.acceptRatio < thisBestHigh.acceptRatio) {  thisBestHigh = it->second; }
907                     acceptToSigmaHigh[it->second] = it->first;
908                 }else { //low
909                     if (it->second.acceptRatio < thisBestLow.acceptRatio) {  thisBestLow = it->second; }
910                     acceptToSigmaLow[it->second] = it->first;
911                 }
912                 acceptToSigma[it->second] = it->first;
913 
914                 //if good value
915                 if (it->second.acceptRatio <= 0.05) { return it->second.pos;  }
916                 else {  ptParams->nIter = 1000;  }
917             }
918         }
919 
920         sigmaToAccept[defaults] = defaultRatio;
921         acceptToSigma[defaultRatio] = defaults;
922 
923         vector<double> factors; factors.resize(3, 0); bool badHigh = false; bool badLow = false; double badFactor = 0.0;
924         vector<double> badFactors; badFactors.resize(3, 0);
925 
926         //find best high and check
927         map<acceptRatioPos, vector<double>>::iterator itFind = acceptToSigma.find(thisBestHigh);
928         if (itFind != acceptToSigma.end()) {
929             if (thisBestHigh.acceptRatio > 0.25) {
930                 badHigh = true; for (int i = 0; i < badFactors.size(); i++) { badFactors[i] += itFind->second[i]; }
931             }else {
932                 for (int i = 0; i < factors.size(); i++) { factors[i] += itFind->second[i]; }
933                 sigmaA = itFind->second[0];
934             }
935         }//else no high values
936 
937         //find best low and check
938         itFind = acceptToSigma.find(thisBestLow);
939         if (itFind != acceptToSigma.end()) {
940             if (thisBestLow.acceptRatio > 0.25) { //below 25% acceptance, lets disregard
941                 badLow = true; badHigh = true; for (int i = 0; i < badFactors.size(); i++) { badFactors[i] += itFind->second[i]; }
942             }else {
943                 for (int i = 0; i < factors.size(); i++) { factors[i] += itFind->second[i]; }
944                 if (badHigh) { sigmaA = itFind->second[0]; }
945                 else { if (sigmaA > itFind->second[0]) { sigmaA = itFind->second[0]; } }
946             }
947         }//no low values
948 
949         if (badHigh && badLow) {
950             double increment = badFactor / (double)(fitIters);
951             sigmaA = acceptToSigma.begin()->second[0]; //sigma for best try
952             sigmaA -= (increment*(fitIters/(double)2.0));
953             for (int i = 0; i < factors.size(); i++) { factors[i] = badFactors[i] / (double)(fitIters); }
954         }else if (badHigh || badLow)  {
955             for (int i = 0; i < factors.size(); i++) {
956                 double increment = factors[i] / (double)(fitIters);
957                 sigmaA -= (increment*(fitIters/(double)2.0));
958             }
959         }else { //good high and low
960             for (int i = 0; i < factors.size(); i++) { factors[i] /= (double)(fitIters); }
961         }
962 
963         for (int i = 0; i < factors.size(); i++) {
964             if (util.isEqual(factors[i], 0)) { factors[i] = 0.1; }
965         }
966 
967         ptParams->dSigmaX = acceptToSigma.begin()->second[0]; ptParams->dSigmaY = acceptToSigma.begin()->second[1]; ptParams->dSigmaN = acceptToSigma.begin()->second[2];
968         ptParams->nIter = savedIters;
969 
970         while ((thisBestLow.acceptRatio > 0.05) && (numTries < fitIters)) {
971             if (m->getControl_pressed()) { break; }
972 
973             m->mothurOut("\nFit try: " + toString(numTries) + "\n");
974 
975             ptParams->dSigmaX += factors[0]; ptParams->dSigmaY += factors[1]; ptParams->dSigmaN += factors[2];
976 
977             vector<double> theseSettings; theseSettings.push_back(ptParams->dSigmaX); theseSettings.push_back(ptParams->dSigmaY); theseSettings.push_back(ptParams->dSigmaN);
978             map<vector<double>, acceptRatioPos>::iterator it = sigmaToAccept.find(theseSettings);
979 
980             if (it == sigmaToAccept.end()) {
981                 acceptanceRates = mcmc(ptParams, ptData, ptX, f);
982 
983                 thisBestLow = findBest(acceptanceRates);
984 
985                 acceptToSigma[thisBestLow] = theseSettings;
986                 sigmaToAccept[theseSettings] = thisBestLow;
987                 numTries++;
988             }
989         }
990 
991         if (numTries == fitIters) {
992             vector<double> theBestSettings = acceptToSigma.begin()->second;
993 
994             ptParams->dSigmaX = theBestSettings[0]; ptParams->dSigmaY = theBestSettings[1]; ptParams->dSigmaN = theBestSettings[2];
995 
996             acceptanceRates = mcmc(ptParams, ptData, ptX, f);
997 
998             thisBestLow = findBest(acceptanceRates);
999         }
1000 
1001         if ((thisBestLow.acceptRatio > 0.05)) { m->mothurOut("\n[ERROR]: Unable to reach acceptable ratio, please review and set sigma parameters manually.\n"); m->setControl_pressed(true); }
1002 
1003         return thisBestLow.pos;
1004 
1005     }
1006     catch(exception& e) {
1007         m->errorOut(e, "DiversityUtils", "fitSigma");
1008         exit(1);
1009     }
1010 }
1011 /***********************************************************************/
mcmc(t_Params * ptParams,t_Data * ptData,gsl_vector * ptX,void * f (void * pvInitMetro))1012 vector<double> DiversityUtils::mcmc(t_Params *ptParams, t_Data *ptData, gsl_vector* ptX, void* f (void * pvInitMetro)){
1013     try {
1014         int ptXSize = 3;
1015         if ((method == "metrols") || (method == "metrosichel")) {  ptXSize = 4;  }
1016 
1017         pthread_t thread1, thread2, thread3;
1018         int       iret1  , iret2  , iret3;
1019         gsl_vector *ptX1 = gsl_vector_alloc(ptXSize),
1020         *ptX2 = gsl_vector_alloc(ptXSize),
1021         *ptX3 = gsl_vector_alloc(ptXSize);
1022         t_MetroInit atMetroInit[3];
1023 
1024         if (method == "metrols") {  m->mothurOut("\nMCMC iter = " + toString(ptParams->nIter) + " sigmaM = " + toString(ptParams->dSigmaX) +  " sigmaV = " + toString(ptParams->dSigmaY) +  " sigmaN = " + toString(ptParams->dSigmaN) +  " sigmaS = " + toString(ptParams->dSigmaS) + "\n"); }
1025         else if (method == "metrosichel") {  m->mothurOut("\nMCMC iter = " + toString(ptParams->nIter) + " sigmaA = " + toString(ptParams->dSigmaX) +  " sigmaB = " + toString(ptParams->dSigmaY) +  " sigmaG = " + toString(ptParams->dSigmaN) +  " sigmaS = " + toString(ptParams->dSigmaS) + "\n"); }
1026         else { m->mothurOut("\nMCMC iter = " + toString(ptParams->nIter) + " sigmaX = " + toString(ptParams->dSigmaX) +  " sigmaY = " + toString(ptParams->dSigmaY) +  " sigmaS = " + toString(ptParams->dSigmaS) + "\n"); }
1027 
1028         gsl_vector_memcpy(ptX1, ptX);
1029 
1030         gsl_vector_set(ptX2, 0, gsl_vector_get(ptX,0) + 2.0*ptParams->dSigmaX);
1031         gsl_vector_set(ptX2, 1, gsl_vector_get(ptX,1) + 2.0*ptParams->dSigmaY);
1032 
1033         if ((method == "metrols") || (method == "metrosichel")) {
1034             gsl_vector_set(ptX2, 2, gsl_vector_get(ptX,2) + 2.0*ptParams->dSigmaN);
1035             gsl_vector_set(ptX2, 3, gsl_vector_get(ptX,3) + 2.0*ptParams->dSigmaS);
1036         }
1037         else { gsl_vector_set(ptX2, 2, gsl_vector_get(ptX,2) + 2.0*ptParams->dSigmaS);  }
1038 
1039 
1040         gsl_vector_set(ptX3, 0, gsl_vector_get(ptX,0) - 2.0*ptParams->dSigmaX);
1041         gsl_vector_set(ptX3, 1, gsl_vector_get(ptX,1) - 2.0*ptParams->dSigmaY);
1042 
1043         if ((method == "metrols") || (method == "metrosichel")) {
1044             gsl_vector_set(ptX3, 2, gsl_vector_get(ptX,2) - 2.0*ptParams->dSigmaN);
1045 
1046             if(gsl_vector_get(ptX,3) - 2.0*ptParams->dSigmaS > (double) ptData->nL){
1047                 gsl_vector_set(ptX3, 3, gsl_vector_get(ptX,3) - 2.0*ptParams->dSigmaS);   }
1048             else{ gsl_vector_set(ptX3, 3, (double) ptData->nL);   }
1049         }else {
1050 
1051             if(gsl_vector_get(ptX,2) - 2.0*ptParams->dSigmaS > (double) ptData->nL){ gsl_vector_set(ptX3, 2, gsl_vector_get(ptX,2) - 2.0*ptParams->dSigmaS); }
1052             else{ gsl_vector_set(ptX3, 2, (double) ptData->nL); }
1053         }
1054 
1055         atMetroInit[0].ptParams = ptParams;
1056         atMetroInit[0].ptData   = ptData;
1057         atMetroInit[0].ptX      = ptX1;
1058         atMetroInit[0].nThread  = 0;
1059         atMetroInit[0].lSeed    = ptParams->lSeed;
1060         atMetroInit[0].nAccepted = 0;
1061 
1062         //write thread 0
1063 
1064         if ((method == "metrols") || (method == "metrosichel")) { m->mothurOut(toString(atMetroInit[0].nThread) + ": a = " + toString(gsl_vector_get(ptX1, 0)) +  " b = " + toString(gsl_vector_get(ptX1, 1)) +  " g = " + toString(gsl_vector_get(ptX1, 2)) +  " S = " + toString(gsl_vector_get(ptX1, 3)) + "\n"); }
1065         else { m->mothurOut(toString(atMetroInit[0].nThread) + ": a = " + toString(gsl_vector_get(ptX1, 0)) +  " b = " + toString(gsl_vector_get(ptX1, 1)) +  " S = " + toString(gsl_vector_get(ptX1, 2)) + "\n"); }
1066 
1067         atMetroInit[1].ptParams = ptParams;
1068         atMetroInit[1].ptData   = ptData;
1069         atMetroInit[1].ptX      = ptX2;
1070         atMetroInit[1].nThread  = 1;
1071         atMetroInit[1].lSeed    = ptParams->lSeed + 1;
1072         atMetroInit[1].nAccepted = 0;
1073 
1074         //write thread 1
1075         if ((method == "metrols") || (method == "metrosichel")) { m->mothurOut(toString(atMetroInit[1].nThread) + ": a = " + toString(gsl_vector_get(ptX2, 0)) +  " b = " + toString(gsl_vector_get(ptX2, 1)) +  " g = " + toString(gsl_vector_get(ptX2, 2)) +  " S = " + toString(gsl_vector_get(ptX2, 3)) + "\n"); }
1076         else { m->mothurOut(toString(atMetroInit[1].nThread) + ": a = " + toString(gsl_vector_get(ptX2, 0)) +  " b = " + toString(gsl_vector_get(ptX2, 1)) +  " S = " + toString(gsl_vector_get(ptX2, 2)) + "\n"); }
1077 
1078         atMetroInit[2].ptParams = ptParams;
1079         atMetroInit[2].ptData   = ptData;
1080         atMetroInit[2].ptX      = ptX3;
1081         atMetroInit[2].nThread  = 2;
1082         atMetroInit[2].lSeed    = ptParams->lSeed + 2;
1083         atMetroInit[2].nAccepted = 0;
1084 
1085         //write thread 2
1086         if ((method == "metrols") || (method == "metrosichel")) { m->mothurOut(toString(atMetroInit[2].nThread) + ": a = " + toString(gsl_vector_get(ptX3, 0)) +  " b = " + toString(gsl_vector_get(ptX3, 1)) +  " g = " + toString(gsl_vector_get(ptX3, 2)) +  " S = " + toString(gsl_vector_get(ptX3, 3)) + "\n"); }
1087         else { m->mothurOut(toString(atMetroInit[2].nThread) + ": a = " + toString(gsl_vector_get(ptX3, 0)) +  " b = " + toString(gsl_vector_get(ptX3, 1)) +  " S = " + toString(gsl_vector_get(ptX3, 2)) + "\n"); }
1088 
1089         iret1 = pthread_create(&thread1, NULL, f, (void*) &atMetroInit[0]);
1090         iret2 = pthread_create(&thread2, NULL, f, (void*) &atMetroInit[1]);
1091         iret3 = pthread_create(&thread3, NULL, f, (void*) &atMetroInit[2]);
1092         pthread_join(thread1, NULL);
1093         pthread_join(thread2, NULL);
1094         pthread_join(thread3, NULL);
1095 
1096         m->mothurOut(toString(atMetroInit[0].nThread) +": accept. ratio " + toString(atMetroInit[0].nAccepted) + "/" + toString(ptParams->nIter) +  " = " + toString(((double) atMetroInit[0].nAccepted)/((double) ptParams->nIter)) +  "\n");
1097         m->mothurOut(toString(atMetroInit[1].nThread) +": accept. ratio " + toString(atMetroInit[1].nAccepted) + "/" + toString(ptParams->nIter) +  " = " + toString(((double) atMetroInit[1].nAccepted)/((double) ptParams->nIter)) +  "\n");
1098         m->mothurOut(toString(atMetroInit[2].nThread) +": accept. ratio " + toString(atMetroInit[2].nAccepted) + "/" + toString(ptParams->nIter) +  " = " + toString(((double) atMetroInit[2].nAccepted)/((double) ptParams->nIter)) +  "\n");
1099 
1100         vector<double> results;
1101         results.push_back(atMetroInit[0].nAccepted/((double) ptParams->nIter));
1102         results.push_back(atMetroInit[1].nAccepted/((double) ptParams->nIter));
1103         results.push_back(atMetroInit[2].nAccepted/((double) ptParams->nIter));
1104 
1105         gsl_vector_free(ptX1); gsl_vector_free(ptX2); gsl_vector_free(ptX3);
1106 
1107         return results;
1108     }
1109     catch(exception& e) {
1110         m->errorOut(e, "DiversityUtils", "mcmc");
1111         exit(1);
1112     }
1113 }
1114 
1115 #endif
1116 
1117 /***********************************************************************/
findBest(vector<double> acceptanceRates)1118 acceptRatioPos DiversityUtils::findBest(vector<double> acceptanceRates){
1119     try {
1120 
1121         double defaultSigmaAcc = fabs(0.5 - acceptanceRates[0]);  //"0" version
1122         bool high = true;
1123         if ((0.5 - acceptanceRates[0]) > 0.0) { high = false; }
1124         acceptRatioPos defaultRatio(defaultSigmaAcc, 0, high);
1125 
1126 
1127         if (defaultRatio.acceptRatio > fabs(0.5 - acceptanceRates[1])) {  //is the "1" version better?
1128             defaultRatio.acceptRatio = fabs(0.5 - acceptanceRates[1]);
1129             defaultRatio.pos = 1;
1130             defaultRatio.high = true;
1131             if ((0.5 - acceptanceRates[1]) > 0.0) { defaultRatio.high = false; }
1132         }
1133         if (defaultRatio.acceptRatio > fabs(0.5 - acceptanceRates[2])) {  //is the "2" version better?
1134             defaultRatio.acceptRatio = fabs(0.5 - acceptanceRates[2]);
1135             defaultRatio.pos = 2;
1136             defaultRatio.high = true;
1137             if ((0.5 - acceptanceRates[2]) > 0.0) { defaultRatio.high = false; }
1138         }
1139 
1140         return defaultRatio;
1141     }
1142     catch(exception& e) {
1143         m->errorOut(e, "DiversityUtils", "findBest");
1144         exit(1);
1145     }
1146 }
1147 
1148 /***********************************************************************/
1149 
1150 
1151 
1152 
1153 
1154 
1155