5 // Created by Patrick Schloss on 11/8/12.
6 // Copyright (c) 2012 University of Michigan. All rights reserved.
9 #include "qFinderDMM.h"
10 #include "linearalgebra.h"
12 #define EPSILON numeric_limits<double>::epsilon()
14 /**************************************************************************************************/
16 qFinderDMM::qFinderDMM(vector<vector<int> > cm, int p): countMatrix(cm), numPartitions(p){
18 m = MothurOut::getInstance();
19 numSamples = (int)countMatrix.size();
20 numOTUs = (int)countMatrix[0].size();
24 // cout << "done kMeans" << endl;
29 // cout << "done optimizeLambda" << endl;
31 double change = 1.0000;
36 while(change > 1.0e-6 && iter < 100){
38 // cout << "Calc_Z: ";
43 // printf("Iter:%d\t",iter);
45 for(int i=0;i<numPartitions;i++){
47 for(int j=0;j<numSamples;j++){
48 weights[i] += zMatrix[i][j];
50 // printf("w_%d=%.3f\t",i,weights[i]);
54 double nLL = getNegativeLogLikelihood();
56 change = abs(nLL - currNLL);
60 // printf("NLL=%.5f\tDelta=%.4e\n",currNLL, change);
65 error.resize(numPartitions);
67 logDeterminant = 0.0000;
71 for(currentPartition=0;currentPartition<numPartitions;currentPartition++){
73 error[currentPartition].assign(numOTUs, 0.0000);
75 if(currentPartition > 0){
76 logDeterminant += (2.0 * log(numSamples) - log(weights[currentPartition]));
78 vector<vector<double> > hessian = getHessian();
79 vector<vector<double> > invHessian = l.getInverse(hessian);
81 for(int i=0;i<numOTUs;i++){
82 logDeterminant += log(abs(hessian[i][i]));
83 error[currentPartition][i] = invHessian[i][i];
87 int numParameters = numPartitions * numOTUs + numPartitions - 1;
88 laplace = currNLL + 0.5 * logDeterminant - 0.5 * numParameters * log(2.0 * 3.14159);
89 bic = currNLL + 0.5 * log(numSamples) * numParameters;
90 aic = currNLL + numParameters;
93 m->errorOut(e, "qFinderDMM", "qFinderDMM");
98 /**************************************************************************************************/
100 void qFinderDMM::printZMatrix(string fileName, vector<string> sampleName){
102 ofstream printMatrix;
103 m->openOutputFile(fileName, printMatrix); //(fileName.c_str());
104 printMatrix.setf(ios::fixed, ios::floatfield);
105 printMatrix.setf(ios::showpoint);
107 for(int i=0;i<numPartitions;i++){ printMatrix << "\tPartition_" << i+1; } printMatrix << endl;
109 for(int i=0;i<numSamples;i++){
110 printMatrix << sampleName[i];
111 for(int j=0;j<numPartitions;j++){
112 printMatrix << setprecision(4) << '\t' << zMatrix[j][i];
118 catch(exception& e) {
119 m->errorOut(e, "qFinderDMM", "printZMatrix");
124 /**************************************************************************************************/
126 void qFinderDMM::printRelAbund(string fileName, vector<string> otuNames){
129 m->openOutputFile(fileName, printRA); //(fileName.c_str());
130 printRA.setf(ios::fixed, ios::floatfield);
131 printRA.setf(ios::showpoint);
133 vector<double> totals(numPartitions, 0.0000);
134 for(int i=0;i<numPartitions;i++){
135 for(int j=0;j<numOTUs;j++){
136 totals[i] += exp(lambdaMatrix[i][j]);
141 for(int i=0;i<numPartitions;i++){
142 printRA << "\tPartition_" << i+1 << '_' << setprecision(4) << totals[i];
143 printRA << "\tPartition_" << i+1 <<"_LCI" << "\tPartition_" << i+1 << "_UCI";
147 for(int i=0;i<numOTUs;i++){
149 if (m->control_pressed) { break; }
151 printRA << otuNames[i];
152 for(int j=0;j<numPartitions;j++){
154 if(error[j][i] >= 0.0000){
155 double std = sqrt(error[j][i]);
156 printRA << '\t' << 100 * exp(lambdaMatrix[j][i]) / totals[j];
157 printRA << '\t' << 100 * exp(lambdaMatrix[j][i] - 2.0 * std) / totals[j];
158 printRA << '\t' << 100 * exp(lambdaMatrix[j][i] + 2.0 * std) / totals[j];
161 printRA << '\t' << 100 * exp(lambdaMatrix[j][i]) / totals[j];
162 printRA << '\t' << "NA";
163 printRA << '\t' << "NA";
171 catch(exception& e) {
172 m->errorOut(e, "qFinderDMM", "printRelAbund");
177 /**************************************************************************************************/
179 // these functions for bfgs2 solver were lifted from the gnu_gsl source code...
181 /* Find a minimum in x=[0,1] of the interpolating quadratic through
182 * (0,f0) (1,f1) with derivative fp0 at x=0. The interpolating
183 * polynomial is q(x) = f0 + fp0 * z + (f1-f0-fp0) * z^2
187 interp_quad (double f0, double fp0, double f1, double zl, double zh)
189 double fl = f0 + zl*(fp0 + zl*(f1 - f0 -fp0));
190 double fh = f0 + zh*(fp0 + zh*(f1 - f0 -fp0));
191 double c = 2 * (f1 - f0 - fp0); /* curvature */
193 double zmin = zl, fmin = fl;
195 if (fh < fmin) { zmin = zh; fmin = fh; }
197 if (c > 0) /* positive curvature required for a minimum */
199 double z = -fp0 / c; /* location of minimum */
200 if (z > zl && z < zh) {
201 double f = f0 + z*(fp0 + z*(f1 - f0 -fp0));
202 if (f < fmin) { zmin = z; fmin = f; };
209 /**************************************************************************************************/
211 /* Find a minimum in x=[0,1] of the interpolating cubic through
212 * (0,f0) (1,f1) with derivatives fp0 at x=0 and fp1 at x=1.
214 * The interpolating polynomial is:
216 * c(x) = f0 + fp0 * z + eta * z^2 + xi * z^3
218 * where eta=3*(f1-f0)-2*fp0-fp1, xi=fp0+fp1-2*(f1-f0).
221 double cubic (double c0, double c1, double c2, double c3, double z){
222 return c0 + z * (c1 + z * (c2 + z * c3));
225 /**************************************************************************************************/
227 void check_extremum (double c0, double c1, double c2, double c3, double z,
228 double *zmin, double *fmin){
229 /* could make an early return by testing curvature >0 for minimum */
231 double y = cubic (c0, c1, c2, c3, z);
235 *zmin = z; /* accepted new point*/
240 /**************************************************************************************************/
242 int gsl_poly_solve_quadratic (double a, double b, double c,
243 double *x0, double *x1)
246 double disc = b * b - 4 * a * c;
248 if (a == 0) /* Handle linear case */
265 double r = fabs (0.5 * sqrt (disc) / a);
271 double sgnb = (b > 0 ? 1 : -1);
272 double temp = -0.5 * (b + sgnb * sqrt (disc));
273 double r1 = temp / a ;
274 double r2 = c / temp ;
302 /**************************************************************************************************/
304 double interp_cubic (double f0, double fp0, double f1, double fp1, double zl, double zh){
305 double eta = 3 * (f1 - f0) - 2 * fp0 - fp1;
306 double xi = fp0 + fp1 - 2 * (f1 - f0);
307 double c0 = f0, c1 = fp0, c2 = eta, c3 = xi;
311 zmin = zl; fmin = cubic(c0, c1, c2, c3, zl);
312 check_extremum (c0, c1, c2, c3, zh, &zmin, &fmin);
315 int n = gsl_poly_solve_quadratic (3 * c3, 2 * c2, c1, &z0, &z1);
317 if (n == 2) /* found 2 roots */
319 if (z0 > zl && z0 < zh)
320 check_extremum (c0, c1, c2, c3, z0, &zmin, &fmin);
321 if (z1 > zl && z1 < zh)
322 check_extremum (c0, c1, c2, c3, z1, &zmin, &fmin);
324 else if (n == 1) /* found 1 root */
326 if (z0 > zl && z0 < zh)
327 check_extremum (c0, c1, c2, c3, z0, &zmin, &fmin);
334 /**************************************************************************************************/
336 double interpolate (double a, double fa, double fpa,
337 double b, double fb, double fpb, double xmin, double xmax){
338 /* Map [a,b] to [0,1] */
339 double z, alpha, zmin, zmax;
341 zmin = (xmin - a) / (b - a);
342 zmax = (xmax - a) / (b - a);
352 z = interp_cubic (fa, fpa * (b - a), fb, fpb * (b - a), zmin, zmax);
355 z = interp_quad(fa, fpa * (b - a), fb, zmin, zmax);
359 alpha = a + z * (b - a);
364 /**************************************************************************************************/
366 int qFinderDMM::lineMinimizeFletcher(vector<double>& x, vector<double>& p, double f0, double df0, double alpha1, double& alphaNew, double& fAlpha, vector<double>& xalpha, vector<double>& gradient ){
375 double alpha = alpha1;
376 double alpha_prev = 0.0000;
378 xalpha.resize(numOTUs, 0.0000);
380 double falpha_prev = f0;
381 double dfalpha_prev = df0;
383 double a = 0.0000; double b = alpha;
384 double fa = f0; double fb = 0.0000;
385 double dfa = df0; double dfb = 0.0/0.0;
389 while(iter++ < maxIters){
390 if (m->control_pressed) { break; }
392 for(int i=0;i<numOTUs;i++){
393 xalpha[i] = x[i] + alpha * p[i];
396 fAlpha = negativeLogEvidenceLambdaPi(xalpha);
398 if(fAlpha > f0 + alpha * rho * df0 || fAlpha >= falpha_prev){
399 a = alpha_prev; b = alpha;
400 fa = falpha_prev; fb = fAlpha;
401 dfa = dfalpha_prev; dfb = 0.0/0.0;
405 negativeLogDerivEvidenceLambdaPi(xalpha, gradient);
406 double dfalpha = 0.0000;
407 for(int i=0;i<numOTUs;i++){ dfalpha += gradient[i] * p[i]; }
409 if(abs(dfalpha) <= -sigma * df0){
415 a = alpha; b = alpha_prev;
416 fa = fAlpha; fb = falpha_prev;
417 dfa = dfalpha; dfb = dfalpha_prev;
421 double delta = alpha - alpha_prev;
423 double lower = alpha + delta;
424 double upper = alpha + tau1 * delta;
426 double alphaNext = interpolate(alpha_prev, falpha_prev, dfalpha_prev, alpha, fAlpha, dfalpha, lower, upper);
429 falpha_prev = fAlpha;
430 dfalpha_prev = dfalpha;
435 while(iter++ < maxIters){
436 if (m->control_pressed) { break; }
438 double delta = b - a;
440 double lower = a + tau2 * delta;
441 double upper = b - tau3 * delta;
443 alpha = interpolate(a, fa, dfa, b, fb, dfb, lower, upper);
445 for(int i=0;i<numOTUs;i++){
446 xalpha[i] = x[i] + alpha * p[i];
449 fAlpha = negativeLogEvidenceLambdaPi(xalpha);
451 if((a - alpha) * dfa <= EPSILON){
455 if(fAlpha > f0 + rho * alpha * df0 || fAlpha >= fa){
461 double dfalpha = 0.0000;
463 negativeLogDerivEvidenceLambdaPi(xalpha, gradient);
465 for(int i=0;i<numOTUs;i++){ dfalpha += gradient[i] * p[i]; }
467 if(abs(dfalpha) <= -sigma * df0){
472 if(((b-a >= 0 && dfalpha >= 0) || ((b-a) <= 0.000 && dfalpha <= 0))){
473 b = a; fb = fa; dfb = dfa;
474 a = alpha; fa = fAlpha; dfa = dfalpha;
488 catch(exception& e) {
489 m->errorOut(e, "qFinderDMM", "lineMinimizeFletcher");
494 /**************************************************************************************************/
496 int qFinderDMM::bfgs2_Solver(vector<double>& x){
498 // cout << "bfgs2_Solver" << endl;
500 double step = 1.0e-6;
501 double delta_f = 0.0000;//f-f0;
503 vector<double> gradient;
504 double f = negativeLogEvidenceLambdaPi(x);
506 // cout << "after negLE" << endl;
508 negativeLogDerivEvidenceLambdaPi(x, gradient);
510 // cout << "after negLDE" << endl;
512 vector<double> x0 = x;
513 vector<double> g0 = gradient;
516 for(int i=0;i<numOTUs;i++){
517 g0norm += g0[i] * g0[i];
519 g0norm = sqrt(g0norm);
521 vector<double> p = gradient;
523 for(int i=0;i<numOTUs;i++){
525 pNorm += p[i] * p[i];
528 double df0 = -g0norm;
532 // cout << "before while" << endl;
534 while(g0norm > 0.001 && bfgsIter++ < maxIter){
535 if (m->control_pressed) { return 0; }
538 vector<double> dx(numOTUs, 0.0000);
540 double alphaOld, alphaNew;
542 if(pNorm == 0 || g0norm == 0 || df0 == 0){
543 dx.assign(numOTUs, 0.0000);
547 double delta = max(-delta_f, 10 * EPSILON * abs(f0));
548 alphaOld = min(1.0, 2.0 * delta / (-df0));
554 int success = lineMinimizeFletcher(x0, p, f0, df0, alphaOld, alphaNew, f, x, gradient);
563 vector<double> dx0(numOTUs);
564 vector<double> dg0(numOTUs);
566 for(int i=0;i<numOTUs;i++){
567 dx0[i] = x[i] - x0[i];
568 dg0[i] = gradient[i] - g0[i];
576 for(int i=0;i<numOTUs;i++){
577 dxg += dx0[i] * gradient[i];
578 dgg += dg0[i] * gradient[i];
579 dxdg += dx0[i] * dg0[i];
580 dgnorm += dg0[i] * dg0[i];
582 dgnorm = sqrt(dgnorm);
588 A = -(1.0 + dgnorm*dgnorm /dxdg) * B + dgg / dxdg;
595 for(int i=0;i<numOTUs;i++){ p[i] = gradient[i] - A * dx0[i] - B * dg0[i]; }
605 for(int i=0;i<numOTUs;i++){
606 pg += p[i] * gradient[i];
607 pNorm += p[i] * p[i];
608 g0norm += g0[i] * g0[i];
611 g0norm = sqrt(g0norm);
613 double dir = (pg >= 0.0) ? -1.0 : +1.0;
615 for(int i=0;i<numOTUs;i++){ p[i] *= dir / pNorm; }
619 for(int i=0;i<numOTUs;i++){
620 pNorm += p[i] * p[i];
627 // cout << "bfgsIter:\t" << bfgsIter << endl;
632 m->errorOut(e, "qFinderDMM", "bfgs2_Solver");
637 /**************************************************************************************************/
639 //can we get these psi/psi1 calculations into their own math class?
640 //psi calcualtions swiped from gsl library...
642 static const double psi_cs[23] = {
643 -.038057080835217922,
645 -.056815747821244730,
647 -.001333232857994342,
649 -.000037040238178456,
651 -.000001071263908506,
653 -.000000031353509361,
655 -.000000000921168141,
657 -.000000000027098646,
659 -.000000000000797527,
661 -.000000000000023475,
663 -.000000000000000691,
668 static double apsi_cs[16] = {
687 /**************************************************************************************************/
689 double qFinderDMM::cheb_eval(const double seriesData[], int order, double xx){
694 double x2 = xx * 2.0000;
696 for(int j=order;j>=1;j--){
697 if (m->control_pressed) { return 0; }
699 d = x2 * d - dd + seriesData[j];
703 d = xx * d - dd + 0.5 * seriesData[0];
707 m->errorOut(e, "qFinderDMM", "cheb_eval");
712 /**************************************************************************************************/
714 double qFinderDMM::psi(double xx){
716 double psiX = 0.0000;
720 double t1 = 1.0 / xx;
721 psiX = cheb_eval(psi_cs, 22, 2.0*xx-1.0);
725 else if(xx < 2.0000){
727 const double v = xx - 1.0;
728 psiX = cheb_eval(psi_cs, 22, 2.0*v-1.0);
732 const double t = 8.0/(xx*xx)-1.0;
733 psiX = cheb_eval(apsi_cs, 15, t);
734 psiX += log(xx) - 0.5/xx;
740 m->errorOut(e, "qFinderDMM", "psi");
745 /**************************************************************************************************/
747 /* coefficients for Maclaurin summation in hzeta()
750 static double hzeta_c[15] = {
751 1.00000000000000000000000000000,
752 0.083333333333333333333333333333,
753 -0.00138888888888888888888888888889,
754 0.000033068783068783068783068783069,
755 -8.2671957671957671957671957672e-07,
756 2.0876756987868098979210090321e-08,
757 -5.2841901386874931848476822022e-10,
758 1.3382536530684678832826980975e-11,
759 -3.3896802963225828668301953912e-13,
760 8.5860620562778445641359054504e-15,
761 -2.1748686985580618730415164239e-16,
762 5.5090028283602295152026526089e-18,
763 -1.3954464685812523340707686264e-19,
764 3.5347070396294674716932299778e-21,
765 -8.9535174270375468504026113181e-23
768 /**************************************************************************************************/
770 double qFinderDMM::psi1(double xx){
773 /* Euler-Maclaurin summation formula
774 * [Moshier, p. 400, with several typo corrections]
781 const double pmax = pow(kmax + xx, -s);
783 double pcp = pmax / (kmax + xx);
784 double value = pmax*((kmax+xx)/(s-1.0) + 0.5);
786 for(k=0; k<kmax; k++) {
787 if (m->control_pressed) { return 0; }
788 value += pow(k + xx, -s);
791 for(j=0; j<=jmax; j++) {
792 if (m->control_pressed) { return 0; }
793 double delta = hzeta_c[j+1] * scp * pcp;
796 if(fabs(delta/value) < 0.5*EPSILON) break;
798 scp *= (s+2*j+1)*(s+2*j+2);
799 pcp /= (kmax + xx)*(kmax + xx);
805 m->errorOut(e, "qFinderDMM", "psi1");
810 /**************************************************************************************************/
812 double qFinderDMM::negativeLogEvidenceLambdaPi(vector<double>& x){
814 vector<double> sumAlphaX(numSamples, 0.0000);
816 double logEAlpha = 0.0000;
817 double sumLambda = 0.0000;
818 double sumAlpha = 0.0000;
819 double logE = 0.0000;
821 double eta = 0.10000;
823 double weight = 0.00000;
824 for(int i=0;i<numSamples;i++){
825 weight += zMatrix[currentPartition][i];
828 for(int i=0;i<numOTUs;i++){
829 if (m->control_pressed) { return 0; }
830 double lambda = x[i];
831 double alpha = exp(x[i]);
832 logEAlpha += lgamma(alpha);
836 for(int j=0;j<numSamples;j++){
837 double X = countMatrix[j][i];
838 double alphaX = alpha + X;
839 sumAlphaX[j] += alphaX;
841 logE -= zMatrix[currentPartition][j] * lgamma(alphaX);
845 logEAlpha -= lgamma(sumAlpha);
847 for(int i=0;i<numSamples;i++){
848 logE += zMatrix[currentPartition][i] * lgamma(sumAlphaX[i]);
851 return logE + weight * logEAlpha + nu * sumAlpha - eta * sumLambda;
854 m->errorOut(e, "qFinderDMM", "negativeLogEvidenceLambdaPi");
859 /**************************************************************************************************/
861 void qFinderDMM::negativeLogDerivEvidenceLambdaPi(vector<double>& x, vector<double>& df){
863 // cout << "\tstart negativeLogDerivEvidenceLambdaPi" << endl;
865 vector<double> storeVector(numSamples, 0.0000);
866 vector<double> derivative(numOTUs, 0.0000);
867 vector<double> alpha(numOTUs, 0.0000);
869 double store = 0.0000;
873 double weight = 0.0000;
874 for(int i=0;i<numSamples;i++){
875 weight += zMatrix[currentPartition][i];
879 for(int i=0;i<numOTUs;i++){
880 if (m->control_pressed) { return; }
881 // cout << "start i loop" << endl;
883 // cout << i << '\t' << alpha[i] << '\t' << x[i] << '\t' << exp(x[i]) << '\t' << store << endl;
885 alpha[i] = exp(x[i]);
888 // cout << "before derivative" << endl;
890 derivative[i] = weight * psi(alpha[i]);
892 // cout << "after derivative" << endl;
894 // cout << i << '\t' << alpha[i] << '\t' << psi(alpha[i]) << '\t' << derivative[i] << endl;
896 for(int j=0;j<numSamples;j++){
897 double X = countMatrix[j][i];
898 double alphaX = X + alpha[i];
900 derivative[i] -= zMatrix[currentPartition][j] * psi(alphaX);
901 storeVector[j] += alphaX;
903 // cout << "end i loop" << endl;
906 double sumStore = 0.0000;
907 for(int i=0;i<numSamples;i++){
908 sumStore += zMatrix[currentPartition][i] * psi(storeVector[i]);
911 store = weight * psi(store);
913 df.resize(numOTUs, 0.0000);
915 for(int i=0;i<numOTUs;i++){
916 df[i] = alpha[i] * (nu + derivative[i] - store + sumStore) - eta;
917 // cout << i << '\t' << df[i] << endl;
919 // cout << df.size() << endl;
920 // cout << "\tend negativeLogDerivEvidenceLambdaPi" << endl;
923 m->errorOut(e, "qFinderDMM", "negativeLogDerivEvidenceLambdaPi");
928 /**************************************************************************************************/
930 double qFinderDMM::getNegativeLogEvidence(vector<double>& lambda, int group){
932 double sumAlpha = 0.0000;
933 double sumAlphaX = 0.0000;
934 double sumLnGamAlpha = 0.0000;
935 double logEvidence = 0.0000;
937 for(int i=0;i<numOTUs;i++){
938 if (m->control_pressed) { return 0; }
939 double alpha = exp(lambda[i]);
940 double X = countMatrix[group][i];
941 double alphaX = alpha + X;
943 sumLnGamAlpha += lgamma(alpha);
947 logEvidence -= lgamma(alphaX);
950 sumLnGamAlpha -= lgamma(sumAlpha);
951 logEvidence += lgamma(sumAlphaX);
953 return logEvidence + sumLnGamAlpha;
956 m->errorOut(e, "qFinderDMM", "getNegativeLogEvidence");
961 /**************************************************************************************************/
963 void qFinderDMM::kMeans(){
966 vector<vector<double> > relativeAbundance(numSamples);
967 vector<vector<double> > alphaMatrix;
969 alphaMatrix.resize(numPartitions);
970 lambdaMatrix.resize(numPartitions);
971 for(int i=0;i<numPartitions;i++){
972 alphaMatrix[i].assign(numOTUs, 0);
973 lambdaMatrix[i].assign(numOTUs, 0);
976 //get relative abundance
977 for(int i=0;i<numSamples;i++){
978 if (m->control_pressed) { return; }
981 relativeAbundance[i].assign(numOTUs, 0.0);
983 for(int j=0;j<numOTUs;j++){
984 groupTotal += countMatrix[i][j];
986 for(int j=0;j<numOTUs;j++){
987 relativeAbundance[i][j] = countMatrix[i][j] / (double)groupTotal;
991 //randomly assign samples into partitions
992 zMatrix.resize(numPartitions);
993 for(int i=0;i<numPartitions;i++){
994 zMatrix[i].assign(numSamples, 0);
997 for(int i=0;i<numSamples;i++){
998 zMatrix[rand()%numPartitions][i] = 1;
1001 double maxChange = 1;
1002 int maxIters = 1000;
1005 weights.assign(numPartitions, 0);
1007 while(maxChange > 1e-6 && iteration < maxIters){
1009 if (m->control_pressed) { return; }
1010 //calcualte average relative abundance
1012 for(int i=0;i<numPartitions;i++){
1014 double normChange = 0.0;
1018 for(int j=0;j<numSamples;j++){
1019 weights[i] += (double)zMatrix[i][j];
1022 vector<double> averageRelativeAbundance(numOTUs, 0);
1023 for(int j=0;j<numOTUs;j++){
1024 for(int k=0;k<numSamples;k++){
1025 averageRelativeAbundance[j] += zMatrix[i][k] * relativeAbundance[k][j];
1029 for(int j=0;j<numOTUs;j++){
1030 averageRelativeAbundance[j] /= weights[i];
1031 double difference = averageRelativeAbundance[j] - alphaMatrix[i][j];
1032 normChange += difference * difference;
1033 alphaMatrix[i][j] = averageRelativeAbundance[j];
1036 normChange = sqrt(normChange);
1038 if(normChange > maxChange){ maxChange = normChange; }
1042 //calcualte distance between each sample in partition adn the average relative abundance
1043 for(int i=0;i<numSamples;i++){
1044 if (m->control_pressed) { return; }
1046 double normalizationFactor = 0;
1047 vector<double> totalDistToPartition(numPartitions, 0);
1049 for(int j=0;j<numPartitions;j++){
1050 for(int k=0;k<numOTUs;k++){
1051 double difference = alphaMatrix[j][k] - relativeAbundance[i][k];
1052 totalDistToPartition[j] += difference * difference;
1054 totalDistToPartition[j] = sqrt(totalDistToPartition[j]);
1055 normalizationFactor += exp(-50.0 * totalDistToPartition[j]);
1059 for(int j=0;j<numPartitions;j++){
1060 zMatrix[j][i] = exp(-50.0 * totalDistToPartition[j]) / normalizationFactor;
1066 // cout << "K means: " << iteration << '\t' << maxChange << endl;
1070 // cout << "Iter:-1";
1071 for(int i=0;i<numPartitions;i++){
1072 weights[i] = 0.0000;
1074 for(int j=0;j<numSamples;j++){
1075 weights[i] += zMatrix[i][j];
1077 // printf("\tw_%d=%.3f", i, weights[i]);
1082 for(int i=0;i<numOTUs;i++){
1083 if (m->control_pressed) { return; }
1084 for(int j=0;j<numPartitions;j++){
1085 if(alphaMatrix[j][i] > 0){
1086 lambdaMatrix[j][i] = log(alphaMatrix[j][i]);
1089 lambdaMatrix[j][i] = -10.0;
1094 catch(exception& e){
1095 m->errorOut(e, "qFinderDMM", "kMeans");
1100 /**************************************************************************************************/
1102 void qFinderDMM::optimizeLambda(){
1104 for(currentPartition=0;currentPartition<numPartitions;currentPartition++){
1105 if (m->control_pressed) { return; }
1106 bfgs2_Solver(lambdaMatrix[currentPartition]);
1109 catch(exception& e){
1110 m->errorOut(e, "qFinderDMM", "optimizeLambda");
1115 /**************************************************************************************************/
1117 void qFinderDMM::calculatePiK(){
1119 vector<double> store(numPartitions);
1121 for(int i=0;i<numSamples;i++){
1122 if (m->control_pressed) { return; }
1123 double sum = 0.0000;
1124 double minNegLogEvidence =numeric_limits<double>::max();
1126 for(int j=0;j<numPartitions;j++){
1127 double negLogEvidenceJ = getNegativeLogEvidence(lambdaMatrix[j], i);
1129 if(negLogEvidenceJ < minNegLogEvidence){
1130 minNegLogEvidence = negLogEvidenceJ;
1132 store[j] = negLogEvidenceJ;
1135 for(int j=0;j<numPartitions;j++){
1136 if (m->control_pressed) { return; }
1137 zMatrix[j][i] = weights[j] * exp(-(store[j] - minNegLogEvidence));
1138 sum += zMatrix[j][i];
1141 for(int j=0;j<numPartitions;j++){
1142 zMatrix[j][i] /= sum;
1147 catch(exception& e){
1148 m->errorOut(e, "qFinderDMM", "calculatePiK");
1154 /**************************************************************************************************/
1156 double qFinderDMM::getNegativeLogLikelihood(){
1158 double eta = 0.10000;
1159 double nu = 0.10000;
1161 vector<double> pi(numPartitions, 0.0000);
1162 vector<double> logBAlpha(numPartitions, 0.0000);
1164 double doubleSum = 0.0000;
1166 for(int i=0;i<numPartitions;i++){
1167 if (m->control_pressed) { return 0; }
1168 double sumAlphaK = 0.0000;
1170 pi[i] = weights[i] / (double)numSamples;
1172 for(int j=0;j<numOTUs;j++){
1173 double alpha = exp(lambdaMatrix[i][j]);
1176 logBAlpha[i] += lgamma(alpha);
1178 logBAlpha[i] -= lgamma(sumAlphaK);
1181 for(int i=0;i<numSamples;i++){
1182 if (m->control_pressed) { return 0; }
1184 double probability = 0.0000;
1185 double factor = 0.0000;
1186 double sum = 0.0000;
1187 vector<double> logStore(numPartitions, 0.0000);
1188 double offset = -numeric_limits<double>::max();
1190 for(int j=0;j<numOTUs;j++){
1191 sum += countMatrix[i][j];
1192 factor += lgamma(countMatrix[i][j] + 1.0000);
1194 factor -= lgamma(sum + 1.0);
1196 for(int k=0;k<numPartitions;k++){
1198 double sumAlphaKX = 0.0000;
1199 double logBAlphaX = 0.0000;
1201 for(int j=0;j<numOTUs;j++){
1202 double alphaX = exp(lambdaMatrix[k][j]) + (double)countMatrix[i][j];
1204 sumAlphaKX += alphaX;
1205 logBAlphaX += lgamma(alphaX);
1208 logBAlphaX -= lgamma(sumAlphaKX);
1210 logStore[k] = logBAlphaX - logBAlpha[k] - factor;
1211 if(logStore[k] > offset){
1212 offset = logStore[k];
1217 for(int k=0;k<numPartitions;k++){
1218 probability += pi[k] * exp(-offset + logStore[k]);
1220 doubleSum += log(probability) + offset;
1224 double L5 = - numOTUs * numPartitions * lgamma(eta);
1225 double L6 = eta * numPartitions * numOTUs * log(nu);
1227 double alphaSum, lambdaSum;
1228 alphaSum = lambdaSum = 0.0000;
1230 for(int i=0;i<numPartitions;i++){
1231 for(int j=0;j<numOTUs;j++){
1232 if (m->control_pressed) { return 0; }
1233 alphaSum += exp(lambdaMatrix[i][j]);
1234 lambdaSum += lambdaMatrix[i][j];
1240 return (-doubleSum - L5 - L6 - alphaSum - lambdaSum);
1242 catch(exception& e){
1243 m->errorOut(e, "qFinderDMM", "getNegativeLogLikelihood");
1250 /**************************************************************************************************/
1252 vector<vector<double> > qFinderDMM::getHessian(){
1254 vector<double> alpha(numOTUs, 0.0000);
1255 double alphaSum = 0.0000;
1257 vector<double> pi = zMatrix[currentPartition];
1258 vector<double> psi_ajk(numOTUs, 0.0000);
1259 vector<double> psi_cjk(numOTUs, 0.0000);
1260 vector<double> psi1_ajk(numOTUs, 0.0000);
1261 vector<double> psi1_cjk(numOTUs, 0.0000);
1263 for(int j=0;j<numOTUs;j++){
1265 if (m->control_pressed) { break; }
1267 alpha[j] = exp(lambdaMatrix[currentPartition][j]);
1268 alphaSum += alpha[j];
1270 for(int i=0;i<numSamples;i++){
1271 double X = (double) countMatrix[i][j];
1273 psi_ajk[j] += pi[i] * psi(alpha[j]);
1274 psi1_ajk[j] += pi[i] * psi1(alpha[j]);
1276 psi_cjk[j] += pi[i] * psi(alpha[j] + X);
1277 psi1_cjk[j] += pi[i] * psi1(alpha[j] + X);
1282 double psi_Ck = 0.0000;
1283 double psi1_Ck = 0.0000;
1285 double weight = 0.0000;
1287 for(int i=0;i<numSamples;i++){
1288 if (m->control_pressed) { break; }
1290 double sum = 0.0000;
1291 for(int j=0;j<numOTUs;j++){ sum += alpha[j] + countMatrix[i][j]; }
1293 psi_Ck += pi[i] * psi(sum);
1294 psi1_Ck += pi[i] * psi1(sum);
1297 double psi_Ak = weight * psi(alphaSum);
1298 double psi1_Ak = weight * psi1(alphaSum);
1300 vector<vector<double> > hessian(numOTUs);
1301 for(int i=0;i<numOTUs;i++){ hessian[i].assign(numOTUs, 0.0000); }
1303 for(int i=0;i<numOTUs;i++){
1304 if (m->control_pressed) { break; }
1305 double term1 = -alpha[i] * (- psi_ajk[i] + psi_Ak + psi_cjk[i] - psi_Ck);
1306 double term2 = -alpha[i] * alpha[i] * (-psi1_ajk[i] + psi1_Ak + psi1_cjk[i] - psi1_Ck);
1307 double term3 = 0.1 * alpha[i];
1309 hessian[i][i] = term1 + term2 + term3;
1311 for(int j=0;j<i;j++){
1312 hessian[i][j] = - alpha[i] * alpha[j] * (psi1_Ak - psi1_Ck);
1313 hessian[j][i] = hessian[i][j];
1319 catch(exception& e){
1320 m->errorOut(e, "qFinderDMM", "getHessian");
1325 /**************************************************************************************************/