X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=linearalgebra.cpp;h=cd2ca0397bb21493761c69f729893103fc314613;hp=5de5a71dc35c1ab8888ff28ab8d7e92107f8cc72;hb=050a3ff02473a3d4c0980964e1a9ebe52e55d6b8;hpb=37eac2026d91179acda0494c4dcca22f176551b9 diff --git a/linearalgebra.cpp b/linearalgebra.cpp index 5de5a71..cd2ca03 100644 --- a/linearalgebra.cpp +++ b/linearalgebra.cpp @@ -8,7 +8,15 @@ */ #include "linearalgebra.h" +#include "wilcox.h" +// This class references functions used from "Numerical Recipes in C++" // + +/*********************************************************************************************************************************/ +inline double SQR(const double a) +{ + return a*a; +} /*********************************************************************************************************************************/ inline double SIGN(const double a, const double b) @@ -16,6 +24,235 @@ inline double SIGN(const double a, const double b) return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a); } /*********************************************************************************************************************************/ +//NUmerical recipes pg. 245 - Returns the complementary error function erfc(x) with fractional error everywhere less than 1.2 × 10−7. +double LinearAlgebra::erfcc(double x){ + try { + double t,z,ans; + z=fabs(x); + t=1.0/(1.0+0.5*z); + + ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+ + t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+ + t*(-0.82215223+t*0.17087277))))))))); + + //cout << "in erfcc " << t << '\t' << ans<< endl; + return (x >= 0.0 ? ans : 2.0 - ans); + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "betai"); + exit(1); + } +} +/*********************************************************************************************************************************/ +//Numerical Recipes pg. 232 +double LinearAlgebra::betai(const double a, const double b, const double x) { + try { + double bt; + double result = 0.0; + + if (x < 0.0 || x > 1.0) { m->mothurOut("[ERROR]: bad x in betai.\n"); m->control_pressed = true; return 0.0; } + + if (x == 0.0 || x == 1.0) { bt = 0.0; } + else { bt = exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x)); } + + if (x < (a+1.0) / (a+b+2.0)) { result = bt*betacf(a,b,x)/a; } + else { result = 1.0-bt*betacf(b,a,1.0-x)/b; } + + return result; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "betai"); + exit(1); + } +} +/*********************************************************************************************************************************/ +//Numerical Recipes pg. 219 +double LinearAlgebra::gammln(const double xx) { + try { + int j; + double x,y,tmp,ser; + static const double cof[6]={76.18009172947146,-86.50532032941677,24.01409824083091, + -1.231739572450155,0.120858003e-2,-0.536382e-5}; + + y=x=xx; + tmp=x+5.5; + tmp -= (x+0.5)*log(tmp); + ser=1.0; + for (j=0;j<6;j++) { + ser += cof[j]/++y; + } + return -tmp+log(2.5066282746310005*ser/x); + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "gammln"); + exit(1); + } +} +/*********************************************************************************************************************************/ +//Numerical Recipes pg. 223 +double LinearAlgebra::gammp(const double a, const double x) { + try { + double gamser,gammcf,gln; + + if (x < 0.0 || a <= 0.0) { m->mothurOut("[ERROR]: Invalid arguments in routine GAMMP\n"); m->control_pressed = true; return 0.0;} + if (x < (a+1.0)) { + gser(gamser,a,x,gln); + return gamser; + } else { + gcf(gammcf,a,x,gln); + return 1.0-gammcf; + } + + return 0; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "gammp"); + exit(1); + } +} +/*********************************************************************************************************************************/ +//Numerical Recipes pg. 223 +double LinearAlgebra::gammq(const double a, const double x) { + try { + double gamser,gammcf,gln; + + if (x < 0.0 || a <= 0.0) { m->mothurOut("[ERROR]: Invalid arguments in routine GAMMQ\n"); m->control_pressed = true; return 0.0; } + if (x < (a+1.0)) { + gser(gamser,a,x,gln); + return 1.0-gamser; + } else { + gcf(gammcf,a,x,gln); + return gammcf; + } + + return 0; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "gammp"); + exit(1); + } +} +/*********************************************************************************************************************************/ +//Numerical Recipes pg. 224 +double LinearAlgebra::gcf(double& gammcf, const double a, const double x, double& gln){ + try { + const int ITMAX=100; + const double EPS=numeric_limits::epsilon(); + const double FPMIN=numeric_limits::min()/EPS; + int i; + double an,b,c,d,del,h; + + gln=gammln(a); + b=x+1.0-a; + c=1.0/FPMIN; + d=1.0/b; + h=d; + for (i=1;i<=ITMAX;i++) { + an = -i*(i-a); + b += 2.0; + d=an*d+b; + if (fabs(d) < FPMIN) { d=FPMIN; } + c=b+an/c; + if (fabs(c) < FPMIN) { c=FPMIN; } + d=1.0/d; + del=d*c; + h *= del; + if (fabs(del-1.0) <= EPS) break; + } + if (i > ITMAX) { m->mothurOut("[ERROR]: a too large, ITMAX too small in gcf\n"); m->control_pressed = true; } + gammcf=exp(-x+a*log(x)-gln)*h; + + return 0.0; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "gcf"); + exit(1); + } + +} +/*********************************************************************************************************************************/ +//Numerical Recipes pg. 223 +double LinearAlgebra::gser(double& gamser, const double a, const double x, double& gln) { + try { + int n; + double sum,del,ap; + const double EPS = numeric_limits::epsilon(); + + gln=gammln(a); + if (x <= 0.0) { + if (x < 0.0) { m->mothurOut("[ERROR]: x less than 0 in routine GSER\n"); m->control_pressed = true; } + gamser=0.0; return 0.0; + } else { + ap=a; + del=sum=1.0/a; + for (n=0;n<100;n++) { + ++ap; + del *= x/ap; + sum += del; + if (fabs(del) < fabs(sum)*EPS) { + gamser=sum*exp(-x+a*log(x)-gln); + return 0.0; + } + } + + m->mothurOut("[ERROR]: a too large, ITMAX too small in routine GSER\n"); + return 0.0; + } + return 0; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "gser"); + exit(1); + } +} +/*********************************************************************************************************************************/ +//Numerical Recipes pg. 233 +double LinearAlgebra::betacf(const double a, const double b, const double x) { + try { + const int MAXIT = 100; + const double EPS = numeric_limits::epsilon(); + const double FPMIN = numeric_limits::min() / EPS; + int m1, m2; + double aa, c, d, del, h, qab, qam, qap; + + qab=a+b; + qap=a+1.0; + qam=a-1.0; + c=1.0; + d=1.0-qab*x/qap; + if (fabs(d) < FPMIN) d=FPMIN; + d=1.0/d; + h=d; + for (m1=1;m1<=MAXIT;m1++) { + m2=2*m1; + aa=m1*(b-m1)*x/((qam+m2)*(a+m2)); + d=1.0+aa*d; + if (fabs(d) < FPMIN) d=FPMIN; + c=1.0+aa/c; + if (fabs(c) < FPMIN) c=FPMIN; + d=1.0/d; + h *= d*c; + aa = -(a+m1)*(qab+m1)*x/((a+m2)*(qap+m2)); + d=1.0+aa*d; + if (fabs(d) < FPMIN) d=FPMIN; + c=1.0+aa/c; + if (fabs(c) < FPMIN) c=FPMIN; + d=1.0/d; + del=d*c; + h *= del; + if (fabs(del-1.0) < EPS) break; + } + + if (m1 > MAXIT) { m->mothurOut("[ERROR]: a or b too big or MAXIT too small in betacf."); m->mothurOutEndLine(); m->control_pressed = true; } + return h; + + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "betacf"); + exit(1); + } +} +/*********************************************************************************************************************************/ vector > LinearAlgebra::matrix_mult(vector > first, vector > second){ try { @@ -27,7 +264,7 @@ vector > LinearAlgebra::matrix_mult(vector > first product.resize(first_rows); for(int i=0;i > LinearAlgebra::matrix_mult(vector > first /*********************************************************************************************************************************/ +void LinearAlgebra::recenter(double offset, vector > D, vector >& G){ + try { + int rank = D.size(); + + vector > A(rank); + vector > C(rank); + for(int i=0;ierrorOut(e, "LinearAlgebra", "recenter"); + exit(1); + } + +} +/*********************************************************************************************************************************/ + // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479 int LinearAlgebra::tred2(vector >& a, vector& d, vector& e){ @@ -172,7 +442,7 @@ int LinearAlgebra::qtli(vector& d, vector& e, vector& d, vector& e, vector > LinearAlgebra::calculateEuclidianDistance(vector< vector >& axes, int dimensions){ try { //make square matrix @@ -252,38 +523,70 @@ vector< vector > LinearAlgebra::calculateEuclidianDistance(vector< vecto } } - }else if (dimensions == 2) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2) + }else if (dimensions > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)... for (int i = 0; i < dists.size(); i++) { if (m->control_pressed) { return dists; } for (int j = 0; j < i; j++) { - double firstDim = ((axes[i][0] - axes[j][0]) * (axes[i][0] - axes[j][0])); - double secondDim = ((axes[i][1] - axes[j][1]) * (axes[i][1] - axes[j][1])); + double sum = 0.0; + for (int k = 0; k < dimensions; k++) { + sum += ((axes[i][k] - axes[j][k]) * (axes[i][k] - axes[j][k])); + } - dists[i][j] = sqrt((firstDim + secondDim)); + dists[i][j] = sqrt(sum); dists[j][i] = dists[i][j]; } } - }else if (dimensions == 3) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2 + (x3 - y3)^2) + } + + return dists; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance"); + exit(1); + } +} +/*********************************************************************************************************************************/ +//returns groups by dimensions from dimensions by groups +vector< vector > LinearAlgebra::calculateEuclidianDistance(vector< vector >& axes){ + try { + //make square matrix + vector< vector > dists; dists.resize(axes[0].size()); + for (int i = 0; i < dists.size(); i++) { dists[i].resize(axes[0].size(), 0.0); } + + if (axes.size() == 1) { //one dimension calc = abs(x-y) for (int i = 0; i < dists.size(); i++) { if (m->control_pressed) { return dists; } for (int j = 0; j < i; j++) { - double firstDim = ((axes[i][0] - axes[j][0]) * (axes[i][0] - axes[j][0])); - double secondDim = ((axes[i][1] - axes[j][1]) * (axes[i][1] - axes[j][1])); - double thirdDim = ((axes[i][2] - axes[j][2]) * (axes[i][2] - axes[j][2])); + dists[i][j] = abs(axes[0][i] - axes[0][j]); + dists[j][i] = dists[i][j]; + } + } + + }else if (axes.size() > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)... + + for (int i = 0; i < dists[0].size(); i++) { + + if (m->control_pressed) { return dists; } + + for (int j = 0; j < i; j++) { + double sum = 0.0; + for (int k = 0; k < axes.size(); k++) { + sum += ((axes[k][i] - axes[k][j]) * (axes[k][i] - axes[k][j])); + } - dists[i][j] = sqrt((firstDim + secondDim + thirdDim)); + dists[i][j] = sqrt(sum); dists[j][i] = dists[i][j]; } } - }else { m->mothurOut("[ERROR]: too many dimensions, aborting."); m->mothurOutEndLine(); m->control_pressed = true; } + } return dists; } @@ -293,54 +596,1421 @@ vector< vector > LinearAlgebra::calculateEuclidianDistance(vector< vecto } } /*********************************************************************************************************************************/ +//assumes both matrices are square and the same size double LinearAlgebra::calcPearson(vector< vector >& euclidDists, vector< vector >& userDists){ try { //find average for - X - vector averageEuclid; averageEuclid.resize(euclidDists.size(), 0.0); + int count = 0; + float averageEuclid = 0.0; for (int i = 0; i < euclidDists.size(); i++) { - for (int j = 0; j < euclidDists[i].size(); j++) { - averageEuclid[i] += euclidDists[i][j]; + for (int j = 0; j < i; j++) { + averageEuclid += euclidDists[i][j]; + count++; } } - for (int i = 0; i < averageEuclid.size(); i++) { averageEuclid[i] = averageEuclid[i] / (float) euclidDists.size(); } - + averageEuclid = averageEuclid / (float) count; + //find average for - Y - vector averageUser; averageUser.resize(userDists.size(), 0.0); + count = 0; + float averageUser = 0.0; for (int i = 0; i < userDists.size(); i++) { - for (int j = 0; j < userDists[i].size(); j++) { - averageUser[i] += userDists[i][j]; + for (int j = 0; j < i; j++) { + averageUser += userDists[i][j]; + count++; } } - for (int i = 0; i < averageUser.size(); i++) { averageUser[i] = averageUser[i] / (float) userDists.size(); } - + averageUser = averageUser / (float) count; + double numerator = 0.0; double denomTerm1 = 0.0; double denomTerm2 = 0.0; for (int i = 0; i < euclidDists.size(); i++) { - for (int k = 0; k < i; k++) { + for (int k = 0; k < i; k++) { //just lt dists float Yi = userDists[i][k]; float Xi = euclidDists[i][k]; - numerator += ((Xi - averageEuclid[k]) * (Yi - averageUser[k])); - denomTerm1 += ((Xi - averageEuclid[k]) * (Xi - averageEuclid[k])); - denomTerm2 += ((Yi - averageUser[k]) * (Yi - averageUser[k])); + numerator += ((Xi - averageEuclid) * (Yi - averageUser)); + denomTerm1 += ((Xi - averageEuclid) * (Xi - averageEuclid)); + denomTerm2 += ((Yi - averageUser) * (Yi - averageUser)); } } double denom = (sqrt(denomTerm1) * sqrt(denomTerm2)); double r = numerator / denom; + //divide by zero error + if (isnan(r) || isinf(r)) { r = 0.0; } + return r; + } catch(exception& e) { - m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance"); + m->errorOut(e, "LinearAlgebra", "calcPearson"); exit(1); } } /*********************************************************************************************************************************/ +//assumes both matrices are square and the same size +double LinearAlgebra::calcSpearman(vector< vector >& euclidDists, vector< vector >& userDists){ + try { + double r; + + //format data + map tableX; + map::iterator itTable; + vector scores; + + for (int i = 0; i < euclidDists.size(); i++) { + for (int j = 0; j < i; j++) { + spearmanRank member(toString(scores.size()), euclidDists[i][j]); + scores.push_back(member); + + //count number of repeats + itTable = tableX.find(euclidDists[i][j]); + if (itTable == tableX.end()) { + tableX[euclidDists[i][j]] = 1; + }else { + tableX[euclidDists[i][j]]++; + } + } + } + + //sort scores + sort(scores.begin(), scores.end(), compareSpearman); + //calc LX + double Lx = 0.0; + for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) { + double tx = (double) itTable->second; + Lx += ((pow(tx, 3.0) - tx) / 12.0); + } + + //find ranks of xi + map rankEuclid; + vector ties; + int rankTotal = 0; + for (int j = 0; j < scores.size(); j++) { + rankTotal += (j+1); + ties.push_back(scores[j]); + + if (j != (scores.size()-1)) { // you are not the last so you can look ahead + if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue + + for (int k = 0; k < ties.size(); k++) { + float thisrank = rankTotal / (float) ties.size(); + rankEuclid[ties[k].name] = thisrank; + } + ties.clear(); + rankTotal = 0; + } + }else { // you are the last one + + for (int k = 0; k < ties.size(); k++) { + float thisrank = rankTotal / (float) ties.size(); + rankEuclid[ties[k].name] = thisrank; + } + } + } + + + //format data + map tableY; + scores.clear(); + + for (int i = 0; i < userDists.size(); i++) { + for (int j = 0; j < i; j++) { + spearmanRank member(toString(scores.size()), userDists[i][j]); + scores.push_back(member); + + //count number of repeats + itTable = tableY.find(userDists[i][j]); + if (itTable == tableY.end()) { + tableY[userDists[i][j]] = 1; + }else { + tableY[userDists[i][j]]++; + } + } + } + + //sort scores + sort(scores.begin(), scores.end(), compareSpearman); + + //calc LX + double Ly = 0.0; + for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) { + double ty = (double) itTable->second; + Ly += ((pow(ty, 3.0) - ty) / 12.0); + } + + //find ranks of yi + map rankUser; + ties.clear(); + rankTotal = 0; + for (int j = 0; j < scores.size(); j++) { + rankTotal += (j+1); + ties.push_back(scores[j]); + + if (j != (scores.size()-1)) { // you are not the last so you can look ahead + if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue + + for (int k = 0; k < ties.size(); k++) { + float thisrank = rankTotal / (float) ties.size(); + rankUser[ties[k].name] = thisrank; + } + ties.clear(); + rankTotal = 0; + } + }else { // you are the last one + + for (int k = 0; k < ties.size(); k++) { + float thisrank = rankTotal / (float) ties.size(); + rankUser[ties[k].name] = thisrank; + } + } + } + + + double di = 0.0; + int count = 0; + for (int i = 0; i < userDists.size(); i++) { + for (int j = 0; j < i; j++) { + + float xi = rankEuclid[toString(count)]; + float yi = rankUser[toString(count)]; + + di += ((xi - yi) * (xi - yi)); + + count++; + } + } + + double n = (double) count; + + double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx; + double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly; + + r = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2))); + + //divide by zero error + if (isnan(r) || isinf(r)) { r = 0.0; } + + return r; + + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "calcSpearman"); + exit(1); + } +} +/*********************************************************************************************************************************/ +//assumes both matrices are square and the same size +double LinearAlgebra::calcKendall(vector< vector >& euclidDists, vector< vector >& userDists){ + try { + double r; + + //format data + vector scores; + for (int i = 0; i < euclidDists.size(); i++) { + for (int j = 0; j < i; j++) { + spearmanRank member(toString(scores.size()), euclidDists[i][j]); + scores.push_back(member); + } + } + + //sort scores + sort(scores.begin(), scores.end(), compareSpearman); + + //find ranks of xi + map rankEuclid; + vector ties; + int rankTotal = 0; + for (int j = 0; j < scores.size(); j++) { + rankTotal += (j+1); + ties.push_back(scores[j]); + + if (j != (scores.size()-1)) { // you are not the last so you can look ahead + if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue + + for (int k = 0; k < ties.size(); k++) { + float thisrank = rankTotal / (float) ties.size(); + rankEuclid[ties[k].name] = thisrank; + } + ties.clear(); + rankTotal = 0; + } + }else { // you are the last one + + for (int k = 0; k < ties.size(); k++) { + float thisrank = rankTotal / (float) ties.size(); + rankEuclid[ties[k].name] = thisrank; + } + } + } + + vector scoresUser; + for (int i = 0; i < userDists.size(); i++) { + for (int j = 0; j < i; j++) { + spearmanRank member(toString(scoresUser.size()), userDists[i][j]); + scoresUser.push_back(member); + } + } + + //sort scores + sort(scoresUser.begin(), scoresUser.end(), compareSpearman); + + //find ranks of yi + map rankUser; + ties.clear(); + rankTotal = 0; + for (int j = 0; j < scoresUser.size(); j++) { + rankTotal += (j+1); + ties.push_back(scoresUser[j]); + + if (j != (scoresUser.size()-1)) { // you are not the last so you can look ahead + if (scoresUser[j].score != scoresUser[j+1].score) { // you are done with ties, rank them and continue + + for (int k = 0; k < ties.size(); k++) { + float thisrank = rankTotal / (float) ties.size(); + rankUser[ties[k].name] = thisrank; + } + ties.clear(); + rankTotal = 0; + } + }else { // you are the last one + + for (int k = 0; k < ties.size(); k++) { + float thisrank = rankTotal / (float) ties.size(); + rankUser[ties[k].name] = thisrank; + } + } + } + + int numCoor = 0; + int numDisCoor = 0; + + //order user ranks + vector user; + for (int l = 0; l < scores.size(); l++) { + spearmanRank member(scores[l].name, rankUser[scores[l].name]); + user.push_back(member); + } + + int count = 0; + for (int l = 0; l < scores.size(); l++) { + + int numWithHigherRank = 0; + int numWithLowerRank = 0; + float thisrank = user[l].score; + + for (int u = l+1; u < scores.size(); u++) { + if (user[u].score > thisrank) { numWithHigherRank++; } + else if (user[u].score < thisrank) { numWithLowerRank++; } + count++; + } + + numCoor += numWithHigherRank; + numDisCoor += numWithLowerRank; + } + + r = (numCoor - numDisCoor) / (float) count; + + //divide by zero error + if (isnan(r) || isinf(r)) { r = 0.0; } + + return r; + + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "calcKendall"); + exit(1); + } +} +/*********************************************************************************************************************************/ +double LinearAlgebra::calcKendall(vector& x, vector& y, double& sig){ + try { + if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; } + + //format data + vector xscores; + for (int i = 0; i < x.size(); i++) { + spearmanRank member(toString(i), x[i]); + xscores.push_back(member); + } + + //sort xscores + sort(xscores.begin(), xscores.end(), compareSpearman); + + //convert scores to ranks of x + vector ties; + int rankTotal = 0; + for (int j = 0; j < xscores.size(); j++) { + rankTotal += (j+1); + ties.push_back(&(xscores[j])); + + if (j != xscores.size()-1) { // you are not the last so you can look ahead + if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue + for (int k = 0; k < ties.size(); k++) { + float thisrank = rankTotal / (float) ties.size(); + (*ties[k]).score = thisrank; + } + ties.clear(); + rankTotal = 0; + } + }else { // you are the last one + for (int k = 0; k < ties.size(); k++) { + float thisrank = rankTotal / (float) ties.size(); + (*ties[k]).score = thisrank; + } + } + } + + + //format data + vector yscores; + for (int j = 0; j < y.size(); j++) { + spearmanRank member(toString(j), y[j]); + yscores.push_back(member); + } + + //sort yscores + sort(yscores.begin(), yscores.end(), compareSpearman); + + //convert to ranks + map rank; + vector yties; + rankTotal = 0; + for (int j = 0; j < yscores.size(); j++) { + rankTotal += (j+1); + yties.push_back(yscores[j]); + + if (j != yscores.size()-1) { // you are not the last so you can look ahead + if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue + for (int k = 0; k < yties.size(); k++) { + float thisrank = rankTotal / (float) yties.size(); + rank[yties[k].name] = thisrank; + } + yties.clear(); + rankTotal = 0; + } + }else { // you are the last one + for (int k = 0; k < yties.size(); k++) { + float thisrank = rankTotal / (float) yties.size(); + rank[yties[k].name] = thisrank; + } + } + } + + + int numCoor = 0; + int numDisCoor = 0; + + //associate x and y + vector otus; + for (int l = 0; l < xscores.size(); l++) { + spearmanRank member(xscores[l].name, rank[xscores[l].name]); + otus.push_back(member); + } + + int count = 0; + for (int l = 0; l < xscores.size(); l++) { + + int numWithHigherRank = 0; + int numWithLowerRank = 0; + float thisrank = otus[l].score; + + for (int u = l+1; u < xscores.size(); u++) { + if (otus[u].score > thisrank) { numWithHigherRank++; } + else if (otus[u].score < thisrank) { numWithLowerRank++; } + count++; + } + + numCoor += numWithHigherRank; + numDisCoor += numWithLowerRank; + } + + double p = (numCoor - numDisCoor) / (float) count; + + sig = calcKendallSig(x.size(), p); + + return p; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "calcKendall"); + exit(1); + } +} +double LinearAlgebra::ran0(int& idum) +{ + const int IA=16807,IM=2147483647,IQ=127773; + const int IR=2836,MASK=123459876; + const double AM=1.0/double(IM); + int k; + double ans; + + idum ^= MASK; + k=idum/IQ; + idum=IA*(idum-k*IQ)-IR*k; + if (idum < 0) idum += IM; + ans=AM*idum; + idum ^= MASK; + return ans; +} + +double LinearAlgebra::ran1(int &idum) +{ + const int IA=16807,IM=2147483647,IQ=127773,IR=2836,NTAB=32; + const int NDIV=(1+(IM-1)/NTAB); + const double EPS=3.0e-16,AM=1.0/IM,RNMX=(1.0-EPS); + static int iy=0; + static vector iv(NTAB); + int j,k; + double temp; + + if (idum <= 0 || !iy) { + if (-idum < 1) idum=1; + else idum = -idum; + for (j=NTAB+7;j>=0;j--) { + k=idum/IQ; + idum=IA*(idum-k*IQ)-IR*k; + if (idum < 0) idum += IM; + if (j < NTAB) iv[j] = idum; + } + iy=iv[0]; + } + k=idum/IQ; + idum=IA*(idum-k*IQ)-IR*k; + if (idum < 0) idum += IM; + j=iy/NDIV; + iy=iv[j]; + iv[j] = idum; + if ((temp=AM*iy) > RNMX) return RNMX; + else return temp; +} + +double LinearAlgebra::ran2(int &idum) +{ + const int IM1=2147483563,IM2=2147483399; + const int IA1=40014,IA2=40692,IQ1=53668,IQ2=52774; + const int IR1=12211,IR2=3791,NTAB=32,IMM1=IM1-1; + const int NDIV=1+IMM1/NTAB; + const double EPS=3.0e-16,RNMX=1.0-EPS,AM=1.0/double(IM1); + static int idum2=123456789,iy=0; + static vector iv(NTAB); + int j,k; + double temp; + + if (idum <= 0) { + idum=(idum==0 ? 1 : -idum); + idum2=idum; + for (j=NTAB+7;j>=0;j--) { + k=idum/IQ1; + idum=IA1*(idum-k*IQ1)-k*IR1; + if (idum < 0) idum += IM1; + if (j < NTAB) iv[j] = idum; + } + iy=iv[0]; + } + k=idum/IQ1; + idum=IA1*(idum-k*IQ1)-k*IR1; + if (idum < 0) idum += IM1; + k=idum2/IQ2; + idum2=IA2*(idum2-k*IQ2)-k*IR2; + if (idum2 < 0) idum2 += IM2; + j=iy/NDIV; + iy=iv[j]-idum2; + iv[j] = idum; + if (iy < 1) iy += IMM1; + if ((temp=AM*iy) > RNMX) return RNMX; + else return temp; +} + +double LinearAlgebra::ran3(int &idum) +{ + static int inext,inextp; + static int iff=0; + const int MBIG=1000000000,MSEED=161803398,MZ=0; + const double FAC=(1.0/MBIG); + static vector ma(56); + int i,ii,k,mj,mk; + + if (idum < 0 || iff == 0) { + iff=1; + mj=labs(MSEED-labs(idum)); + mj %= MBIG; + ma[55]=mj; + mk=1; + for (i=1;i<=54;i++) { + ii=(21*i) % 55; + ma[ii]=mk; + mk=mj-mk; + if (mk < int(MZ)) mk += MBIG; + mj=ma[ii]; + } + for (k=0;k<4;k++) + for (i=1;i<=55;i++) { + ma[i] -= ma[1+(i+30) % 55]; + if (ma[i] < int(MZ)) ma[i] += MBIG; + } + inext=0; + inextp=31; + idum=1; + } + if (++inext == 56) inext=1; + if (++inextp == 56) inextp=1; + mj=ma[inext]-ma[inextp]; + if (mj < int(MZ)) mj += MBIG; + ma[inext]=mj; + return mj*FAC; +} + +double LinearAlgebra::ran4(int &idum) +{ +#if defined(vax) || defined(_vax_) || defined(__vax__) || defined(VAX) + static const unsigned long jflone = 0x00004080; + static const unsigned long jflmsk = 0xffff007f; +#else + static const unsigned long jflone = 0x3f800000; + static const unsigned long jflmsk = 0x007fffff; +#endif + unsigned long irword,itemp,lword; + static int idums = 0; + + if (idum < 0) { + idums = -idum; + idum=1; + } + irword=idum; + lword=idums; + psdes(lword,irword); + itemp=jflone | (jflmsk & irword); + ++idum; + return (*(float *)&itemp)-1.0; +} + +void LinearAlgebra::psdes(unsigned long &lword, unsigned long &irword) +{ + const int NITER=4; + static const unsigned long c1[NITER]={ + 0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L}; + static const unsigned long c2[NITER]={ + 0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L}; + unsigned long i,ia,ib,iswap,itmph=0,itmpl=0; + + for (i=0;i> 16; + ib=itmpl*itmpl+ ~(itmph*itmph); + irword=lword ^ (((ia = (ib >> 16) | + ((ib & 0xffff) << 16)) ^ c2[i])+itmpl*itmph); + lword=iswap; + } +} +/*********************************************************************************************************************************/ +double LinearAlgebra::calcKendallSig(double n, double r){ + try { + + double sig = 0.0; + double svar=(4.0*n+10.0)/(9.0*n*(n-1.0)); + double z= r/sqrt(svar); + sig=erfcc(fabs(z)/1.4142136); + + if (isnan(sig) || isinf(sig)) { sig = 0.0; } + + return sig; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "calcKendallSig"); + exit(1); + } +} +/*********************************************************************************************************************************/ +double LinearAlgebra::calcKruskalWallis(vector& values, double& pValue){ + try { + double H; + set treatments; + + //rank values + sort(values.begin(), values.end(), compareSpearman); + vector ties; + int rankTotal = 0; + vector TIES; + for (int j = 0; j < values.size(); j++) { + treatments.insert(values[j].name); + rankTotal += (j+1); + ties.push_back(&(values[j])); + + if (j != values.size()-1) { // you are not the last so you can look ahead + if (values[j].score != values[j+1].score) { // you are done with ties, rank them and continue + if (ties.size() > 1) { TIES.push_back(ties.size()); } + for (int k = 0; k < ties.size(); k++) { + double thisrank = rankTotal / (double) ties.size(); + (*ties[k]).score = thisrank; + } + ties.clear(); + rankTotal = 0; + } + }else { // you are the last one + if (ties.size() > 1) { TIES.push_back(ties.size()); } + for (int k = 0; k < ties.size(); k++) { + double thisrank = rankTotal / (double) ties.size(); + (*ties[k]).score = thisrank; + } + } + } + + + // H = 12/(N*(N+1)) * (sum Ti^2/n) - 3(N+1) + map sums; + map counts; + for (set::iterator it = treatments.begin(); it != treatments.end(); it++) { sums[*it] = 0.0; counts[*it] = 0; } + + for (int j = 0; j < values.size(); j++) { + sums[values[j].name] += values[j].score; + counts[values[j].name]+= 1.0; + } + + double middleTerm = 0.0; + for (set::iterator it = treatments.begin(); it != treatments.end(); it++) { + middleTerm += ((sums[*it]*sums[*it])/counts[*it]); + } + + double firstTerm = 12 / (double) (values.size()*(values.size()+1)); + double lastTerm = 3 * (values.size()+1); + + H = firstTerm * middleTerm - lastTerm; + + //adjust for ties + if (TIES.size() != 0) { + double sum = 0.0; + for (int j = 0; j < TIES.size(); j++) { sum += ((TIES[j]*TIES[j]*TIES[j])-TIES[j]); } + double result = 1.0 - (sum / (double) ((values.size()*values.size()*values.size())-values.size())); + H /= result; + } + + //Numerical Recipes pg221 + pValue = 1.0 - (gammp(((treatments.size()-1)/(double)2.0), H/2.0)); + + return H; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "calcKruskalWallis"); + exit(1); + } +} +/*********************************************************************************************************************************/ +//thanks http://www.johndcook.com/cpp_phi.html +double LinearAlgebra::pnorm(double x){ + try { + // constants + double a1 = 0.254829592; + double a2 = -0.284496736; + double a3 = 1.421413741; + double a4 = -1.453152027; + double a5 = 1.061405429; + double p = 0.3275911; + + // Save the sign of x + int sign = 1; + if (x < 0) + sign = -1; + x = fabs(x)/sqrt(2.0); + + // A&S formula 7.1.26 + double t = 1.0/(1.0 + p*x); + double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x); + + return 0.5*(1.0 + sign*y); + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "pnorm"); + exit(1); + } +} + +/*********************************************************************************************************************************/ +double LinearAlgebra::calcWilcoxon(vector& x, vector& y, double& sig){ + try { + double W = 0.0; + sig = 0.0; + + vector ranks; + for (int i = 0; i < x.size(); i++) { + if (m->control_pressed) { return W; } + spearmanRank member("x", x[i]); + ranks.push_back(member); + } + + for (int i = 0; i < y.size(); i++) { + if (m->control_pressed) { return W; } + spearmanRank member("y", y[i]); + ranks.push_back(member); + } + + //sort values + sort(ranks.begin(), ranks.end(), compareSpearman); + + //convert scores to ranks of x + vector ties; + int rankTotal = 0; + vector TIES; + for (int j = 0; j < ranks.size(); j++) { + if (m->control_pressed) { return W; } + rankTotal += (j+1); + ties.push_back(&(ranks[j])); + + if (j != ranks.size()-1) { // you are not the last so you can look ahead + if (ranks[j].score != ranks[j+1].score) { // you are done with ties, rank them and continue + if (ties.size() > 1) { TIES.push_back(ties.size()); } + for (int k = 0; k < ties.size(); k++) { + float thisrank = rankTotal / (float) ties.size(); + (*ties[k]).score = thisrank; + } + ties.clear(); + rankTotal = 0; + } + }else { // you are the last one + if (ties.size() > 1) { TIES.push_back(ties.size()); } + for (int k = 0; k < ties.size(); k++) { + float thisrank = rankTotal / (float) ties.size(); + (*ties[k]).score = thisrank; + } + } + } + + //from R wilcox.test function + //STATISTIC <- sum(r[seq_along(x)]) - n.x * (n.x + 1)/2 + double sumRanks = 0.0; + for (int i = 0; i < ranks.size(); i++) { + if (m->control_pressed) { return W; } + if (ranks[i].name == "x") { sumRanks += ranks[i].score; } + } + + W = sumRanks - x.size() * ((double)(x.size() + 1)) / 2.0; + + //exact <- (n.x < 50) && (n.y < 50) + bool findExact = false; + if ((x.size() < 50) && (y.size() < 50)) { findExact = true; } + + + if (findExact && (TIES.size() == 0)) { //find exact and no ties + //PVAL <- switch(alternative, two.sided = { + //p <- if (STATISTIC > (n.x * n.y/2)) + PWilcox wilcox; + double pval = 0.0; + if (W > ((double)x.size()*y.size()/2.0)) { + //pwilcox(STATISTIC-1, n.x, n.y, lower.tail = FALSE) + pval = wilcox.pwilcox(W-1, x.size(), y.size(), false); + }else { + //pwilcox(STATISTIC,n.x, n.y) + pval = wilcox.pwilcox(W, x.size(), y.size(), true); + } + sig = 2.0 * pval; + if (1.0 < sig) { sig = 1.0; } + }else { + //z <- STATISTIC - n.x * n.y/2 + double z = W - (double)(x.size() * y.size()/2.0); + //NTIES <- table(r) + double sum = 0.0; + for (int j = 0; j < TIES.size(); j++) { sum += ((TIES[j]*TIES[j]*TIES[j])-TIES[j]); } + + //SIGMA <- sqrt((n.x * n.y/12) * ((n.x + n.y + 1) - + //sum(NTIES^3 - NTIES)/((n.x + n.y) * (n.x + n.y - + //1)))) + double sigma = 0.0; + double firstTerm = (double)(x.size() * y.size()/12.0); + double secondTerm = (double)(x.size() + y.size() + 1) - sum / (double)((x.size() + y.size()) * (x.size() + y.size() - 1)); + sigma = sqrt(firstTerm * secondTerm); + + //CORRECTION <- switch(alternative, two.sided = sign(z) * 0.5, greater = 0.5, less = -0.5) + double CORRECTION = 0.0; + if (z < 0) { CORRECTION = -1.0; } + else if (z > 0) { CORRECTION = 1.0; } + CORRECTION *= 0.5; + + z = (z - CORRECTION)/sigma; + + //PVAL <- switch(alternative, two.sided = 2 * min(pnorm(z), pnorm(z, lower.tail = FALSE))) + sig = pnorm(z); + if ((1.0-sig) < sig) { sig = 1.0 - sig; } + sig *= 2; + } + + return W; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "calcWilcoxon"); + exit(1); + } +} + +/*********************************************************************************************************************************/ +double LinearAlgebra::choose(double n, double k){ + try { + n = floor(n + 0.5); + k = floor(k + 0.5); + + double lchoose = gammln(n + 1.0) - gammln(k + 1.0) - gammln(n - k + 1.0); + + return (floor(exp(lchoose) + 0.5)); + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "choose"); + exit(1); + } +} +/*********************************************************************************************************************************/ +double LinearAlgebra::calcSpearman(vector& x, vector& y, double& sig){ + try { + if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; } + + //format data + double sf = 0.0; //f^3 - f where f is the number of ties in x; + double sg = 0.0; //f^3 - f where f is the number of ties in y; + map tableX; + map::iterator itTable; + vector xscores; + + for (int i = 0; i < x.size(); i++) { + spearmanRank member(toString(i), x[i]); + xscores.push_back(member); + + //count number of repeats + itTable = tableX.find(x[i]); + if (itTable == tableX.end()) { + tableX[x[i]] = 1; + }else { + tableX[x[i]]++; + } + } + + + //calc LX + double Lx = 0.0; + for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) { + double tx = (double) itTable->second; + Lx += ((pow(tx, 3.0) - tx) / 12.0); + } + + + //sort x + sort(xscores.begin(), xscores.end(), compareSpearman); + + //convert scores to ranks of x + //convert to ranks + map rankx; + vector xties; + int rankTotal = 0; + for (int j = 0; j < xscores.size(); j++) { + rankTotal += (j+1); + xties.push_back(xscores[j]); + + if (j != xscores.size()-1) { // you are not the last so you can look ahead + if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue + for (int k = 0; k < xties.size(); k++) { + float thisrank = rankTotal / (float) xties.size(); + rankx[xties[k].name] = thisrank; + } + int t = xties.size(); + sf += (t*t*t-t); + xties.clear(); + rankTotal = 0; + } + }else { // you are the last one + for (int k = 0; k < xties.size(); k++) { + float thisrank = rankTotal / (float) xties.size(); + rankx[xties[k].name] = thisrank; + } + } + } + + //format x + vector yscores; + map tableY; + for (int j = 0; j < y.size(); j++) { + spearmanRank member(toString(j), y[j]); + yscores.push_back(member); + + itTable = tableY.find(member.score); + if (itTable == tableY.end()) { + tableY[member.score] = 1; + }else { + tableY[member.score]++; + } + + } + + //calc Ly + double Ly = 0.0; + for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) { + double ty = (double) itTable->second; + Ly += ((pow(ty, 3.0) - ty) / 12.0); + } + + sort(yscores.begin(), yscores.end(), compareSpearman); + + //convert to ranks + map rank; + vector yties; + rankTotal = 0; + for (int j = 0; j < yscores.size(); j++) { + rankTotal += (j+1); + yties.push_back(yscores[j]); + + if (j != yscores.size()-1) { // you are not the last so you can look ahead + if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue + for (int k = 0; k < yties.size(); k++) { + float thisrank = rankTotal / (float) yties.size(); + rank[yties[k].name] = thisrank; + } + int t = yties.size(); + sg += (t*t*t-t); + yties.clear(); + rankTotal = 0; + } + }else { // you are the last one + for (int k = 0; k < yties.size(); k++) { + float thisrank = rankTotal / (float) yties.size(); + rank[yties[k].name] = thisrank; + } + } + } + + double di = 0.0; + for (int k = 0; k < x.size(); k++) { + + float xi = rankx[toString(k)]; + float yi = rank[toString(k)]; + + di += ((xi - yi) * (xi - yi)); + } + + double p = 0.0; + + double n = (double) x.size(); + double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx; + double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly; + + p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2))); + + //Numerical Recipes 646 + sig = calcSpearmanSig(n, sf, sg, di); + + return p; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "calcSpearman"); + exit(1); + } +} +/*********************************************************************************************************************************/ +double LinearAlgebra::calcSpearmanSig(double n, double sf, double sg, double d){ + try { + + double sig = 0.0; + double probrs = 0.0; + double en=n; + double en3n=en*en*en-en; + double aved=en3n/6.0-(sf+sg)/12.0; + double fac=(1.0-sf/en3n)*(1.0-sg/en3n); + double vard=((en-1.0)*en*en*SQR(en+1.0)/36.0)*fac; + double zd=(d-aved)/sqrt(vard); + double probd=erfcc(fabs(zd)/1.4142136); + double rs=(1.0-(6.0/en3n)*(d+(sf+sg)/12.0))/sqrt(fac); + fac=(rs+1.0)*(1.0-rs); + if (fac > 0.0) { + double t=rs*sqrt((en-2.0)/fac); + double df=en-2.0; + probrs=betai(0.5*df,0.5,df/(df+t*t)); + }else { + probrs = 0.0; + } + + //smaller of probd and probrs is sig + sig = probrs; + if (probd < probrs) { sig = probd; } + + if (isnan(sig) || isinf(sig)) { sig = 0.0; } + + return sig; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "calcSpearmanSig"); + exit(1); + } +} +/*********************************************************************************************************************************/ +double LinearAlgebra::calcPearson(vector& x, vector& y, double& sig){ + try { + if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; } + + //find average X + float averageX = 0.0; + for (int i = 0; i < x.size(); i++) { averageX += x[i]; } + averageX = averageX / (float) x.size(); + + //find average Y + float sumY = 0.0; + for (int j = 0; j < y.size(); j++) { sumY += y[j]; } + float Ybar = sumY / (float) y.size(); + + double r = 0.0; + double numerator = 0.0; + double denomTerm1 = 0.0; + double denomTerm2 = 0.0; + + for (int j = 0; j < x.size(); j++) { + float Yi = y[j]; + float Xi = x[j]; + + numerator += ((Xi - averageX) * (Yi - Ybar)); + denomTerm1 += ((Xi - averageX) * (Xi - averageX)); + denomTerm2 += ((Yi - Ybar) * (Yi - Ybar)); + } + + double denom = (sqrt(denomTerm1) * sqrt(denomTerm2)); + + r = numerator / denom; + + //Numerical Recipes pg.644 + sig = calcPearsonSig(x.size(), r); + + return r; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "calcPearson"); + exit(1); + } +} +/*********************************************************************************************************************************/ +double LinearAlgebra::calcPearsonSig(double n, double r){ + try { + + double sig = 0.0; + const double TINY = 1.0e-20; + double z = 0.5*log((1.0+r+TINY)/(1.0-r+TINY)); //Fisher's z transformation + + //code below was giving an error in betacf with sop files + //int df = n-2; + //double t = r*sqrt(df/((1.0-r+TINY)*(1.0+r+TINY))); + //sig = betai(0.5+df, 0.5, df/(df+t*t)); + + //Numerical Recipes says code below gives approximately the same result + sig = erfcc(fabs(z*sqrt(n-1.0))/1.4142136); + if (isnan(sig) || isinf(sig)) { sig = 0.0; } + + return sig; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "calcPearsonSig"); + exit(1); + } +} +/*********************************************************************************************************************************/ + +vector > LinearAlgebra::getObservedEuclideanDistance(vector >& relAbundData){ + try { + + int numSamples = relAbundData.size(); + int numOTUs = relAbundData[0].size(); + + vector > dMatrix(numSamples); + for(int i=0;icontrol_pressed) { return dMatrix; } + + double d = 0; + for(int k=0;kerrorOut(e, "LinearAlgebra", "getObservedEuclideanDistance"); + exit(1); + } +} + +/*********************************************************************************************************************************/ +vector LinearAlgebra::solveEquations(vector > A, vector b){ + try { + int length = (int)b.size(); + vector x(length, 0); + vector index(length); + for(int i=0;icontrol_pressed) { return b; } + lubksb(A, index, b); + + return b; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "solveEquations"); + exit(1); + } +} +/*********************************************************************************************************************************/ +vector LinearAlgebra::solveEquations(vector > A, vector b){ + try { + int length = (int)b.size(); + vector x(length, 0); + vector index(length); + for(int i=0;icontrol_pressed) { return b; } + lubksb(A, index, b); + + return b; + } + catch(exception& e) { + m->errorOut(e, "LinearAlgebra", "solveEquations"); + exit(1); + } +} + +/*********************************************************************************************************************************/ + +void LinearAlgebra::ludcmp(vector >& A, vector& index, double& d){ + try { + double tiny = 1e-20; + + int n = (int)A.size(); + vector vv(n, 0.0); + double temp; + int imax; + + d = 1.0; + + for(int i=0;i big ) big=temp; } + if(big==0.0){ m->mothurOut("Singular matrix in routine ludcmp\n"); } + vv[i] = 1.0/big; + } + + for(int j=0;jcontrol_pressed) { break; } + for(int i=0;i= big){ + big = dum; + imax = i; + } + } + if(j != imax){ + for(int k=0;kerrorOut(e, "LinearAlgebra", "ludcmp"); + exit(1); + } +} + +/*********************************************************************************************************************************/ + +void LinearAlgebra::lubksb(vector >& A, vector& index, vector& b){ + try { + double total; + int n = (int)A.size(); + int ii = 0; + + for(int i=0;icontrol_pressed) { break; } + int ip = index[i]; + total = b[ip]; + b[ip] = b[i]; + + if (ii != 0) { + for(int j=ii-1;j=0;i--){ + total = b[i]; + for(int j=i+1;jerrorOut(e, "LinearAlgebra", "lubksb"); + exit(1); + } +} +/*********************************************************************************************************************************/ + +void LinearAlgebra::ludcmp(vector >& A, vector& index, float& d){ + try { + double tiny = 1e-20; + + int n = (int)A.size(); + vector vv(n, 0.0); + double temp; + int imax; + + d = 1.0; + + for(int i=0;i big ) big=temp; } + if(big==0.0){ m->mothurOut("Singular matrix in routine ludcmp\n"); } + vv[i] = 1.0/big; + } + + for(int j=0;jcontrol_pressed) { break; } + for(int i=0;i= big){ + big = dum; + imax = i; + } + } + if(j != imax){ + for(int k=0;kerrorOut(e, "LinearAlgebra", "ludcmp"); + exit(1); + } +} + +/*********************************************************************************************************************************/ + +void LinearAlgebra::lubksb(vector >& A, vector& index, vector& b){ + try { + float total; + int n = (int)A.size(); + int ii = 0; + + for(int i=0;icontrol_pressed) { break; } + int ip = index[i]; + total = b[ip]; + b[ip] = b[i]; + + if (ii != 0) { + for(int j=ii-1;j=0;i--){ + total = b[i]; + for(int j=i+1;jerrorOut(e, "LinearAlgebra", "lubksb"); + exit(1); + } +} + +/*********************************************************************************************************************************/ + +vector > LinearAlgebra::getInverse(vector > matrix){ + try { + int n = (int)matrix.size(); + + vector > inverse(n); + for(int i=0;i column(n, 0.0000); + vector index(n, 0); + double dummy; + + ludcmp(matrix, index, dummy); + + for(int j=0;jcontrol_pressed) { break; } + + column.assign(n, 0); + + column[j] = 1.0000; + + lubksb(matrix, index, column); + + for(int i=0;ierrorOut(e, "LinearAlgebra", "getInverse"); + exit(1); + } +}