5 * Created by westcott on 1/7/11.
6 * Copyright 2011 Schloss Lab. All rights reserved.
10 #include "linearalgebra.h"
13 // This class references functions used from "Numerical Recipes in C++" //
15 /*********************************************************************************************************************************/
16 inline double SQR(const double a)
20 /*********************************************************************************************************************************/
22 inline double SIGN(const double a, const double b)
24 return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a);
26 /*********************************************************************************************************************************/
27 //NUmerical recipes pg. 245 - Returns the complementary error function erfc(x) with fractional error everywhere less than 1.2 × 10−7.
28 double LinearAlgebra::erfcc(double x){
34 ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
35 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
36 t*(-0.82215223+t*0.17087277)))))))));
38 //cout << "in erfcc " << t << '\t' << ans<< endl;
39 return (x >= 0.0 ? ans : 2.0 - ans);
42 m->errorOut(e, "LinearAlgebra", "betai");
46 /*********************************************************************************************************************************/
47 //Numerical Recipes pg. 232
48 double LinearAlgebra::betai(const double a, const double b, const double x) {
53 if (x < 0.0 || x > 1.0) { m->mothurOut("[ERROR]: bad x in betai.\n"); m->control_pressed = true; return 0.0; }
55 if (x == 0.0 || x == 1.0) { bt = 0.0; }
56 else { bt = exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x)); }
58 if (x < (a+1.0) / (a+b+2.0)) { result = bt*betacf(a,b,x)/a; }
59 else { result = 1.0-bt*betacf(b,a,1.0-x)/b; }
64 m->errorOut(e, "LinearAlgebra", "betai");
68 /*********************************************************************************************************************************/
69 //Numerical Recipes pg. 219
70 double LinearAlgebra::gammln(const double xx) {
74 static const double cof[6]={76.18009172947146,-86.50532032941677,24.01409824083091,
75 -1.231739572450155,0.120858003e-2,-0.536382e-5};
79 tmp -= (x+0.5)*log(tmp);
84 return -tmp+log(2.5066282746310005*ser/x);
87 m->errorOut(e, "LinearAlgebra", "gammln");
91 /*********************************************************************************************************************************/
92 //Numerical Recipes pg. 223
93 double LinearAlgebra::gammp(const double a, const double x) {
95 double gamser,gammcf,gln;
97 if (x < 0.0 || a <= 0.0) { m->mothurOut("[ERROR]: Invalid arguments in routine GAMMP\n"); m->control_pressed = true; return 0.0;}
108 catch(exception& e) {
109 m->errorOut(e, "LinearAlgebra", "gammp");
113 /*********************************************************************************************************************************/
114 //Numerical Recipes pg. 223
115 double LinearAlgebra::gammq(const double a, const double x) {
117 double gamser,gammcf,gln;
119 if (x < 0.0 || a <= 0.0) { m->mothurOut("[ERROR]: Invalid arguments in routine GAMMQ\n"); m->control_pressed = true; return 0.0; }
121 gser(gamser,a,x,gln);
130 catch(exception& e) {
131 m->errorOut(e, "LinearAlgebra", "gammp");
135 /*********************************************************************************************************************************/
136 //Numerical Recipes pg. 224
137 double LinearAlgebra::gcf(double& gammcf, const double a, const double x, double& gln){
140 const double EPS=numeric_limits<double>::epsilon();
141 const double FPMIN=numeric_limits<double>::min()/EPS;
143 double an,b,c,d,del,h;
150 for (i=1;i<=ITMAX;i++) {
154 if (fabs(d) < FPMIN) { d=FPMIN; }
156 if (fabs(c) < FPMIN) { c=FPMIN; }
160 if (fabs(del-1.0) <= EPS) break;
162 if (i > ITMAX) { m->mothurOut("[ERROR]: a too large, ITMAX too small in gcf\n"); m->control_pressed = true; }
163 gammcf=exp(-x+a*log(x)-gln)*h;
167 catch(exception& e) {
168 m->errorOut(e, "LinearAlgebra", "gcf");
173 /*********************************************************************************************************************************/
174 //Numerical Recipes pg. 223
175 double LinearAlgebra::gser(double& gamser, const double a, const double x, double& gln) {
179 const double EPS = numeric_limits<double>::epsilon();
183 if (x < 0.0) { m->mothurOut("[ERROR]: x less than 0 in routine GSER\n"); m->control_pressed = true; }
184 gamser=0.0; return 0.0;
188 for (n=0;n<100;n++) {
192 if (fabs(del) < fabs(sum)*EPS) {
193 gamser=sum*exp(-x+a*log(x)-gln);
198 m->mothurOut("[ERROR]: a too large, ITMAX too small in routine GSER\n");
203 catch(exception& e) {
204 m->errorOut(e, "LinearAlgebra", "gser");
208 /*********************************************************************************************************************************/
209 //Numerical Recipes pg. 233
210 double LinearAlgebra::betacf(const double a, const double b, const double x) {
212 const int MAXIT = 100;
213 const double EPS = numeric_limits<double>::epsilon();
214 const double FPMIN = numeric_limits<double>::min() / EPS;
216 double aa, c, d, del, h, qab, qam, qap;
223 if (fabs(d) < FPMIN) d=FPMIN;
226 for (m1=1;m1<=MAXIT;m1++) {
228 aa=m1*(b-m1)*x/((qam+m2)*(a+m2));
230 if (fabs(d) < FPMIN) d=FPMIN;
232 if (fabs(c) < FPMIN) c=FPMIN;
235 aa = -(a+m1)*(qab+m1)*x/((a+m2)*(qap+m2));
237 if (fabs(d) < FPMIN) d=FPMIN;
239 if (fabs(c) < FPMIN) c=FPMIN;
243 if (fabs(del-1.0) < EPS) break;
246 if (m1 > MAXIT) { m->mothurOut("[ERROR]: a or b too big or MAXIT too small in betacf."); m->mothurOutEndLine(); m->control_pressed = true; }
250 catch(exception& e) {
251 m->errorOut(e, "LinearAlgebra", "betacf");
255 /*********************************************************************************************************************************/
257 vector<vector<double> > LinearAlgebra::matrix_mult(vector<vector<double> > first, vector<vector<double> > second){
259 vector<vector<double> > product;
261 int first_rows = first.size();
262 int first_cols = first[0].size();
263 int second_cols = second[0].size();
265 product.resize(first_rows);
266 for(int i=0;i<first_rows;i++){
267 product[i].resize(second_cols);
270 for(int i=0;i<first_rows;i++){
271 for(int j=0;j<second_cols;j++){
273 if (m->control_pressed) { return product; }
276 for(int k=0;k<first_cols;k++){
277 product[i][j] += first[i][k] * second[k][j];
284 catch(exception& e) {
285 m->errorOut(e, "LinearAlgebra", "matrix_mult");
291 /*********************************************************************************************************************************/
293 void LinearAlgebra::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
297 vector<vector<double> > A(rank);
298 vector<vector<double> > C(rank);
299 for(int i=0;i<rank;i++){
304 double scale = -1.0000 / (double) rank;
306 for(int i=0;i<rank;i++){
308 C[i][i] = 1.0000 + scale;
309 for(int j=i+1;j<rank;j++){
310 A[i][j] = A[j][i] = -0.5 * D[i][j] * D[i][j] + offset;
311 C[i][j] = C[j][i] = scale;
315 A = matrix_mult(C,A);
316 G = matrix_mult(A,C);
318 catch(exception& e) {
319 m->errorOut(e, "LinearAlgebra", "recenter");
324 /*********************************************************************************************************************************/
326 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
328 int LinearAlgebra::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
330 double scale, hh, h, g, f;
337 for(int i=n-1;i>0;i--){
341 for(int k=0;k<l+1;k++){
342 scale += fabs(a[i][k]);
348 for(int k=0;k<l+1;k++){
350 h += a[i][k] * a[i][k];
353 g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
358 for(int j=0;j<l+1;j++){
359 a[j][i] = a[i][j] / h;
361 for(int k=0;k<j+1;k++){
362 g += a[j][k] * a[i][k];
364 for(int k=j+1;k<l+1;k++){
365 g += a[k][j] * a[i][k];
371 for(int j=0;j<l+1;j++){
373 e[j] = g = e[j] - hh * f;
374 for(int k=0;k<j+1;k++){
375 a[j][k] -= (f * e[k] + g * a[i][k]);
390 for(int i=0;i<n;i++){
393 for(int j=0;j<l;j++){
395 for(int k=0;k<l;k++){
396 g += a[i][k] * a[k][j];
398 for(int k=0;k<l;k++){
399 a[k][j] -= g * a[k][i];
405 for(int j=0;j<l;j++){
406 a[j][i] = a[i][j] = 0.0;
412 catch(exception& e) {
413 m->errorOut(e, "LinearAlgebra", "tred2");
418 /*********************************************************************************************************************************/
420 double LinearAlgebra::pythag(double a, double b) { return(pow(a*a+b*b,0.5)); }
422 /*********************************************************************************************************************************/
424 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
426 int LinearAlgebra::qtli(vector<double>& d, vector<double>& e, vector<vector<double> >& z) {
429 double s, r, p, g, f, dd, c, b;
432 for(int i=1;i<=n;i++){
437 for(int l=0;l<n;l++){
440 for(myM=l;myM<n-1;myM++){
441 dd = fabs(d[myM]) + fabs(d[myM+1]);
442 if(fabs(e[myM])+dd == dd) break;
445 if(iter++ == 3000) cerr << "Too many iterations in tqli\n";
446 g = (d[l+1]-d[l]) / (2.0 * e[l]);
448 g = d[myM] - d[l] + e[l] / (g + SIGN(r,g));
451 for(i=myM-1;i>=l;i--){
454 e[i+1] = (r=pythag(f,g));
463 r = (d[i] - g) * s + 2.0 * c * b;
464 d[i+1] = g + ( p = s * r);
466 for(int k=0;k<n;k++){
468 z[k][i+1] = s * z[k][i] + c * f;
469 z[k][i] = c * z[k][i] - s * f;
472 if(r == 0.00 && i >= l) continue;
481 for(int i=0;i<n;i++){
483 for(int j=i;j<n;j++){
491 for(int j=0;j<n;j++){
501 catch(exception& e) {
502 m->errorOut(e, "LinearAlgebra", "qtli");
506 /*********************************************************************************************************************************/
507 //groups by dimension
508 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes, int dimensions){
511 vector< vector<double> > dists; dists.resize(axes.size());
512 for (int i = 0; i < dists.size(); i++) { dists[i].resize(axes.size(), 0.0); }
514 if (dimensions == 1) { //one dimension calc = abs(x-y)
516 for (int i = 0; i < dists.size(); i++) {
518 if (m->control_pressed) { return dists; }
520 for (int j = 0; j < i; j++) {
521 dists[i][j] = abs(axes[i][0] - axes[j][0]);
522 dists[j][i] = dists[i][j];
526 }else if (dimensions > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
528 for (int i = 0; i < dists.size(); i++) {
530 if (m->control_pressed) { return dists; }
532 for (int j = 0; j < i; j++) {
534 for (int k = 0; k < dimensions; k++) {
535 sum += ((axes[i][k] - axes[j][k]) * (axes[i][k] - axes[j][k]));
538 dists[i][j] = sqrt(sum);
539 dists[j][i] = dists[i][j];
547 catch(exception& e) {
548 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
552 /*********************************************************************************************************************************/
553 //returns groups by dimensions from dimensions by groups
554 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes){
557 vector< vector<double> > dists; dists.resize(axes[0].size());
558 for (int i = 0; i < dists.size(); i++) { dists[i].resize(axes[0].size(), 0.0); }
560 if (axes.size() == 1) { //one dimension calc = abs(x-y)
562 for (int i = 0; i < dists.size(); i++) {
564 if (m->control_pressed) { return dists; }
566 for (int j = 0; j < i; j++) {
567 dists[i][j] = abs(axes[0][i] - axes[0][j]);
568 dists[j][i] = dists[i][j];
572 }else if (axes.size() > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
574 for (int i = 0; i < dists[0].size(); i++) {
576 if (m->control_pressed) { return dists; }
578 for (int j = 0; j < i; j++) {
580 for (int k = 0; k < axes.size(); k++) {
581 sum += ((axes[k][i] - axes[k][j]) * (axes[k][i] - axes[k][j]));
584 dists[i][j] = sqrt(sum);
585 dists[j][i] = dists[i][j];
593 catch(exception& e) {
594 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
598 /*********************************************************************************************************************************/
599 //assumes both matrices are square and the same size
600 double LinearAlgebra::calcPearson(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
603 //find average for - X
605 float averageEuclid = 0.0;
606 for (int i = 0; i < euclidDists.size(); i++) {
607 for (int j = 0; j < i; j++) {
608 averageEuclid += euclidDists[i][j];
612 averageEuclid = averageEuclid / (float) count;
614 //find average for - Y
616 float averageUser = 0.0;
617 for (int i = 0; i < userDists.size(); i++) {
618 for (int j = 0; j < i; j++) {
619 averageUser += userDists[i][j];
623 averageUser = averageUser / (float) count;
625 double numerator = 0.0;
626 double denomTerm1 = 0.0;
627 double denomTerm2 = 0.0;
629 for (int i = 0; i < euclidDists.size(); i++) {
631 for (int k = 0; k < i; k++) { //just lt dists
633 float Yi = userDists[i][k];
634 float Xi = euclidDists[i][k];
636 numerator += ((Xi - averageEuclid) * (Yi - averageUser));
637 denomTerm1 += ((Xi - averageEuclid) * (Xi - averageEuclid));
638 denomTerm2 += ((Yi - averageUser) * (Yi - averageUser));
642 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
643 double r = numerator / denom;
645 //divide by zero error
646 if (isnan(r) || isinf(r)) { r = 0.0; }
651 catch(exception& e) {
652 m->errorOut(e, "LinearAlgebra", "calcPearson");
656 /*********************************************************************************************************************************/
657 //assumes both matrices are square and the same size
658 double LinearAlgebra::calcSpearman(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
663 map<float, int> tableX;
664 map<float, int>::iterator itTable;
665 vector<spearmanRank> scores;
667 for (int i = 0; i < euclidDists.size(); i++) {
668 for (int j = 0; j < i; j++) {
669 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
670 scores.push_back(member);
672 //count number of repeats
673 itTable = tableX.find(euclidDists[i][j]);
674 if (itTable == tableX.end()) {
675 tableX[euclidDists[i][j]] = 1;
677 tableX[euclidDists[i][j]]++;
683 sort(scores.begin(), scores.end(), compareSpearman);
687 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
688 double tx = (double) itTable->second;
689 Lx += ((pow(tx, 3.0) - tx) / 12.0);
693 map<string, float> rankEuclid;
694 vector<spearmanRank> ties;
696 for (int j = 0; j < scores.size(); j++) {
698 ties.push_back(scores[j]);
700 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
701 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
703 for (int k = 0; k < ties.size(); k++) {
704 float thisrank = rankTotal / (float) ties.size();
705 rankEuclid[ties[k].name] = thisrank;
710 }else { // you are the last one
712 for (int k = 0; k < ties.size(); k++) {
713 float thisrank = rankTotal / (float) ties.size();
714 rankEuclid[ties[k].name] = thisrank;
721 map<float, int> tableY;
724 for (int i = 0; i < userDists.size(); i++) {
725 for (int j = 0; j < i; j++) {
726 spearmanRank member(toString(scores.size()), userDists[i][j]);
727 scores.push_back(member);
729 //count number of repeats
730 itTable = tableY.find(userDists[i][j]);
731 if (itTable == tableY.end()) {
732 tableY[userDists[i][j]] = 1;
734 tableY[userDists[i][j]]++;
740 sort(scores.begin(), scores.end(), compareSpearman);
744 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
745 double ty = (double) itTable->second;
746 Ly += ((pow(ty, 3.0) - ty) / 12.0);
750 map<string, float> rankUser;
753 for (int j = 0; j < scores.size(); j++) {
755 ties.push_back(scores[j]);
757 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
758 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
760 for (int k = 0; k < ties.size(); k++) {
761 float thisrank = rankTotal / (float) ties.size();
762 rankUser[ties[k].name] = thisrank;
767 }else { // you are the last one
769 for (int k = 0; k < ties.size(); k++) {
770 float thisrank = rankTotal / (float) ties.size();
771 rankUser[ties[k].name] = thisrank;
779 for (int i = 0; i < userDists.size(); i++) {
780 for (int j = 0; j < i; j++) {
782 float xi = rankEuclid[toString(count)];
783 float yi = rankUser[toString(count)];
785 di += ((xi - yi) * (xi - yi));
791 double n = (double) count;
793 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
794 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
796 r = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
798 //divide by zero error
799 if (isnan(r) || isinf(r)) { r = 0.0; }
804 catch(exception& e) {
805 m->errorOut(e, "LinearAlgebra", "calcSpearman");
810 /*********************************************************************************************************************************/
811 //assumes both matrices are square and the same size
812 double LinearAlgebra::calcKendall(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
817 vector<spearmanRank> scores;
818 for (int i = 0; i < euclidDists.size(); i++) {
819 for (int j = 0; j < i; j++) {
820 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
821 scores.push_back(member);
826 sort(scores.begin(), scores.end(), compareSpearman);
829 map<string, float> rankEuclid;
830 vector<spearmanRank> ties;
832 for (int j = 0; j < scores.size(); j++) {
834 ties.push_back(scores[j]);
836 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
837 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
839 for (int k = 0; k < ties.size(); k++) {
840 float thisrank = rankTotal / (float) ties.size();
841 rankEuclid[ties[k].name] = thisrank;
846 }else { // you are the last one
848 for (int k = 0; k < ties.size(); k++) {
849 float thisrank = rankTotal / (float) ties.size();
850 rankEuclid[ties[k].name] = thisrank;
855 vector<spearmanRank> scoresUser;
856 for (int i = 0; i < userDists.size(); i++) {
857 for (int j = 0; j < i; j++) {
858 spearmanRank member(toString(scoresUser.size()), userDists[i][j]);
859 scoresUser.push_back(member);
864 sort(scoresUser.begin(), scoresUser.end(), compareSpearman);
867 map<string, float> rankUser;
870 for (int j = 0; j < scoresUser.size(); j++) {
872 ties.push_back(scoresUser[j]);
874 if (j != (scoresUser.size()-1)) { // you are not the last so you can look ahead
875 if (scoresUser[j].score != scoresUser[j+1].score) { // you are done with ties, rank them and continue
877 for (int k = 0; k < ties.size(); k++) {
878 float thisrank = rankTotal / (float) ties.size();
879 rankUser[ties[k].name] = thisrank;
884 }else { // you are the last one
886 for (int k = 0; k < ties.size(); k++) {
887 float thisrank = rankTotal / (float) ties.size();
888 rankUser[ties[k].name] = thisrank;
897 vector<spearmanRank> user;
898 for (int l = 0; l < scores.size(); l++) {
899 spearmanRank member(scores[l].name, rankUser[scores[l].name]);
900 user.push_back(member);
904 for (int l = 0; l < scores.size(); l++) {
906 int numWithHigherRank = 0;
907 int numWithLowerRank = 0;
908 float thisrank = user[l].score;
910 for (int u = l+1; u < scores.size(); u++) {
911 if (user[u].score > thisrank) { numWithHigherRank++; }
912 else if (user[u].score < thisrank) { numWithLowerRank++; }
916 numCoor += numWithHigherRank;
917 numDisCoor += numWithLowerRank;
920 r = (numCoor - numDisCoor) / (float) count;
922 //divide by zero error
923 if (isnan(r) || isinf(r)) { r = 0.0; }
928 catch(exception& e) {
929 m->errorOut(e, "LinearAlgebra", "calcKendall");
933 /*********************************************************************************************************************************/
934 double LinearAlgebra::calcKendall(vector<double>& x, vector<double>& y, double& sig){
936 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
939 vector<spearmanRank> xscores;
940 for (int i = 0; i < x.size(); i++) {
941 spearmanRank member(toString(i), x[i]);
942 xscores.push_back(member);
946 sort(xscores.begin(), xscores.end(), compareSpearman);
948 //convert scores to ranks of x
949 vector<spearmanRank*> ties;
951 for (int j = 0; j < xscores.size(); j++) {
953 ties.push_back(&(xscores[j]));
955 if (j != xscores.size()-1) { // you are not the last so you can look ahead
956 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
957 for (int k = 0; k < ties.size(); k++) {
958 float thisrank = rankTotal / (float) ties.size();
959 (*ties[k]).score = thisrank;
964 }else { // you are the last one
965 for (int k = 0; k < ties.size(); k++) {
966 float thisrank = rankTotal / (float) ties.size();
967 (*ties[k]).score = thisrank;
974 vector<spearmanRank> yscores;
975 for (int j = 0; j < y.size(); j++) {
976 spearmanRank member(toString(j), y[j]);
977 yscores.push_back(member);
981 sort(yscores.begin(), yscores.end(), compareSpearman);
984 map<string, float> rank;
985 vector<spearmanRank> yties;
987 for (int j = 0; j < yscores.size(); j++) {
989 yties.push_back(yscores[j]);
991 if (j != yscores.size()-1) { // you are not the last so you can look ahead
992 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
993 for (int k = 0; k < yties.size(); k++) {
994 float thisrank = rankTotal / (float) yties.size();
995 rank[yties[k].name] = thisrank;
1000 }else { // you are the last one
1001 for (int k = 0; k < yties.size(); k++) {
1002 float thisrank = rankTotal / (float) yties.size();
1003 rank[yties[k].name] = thisrank;
1013 vector<spearmanRank> otus;
1014 for (int l = 0; l < xscores.size(); l++) {
1015 spearmanRank member(xscores[l].name, rank[xscores[l].name]);
1016 otus.push_back(member);
1020 for (int l = 0; l < xscores.size(); l++) {
1022 int numWithHigherRank = 0;
1023 int numWithLowerRank = 0;
1024 float thisrank = otus[l].score;
1026 for (int u = l+1; u < xscores.size(); u++) {
1027 if (otus[u].score > thisrank) { numWithHigherRank++; }
1028 else if (otus[u].score < thisrank) { numWithLowerRank++; }
1032 numCoor += numWithHigherRank;
1033 numDisCoor += numWithLowerRank;
1036 double p = (numCoor - numDisCoor) / (float) count;
1038 sig = calcKendallSig(x.size(), p);
1042 catch(exception& e) {
1043 m->errorOut(e, "LinearAlgebra", "calcKendall");
1047 double LinearAlgebra::ran0(int& idum)
1049 const int IA=16807,IM=2147483647,IQ=127773;
1050 const int IR=2836,MASK=123459876;
1051 const double AM=1.0/double(IM);
1057 idum=IA*(idum-k*IQ)-IR*k;
1058 if (idum < 0) idum += IM;
1064 double LinearAlgebra::ran1(int &idum)
1066 const int IA=16807,IM=2147483647,IQ=127773,IR=2836,NTAB=32;
1067 const int NDIV=(1+(IM-1)/NTAB);
1068 const double EPS=3.0e-16,AM=1.0/IM,RNMX=(1.0-EPS);
1070 static vector<int> iv(NTAB);
1074 if (idum <= 0 || !iy) {
1075 if (-idum < 1) idum=1;
1077 for (j=NTAB+7;j>=0;j--) {
1079 idum=IA*(idum-k*IQ)-IR*k;
1080 if (idum < 0) idum += IM;
1081 if (j < NTAB) iv[j] = idum;
1086 idum=IA*(idum-k*IQ)-IR*k;
1087 if (idum < 0) idum += IM;
1091 if ((temp=AM*iy) > RNMX) return RNMX;
1095 double LinearAlgebra::ran2(int &idum)
1097 const int IM1=2147483563,IM2=2147483399;
1098 const int IA1=40014,IA2=40692,IQ1=53668,IQ2=52774;
1099 const int IR1=12211,IR2=3791,NTAB=32,IMM1=IM1-1;
1100 const int NDIV=1+IMM1/NTAB;
1101 const double EPS=3.0e-16,RNMX=1.0-EPS,AM=1.0/double(IM1);
1102 static int idum2=123456789,iy=0;
1103 static vector<int> iv(NTAB);
1108 idum=(idum==0 ? 1 : -idum);
1110 for (j=NTAB+7;j>=0;j--) {
1112 idum=IA1*(idum-k*IQ1)-k*IR1;
1113 if (idum < 0) idum += IM1;
1114 if (j < NTAB) iv[j] = idum;
1119 idum=IA1*(idum-k*IQ1)-k*IR1;
1120 if (idum < 0) idum += IM1;
1122 idum2=IA2*(idum2-k*IQ2)-k*IR2;
1123 if (idum2 < 0) idum2 += IM2;
1127 if (iy < 1) iy += IMM1;
1128 if ((temp=AM*iy) > RNMX) return RNMX;
1132 double LinearAlgebra::ran3(int &idum)
1134 static int inext,inextp;
1136 const int MBIG=1000000000,MSEED=161803398,MZ=0;
1137 const double FAC=(1.0/MBIG);
1138 static vector<int> ma(56);
1141 if (idum < 0 || iff == 0) {
1143 mj=labs(MSEED-labs(idum));
1147 for (i=1;i<=54;i++) {
1151 if (mk < int(MZ)) mk += MBIG;
1155 for (i=1;i<=55;i++) {
1156 ma[i] -= ma[1+(i+30) % 55];
1157 if (ma[i] < int(MZ)) ma[i] += MBIG;
1163 if (++inext == 56) inext=1;
1164 if (++inextp == 56) inextp=1;
1165 mj=ma[inext]-ma[inextp];
1166 if (mj < int(MZ)) mj += MBIG;
1171 double LinearAlgebra::ran4(int &idum)
1173 #if defined(vax) || defined(_vax_) || defined(__vax__) || defined(VAX)
1174 static const unsigned long jflone = 0x00004080;
1175 static const unsigned long jflmsk = 0xffff007f;
1177 static const unsigned long jflone = 0x3f800000;
1178 static const unsigned long jflmsk = 0x007fffff;
1180 unsigned long irword,itemp,lword;
1181 static int idums = 0;
1189 psdes(lword,irword);
1190 itemp=jflone | (jflmsk & irword);
1192 return (*(float *)&itemp)-1.0;
1195 void LinearAlgebra::psdes(unsigned long &lword, unsigned long &irword)
1198 static const unsigned long c1[NITER]={
1199 0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L};
1200 static const unsigned long c2[NITER]={
1201 0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L};
1202 unsigned long i,ia,ib,iswap,itmph=0,itmpl=0;
1204 for (i=0;i<NITER;i++) {
1205 ia=(iswap=irword) ^ c1[i];
1206 itmpl = ia & 0xffff;
1208 ib=itmpl*itmpl+ ~(itmph*itmph);
1209 irword=lword ^ (((ia = (ib >> 16) |
1210 ((ib & 0xffff) << 16)) ^ c2[i])+itmpl*itmph);
1214 /*********************************************************************************************************************************/
1215 double LinearAlgebra::calcKendallSig(double n, double r){
1219 double svar=(4.0*n+10.0)/(9.0*n*(n-1.0));
1220 double z= r/sqrt(svar);
1221 sig=erfcc(fabs(z)/1.4142136);
1223 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1227 catch(exception& e) {
1228 m->errorOut(e, "LinearAlgebra", "calcKendallSig");
1232 /*********************************************************************************************************************************/
1233 double LinearAlgebra::calcKruskalWallis(vector<spearmanRank>& values, double& pValue){
1236 set<string> treatments;
1239 sort(values.begin(), values.end(), compareSpearman);
1240 vector<spearmanRank*> ties;
1243 for (int j = 0; j < values.size(); j++) {
1244 treatments.insert(values[j].name);
1246 ties.push_back(&(values[j]));
1248 if (j != values.size()-1) { // you are not the last so you can look ahead
1249 if (values[j].score != values[j+1].score) { // you are done with ties, rank them and continue
1250 if (ties.size() > 1) { TIES.push_back(ties.size()); }
1251 for (int k = 0; k < ties.size(); k++) {
1252 double thisrank = rankTotal / (double) ties.size();
1253 (*ties[k]).score = thisrank;
1258 }else { // you are the last one
1259 if (ties.size() > 1) { TIES.push_back(ties.size()); }
1260 for (int k = 0; k < ties.size(); k++) {
1261 double thisrank = rankTotal / (double) ties.size();
1262 (*ties[k]).score = thisrank;
1268 // H = 12/(N*(N+1)) * (sum Ti^2/n) - 3(N+1)
1269 map<string, double> sums;
1270 map<string, double> counts;
1271 for (set<string>::iterator it = treatments.begin(); it != treatments.end(); it++) { sums[*it] = 0.0; counts[*it] = 0; }
1273 for (int j = 0; j < values.size(); j++) {
1274 sums[values[j].name] += values[j].score;
1275 counts[values[j].name]+= 1.0;
1278 double middleTerm = 0.0;
1279 for (set<string>::iterator it = treatments.begin(); it != treatments.end(); it++) {
1280 middleTerm += ((sums[*it]*sums[*it])/counts[*it]);
1283 double firstTerm = 12 / (double) (values.size()*(values.size()+1));
1284 double lastTerm = 3 * (values.size()+1);
1286 H = firstTerm * middleTerm - lastTerm;
1289 if (TIES.size() != 0) {
1291 for (int j = 0; j < TIES.size(); j++) { sum += ((TIES[j]*TIES[j]*TIES[j])-TIES[j]); }
1292 double result = 1.0 - (sum / (double) ((values.size()*values.size()*values.size())-values.size()));
1296 //Numerical Recipes pg221
1297 pValue = 1.0 - (gammp(((treatments.size()-1)/(double)2.0), H/2.0));
1301 catch(exception& e) {
1302 m->errorOut(e, "LinearAlgebra", "calcKruskalWallis");
1306 /*********************************************************************************************************************************/
1307 //thanks http://www.johndcook.com/cpp_phi.html
1308 double LinearAlgebra::pnorm(double x){
1311 double a1 = 0.254829592;
1312 double a2 = -0.284496736;
1313 double a3 = 1.421413741;
1314 double a4 = -1.453152027;
1315 double a5 = 1.061405429;
1316 double p = 0.3275911;
1318 // Save the sign of x
1322 x = fabs(x)/sqrt(2.0);
1324 // A&S formula 7.1.26
1325 double t = 1.0/(1.0 + p*x);
1326 double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);
1328 return 0.5*(1.0 + sign*y);
1330 catch(exception& e) {
1331 m->errorOut(e, "LinearAlgebra", "pnorm");
1336 /*********************************************************************************************************************************/
1337 double LinearAlgebra::calcWilcoxon(vector<double>& x, vector<double>& y, double& sig){
1342 vector<spearmanRank> ranks;
1343 for (int i = 0; i < x.size(); i++) {
1344 if (m->control_pressed) { return W; }
1345 spearmanRank member("x", x[i]);
1346 ranks.push_back(member);
1349 for (int i = 0; i < y.size(); i++) {
1350 if (m->control_pressed) { return W; }
1351 spearmanRank member("y", y[i]);
1352 ranks.push_back(member);
1356 sort(ranks.begin(), ranks.end(), compareSpearman);
1358 //convert scores to ranks of x
1359 vector<spearmanRank*> ties;
1362 for (int j = 0; j < ranks.size(); j++) {
1363 if (m->control_pressed) { return W; }
1365 ties.push_back(&(ranks[j]));
1367 if (j != ranks.size()-1) { // you are not the last so you can look ahead
1368 if (ranks[j].score != ranks[j+1].score) { // you are done with ties, rank them and continue
1369 if (ties.size() > 1) { TIES.push_back(ties.size()); }
1370 for (int k = 0; k < ties.size(); k++) {
1371 float thisrank = rankTotal / (float) ties.size();
1372 (*ties[k]).score = thisrank;
1377 }else { // you are the last one
1378 if (ties.size() > 1) { TIES.push_back(ties.size()); }
1379 for (int k = 0; k < ties.size(); k++) {
1380 float thisrank = rankTotal / (float) ties.size();
1381 (*ties[k]).score = thisrank;
1386 //from R wilcox.test function
1387 //STATISTIC <- sum(r[seq_along(x)]) - n.x * (n.x + 1)/2
1388 double sumRanks = 0.0;
1389 for (int i = 0; i < ranks.size(); i++) {
1390 if (m->control_pressed) { return W; }
1391 if (ranks[i].name == "x") { sumRanks += ranks[i].score; }
1394 W = sumRanks - x.size() * ((double)(x.size() + 1)) / 2.0;
1396 //exact <- (n.x < 50) && (n.y < 50)
1397 bool findExact = false;
1398 if ((x.size() < 50) && (y.size() < 50)) { findExact = true; }
1401 if (findExact && (TIES.size() == 0)) { //find exact and no ties
1402 //PVAL <- switch(alternative, two.sided = {
1403 //p <- if (STATISTIC > (n.x * n.y/2))
1406 if (W > ((double)x.size()*y.size()/2.0)) {
1407 //pwilcox(STATISTIC-1, n.x, n.y, lower.tail = FALSE)
1408 pval = wilcox.pwilcox(W-1, x.size(), y.size(), false);
1410 //pwilcox(STATISTIC,n.x, n.y)
1411 pval = wilcox.pwilcox(W, x.size(), y.size(), true);
1414 if (1.0 < sig) { sig = 1.0; }
1416 //z <- STATISTIC - n.x * n.y/2
1417 double z = W - (double)(x.size() * y.size()/2.0);
1420 for (int j = 0; j < TIES.size(); j++) { sum += ((TIES[j]*TIES[j]*TIES[j])-TIES[j]); }
1422 //SIGMA <- sqrt((n.x * n.y/12) * ((n.x + n.y + 1) -
1423 //sum(NTIES^3 - NTIES)/((n.x + n.y) * (n.x + n.y -
1426 double firstTerm = (double)(x.size() * y.size()/12.0);
1427 double secondTerm = (double)(x.size() + y.size() + 1) - sum / (double)((x.size() + y.size()) * (x.size() + y.size() - 1));
1428 sigma = sqrt(firstTerm * secondTerm);
1430 //CORRECTION <- switch(alternative, two.sided = sign(z) * 0.5, greater = 0.5, less = -0.5)
1431 double CORRECTION = 0.0;
1432 if (z < 0) { CORRECTION = -1.0; }
1433 else if (z > 0) { CORRECTION = 1.0; }
1436 z = (z - CORRECTION)/sigma;
1438 //PVAL <- switch(alternative, two.sided = 2 * min(pnorm(z), pnorm(z, lower.tail = FALSE)))
1440 if ((1.0-sig) < sig) { sig = 1.0 - sig; }
1446 catch(exception& e) {
1447 m->errorOut(e, "LinearAlgebra", "calcWilcoxon");
1452 /*********************************************************************************************************************************/
1453 double LinearAlgebra::choose(double n, double k){
1458 double lchoose = gammln(n + 1.0) - gammln(k + 1.0) - gammln(n - k + 1.0);
1460 return (floor(exp(lchoose) + 0.5));
1462 catch(exception& e) {
1463 m->errorOut(e, "LinearAlgebra", "choose");
1467 /*********************************************************************************************************************************/
1468 double LinearAlgebra::calcSpearman(vector<double>& x, vector<double>& y, double& sig){
1470 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1473 double sf = 0.0; //f^3 - f where f is the number of ties in x;
1474 double sg = 0.0; //f^3 - f where f is the number of ties in y;
1475 map<float, int> tableX;
1476 map<float, int>::iterator itTable;
1477 vector<spearmanRank> xscores;
1479 for (int i = 0; i < x.size(); i++) {
1480 spearmanRank member(toString(i), x[i]);
1481 xscores.push_back(member);
1483 //count number of repeats
1484 itTable = tableX.find(x[i]);
1485 if (itTable == tableX.end()) {
1495 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
1496 double tx = (double) itTable->second;
1497 Lx += ((pow(tx, 3.0) - tx) / 12.0);
1502 sort(xscores.begin(), xscores.end(), compareSpearman);
1504 //convert scores to ranks of x
1506 map<string, float> rankx;
1507 vector<spearmanRank> xties;
1509 for (int j = 0; j < xscores.size(); j++) {
1511 xties.push_back(xscores[j]);
1513 if (j != xscores.size()-1) { // you are not the last so you can look ahead
1514 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
1515 for (int k = 0; k < xties.size(); k++) {
1516 float thisrank = rankTotal / (float) xties.size();
1517 rankx[xties[k].name] = thisrank;
1519 int t = xties.size();
1524 }else { // you are the last one
1525 for (int k = 0; k < xties.size(); k++) {
1526 float thisrank = rankTotal / (float) xties.size();
1527 rankx[xties[k].name] = thisrank;
1533 vector<spearmanRank> yscores;
1534 map<float, int> tableY;
1535 for (int j = 0; j < y.size(); j++) {
1536 spearmanRank member(toString(j), y[j]);
1537 yscores.push_back(member);
1539 itTable = tableY.find(member.score);
1540 if (itTable == tableY.end()) {
1541 tableY[member.score] = 1;
1543 tableY[member.score]++;
1550 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
1551 double ty = (double) itTable->second;
1552 Ly += ((pow(ty, 3.0) - ty) / 12.0);
1555 sort(yscores.begin(), yscores.end(), compareSpearman);
1558 map<string, float> rank;
1559 vector<spearmanRank> yties;
1561 for (int j = 0; j < yscores.size(); j++) {
1563 yties.push_back(yscores[j]);
1565 if (j != yscores.size()-1) { // you are not the last so you can look ahead
1566 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
1567 for (int k = 0; k < yties.size(); k++) {
1568 float thisrank = rankTotal / (float) yties.size();
1569 rank[yties[k].name] = thisrank;
1571 int t = yties.size();
1576 }else { // you are the last one
1577 for (int k = 0; k < yties.size(); k++) {
1578 float thisrank = rankTotal / (float) yties.size();
1579 rank[yties[k].name] = thisrank;
1585 for (int k = 0; k < x.size(); k++) {
1587 float xi = rankx[toString(k)];
1588 float yi = rank[toString(k)];
1590 di += ((xi - yi) * (xi - yi));
1595 double n = (double) x.size();
1596 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
1597 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
1599 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
1601 //Numerical Recipes 646
1602 sig = calcSpearmanSig(n, sf, sg, di);
1606 catch(exception& e) {
1607 m->errorOut(e, "LinearAlgebra", "calcSpearman");
1611 /*********************************************************************************************************************************/
1612 double LinearAlgebra::calcSpearmanSig(double n, double sf, double sg, double d){
1616 double probrs = 0.0;
1618 double en3n=en*en*en-en;
1619 double aved=en3n/6.0-(sf+sg)/12.0;
1620 double fac=(1.0-sf/en3n)*(1.0-sg/en3n);
1621 double vard=((en-1.0)*en*en*SQR(en+1.0)/36.0)*fac;
1622 double zd=(d-aved)/sqrt(vard);
1623 double probd=erfcc(fabs(zd)/1.4142136);
1624 double rs=(1.0-(6.0/en3n)*(d+(sf+sg)/12.0))/sqrt(fac);
1625 fac=(rs+1.0)*(1.0-rs);
1627 double t=rs*sqrt((en-2.0)/fac);
1629 probrs=betai(0.5*df,0.5,df/(df+t*t));
1634 //smaller of probd and probrs is sig
1636 if (probd < probrs) { sig = probd; }
1638 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1642 catch(exception& e) {
1643 m->errorOut(e, "LinearAlgebra", "calcSpearmanSig");
1647 /*********************************************************************************************************************************/
1648 double LinearAlgebra::calcPearson(vector<double>& x, vector<double>& y, double& sig){
1650 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1653 float averageX = 0.0;
1654 for (int i = 0; i < x.size(); i++) { averageX += x[i]; }
1655 averageX = averageX / (float) x.size();
1659 for (int j = 0; j < y.size(); j++) { sumY += y[j]; }
1660 float Ybar = sumY / (float) y.size();
1663 double numerator = 0.0;
1664 double denomTerm1 = 0.0;
1665 double denomTerm2 = 0.0;
1667 for (int j = 0; j < x.size(); j++) {
1671 numerator += ((Xi - averageX) * (Yi - Ybar));
1672 denomTerm1 += ((Xi - averageX) * (Xi - averageX));
1673 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
1676 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
1678 r = numerator / denom;
1680 //Numerical Recipes pg.644
1681 sig = calcPearsonSig(x.size(), r);
1685 catch(exception& e) {
1686 m->errorOut(e, "LinearAlgebra", "calcPearson");
1690 /*********************************************************************************************************************************/
1691 double LinearAlgebra::calcPearsonSig(double n, double r){
1695 const double TINY = 1.0e-20;
1696 double z = 0.5*log((1.0+r+TINY)/(1.0-r+TINY)); //Fisher's z transformation
1698 //code below was giving an error in betacf with sop files
1700 //double t = r*sqrt(df/((1.0-r+TINY)*(1.0+r+TINY)));
1701 //sig = betai(0.5+df, 0.5, df/(df+t*t));
1703 //Numerical Recipes says code below gives approximately the same result
1704 sig = erfcc(fabs(z*sqrt(n-1.0))/1.4142136);
1705 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1709 catch(exception& e) {
1710 m->errorOut(e, "LinearAlgebra", "calcPearsonSig");
1714 /*********************************************************************************************************************************/
1716 vector<vector<double> > LinearAlgebra::getObservedEuclideanDistance(vector<vector<double> >& relAbundData){
1719 int numSamples = relAbundData.size();
1720 int numOTUs = relAbundData[0].size();
1722 vector<vector<double> > dMatrix(numSamples);
1723 for(int i=0;i<numSamples;i++){
1724 dMatrix[i].resize(numSamples);
1727 for(int i=0;i<numSamples;i++){
1728 for(int j=0;j<numSamples;j++){
1730 if (m->control_pressed) { return dMatrix; }
1733 for(int k=0;k<numOTUs;k++){
1734 d += pow((relAbundData[i][k] - relAbundData[j][k]), 2.0000);
1736 dMatrix[i][j] = pow(d, 0.50000);
1737 dMatrix[j][i] = dMatrix[i][j];
1743 catch(exception& e) {
1744 m->errorOut(e, "LinearAlgebra", "getObservedEuclideanDistance");
1749 /*********************************************************************************************************************************/
1750 vector<double> LinearAlgebra::solveEquations(vector<vector<double> > A, vector<double> b){
1752 int length = (int)b.size();
1753 vector<double> x(length, 0);
1754 vector<int> index(length);
1755 for(int i=0;i<length;i++){ index[i] = i; }
1758 ludcmp(A, index, d); if (m->control_pressed) { return b; }
1759 lubksb(A, index, b);
1763 catch(exception& e) {
1764 m->errorOut(e, "LinearAlgebra", "solveEquations");
1768 /*********************************************************************************************************************************/
1769 vector<float> LinearAlgebra::solveEquations(vector<vector<float> > A, vector<float> b){
1771 int length = (int)b.size();
1772 vector<double> x(length, 0);
1773 vector<int> index(length);
1774 for(int i=0;i<length;i++){ index[i] = i; }
1777 ludcmp(A, index, d); if (m->control_pressed) { return b; }
1778 lubksb(A, index, b);
1782 catch(exception& e) {
1783 m->errorOut(e, "LinearAlgebra", "solveEquations");
1788 /*********************************************************************************************************************************/
1790 void LinearAlgebra::ludcmp(vector<vector<double> >& A, vector<int>& index, double& d){
1792 double tiny = 1e-20;
1794 int n = (int)A.size();
1795 vector<double> vv(n, 0.0);
1801 for(int i=0;i<n;i++){
1803 for(int j=0;j<n;j++){ if((temp=fabs(A[i][j])) > big ) big=temp; }
1804 if(big==0.0){ m->mothurOut("Singular matrix in routine ludcmp\n"); }
1808 for(int j=0;j<n;j++){
1809 if (m->control_pressed) { break; }
1810 for(int i=0;i<j;i++){
1811 double sum = A[i][j];
1812 for(int k=0;k<i;k++){ sum -= A[i][k] * A[k][j]; }
1817 for(int i=j;i<n;i++){
1818 double sum = A[i][j];
1819 for(int k=0;k<j;k++){ sum -= A[i][k] * A[k][j]; }
1822 if((dum = vv[i] * fabs(sum)) >= big){
1828 for(int k=0;k<n;k++){
1829 double dum = A[imax][k];
1830 A[imax][k] = A[j][k];
1838 if(A[j][j] == 0.0){ A[j][j] = tiny; }
1841 double dum = 1.0/A[j][j];
1842 for(int i=j+1;i<n;i++){ A[i][j] *= dum; }
1846 catch(exception& e) {
1847 m->errorOut(e, "LinearAlgebra", "ludcmp");
1852 /*********************************************************************************************************************************/
1854 void LinearAlgebra::lubksb(vector<vector<double> >& A, vector<int>& index, vector<double>& b){
1857 int n = (int)A.size();
1860 for(int i=0;i<n;i++){
1861 if (m->control_pressed) { break; }
1867 for(int j=ii-1;j<i;j++){
1868 total -= A[i][j] * b[j];
1871 else if(total != 0){ ii = i+1; }
1874 for(int i=n-1;i>=0;i--){
1876 for(int j=i+1;j<n;j++){ total -= A[i][j] * b[j]; }
1877 b[i] = total / A[i][i];
1880 catch(exception& e) {
1881 m->errorOut(e, "LinearAlgebra", "lubksb");
1885 /*********************************************************************************************************************************/
1887 void LinearAlgebra::ludcmp(vector<vector<float> >& A, vector<int>& index, float& d){
1889 double tiny = 1e-20;
1891 int n = (int)A.size();
1892 vector<float> vv(n, 0.0);
1898 for(int i=0;i<n;i++){
1900 for(int j=0;j<n;j++){ if((temp=fabs(A[i][j])) > big ) big=temp; }
1901 if(big==0.0){ m->mothurOut("Singular matrix in routine ludcmp\n"); }
1905 for(int j=0;j<n;j++){
1906 if (m->control_pressed) { break; }
1907 for(int i=0;i<j;i++){
1908 float sum = A[i][j];
1909 for(int k=0;k<i;k++){ sum -= A[i][k] * A[k][j]; }
1914 for(int i=j;i<n;i++){
1915 float sum = A[i][j];
1916 for(int k=0;k<j;k++){ sum -= A[i][k] * A[k][j]; }
1919 if((dum = vv[i] * fabs(sum)) >= big){
1925 for(int k=0;k<n;k++){
1926 float dum = A[imax][k];
1927 A[imax][k] = A[j][k];
1935 if(A[j][j] == 0.0){ A[j][j] = tiny; }
1938 float dum = 1.0/A[j][j];
1939 for(int i=j+1;i<n;i++){ A[i][j] *= dum; }
1943 catch(exception& e) {
1944 m->errorOut(e, "LinearAlgebra", "ludcmp");
1949 /*********************************************************************************************************************************/
1951 void LinearAlgebra::lubksb(vector<vector<float> >& A, vector<int>& index, vector<float>& b){
1954 int n = (int)A.size();
1957 for(int i=0;i<n;i++){
1958 if (m->control_pressed) { break; }
1964 for(int j=ii-1;j<i;j++){
1965 total -= A[i][j] * b[j];
1968 else if(total != 0){ ii = i+1; }
1971 for(int i=n-1;i>=0;i--){
1973 for(int j=i+1;j<n;j++){ total -= A[i][j] * b[j]; }
1974 b[i] = total / A[i][i];
1977 catch(exception& e) {
1978 m->errorOut(e, "LinearAlgebra", "lubksb");
1983 /*********************************************************************************************************************************/
1985 vector<vector<double> > LinearAlgebra::getInverse(vector<vector<double> > matrix){
1987 int n = (int)matrix.size();
1989 vector<vector<double> > inverse(n);
1990 for(int i=0;i<n;i++){ inverse[i].assign(n, 0.0000); }
1992 vector<double> column(n, 0.0000);
1993 vector<int> index(n, 0);
1996 ludcmp(matrix, index, dummy);
1998 for(int j=0;j<n;j++){
1999 if (m->control_pressed) { break; }
2001 column.assign(n, 0);
2005 lubksb(matrix, index, column);
2007 for(int i=0;i<n;i++){ inverse[i][j] = column[i]; }
2012 catch(exception& e) {
2013 m->errorOut(e, "LinearAlgebra", "getInverse");