5 * Created by westcott on 1/7/11.
6 * Copyright 2011 Schloss Lab. All rights reserved.
10 #include "linearalgebra.h"
12 // This class references functions used from "Numerical Recipes in C++" //
14 /*********************************************************************************************************************************/
15 inline double SQR(const double a)
19 /*********************************************************************************************************************************/
21 inline double SIGN(const double a, const double b)
23 return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a);
25 /*********************************************************************************************************************************/
26 //NUmerical recipes pg. 245 - Returns the complementary error function erfc(x) with fractional error everywhere less than 1.2 × 10−7.
27 double LinearAlgebra::erfcc(double x){
33 ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
34 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
35 t*(-0.82215223+t*0.17087277)))))))));
37 //cout << "in erfcc " << t << '\t' << ans<< endl;
38 return (x >= 0.0 ? ans : 2.0 - ans);
41 m->errorOut(e, "LinearAlgebra", "betai");
45 /*********************************************************************************************************************************/
46 //Numerical Recipes pg. 232
47 double LinearAlgebra::betai(const double a, const double b, const double x) {
52 if (x < 0.0 || x > 1.0) { m->mothurOut("[ERROR]: bad x in betai.\n"); m->control_pressed = true; return 0.0; }
54 if (x == 0.0 || x == 1.0) { bt = 0.0; }
55 else { bt = exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x)); }
57 if (x < (a+1.0) / (a+b+2.0)) { result = bt*betacf(a,b,x)/a; }
58 else { result = 1.0-bt*betacf(b,a,1.0-x)/b; }
63 m->errorOut(e, "LinearAlgebra", "betai");
67 /*********************************************************************************************************************************/
68 //Numerical Recipes pg. 219
69 double LinearAlgebra::gammln(const double xx) {
73 static const double cof[6]={76.18009172947146,-86.50532032941677,24.01409824083091,
74 -1.231739572450155,0.120858003e-2,-0.536382e-5};
78 tmp -= (x+0.5)*log(tmp);
83 return -tmp+log(2.5066282746310005*ser/x);
86 m->errorOut(e, "LinearAlgebra", "gammln");
90 /*********************************************************************************************************************************/
91 //Numerical Recipes pg. 223
92 double LinearAlgebra::gammp(const double a, const double x) {
94 double gamser,gammcf,gln;
96 if (x < 0.0 || a <= 0.0) { m->mothurOut("[ERROR]: Invalid arguments in routine GAMMP\n"); m->control_pressed = true; return 0.0;}
107 catch(exception& e) {
108 m->errorOut(e, "LinearAlgebra", "gammp");
112 /*********************************************************************************************************************************/
113 //Numerical Recipes pg. 223
114 double LinearAlgebra::gammq(const double a, const double x) {
116 double gamser,gammcf,gln;
118 if (x < 0.0 || a <= 0.0) { m->mothurOut("[ERROR]: Invalid arguments in routine GAMMQ\n"); m->control_pressed = true; return 0.0; }
120 gser(gamser,a,x,gln);
129 catch(exception& e) {
130 m->errorOut(e, "LinearAlgebra", "gammp");
134 /*********************************************************************************************************************************/
135 //Numerical Recipes pg. 224
136 double LinearAlgebra::gcf(double& gammcf, const double a, const double x, double& gln){
139 const double EPS=numeric_limits<double>::epsilon();
140 const double FPMIN=numeric_limits<double>::min()/EPS;
142 double an,b,c,d,del,h;
149 for (i=1;i<=ITMAX;i++) {
153 if (fabs(d) < FPMIN) { d=FPMIN; }
155 if (fabs(c) < FPMIN) { c=FPMIN; }
159 if (fabs(del-1.0) <= EPS) break;
161 if (i > ITMAX) { m->mothurOut("[ERROR]: a too large, ITMAX too small in gcf\n"); m->control_pressed = true; }
162 gammcf=exp(-x+a*log(x)-gln)*h;
166 catch(exception& e) {
167 m->errorOut(e, "LinearAlgebra", "gcf");
172 /*********************************************************************************************************************************/
173 //Numerical Recipes pg. 223
174 double LinearAlgebra::gser(double& gamser, const double a, const double x, double& gln) {
178 const double EPS = numeric_limits<double>::epsilon();
182 if (x < 0.0) { m->mothurOut("[ERROR]: x less than 0 in routine GSER\n"); m->control_pressed = true; }
183 gamser=0.0; return 0.0;
187 for (n=0;n<100;n++) {
191 if (fabs(del) < fabs(sum)*EPS) {
192 gamser=sum*exp(-x+a*log(x)-gln);
197 m->mothurOut("[ERROR]: a too large, ITMAX too small in routine GSER\n");
202 catch(exception& e) {
203 m->errorOut(e, "LinearAlgebra", "gser");
207 /*********************************************************************************************************************************/
208 //Numerical Recipes pg. 233
209 double LinearAlgebra::betacf(const double a, const double b, const double x) {
211 const int MAXIT = 100;
212 const double EPS = numeric_limits<double>::epsilon();
213 const double FPMIN = numeric_limits<double>::min() / EPS;
215 double aa, c, d, del, h, qab, qam, qap;
222 if (fabs(d) < FPMIN) d=FPMIN;
225 for (m1=1;m1<=MAXIT;m1++) {
227 aa=m1*(b-m1)*x/((qam+m2)*(a+m2));
229 if (fabs(d) < FPMIN) d=FPMIN;
231 if (fabs(c) < FPMIN) c=FPMIN;
234 aa = -(a+m1)*(qab+m1)*x/((a+m2)*(qap+m2));
236 if (fabs(d) < FPMIN) d=FPMIN;
238 if (fabs(c) < FPMIN) c=FPMIN;
242 if (fabs(del-1.0) < EPS) break;
245 if (m1 > MAXIT) { m->mothurOut("[ERROR]: a or b too big or MAXIT too small in betacf."); m->mothurOutEndLine(); m->control_pressed = true; }
249 catch(exception& e) {
250 m->errorOut(e, "LinearAlgebra", "betacf");
254 /*********************************************************************************************************************************/
256 vector<vector<double> > LinearAlgebra::matrix_mult(vector<vector<double> > first, vector<vector<double> > second){
258 vector<vector<double> > product;
260 int first_rows = first.size();
261 int first_cols = first[0].size();
262 int second_cols = second[0].size();
264 product.resize(first_rows);
265 for(int i=0;i<first_rows;i++){
266 product[i].resize(second_cols);
269 for(int i=0;i<first_rows;i++){
270 for(int j=0;j<second_cols;j++){
272 if (m->control_pressed) { return product; }
275 for(int k=0;k<first_cols;k++){
276 product[i][j] += first[i][k] * second[k][j];
283 catch(exception& e) {
284 m->errorOut(e, "LinearAlgebra", "matrix_mult");
290 /*********************************************************************************************************************************/
292 void LinearAlgebra::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
296 vector<vector<double> > A(rank);
297 vector<vector<double> > C(rank);
298 for(int i=0;i<rank;i++){
303 double scale = -1.0000 / (double) rank;
305 for(int i=0;i<rank;i++){
307 C[i][i] = 1.0000 + scale;
308 for(int j=i+1;j<rank;j++){
309 A[i][j] = A[j][i] = -0.5 * D[i][j] * D[i][j] + offset;
310 C[i][j] = C[j][i] = scale;
314 A = matrix_mult(C,A);
315 G = matrix_mult(A,C);
317 catch(exception& e) {
318 m->errorOut(e, "LinearAlgebra", "recenter");
323 /*********************************************************************************************************************************/
325 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
327 int LinearAlgebra::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
329 double scale, hh, h, g, f;
336 for(int i=n-1;i>0;i--){
340 for(int k=0;k<l+1;k++){
341 scale += fabs(a[i][k]);
347 for(int k=0;k<l+1;k++){
349 h += a[i][k] * a[i][k];
352 g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
357 for(int j=0;j<l+1;j++){
358 a[j][i] = a[i][j] / h;
360 for(int k=0;k<j+1;k++){
361 g += a[j][k] * a[i][k];
363 for(int k=j+1;k<l+1;k++){
364 g += a[k][j] * a[i][k];
370 for(int j=0;j<l+1;j++){
372 e[j] = g = e[j] - hh * f;
373 for(int k=0;k<j+1;k++){
374 a[j][k] -= (f * e[k] + g * a[i][k]);
389 for(int i=0;i<n;i++){
392 for(int j=0;j<l;j++){
394 for(int k=0;k<l;k++){
395 g += a[i][k] * a[k][j];
397 for(int k=0;k<l;k++){
398 a[k][j] -= g * a[k][i];
404 for(int j=0;j<l;j++){
405 a[j][i] = a[i][j] = 0.0;
411 catch(exception& e) {
412 m->errorOut(e, "LinearAlgebra", "tred2");
417 /*********************************************************************************************************************************/
419 double LinearAlgebra::pythag(double a, double b) { return(pow(a*a+b*b,0.5)); }
421 /*********************************************************************************************************************************/
423 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
425 int LinearAlgebra::qtli(vector<double>& d, vector<double>& e, vector<vector<double> >& z) {
428 double s, r, p, g, f, dd, c, b;
431 for(int i=1;i<=n;i++){
436 for(int l=0;l<n;l++){
439 for(myM=l;myM<n-1;myM++){
440 dd = fabs(d[myM]) + fabs(d[myM+1]);
441 if(fabs(e[myM])+dd == dd) break;
444 if(iter++ == 3000) cerr << "Too many iterations in tqli\n";
445 g = (d[l+1]-d[l]) / (2.0 * e[l]);
447 g = d[myM] - d[l] + e[l] / (g + SIGN(r,g));
450 for(i=myM-1;i>=l;i--){
453 e[i+1] = (r=pythag(f,g));
462 r = (d[i] - g) * s + 2.0 * c * b;
463 d[i+1] = g + ( p = s * r);
465 for(int k=0;k<n;k++){
467 z[k][i+1] = s * z[k][i] + c * f;
468 z[k][i] = c * z[k][i] - s * f;
471 if(r == 0.00 && i >= l) continue;
480 for(int i=0;i<n;i++){
482 for(int j=i;j<n;j++){
490 for(int j=0;j<n;j++){
500 catch(exception& e) {
501 m->errorOut(e, "LinearAlgebra", "qtli");
505 /*********************************************************************************************************************************/
506 //groups by dimension
507 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes, int dimensions){
510 vector< vector<double> > dists; dists.resize(axes.size());
511 for (int i = 0; i < dists.size(); i++) { dists[i].resize(axes.size(), 0.0); }
513 if (dimensions == 1) { //one dimension calc = abs(x-y)
515 for (int i = 0; i < dists.size(); i++) {
517 if (m->control_pressed) { return dists; }
519 for (int j = 0; j < i; j++) {
520 dists[i][j] = abs(axes[i][0] - axes[j][0]);
521 dists[j][i] = dists[i][j];
525 }else if (dimensions > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
527 for (int i = 0; i < dists.size(); i++) {
529 if (m->control_pressed) { return dists; }
531 for (int j = 0; j < i; j++) {
533 for (int k = 0; k < dimensions; k++) {
534 sum += ((axes[i][k] - axes[j][k]) * (axes[i][k] - axes[j][k]));
537 dists[i][j] = sqrt(sum);
538 dists[j][i] = dists[i][j];
546 catch(exception& e) {
547 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
551 /*********************************************************************************************************************************/
552 //returns groups by dimensions from dimensions by groups
553 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes){
556 vector< vector<double> > dists; dists.resize(axes[0].size());
557 for (int i = 0; i < dists.size(); i++) { dists[i].resize(axes[0].size(), 0.0); }
559 if (axes.size() == 1) { //one dimension calc = abs(x-y)
561 for (int i = 0; i < dists.size(); i++) {
563 if (m->control_pressed) { return dists; }
565 for (int j = 0; j < i; j++) {
566 dists[i][j] = abs(axes[0][i] - axes[0][j]);
567 dists[j][i] = dists[i][j];
571 }else if (axes.size() > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
573 for (int i = 0; i < dists[0].size(); i++) {
575 if (m->control_pressed) { return dists; }
577 for (int j = 0; j < i; j++) {
579 for (int k = 0; k < axes.size(); k++) {
580 sum += ((axes[k][i] - axes[k][j]) * (axes[k][i] - axes[k][j]));
583 dists[i][j] = sqrt(sum);
584 dists[j][i] = dists[i][j];
592 catch(exception& e) {
593 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
597 /*********************************************************************************************************************************/
598 //assumes both matrices are square and the same size
599 double LinearAlgebra::calcPearson(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
602 //find average for - X
604 float averageEuclid = 0.0;
605 for (int i = 0; i < euclidDists.size(); i++) {
606 for (int j = 0; j < i; j++) {
607 averageEuclid += euclidDists[i][j];
611 averageEuclid = averageEuclid / (float) count;
613 //find average for - Y
615 float averageUser = 0.0;
616 for (int i = 0; i < userDists.size(); i++) {
617 for (int j = 0; j < i; j++) {
618 averageUser += userDists[i][j];
622 averageUser = averageUser / (float) count;
624 double numerator = 0.0;
625 double denomTerm1 = 0.0;
626 double denomTerm2 = 0.0;
628 for (int i = 0; i < euclidDists.size(); i++) {
630 for (int k = 0; k < i; k++) { //just lt dists
632 float Yi = userDists[i][k];
633 float Xi = euclidDists[i][k];
635 numerator += ((Xi - averageEuclid) * (Yi - averageUser));
636 denomTerm1 += ((Xi - averageEuclid) * (Xi - averageEuclid));
637 denomTerm2 += ((Yi - averageUser) * (Yi - averageUser));
641 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
642 double r = numerator / denom;
644 //divide by zero error
645 if (isnan(r) || isinf(r)) { r = 0.0; }
650 catch(exception& e) {
651 m->errorOut(e, "LinearAlgebra", "calcPearson");
655 /*********************************************************************************************************************************/
656 //assumes both matrices are square and the same size
657 double LinearAlgebra::calcSpearman(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
662 map<float, int> tableX;
663 map<float, int>::iterator itTable;
664 vector<spearmanRank> scores;
666 for (int i = 0; i < euclidDists.size(); i++) {
667 for (int j = 0; j < i; j++) {
668 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
669 scores.push_back(member);
671 //count number of repeats
672 itTable = tableX.find(euclidDists[i][j]);
673 if (itTable == tableX.end()) {
674 tableX[euclidDists[i][j]] = 1;
676 tableX[euclidDists[i][j]]++;
682 sort(scores.begin(), scores.end(), compareSpearman);
686 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
687 double tx = (double) itTable->second;
688 Lx += ((pow(tx, 3.0) - tx) / 12.0);
692 map<string, float> rankEuclid;
693 vector<spearmanRank> ties;
695 for (int j = 0; j < scores.size(); j++) {
697 ties.push_back(scores[j]);
699 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
700 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
702 for (int k = 0; k < ties.size(); k++) {
703 float thisrank = rankTotal / (float) ties.size();
704 rankEuclid[ties[k].name] = thisrank;
709 }else { // you are the last one
711 for (int k = 0; k < ties.size(); k++) {
712 float thisrank = rankTotal / (float) ties.size();
713 rankEuclid[ties[k].name] = thisrank;
720 map<float, int> tableY;
723 for (int i = 0; i < userDists.size(); i++) {
724 for (int j = 0; j < i; j++) {
725 spearmanRank member(toString(scores.size()), userDists[i][j]);
726 scores.push_back(member);
728 //count number of repeats
729 itTable = tableY.find(userDists[i][j]);
730 if (itTable == tableY.end()) {
731 tableY[userDists[i][j]] = 1;
733 tableY[userDists[i][j]]++;
739 sort(scores.begin(), scores.end(), compareSpearman);
743 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
744 double ty = (double) itTable->second;
745 Ly += ((pow(ty, 3.0) - ty) / 12.0);
749 map<string, float> rankUser;
752 for (int j = 0; j < scores.size(); j++) {
754 ties.push_back(scores[j]);
756 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
757 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
759 for (int k = 0; k < ties.size(); k++) {
760 float thisrank = rankTotal / (float) ties.size();
761 rankUser[ties[k].name] = thisrank;
766 }else { // you are the last one
768 for (int k = 0; k < ties.size(); k++) {
769 float thisrank = rankTotal / (float) ties.size();
770 rankUser[ties[k].name] = thisrank;
778 for (int i = 0; i < userDists.size(); i++) {
779 for (int j = 0; j < i; j++) {
781 float xi = rankEuclid[toString(count)];
782 float yi = rankUser[toString(count)];
784 di += ((xi - yi) * (xi - yi));
790 double n = (double) count;
792 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
793 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
795 r = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
797 //divide by zero error
798 if (isnan(r) || isinf(r)) { r = 0.0; }
803 catch(exception& e) {
804 m->errorOut(e, "LinearAlgebra", "calcSpearman");
809 /*********************************************************************************************************************************/
810 //assumes both matrices are square and the same size
811 double LinearAlgebra::calcKendall(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
816 vector<spearmanRank> scores;
817 for (int i = 0; i < euclidDists.size(); i++) {
818 for (int j = 0; j < i; j++) {
819 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
820 scores.push_back(member);
825 sort(scores.begin(), scores.end(), compareSpearman);
828 map<string, float> rankEuclid;
829 vector<spearmanRank> ties;
831 for (int j = 0; j < scores.size(); j++) {
833 ties.push_back(scores[j]);
835 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
836 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
838 for (int k = 0; k < ties.size(); k++) {
839 float thisrank = rankTotal / (float) ties.size();
840 rankEuclid[ties[k].name] = thisrank;
845 }else { // you are the last one
847 for (int k = 0; k < ties.size(); k++) {
848 float thisrank = rankTotal / (float) ties.size();
849 rankEuclid[ties[k].name] = thisrank;
854 vector<spearmanRank> scoresUser;
855 for (int i = 0; i < userDists.size(); i++) {
856 for (int j = 0; j < i; j++) {
857 spearmanRank member(toString(scoresUser.size()), userDists[i][j]);
858 scoresUser.push_back(member);
863 sort(scoresUser.begin(), scoresUser.end(), compareSpearman);
866 map<string, float> rankUser;
869 for (int j = 0; j < scoresUser.size(); j++) {
871 ties.push_back(scoresUser[j]);
873 if (j != (scoresUser.size()-1)) { // you are not the last so you can look ahead
874 if (scoresUser[j].score != scoresUser[j+1].score) { // you are done with ties, rank them and continue
876 for (int k = 0; k < ties.size(); k++) {
877 float thisrank = rankTotal / (float) ties.size();
878 rankUser[ties[k].name] = thisrank;
883 }else { // you are the last one
885 for (int k = 0; k < ties.size(); k++) {
886 float thisrank = rankTotal / (float) ties.size();
887 rankUser[ties[k].name] = thisrank;
896 vector<spearmanRank> user;
897 for (int l = 0; l < scores.size(); l++) {
898 spearmanRank member(scores[l].name, rankUser[scores[l].name]);
899 user.push_back(member);
903 for (int l = 0; l < scores.size(); l++) {
905 int numWithHigherRank = 0;
906 int numWithLowerRank = 0;
907 float thisrank = user[l].score;
909 for (int u = l+1; u < scores.size(); u++) {
910 if (user[u].score > thisrank) { numWithHigherRank++; }
911 else if (user[u].score < thisrank) { numWithLowerRank++; }
915 numCoor += numWithHigherRank;
916 numDisCoor += numWithLowerRank;
919 r = (numCoor - numDisCoor) / (float) count;
921 //divide by zero error
922 if (isnan(r) || isinf(r)) { r = 0.0; }
927 catch(exception& e) {
928 m->errorOut(e, "LinearAlgebra", "calcKendall");
932 /*********************************************************************************************************************************/
933 double LinearAlgebra::calcKendall(vector<double>& x, vector<double>& y, double& sig){
935 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
938 vector<spearmanRank> xscores;
939 for (int i = 0; i < x.size(); i++) {
940 spearmanRank member(toString(i), x[i]);
941 xscores.push_back(member);
945 sort(xscores.begin(), xscores.end(), compareSpearman);
947 //convert scores to ranks of x
948 vector<spearmanRank*> ties;
950 for (int j = 0; j < xscores.size(); j++) {
952 ties.push_back(&(xscores[j]));
954 if (j != xscores.size()-1) { // you are not the last so you can look ahead
955 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
956 for (int k = 0; k < ties.size(); k++) {
957 float thisrank = rankTotal / (float) ties.size();
958 (*ties[k]).score = thisrank;
963 }else { // you are the last one
964 for (int k = 0; k < ties.size(); k++) {
965 float thisrank = rankTotal / (float) ties.size();
966 (*ties[k]).score = thisrank;
973 vector<spearmanRank> yscores;
974 for (int j = 0; j < y.size(); j++) {
975 spearmanRank member(toString(j), y[j]);
976 yscores.push_back(member);
980 sort(yscores.begin(), yscores.end(), compareSpearman);
983 map<string, float> rank;
984 vector<spearmanRank> yties;
986 for (int j = 0; j < yscores.size(); j++) {
988 yties.push_back(yscores[j]);
990 if (j != yscores.size()-1) { // you are not the last so you can look ahead
991 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
992 for (int k = 0; k < yties.size(); k++) {
993 float thisrank = rankTotal / (float) yties.size();
994 rank[yties[k].name] = thisrank;
999 }else { // you are the last one
1000 for (int k = 0; k < yties.size(); k++) {
1001 float thisrank = rankTotal / (float) yties.size();
1002 rank[yties[k].name] = thisrank;
1012 vector<spearmanRank> otus;
1013 for (int l = 0; l < xscores.size(); l++) {
1014 spearmanRank member(xscores[l].name, rank[xscores[l].name]);
1015 otus.push_back(member);
1019 for (int l = 0; l < xscores.size(); l++) {
1021 int numWithHigherRank = 0;
1022 int numWithLowerRank = 0;
1023 float thisrank = otus[l].score;
1025 for (int u = l+1; u < xscores.size(); u++) {
1026 if (otus[u].score > thisrank) { numWithHigherRank++; }
1027 else if (otus[u].score < thisrank) { numWithLowerRank++; }
1031 numCoor += numWithHigherRank;
1032 numDisCoor += numWithLowerRank;
1035 double p = (numCoor - numDisCoor) / (float) count;
1037 sig = calcKendallSig(x.size(), p);
1041 catch(exception& e) {
1042 m->errorOut(e, "LinearAlgebra", "calcKendall");
1046 double LinearAlgebra::ran0(int& idum)
1048 const int IA=16807,IM=2147483647,IQ=127773;
1049 const int IR=2836,MASK=123459876;
1050 const double AM=1.0/double(IM);
1056 idum=IA*(idum-k*IQ)-IR*k;
1057 if (idum < 0) idum += IM;
1063 double LinearAlgebra::ran1(int &idum)
1065 const int IA=16807,IM=2147483647,IQ=127773,IR=2836,NTAB=32;
1066 const int NDIV=(1+(IM-1)/NTAB);
1067 const double EPS=3.0e-16,AM=1.0/IM,RNMX=(1.0-EPS);
1069 static vector<int> iv(NTAB);
1073 if (idum <= 0 || !iy) {
1074 if (-idum < 1) idum=1;
1076 for (j=NTAB+7;j>=0;j--) {
1078 idum=IA*(idum-k*IQ)-IR*k;
1079 if (idum < 0) idum += IM;
1080 if (j < NTAB) iv[j] = idum;
1085 idum=IA*(idum-k*IQ)-IR*k;
1086 if (idum < 0) idum += IM;
1090 if ((temp=AM*iy) > RNMX) return RNMX;
1094 double LinearAlgebra::ran2(int &idum)
1096 const int IM1=2147483563,IM2=2147483399;
1097 const int IA1=40014,IA2=40692,IQ1=53668,IQ2=52774;
1098 const int IR1=12211,IR2=3791,NTAB=32,IMM1=IM1-1;
1099 const int NDIV=1+IMM1/NTAB;
1100 const double EPS=3.0e-16,RNMX=1.0-EPS,AM=1.0/double(IM1);
1101 static int idum2=123456789,iy=0;
1102 static vector<int> iv(NTAB);
1107 idum=(idum==0 ? 1 : -idum);
1109 for (j=NTAB+7;j>=0;j--) {
1111 idum=IA1*(idum-k*IQ1)-k*IR1;
1112 if (idum < 0) idum += IM1;
1113 if (j < NTAB) iv[j] = idum;
1118 idum=IA1*(idum-k*IQ1)-k*IR1;
1119 if (idum < 0) idum += IM1;
1121 idum2=IA2*(idum2-k*IQ2)-k*IR2;
1122 if (idum2 < 0) idum2 += IM2;
1126 if (iy < 1) iy += IMM1;
1127 if ((temp=AM*iy) > RNMX) return RNMX;
1131 double LinearAlgebra::ran3(int &idum)
1133 static int inext,inextp;
1135 const int MBIG=1000000000,MSEED=161803398,MZ=0;
1136 const double FAC=(1.0/MBIG);
1137 static vector<int> ma(56);
1140 if (idum < 0 || iff == 0) {
1142 mj=labs(MSEED-labs(idum));
1146 for (i=1;i<=54;i++) {
1150 if (mk < int(MZ)) mk += MBIG;
1154 for (i=1;i<=55;i++) {
1155 ma[i] -= ma[1+(i+30) % 55];
1156 if (ma[i] < int(MZ)) ma[i] += MBIG;
1162 if (++inext == 56) inext=1;
1163 if (++inextp == 56) inextp=1;
1164 mj=ma[inext]-ma[inextp];
1165 if (mj < int(MZ)) mj += MBIG;
1170 double LinearAlgebra::ran4(int &idum)
1172 #if defined(vax) || defined(_vax_) || defined(__vax__) || defined(VAX)
1173 static const unsigned long jflone = 0x00004080;
1174 static const unsigned long jflmsk = 0xffff007f;
1176 static const unsigned long jflone = 0x3f800000;
1177 static const unsigned long jflmsk = 0x007fffff;
1179 unsigned long irword,itemp,lword;
1180 static int idums = 0;
1188 psdes(lword,irword);
1189 itemp=jflone | (jflmsk & irword);
1191 return (*(float *)&itemp)-1.0;
1194 void LinearAlgebra::psdes(unsigned long &lword, unsigned long &irword)
1197 static const unsigned long c1[NITER]={
1198 0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L};
1199 static const unsigned long c2[NITER]={
1200 0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L};
1201 unsigned long i,ia,ib,iswap,itmph=0,itmpl=0;
1203 for (i=0;i<NITER;i++) {
1204 ia=(iswap=irword) ^ c1[i];
1205 itmpl = ia & 0xffff;
1207 ib=itmpl*itmpl+ ~(itmph*itmph);
1208 irword=lword ^ (((ia = (ib >> 16) |
1209 ((ib & 0xffff) << 16)) ^ c2[i])+itmpl*itmph);
1213 /*********************************************************************************************************************************/
1214 double LinearAlgebra::calcKendallSig(double n, double r){
1218 double svar=(4.0*n+10.0)/(9.0*n*(n-1.0));
1219 double z= r/sqrt(svar);
1220 sig=erfcc(fabs(z)/1.4142136);
1222 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1226 catch(exception& e) {
1227 m->errorOut(e, "LinearAlgebra", "calcKendallSig");
1232 /*********************************************************************************************************************************/
1233 double LinearAlgebra::calcSpearman(vector<double>& x, vector<double>& y, double& sig){
1235 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1238 double sf = 0.0; //f^3 - f where f is the number of ties in x;
1239 double sg = 0.0; //f^3 - f where f is the number of ties in y;
1240 map<float, int> tableX;
1241 map<float, int>::iterator itTable;
1242 vector<spearmanRank> xscores;
1244 for (int i = 0; i < x.size(); i++) {
1245 spearmanRank member(toString(i), x[i]);
1246 xscores.push_back(member);
1248 //count number of repeats
1249 itTable = tableX.find(x[i]);
1250 if (itTable == tableX.end()) {
1260 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
1261 double tx = (double) itTable->second;
1262 Lx += ((pow(tx, 3.0) - tx) / 12.0);
1267 sort(xscores.begin(), xscores.end(), compareSpearman);
1269 //convert scores to ranks of x
1271 map<string, float> rankx;
1272 vector<spearmanRank> xties;
1274 for (int j = 0; j < xscores.size(); j++) {
1276 xties.push_back(xscores[j]);
1278 if (j != xscores.size()-1) { // you are not the last so you can look ahead
1279 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
1280 for (int k = 0; k < xties.size(); k++) {
1281 float thisrank = rankTotal / (float) xties.size();
1282 rankx[xties[k].name] = thisrank;
1284 int t = xties.size();
1289 }else { // you are the last one
1290 for (int k = 0; k < xties.size(); k++) {
1291 float thisrank = rankTotal / (float) xties.size();
1292 rankx[xties[k].name] = thisrank;
1298 vector<spearmanRank> yscores;
1299 map<float, int> tableY;
1300 for (int j = 0; j < y.size(); j++) {
1301 spearmanRank member(toString(j), y[j]);
1302 yscores.push_back(member);
1304 itTable = tableY.find(member.score);
1305 if (itTable == tableY.end()) {
1306 tableY[member.score] = 1;
1308 tableY[member.score]++;
1315 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
1316 double ty = (double) itTable->second;
1317 Ly += ((pow(ty, 3.0) - ty) / 12.0);
1320 sort(yscores.begin(), yscores.end(), compareSpearman);
1323 map<string, float> rank;
1324 vector<spearmanRank> yties;
1326 for (int j = 0; j < yscores.size(); j++) {
1328 yties.push_back(yscores[j]);
1330 if (j != yscores.size()-1) { // you are not the last so you can look ahead
1331 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
1332 for (int k = 0; k < yties.size(); k++) {
1333 float thisrank = rankTotal / (float) yties.size();
1334 rank[yties[k].name] = thisrank;
1336 int t = yties.size();
1341 }else { // you are the last one
1342 for (int k = 0; k < yties.size(); k++) {
1343 float thisrank = rankTotal / (float) yties.size();
1344 rank[yties[k].name] = thisrank;
1350 for (int k = 0; k < x.size(); k++) {
1352 float xi = rankx[toString(k)];
1353 float yi = rank[toString(k)];
1355 di += ((xi - yi) * (xi - yi));
1360 double n = (double) x.size();
1361 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
1362 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
1364 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
1366 //Numerical Recipes 646
1367 sig = calcSpearmanSig(n, sf, sg, di);
1371 catch(exception& e) {
1372 m->errorOut(e, "LinearAlgebra", "calcSpearman");
1376 /*********************************************************************************************************************************/
1377 double LinearAlgebra::calcSpearmanSig(double n, double sf, double sg, double d){
1381 double probrs = 0.0;
1383 double en3n=en*en*en-en;
1384 double aved=en3n/6.0-(sf+sg)/12.0;
1385 double fac=(1.0-sf/en3n)*(1.0-sg/en3n);
1386 double vard=((en-1.0)*en*en*SQR(en+1.0)/36.0)*fac;
1387 double zd=(d-aved)/sqrt(vard);
1388 double probd=erfcc(fabs(zd)/1.4142136);
1389 double rs=(1.0-(6.0/en3n)*(d+(sf+sg)/12.0))/sqrt(fac);
1390 fac=(rs+1.0)*(1.0-rs);
1392 double t=rs*sqrt((en-2.0)/fac);
1394 probrs=betai(0.5*df,0.5,df/(df+t*t));
1399 //smaller of probd and probrs is sig
1401 if (probd < probrs) { sig = probd; }
1403 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1407 catch(exception& e) {
1408 m->errorOut(e, "LinearAlgebra", "calcSpearmanSig");
1412 /*********************************************************************************************************************************/
1413 double LinearAlgebra::calcPearson(vector<double>& x, vector<double>& y, double& sig){
1415 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1418 float averageX = 0.0;
1419 for (int i = 0; i < x.size(); i++) { averageX += x[i]; }
1420 averageX = averageX / (float) x.size();
1424 for (int j = 0; j < y.size(); j++) { sumY += y[j]; }
1425 float Ybar = sumY / (float) y.size();
1428 double numerator = 0.0;
1429 double denomTerm1 = 0.0;
1430 double denomTerm2 = 0.0;
1432 for (int j = 0; j < x.size(); j++) {
1436 numerator += ((Xi - averageX) * (Yi - Ybar));
1437 denomTerm1 += ((Xi - averageX) * (Xi - averageX));
1438 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
1441 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
1443 r = numerator / denom;
1445 //Numerical Recipes pg.644
1446 sig = calcPearsonSig(x.size(), r);
1450 catch(exception& e) {
1451 m->errorOut(e, "LinearAlgebra", "calcPearson");
1455 /*********************************************************************************************************************************/
1456 double LinearAlgebra::calcPearsonSig(double n, double r){
1460 const double TINY = 1.0e-20;
1461 double z = 0.5*log((1.0+r+TINY)/(1.0-r+TINY)); //Fisher's z transformation
1463 //code below was giving an error in betacf with sop files
1465 //double t = r*sqrt(df/((1.0-r+TINY)*(1.0+r+TINY)));
1466 //sig = betai(0.5+df, 0.5, df/(df+t*t));
1468 //Numerical Recipes says code below gives approximately the same result
1469 sig = erfcc(fabs(z*sqrt(n-1.0))/1.4142136);
1470 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1474 catch(exception& e) {
1475 m->errorOut(e, "LinearAlgebra", "calcPearsonSig");
1479 /*********************************************************************************************************************************/
1481 vector<vector<double> > LinearAlgebra::getObservedEuclideanDistance(vector<vector<double> >& relAbundData){
1484 int numSamples = relAbundData.size();
1485 int numOTUs = relAbundData[0].size();
1487 vector<vector<double> > dMatrix(numSamples);
1488 for(int i=0;i<numSamples;i++){
1489 dMatrix[i].resize(numSamples);
1492 for(int i=0;i<numSamples;i++){
1493 for(int j=0;j<numSamples;j++){
1495 if (m->control_pressed) { return dMatrix; }
1498 for(int k=0;k<numOTUs;k++){
1499 d += pow((relAbundData[i][k] - relAbundData[j][k]), 2.0000);
1501 dMatrix[i][j] = pow(d, 0.50000);
1502 dMatrix[j][i] = dMatrix[i][j];
1508 catch(exception& e) {
1509 m->errorOut(e, "LinearAlgebra", "getObservedEuclideanDistance");
1514 /*********************************************************************************************************************************/
1515 vector<double> LinearAlgebra::solveEquations(vector<vector<double> > A, vector<double> b){
1517 int length = (int)b.size();
1518 vector<double> x(length, 0);
1519 vector<int> index(length);
1520 for(int i=0;i<length;i++){ index[i] = i; }
1523 ludcmp(A, index, d); if (m->control_pressed) { return b; }
1524 lubksb(A, index, b);
1528 catch(exception& e) {
1529 m->errorOut(e, "LinearAlgebra", "solveEquations");
1533 /*********************************************************************************************************************************/
1534 vector<float> LinearAlgebra::solveEquations(vector<vector<float> > A, vector<float> b){
1536 int length = (int)b.size();
1537 vector<double> x(length, 0);
1538 vector<int> index(length);
1539 for(int i=0;i<length;i++){ index[i] = i; }
1542 ludcmp(A, index, d); if (m->control_pressed) { return b; }
1543 lubksb(A, index, b);
1547 catch(exception& e) {
1548 m->errorOut(e, "LinearAlgebra", "solveEquations");
1553 /*********************************************************************************************************************************/
1555 void LinearAlgebra::ludcmp(vector<vector<double> >& A, vector<int>& index, double& d){
1557 double tiny = 1e-20;
1559 int n = (int)A.size();
1560 vector<double> vv(n, 0.0);
1566 for(int i=0;i<n;i++){
1568 for(int j=0;j<n;j++){ if((temp=fabs(A[i][j])) > big ) big=temp; }
1569 if(big==0.0){ m->mothurOut("Singular matrix in routine ludcmp\n"); }
1573 for(int j=0;j<n;j++){
1574 if (m->control_pressed) { break; }
1575 for(int i=0;i<j;i++){
1576 double sum = A[i][j];
1577 for(int k=0;k<i;k++){ sum -= A[i][k] * A[k][j]; }
1582 for(int i=j;i<n;i++){
1583 double sum = A[i][j];
1584 for(int k=0;k<j;k++){ sum -= A[i][k] * A[k][j]; }
1587 if((dum = vv[i] * fabs(sum)) >= big){
1593 for(int k=0;k<n;k++){
1594 double dum = A[imax][k];
1595 A[imax][k] = A[j][k];
1603 if(A[j][j] == 0.0){ A[j][j] = tiny; }
1606 double dum = 1.0/A[j][j];
1607 for(int i=j+1;i<n;i++){ A[i][j] *= dum; }
1611 catch(exception& e) {
1612 m->errorOut(e, "LinearAlgebra", "ludcmp");
1617 /*********************************************************************************************************************************/
1619 void LinearAlgebra::lubksb(vector<vector<double> >& A, vector<int>& index, vector<double>& b){
1622 int n = (int)A.size();
1625 for(int i=0;i<n;i++){
1626 if (m->control_pressed) { break; }
1632 for(int j=ii-1;j<i;j++){
1633 total -= A[i][j] * b[j];
1636 else if(total != 0){ ii = i+1; }
1639 for(int i=n-1;i>=0;i--){
1641 for(int j=i+1;j<n;j++){ total -= A[i][j] * b[j]; }
1642 b[i] = total / A[i][i];
1645 catch(exception& e) {
1646 m->errorOut(e, "LinearAlgebra", "lubksb");
1650 /*********************************************************************************************************************************/
1652 void LinearAlgebra::ludcmp(vector<vector<float> >& A, vector<int>& index, float& d){
1654 double tiny = 1e-20;
1656 int n = (int)A.size();
1657 vector<float> vv(n, 0.0);
1663 for(int i=0;i<n;i++){
1665 for(int j=0;j<n;j++){ if((temp=fabs(A[i][j])) > big ) big=temp; }
1666 if(big==0.0){ m->mothurOut("Singular matrix in routine ludcmp\n"); }
1670 for(int j=0;j<n;j++){
1671 if (m->control_pressed) { break; }
1672 for(int i=0;i<j;i++){
1673 float sum = A[i][j];
1674 for(int k=0;k<i;k++){ sum -= A[i][k] * A[k][j]; }
1679 for(int i=j;i<n;i++){
1680 float sum = A[i][j];
1681 for(int k=0;k<j;k++){ sum -= A[i][k] * A[k][j]; }
1684 if((dum = vv[i] * fabs(sum)) >= big){
1690 for(int k=0;k<n;k++){
1691 float dum = A[imax][k];
1692 A[imax][k] = A[j][k];
1700 if(A[j][j] == 0.0){ A[j][j] = tiny; }
1703 float dum = 1.0/A[j][j];
1704 for(int i=j+1;i<n;i++){ A[i][j] *= dum; }
1708 catch(exception& e) {
1709 m->errorOut(e, "LinearAlgebra", "ludcmp");
1714 /*********************************************************************************************************************************/
1716 void LinearAlgebra::lubksb(vector<vector<float> >& A, vector<int>& index, vector<float>& b){
1719 int n = (int)A.size();
1722 for(int i=0;i<n;i++){
1723 if (m->control_pressed) { break; }
1729 for(int j=ii-1;j<i;j++){
1730 total -= A[i][j] * b[j];
1733 else if(total != 0){ ii = i+1; }
1736 for(int i=n-1;i>=0;i--){
1738 for(int j=i+1;j<n;j++){ total -= A[i][j] * b[j]; }
1739 b[i] = total / A[i][i];
1742 catch(exception& e) {
1743 m->errorOut(e, "LinearAlgebra", "lubksb");
1749 /*********************************************************************************************************************************/
1751 vector<vector<double> > LinearAlgebra::getInverse(vector<vector<double> > matrix){
1753 int n = (int)matrix.size();
1755 vector<vector<double> > inverse(n);
1756 for(int i=0;i<n;i++){ inverse[i].assign(n, 0.0000); }
1758 vector<double> column(n, 0.0000);
1759 vector<int> index(n, 0);
1762 ludcmp(matrix, index, dummy);
1764 for(int j=0;j<n;j++){
1765 if (m->control_pressed) { break; }
1767 column.assign(n, 0);
1771 lubksb(matrix, index, column);
1773 for(int i=0;i<n;i++){ inverse[i][j] = column[i]; }
1778 catch(exception& e) {
1779 m->errorOut(e, "LinearAlgebra", "getInverse");
1784 /*********************************************************************************************************************************/