+ if (isnan(sig) || isinf(sig)) { sig = 0.0; }
+
+ return sig;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "LinearAlgebra", "calcKendallSig");
+ exit(1);
+ }
+}
+/*********************************************************************************************************************************/
+double LinearAlgebra::calcWilcoxon(vector<double>& x, vector<double>& y, double& sig){
+ try {
+ if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
+
+ double W = 0.0;
+ sig = 0.0;
+
+ vector<double> signPairs;
+ vector<spearmanRank> absV;
+ for (int i = 0; i < x.size(); i++) {
+ if (m->control_pressed) { return W; }
+ double temp = y[i]-x[i];
+ double signV = sign(temp);
+ if (signV != 0) { // exclude zeros
+ spearmanRank member(toString(i), abs(temp));
+ absV.push_back(member);
+ signPairs.push_back(signV);
+ }
+ }
+
+ //rank absV
+ //sort xscores
+ sort(absV.begin(), absV.end(), compareSpearman);
+
+ //convert scores to ranks of x
+ vector<spearmanRank*> ties;
+ int rankTotal = 0;
+ for (int j = 0; j < absV.size(); j++) {
+ if (m->control_pressed) { return W; }
+ rankTotal += (j+1);
+ ties.push_back(&(absV[j]));
+
+ if (j != absV.size()-1) { // you are not the last so you can look ahead
+ if (absV[j].score != absV[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;
+ }
+ }
+ }
+
+ //sum ranks times sign
+ for (int i = 0; i < absV.size(); i++) {
+ if (m->control_pressed) { return W; }
+ W += (absV[i].score*signPairs[i]);
+ }
+ W = abs(W);
+
+ //find zScore
+ cout << "still need to find sig!!" << endl;
+
+ return W;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "LinearAlgebra", "calcWilcoxon");
+ exit(1);
+ }
+}
+/*********************************************************************************************************************************/
+double LinearAlgebra::calcSpearman(vector<double>& x, vector<double>& 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<float, int> tableX;
+ map<float, int>::iterator itTable;
+ vector<spearmanRank> 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<string, float> rankx;
+ vector<spearmanRank> 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<spearmanRank> yscores;
+ map<float, int> 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<string, float> rank;
+ vector<spearmanRank> 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<double>& x, vector<double>& 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<vector<double> > LinearAlgebra::getObservedEuclideanDistance(vector<vector<double> >& relAbundData){
+ try {
+
+ int numSamples = relAbundData.size();
+ int numOTUs = relAbundData[0].size();
+
+ vector<vector<double> > dMatrix(numSamples);
+ for(int i=0;i<numSamples;i++){
+ dMatrix[i].resize(numSamples);
+ }
+
+ for(int i=0;i<numSamples;i++){
+ for(int j=0;j<numSamples;j++){
+
+ if (m->control_pressed) { return dMatrix; }
+
+ double d = 0;
+ for(int k=0;k<numOTUs;k++){
+ d += pow((relAbundData[i][k] - relAbundData[j][k]), 2.0000);
+ }
+ dMatrix[i][j] = pow(d, 0.50000);
+ dMatrix[j][i] = dMatrix[i][j];
+
+ }
+ }
+ return dMatrix;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "LinearAlgebra", "getObservedEuclideanDistance");
+ exit(1);
+ }
+}
+
+/*********************************************************************************************************************************/
+vector<double> LinearAlgebra::solveEquations(vector<vector<double> > A, vector<double> b){
+ try {
+ int length = (int)b.size();
+ vector<double> x(length, 0);
+ vector<int> index(length);
+ for(int i=0;i<length;i++){ index[i] = i; }
+ double d;
+
+ ludcmp(A, index, d); if (m->control_pressed) { return b; }
+ lubksb(A, index, b);
+
+ return b;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "LinearAlgebra", "solveEquations");
+ exit(1);
+ }
+}
+/*********************************************************************************************************************************/
+vector<float> LinearAlgebra::solveEquations(vector<vector<float> > A, vector<float> b){
+ try {
+ int length = (int)b.size();
+ vector<double> x(length, 0);
+ vector<int> index(length);
+ for(int i=0;i<length;i++){ index[i] = i; }
+ float d;
+
+ ludcmp(A, index, d); if (m->control_pressed) { return b; }
+ lubksb(A, index, b);
+
+ return b;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "LinearAlgebra", "solveEquations");
+ exit(1);
+ }
+}
+
+/*********************************************************************************************************************************/
+
+void LinearAlgebra::ludcmp(vector<vector<double> >& A, vector<int>& index, double& d){
+ try {
+ double tiny = 1e-20;
+
+ int n = (int)A.size();
+ vector<double> vv(n, 0.0);
+ double temp;
+ int imax;
+
+ d = 1.0;
+
+ for(int i=0;i<n;i++){
+ double big = 0.0;
+ for(int j=0;j<n;j++){ if((temp=fabs(A[i][j])) > big ) big=temp; }
+ if(big==0.0){ m->mothurOut("Singular matrix in routine ludcmp\n"); }
+ vv[i] = 1.0/big;
+ }
+
+ for(int j=0;j<n;j++){
+ if (m->control_pressed) { break; }
+ for(int i=0;i<j;i++){
+ double sum = A[i][j];
+ for(int k=0;k<i;k++){ sum -= A[i][k] * A[k][j]; }
+ A[i][j] = sum;
+ }
+
+ double big = 0.0;
+ for(int i=j;i<n;i++){
+ double sum = A[i][j];
+ for(int k=0;k<j;k++){ sum -= A[i][k] * A[k][j]; }
+ A[i][j] = sum;
+ double dum;
+ if((dum = vv[i] * fabs(sum)) >= big){
+ big = dum;
+ imax = i;
+ }
+ }
+ if(j != imax){
+ for(int k=0;k<n;k++){
+ double dum = A[imax][k];
+ A[imax][k] = A[j][k];
+ A[j][k] = dum;
+ }
+ d = -d;
+ vv[imax] = vv[j];
+ }
+ index[j] = imax;
+
+ if(A[j][j] == 0.0){ A[j][j] = tiny; }
+
+ if(j != n-1){
+ double dum = 1.0/A[j][j];
+ for(int i=j+1;i<n;i++){ A[i][j] *= dum; }
+ }
+ }
+ }
+ catch(exception& e) {
+ m->errorOut(e, "LinearAlgebra", "ludcmp");
+ exit(1);
+ }
+}
+
+/*********************************************************************************************************************************/
+
+void LinearAlgebra::lubksb(vector<vector<double> >& A, vector<int>& index, vector<double>& b){
+ try {
+ double total;
+ int n = (int)A.size();
+ int ii = 0;
+
+ for(int i=0;i<n;i++){
+ if (m->control_pressed) { break; }
+ int ip = index[i];
+ total = b[ip];
+ b[ip] = b[i];
+
+ if (ii != 0) {
+ for(int j=ii-1;j<i;j++){
+ total -= A[i][j] * b[j];
+ }
+ }
+ else if(total != 0){ ii = i+1; }
+ b[i] = total;
+ }
+ for(int i=n-1;i>=0;i--){
+ total = b[i];
+ for(int j=i+1;j<n;j++){ total -= A[i][j] * b[j]; }
+ b[i] = total / A[i][i];
+ }
+ }
+ catch(exception& e) {
+ m->errorOut(e, "LinearAlgebra", "lubksb");
+ exit(1);
+ }
+}
+/*********************************************************************************************************************************/