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");
808 /*********************************************************************************************************************************/
809 double LinearAlgebra::calcKruskalWallis(vector<spearmanRank>& values, double& pValue){
812 set<string> treatments;
815 sort(values.begin(), values.end(), compareSpearman);
816 vector<spearmanRank*> ties;
819 for (int j = 0; j < values.size(); j++) {
820 treatments.insert(values[j].name);
822 ties.push_back(&(values[j]));
824 if (j != values.size()-1) { // you are not the last so you can look ahead
825 if (values[j].score != values[j+1].score) { // you are done with ties, rank them and continue
826 if (ties.size() > 1) { TIES.push_back(ties.size()); }
827 for (int k = 0; k < ties.size(); k++) {
828 double thisrank = rankTotal / (double) ties.size();
829 (*ties[k]).score = thisrank;
834 }else { // you are the last one
835 if (ties.size() > 1) { TIES.push_back(ties.size()); }
836 for (int k = 0; k < ties.size(); k++) {
837 double thisrank = rankTotal / (double) ties.size();
838 (*ties[k]).score = thisrank;
844 // H = 12/(N*(N+1)) * (sum Ti^2/n) - 3(N+1)
845 map<string, double> sums;
846 map<string, double> counts;
847 for (set<string>::iterator it = treatments.begin(); it != treatments.end(); it++) { sums[*it] = 0.0; counts[*it] = 0; }
849 for (int j = 0; j < values.size(); j++) {
850 sums[values[j].name] += values[j].score;
851 counts[values[j].name]+= 1.0;
854 double middleTerm = 0.0;
855 for (set<string>::iterator it = treatments.begin(); it != treatments.end(); it++) {
856 middleTerm += ((sums[*it]*sums[*it])/counts[*it]);
859 double firstTerm = 12 / (double) (values.size()*(values.size()+1));
860 double lastTerm = 3 * (values.size()+1);
862 H = firstTerm * middleTerm - lastTerm;
865 if (TIES.size() != 0) {
867 for (int j = 0; j < TIES.size(); j++) { sum += ((TIES[j]*TIES[j]*TIES[j])-TIES[j]); }
868 double result = 1.0 - (sum / (double) ((values.size()*values.size()*values.size())-values.size()));
872 //Numerical Recipes pg221
873 pValue = 1.0 - (gammp(((treatments.size()-1)/(double)2.0), H/2.0));
877 catch(exception& e) {
878 m->errorOut(e, "LinearAlgebra", "calcKruskalWallis");
882 /*********************************************************************************************************************************/
883 //assumes both matrices are square and the same size
884 double LinearAlgebra::calcKendall(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
889 vector<spearmanRank> scores;
890 for (int i = 0; i < euclidDists.size(); i++) {
891 for (int j = 0; j < i; j++) {
892 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
893 scores.push_back(member);
898 sort(scores.begin(), scores.end(), compareSpearman);
901 map<string, float> rankEuclid;
902 vector<spearmanRank> ties;
904 for (int j = 0; j < scores.size(); j++) {
906 ties.push_back(scores[j]);
908 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
909 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
911 for (int k = 0; k < ties.size(); k++) {
912 float thisrank = rankTotal / (float) ties.size();
913 rankEuclid[ties[k].name] = thisrank;
918 }else { // you are the last one
920 for (int k = 0; k < ties.size(); k++) {
921 float thisrank = rankTotal / (float) ties.size();
922 rankEuclid[ties[k].name] = thisrank;
927 vector<spearmanRank> scoresUser;
928 for (int i = 0; i < userDists.size(); i++) {
929 for (int j = 0; j < i; j++) {
930 spearmanRank member(toString(scoresUser.size()), userDists[i][j]);
931 scoresUser.push_back(member);
936 sort(scoresUser.begin(), scoresUser.end(), compareSpearman);
939 map<string, float> rankUser;
942 for (int j = 0; j < scoresUser.size(); j++) {
944 ties.push_back(scoresUser[j]);
946 if (j != (scoresUser.size()-1)) { // you are not the last so you can look ahead
947 if (scoresUser[j].score != scoresUser[j+1].score) { // you are done with ties, rank them and continue
949 for (int k = 0; k < ties.size(); k++) {
950 float thisrank = rankTotal / (float) ties.size();
951 rankUser[ties[k].name] = thisrank;
956 }else { // you are the last one
958 for (int k = 0; k < ties.size(); k++) {
959 float thisrank = rankTotal / (float) ties.size();
960 rankUser[ties[k].name] = thisrank;
969 vector<spearmanRank> user;
970 for (int l = 0; l < scores.size(); l++) {
971 spearmanRank member(scores[l].name, rankUser[scores[l].name]);
972 user.push_back(member);
976 for (int l = 0; l < scores.size(); l++) {
978 int numWithHigherRank = 0;
979 int numWithLowerRank = 0;
980 float thisrank = user[l].score;
982 for (int u = l+1; u < scores.size(); u++) {
983 if (user[u].score > thisrank) { numWithHigherRank++; }
984 else if (user[u].score < thisrank) { numWithLowerRank++; }
988 numCoor += numWithHigherRank;
989 numDisCoor += numWithLowerRank;
992 r = (numCoor - numDisCoor) / (float) count;
994 //divide by zero error
995 if (isnan(r) || isinf(r)) { r = 0.0; }
1000 catch(exception& e) {
1001 m->errorOut(e, "LinearAlgebra", "calcKendall");
1005 /*********************************************************************************************************************************/
1006 double LinearAlgebra::calcKendall(vector<double>& x, vector<double>& y, double& sig){
1008 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1011 vector<spearmanRank> xscores;
1012 for (int i = 0; i < x.size(); i++) {
1013 spearmanRank member(toString(i), x[i]);
1014 xscores.push_back(member);
1018 sort(xscores.begin(), xscores.end(), compareSpearman);
1020 //convert scores to ranks of x
1021 vector<spearmanRank*> ties;
1023 for (int j = 0; j < xscores.size(); j++) {
1025 ties.push_back(&(xscores[j]));
1027 if (j != xscores.size()-1) { // you are not the last so you can look ahead
1028 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
1029 for (int k = 0; k < ties.size(); k++) {
1030 float thisrank = rankTotal / (float) ties.size();
1031 (*ties[k]).score = thisrank;
1036 }else { // you are the last one
1037 for (int k = 0; k < ties.size(); k++) {
1038 float thisrank = rankTotal / (float) ties.size();
1039 (*ties[k]).score = thisrank;
1046 vector<spearmanRank> yscores;
1047 for (int j = 0; j < y.size(); j++) {
1048 spearmanRank member(toString(j), y[j]);
1049 yscores.push_back(member);
1053 sort(yscores.begin(), yscores.end(), compareSpearman);
1056 map<string, float> rank;
1057 vector<spearmanRank> yties;
1059 for (int j = 0; j < yscores.size(); j++) {
1061 yties.push_back(yscores[j]);
1063 if (j != yscores.size()-1) { // you are not the last so you can look ahead
1064 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
1065 for (int k = 0; k < yties.size(); k++) {
1066 float thisrank = rankTotal / (float) yties.size();
1067 rank[yties[k].name] = thisrank;
1072 }else { // you are the last one
1073 for (int k = 0; k < yties.size(); k++) {
1074 float thisrank = rankTotal / (float) yties.size();
1075 rank[yties[k].name] = thisrank;
1085 vector<spearmanRank> otus;
1086 for (int l = 0; l < xscores.size(); l++) {
1087 spearmanRank member(xscores[l].name, rank[xscores[l].name]);
1088 otus.push_back(member);
1092 for (int l = 0; l < xscores.size(); l++) {
1094 int numWithHigherRank = 0;
1095 int numWithLowerRank = 0;
1096 float thisrank = otus[l].score;
1098 for (int u = l+1; u < xscores.size(); u++) {
1099 if (otus[u].score > thisrank) { numWithHigherRank++; }
1100 else if (otus[u].score < thisrank) { numWithLowerRank++; }
1104 numCoor += numWithHigherRank;
1105 numDisCoor += numWithLowerRank;
1108 double p = (numCoor - numDisCoor) / (float) count;
1110 sig = calcKendallSig(x.size(), p);
1114 catch(exception& e) {
1115 m->errorOut(e, "LinearAlgebra", "calcKendall");
1119 double LinearAlgebra::ran0(int& idum)
1121 const int IA=16807,IM=2147483647,IQ=127773;
1122 const int IR=2836,MASK=123459876;
1123 const double AM=1.0/double(IM);
1129 idum=IA*(idum-k*IQ)-IR*k;
1130 if (idum < 0) idum += IM;
1136 double LinearAlgebra::ran1(int &idum)
1138 const int IA=16807,IM=2147483647,IQ=127773,IR=2836,NTAB=32;
1139 const int NDIV=(1+(IM-1)/NTAB);
1140 const double EPS=3.0e-16,AM=1.0/IM,RNMX=(1.0-EPS);
1142 static vector<int> iv(NTAB);
1146 if (idum <= 0 || !iy) {
1147 if (-idum < 1) idum=1;
1149 for (j=NTAB+7;j>=0;j--) {
1151 idum=IA*(idum-k*IQ)-IR*k;
1152 if (idum < 0) idum += IM;
1153 if (j < NTAB) iv[j] = idum;
1158 idum=IA*(idum-k*IQ)-IR*k;
1159 if (idum < 0) idum += IM;
1163 if ((temp=AM*iy) > RNMX) return RNMX;
1167 double LinearAlgebra::ran2(int &idum)
1169 const int IM1=2147483563,IM2=2147483399;
1170 const int IA1=40014,IA2=40692,IQ1=53668,IQ2=52774;
1171 const int IR1=12211,IR2=3791,NTAB=32,IMM1=IM1-1;
1172 const int NDIV=1+IMM1/NTAB;
1173 const double EPS=3.0e-16,RNMX=1.0-EPS,AM=1.0/double(IM1);
1174 static int idum2=123456789,iy=0;
1175 static vector<int> iv(NTAB);
1180 idum=(idum==0 ? 1 : -idum);
1182 for (j=NTAB+7;j>=0;j--) {
1184 idum=IA1*(idum-k*IQ1)-k*IR1;
1185 if (idum < 0) idum += IM1;
1186 if (j < NTAB) iv[j] = idum;
1191 idum=IA1*(idum-k*IQ1)-k*IR1;
1192 if (idum < 0) idum += IM1;
1194 idum2=IA2*(idum2-k*IQ2)-k*IR2;
1195 if (idum2 < 0) idum2 += IM2;
1199 if (iy < 1) iy += IMM1;
1200 if ((temp=AM*iy) > RNMX) return RNMX;
1204 double LinearAlgebra::ran3(int &idum)
1206 static int inext,inextp;
1208 const int MBIG=1000000000,MSEED=161803398,MZ=0;
1209 const double FAC=(1.0/MBIG);
1210 static vector<int> ma(56);
1213 if (idum < 0 || iff == 0) {
1215 mj=labs(MSEED-labs(idum));
1219 for (i=1;i<=54;i++) {
1223 if (mk < int(MZ)) mk += MBIG;
1227 for (i=1;i<=55;i++) {
1228 ma[i] -= ma[1+(i+30) % 55];
1229 if (ma[i] < int(MZ)) ma[i] += MBIG;
1235 if (++inext == 56) inext=1;
1236 if (++inextp == 56) inextp=1;
1237 mj=ma[inext]-ma[inextp];
1238 if (mj < int(MZ)) mj += MBIG;
1243 double LinearAlgebra::ran4(int &idum)
1245 #if defined(vax) || defined(_vax_) || defined(__vax__) || defined(VAX)
1246 static const unsigned long jflone = 0x00004080;
1247 static const unsigned long jflmsk = 0xffff007f;
1249 static const unsigned long jflone = 0x3f800000;
1250 static const unsigned long jflmsk = 0x007fffff;
1252 unsigned long irword,itemp,lword;
1253 static int idums = 0;
1261 psdes(lword,irword);
1262 itemp=jflone | (jflmsk & irword);
1264 return (*(float *)&itemp)-1.0;
1267 void LinearAlgebra::psdes(unsigned long &lword, unsigned long &irword)
1270 static const unsigned long c1[NITER]={
1271 0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L};
1272 static const unsigned long c2[NITER]={
1273 0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L};
1274 unsigned long i,ia,ib,iswap,itmph=0,itmpl=0;
1276 for (i=0;i<NITER;i++) {
1277 ia=(iswap=irword) ^ c1[i];
1278 itmpl = ia & 0xffff;
1280 ib=itmpl*itmpl+ ~(itmph*itmph);
1281 irword=lword ^ (((ia = (ib >> 16) |
1282 ((ib & 0xffff) << 16)) ^ c2[i])+itmpl*itmph);
1286 /*********************************************************************************************************************************/
1287 double LinearAlgebra::calcKendallSig(double n, double r){
1291 double svar=(4.0*n+10.0)/(9.0*n*(n-1.0));
1292 double z= r/sqrt(svar);
1293 sig=erfcc(fabs(z)/1.4142136);
1295 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1299 catch(exception& e) {
1300 m->errorOut(e, "LinearAlgebra", "calcKendallSig");
1304 /*********************************************************************************************************************************/
1305 double LinearAlgebra::calcWilcoxon(vector<double>& x, vector<double>& y, double& sig){
1307 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1312 vector<double> signPairs;
1313 vector<spearmanRank> absV;
1314 for (int i = 0; i < x.size(); i++) {
1315 if (m->control_pressed) { return W; }
1316 double temp = y[i]-x[i];
1317 double signV = sign(temp);
1318 if (signV != 0) { // exclude zeros
1319 spearmanRank member(toString(i), abs(temp));
1320 absV.push_back(member);
1321 signPairs.push_back(signV);
1327 sort(absV.begin(), absV.end(), compareSpearman);
1329 //convert scores to ranks of x
1330 vector<spearmanRank*> ties;
1332 for (int j = 0; j < absV.size(); j++) {
1333 if (m->control_pressed) { return W; }
1335 ties.push_back(&(absV[j]));
1337 if (j != absV.size()-1) { // you are not the last so you can look ahead
1338 if (absV[j].score != absV[j+1].score) { // you are done with ties, rank them and continue
1339 for (int k = 0; k < ties.size(); k++) {
1340 float thisrank = rankTotal / (float) ties.size();
1341 (*ties[k]).score = thisrank;
1346 }else { // you are the last one
1347 for (int k = 0; k < ties.size(); k++) {
1348 float thisrank = rankTotal / (float) ties.size();
1349 (*ties[k]).score = thisrank;
1354 //sum ranks times sign
1355 for (int i = 0; i < absV.size(); i++) {
1356 if (m->control_pressed) { return W; }
1357 W += (absV[i].score*signPairs[i]);
1362 cout << "still need to find sig!!" << endl;
1366 catch(exception& e) {
1367 m->errorOut(e, "LinearAlgebra", "calcWilcoxon");
1371 /*********************************************************************************************************************************/
1372 double LinearAlgebra::calcSpearman(vector<double>& x, vector<double>& y, double& sig){
1374 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1377 double sf = 0.0; //f^3 - f where f is the number of ties in x;
1378 double sg = 0.0; //f^3 - f where f is the number of ties in y;
1379 map<float, int> tableX;
1380 map<float, int>::iterator itTable;
1381 vector<spearmanRank> xscores;
1383 for (int i = 0; i < x.size(); i++) {
1384 spearmanRank member(toString(i), x[i]);
1385 xscores.push_back(member);
1387 //count number of repeats
1388 itTable = tableX.find(x[i]);
1389 if (itTable == tableX.end()) {
1399 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
1400 double tx = (double) itTable->second;
1401 Lx += ((pow(tx, 3.0) - tx) / 12.0);
1406 sort(xscores.begin(), xscores.end(), compareSpearman);
1408 //convert scores to ranks of x
1410 map<string, float> rankx;
1411 vector<spearmanRank> xties;
1413 for (int j = 0; j < xscores.size(); j++) {
1415 xties.push_back(xscores[j]);
1417 if (j != xscores.size()-1) { // you are not the last so you can look ahead
1418 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
1419 for (int k = 0; k < xties.size(); k++) {
1420 float thisrank = rankTotal / (float) xties.size();
1421 rankx[xties[k].name] = thisrank;
1423 int t = xties.size();
1428 }else { // you are the last one
1429 for (int k = 0; k < xties.size(); k++) {
1430 float thisrank = rankTotal / (float) xties.size();
1431 rankx[xties[k].name] = thisrank;
1437 vector<spearmanRank> yscores;
1438 map<float, int> tableY;
1439 for (int j = 0; j < y.size(); j++) {
1440 spearmanRank member(toString(j), y[j]);
1441 yscores.push_back(member);
1443 itTable = tableY.find(member.score);
1444 if (itTable == tableY.end()) {
1445 tableY[member.score] = 1;
1447 tableY[member.score]++;
1454 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
1455 double ty = (double) itTable->second;
1456 Ly += ((pow(ty, 3.0) - ty) / 12.0);
1459 sort(yscores.begin(), yscores.end(), compareSpearman);
1462 map<string, float> rank;
1463 vector<spearmanRank> yties;
1465 for (int j = 0; j < yscores.size(); j++) {
1467 yties.push_back(yscores[j]);
1469 if (j != yscores.size()-1) { // you are not the last so you can look ahead
1470 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
1471 for (int k = 0; k < yties.size(); k++) {
1472 float thisrank = rankTotal / (float) yties.size();
1473 rank[yties[k].name] = thisrank;
1475 int t = yties.size();
1480 }else { // you are the last one
1481 for (int k = 0; k < yties.size(); k++) {
1482 float thisrank = rankTotal / (float) yties.size();
1483 rank[yties[k].name] = thisrank;
1489 for (int k = 0; k < x.size(); k++) {
1491 float xi = rankx[toString(k)];
1492 float yi = rank[toString(k)];
1494 di += ((xi - yi) * (xi - yi));
1499 double n = (double) x.size();
1500 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
1501 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
1503 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
1505 //Numerical Recipes 646
1506 sig = calcSpearmanSig(n, sf, sg, di);
1510 catch(exception& e) {
1511 m->errorOut(e, "LinearAlgebra", "calcSpearman");
1515 /*********************************************************************************************************************************/
1516 double LinearAlgebra::calcSpearmanSig(double n, double sf, double sg, double d){
1520 double probrs = 0.0;
1522 double en3n=en*en*en-en;
1523 double aved=en3n/6.0-(sf+sg)/12.0;
1524 double fac=(1.0-sf/en3n)*(1.0-sg/en3n);
1525 double vard=((en-1.0)*en*en*SQR(en+1.0)/36.0)*fac;
1526 double zd=(d-aved)/sqrt(vard);
1527 double probd=erfcc(fabs(zd)/1.4142136);
1528 double rs=(1.0-(6.0/en3n)*(d+(sf+sg)/12.0))/sqrt(fac);
1529 fac=(rs+1.0)*(1.0-rs);
1531 double t=rs*sqrt((en-2.0)/fac);
1533 probrs=betai(0.5*df,0.5,df/(df+t*t));
1538 //smaller of probd and probrs is sig
1540 if (probd < probrs) { sig = probd; }
1542 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1546 catch(exception& e) {
1547 m->errorOut(e, "LinearAlgebra", "calcSpearmanSig");
1551 /*********************************************************************************************************************************/
1552 double LinearAlgebra::calcPearson(vector<double>& x, vector<double>& y, double& sig){
1554 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1557 float averageX = 0.0;
1558 for (int i = 0; i < x.size(); i++) { averageX += x[i]; }
1559 averageX = averageX / (float) x.size();
1563 for (int j = 0; j < y.size(); j++) { sumY += y[j]; }
1564 float Ybar = sumY / (float) y.size();
1567 double numerator = 0.0;
1568 double denomTerm1 = 0.0;
1569 double denomTerm2 = 0.0;
1571 for (int j = 0; j < x.size(); j++) {
1575 numerator += ((Xi - averageX) * (Yi - Ybar));
1576 denomTerm1 += ((Xi - averageX) * (Xi - averageX));
1577 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
1580 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
1582 r = numerator / denom;
1584 //Numerical Recipes pg.644
1585 sig = calcPearsonSig(x.size(), r);
1589 catch(exception& e) {
1590 m->errorOut(e, "LinearAlgebra", "calcPearson");
1594 /*********************************************************************************************************************************/
1595 double LinearAlgebra::calcPearsonSig(double n, double r){
1599 const double TINY = 1.0e-20;
1600 double z = 0.5*log((1.0+r+TINY)/(1.0-r+TINY)); //Fisher's z transformation
1602 //code below was giving an error in betacf with sop files
1604 //double t = r*sqrt(df/((1.0-r+TINY)*(1.0+r+TINY)));
1605 //sig = betai(0.5+df, 0.5, df/(df+t*t));
1607 //Numerical Recipes says code below gives approximately the same result
1608 sig = erfcc(fabs(z*sqrt(n-1.0))/1.4142136);
1609 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1613 catch(exception& e) {
1614 m->errorOut(e, "LinearAlgebra", "calcPearsonSig");
1618 /*********************************************************************************************************************************/
1620 vector<vector<double> > LinearAlgebra::getObservedEuclideanDistance(vector<vector<double> >& relAbundData){
1623 int numSamples = relAbundData.size();
1624 int numOTUs = relAbundData[0].size();
1626 vector<vector<double> > dMatrix(numSamples);
1627 for(int i=0;i<numSamples;i++){
1628 dMatrix[i].resize(numSamples);
1631 for(int i=0;i<numSamples;i++){
1632 for(int j=0;j<numSamples;j++){
1634 if (m->control_pressed) { return dMatrix; }
1637 for(int k=0;k<numOTUs;k++){
1638 d += pow((relAbundData[i][k] - relAbundData[j][k]), 2.0000);
1640 dMatrix[i][j] = pow(d, 0.50000);
1641 dMatrix[j][i] = dMatrix[i][j];
1647 catch(exception& e) {
1648 m->errorOut(e, "LinearAlgebra", "getObservedEuclideanDistance");
1653 /*********************************************************************************************************************************/
1654 vector<double> LinearAlgebra::solveEquations(vector<vector<double> > A, vector<double> b){
1656 int length = (int)b.size();
1657 vector<double> x(length, 0);
1658 vector<int> index(length);
1659 for(int i=0;i<length;i++){ index[i] = i; }
1662 ludcmp(A, index, d); if (m->control_pressed) { return b; }
1663 lubksb(A, index, b);
1667 catch(exception& e) {
1668 m->errorOut(e, "LinearAlgebra", "solveEquations");
1672 /*********************************************************************************************************************************/
1673 vector<float> LinearAlgebra::solveEquations(vector<vector<float> > A, vector<float> b){
1675 int length = (int)b.size();
1676 vector<double> x(length, 0);
1677 vector<int> index(length);
1678 for(int i=0;i<length;i++){ index[i] = i; }
1681 ludcmp(A, index, d); if (m->control_pressed) { return b; }
1682 lubksb(A, index, b);
1686 catch(exception& e) {
1687 m->errorOut(e, "LinearAlgebra", "solveEquations");
1692 /*********************************************************************************************************************************/
1694 void LinearAlgebra::ludcmp(vector<vector<double> >& A, vector<int>& index, double& d){
1696 double tiny = 1e-20;
1698 int n = (int)A.size();
1699 vector<double> vv(n, 0.0);
1705 for(int i=0;i<n;i++){
1707 for(int j=0;j<n;j++){ if((temp=fabs(A[i][j])) > big ) big=temp; }
1708 if(big==0.0){ m->mothurOut("Singular matrix in routine ludcmp\n"); }
1712 for(int j=0;j<n;j++){
1713 if (m->control_pressed) { break; }
1714 for(int i=0;i<j;i++){
1715 double sum = A[i][j];
1716 for(int k=0;k<i;k++){ sum -= A[i][k] * A[k][j]; }
1721 for(int i=j;i<n;i++){
1722 double sum = A[i][j];
1723 for(int k=0;k<j;k++){ sum -= A[i][k] * A[k][j]; }
1726 if((dum = vv[i] * fabs(sum)) >= big){
1732 for(int k=0;k<n;k++){
1733 double dum = A[imax][k];
1734 A[imax][k] = A[j][k];
1742 if(A[j][j] == 0.0){ A[j][j] = tiny; }
1745 double dum = 1.0/A[j][j];
1746 for(int i=j+1;i<n;i++){ A[i][j] *= dum; }
1750 catch(exception& e) {
1751 m->errorOut(e, "LinearAlgebra", "ludcmp");
1756 /*********************************************************************************************************************************/
1758 void LinearAlgebra::lubksb(vector<vector<double> >& A, vector<int>& index, vector<double>& b){
1761 int n = (int)A.size();
1764 for(int i=0;i<n;i++){
1765 if (m->control_pressed) { break; }
1771 for(int j=ii-1;j<i;j++){
1772 total -= A[i][j] * b[j];
1775 else if(total != 0){ ii = i+1; }
1778 for(int i=n-1;i>=0;i--){
1780 for(int j=i+1;j<n;j++){ total -= A[i][j] * b[j]; }
1781 b[i] = total / A[i][i];
1784 catch(exception& e) {
1785 m->errorOut(e, "LinearAlgebra", "lubksb");
1789 /*********************************************************************************************************************************/
1791 void LinearAlgebra::ludcmp(vector<vector<float> >& A, vector<int>& index, float& d){
1793 double tiny = 1e-20;
1795 int n = (int)A.size();
1796 vector<float> vv(n, 0.0);
1802 for(int i=0;i<n;i++){
1804 for(int j=0;j<n;j++){ if((temp=fabs(A[i][j])) > big ) big=temp; }
1805 if(big==0.0){ m->mothurOut("Singular matrix in routine ludcmp\n"); }
1809 for(int j=0;j<n;j++){
1810 if (m->control_pressed) { break; }
1811 for(int i=0;i<j;i++){
1812 float sum = A[i][j];
1813 for(int k=0;k<i;k++){ sum -= A[i][k] * A[k][j]; }
1818 for(int i=j;i<n;i++){
1819 float sum = A[i][j];
1820 for(int k=0;k<j;k++){ sum -= A[i][k] * A[k][j]; }
1823 if((dum = vv[i] * fabs(sum)) >= big){
1829 for(int k=0;k<n;k++){
1830 float dum = A[imax][k];
1831 A[imax][k] = A[j][k];
1839 if(A[j][j] == 0.0){ A[j][j] = tiny; }
1842 float dum = 1.0/A[j][j];
1843 for(int i=j+1;i<n;i++){ A[i][j] *= dum; }
1847 catch(exception& e) {
1848 m->errorOut(e, "LinearAlgebra", "ludcmp");
1853 /*********************************************************************************************************************************/
1855 void LinearAlgebra::lubksb(vector<vector<float> >& A, vector<int>& index, vector<float>& b){
1858 int n = (int)A.size();
1861 for(int i=0;i<n;i++){
1862 if (m->control_pressed) { break; }
1868 for(int j=ii-1;j<i;j++){
1869 total -= A[i][j] * b[j];
1872 else if(total != 0){ ii = i+1; }
1875 for(int i=n-1;i>=0;i--){
1877 for(int j=i+1;j<n;j++){ total -= A[i][j] * b[j]; }
1878 b[i] = total / A[i][i];
1881 catch(exception& e) {
1882 m->errorOut(e, "LinearAlgebra", "lubksb");
1888 /*********************************************************************************************************************************/
1890 vector<vector<double> > LinearAlgebra::getInverse(vector<vector<double> > matrix){
1892 int n = (int)matrix.size();
1894 vector<vector<double> > inverse(n);
1895 for(int i=0;i<n;i++){ inverse[i].assign(n, 0.0000); }
1897 vector<double> column(n, 0.0000);
1898 vector<int> index(n, 0);
1901 ludcmp(matrix, index, dummy);
1903 for(int j=0;j<n;j++){
1904 if (m->control_pressed) { break; }
1906 column.assign(n, 0);
1910 lubksb(matrix, index, column);
1912 for(int i=0;i<n;i++){ inverse[i][j] = column[i]; }
1917 catch(exception& e) {
1918 m->errorOut(e, "LinearAlgebra", "getInverse");
1923 /*********************************************************************************************************************************/