5 * Created by westcott on 1/7/11.
6 * Copyright 2011 Schloss Lab. All rights reserved.
10 #include "linearalgebra.h"
12 /*********************************************************************************************************************************/
14 inline double SIGN(const double a, const double b)
16 return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a);
18 /*********************************************************************************************************************************/
20 vector<vector<double> > LinearAlgebra::matrix_mult(vector<vector<double> > first, vector<vector<double> > second){
22 vector<vector<double> > product;
24 int first_rows = first.size();
25 int first_cols = first[0].size();
26 int second_cols = second[0].size();
28 product.resize(first_rows);
29 for(int i=0;i<first_rows;i++){
30 product[i].resize(second_cols);
33 for(int i=0;i<first_rows;i++){
34 for(int j=0;j<second_cols;j++){
36 if (m->control_pressed) { return product; }
39 for(int k=0;k<first_cols;k++){
40 product[i][j] += first[i][k] * second[k][j];
48 m->errorOut(e, "LinearAlgebra", "matrix_mult");
54 /*********************************************************************************************************************************/
56 void LinearAlgebra::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
60 vector<vector<double> > A(rank);
61 vector<vector<double> > C(rank);
62 for(int i=0;i<rank;i++){
67 double scale = -1.0000 / (double) rank;
69 for(int i=0;i<rank;i++){
71 C[i][i] = 1.0000 + scale;
72 for(int j=i+1;j<rank;j++){
73 A[i][j] = A[j][i] = -0.5 * D[i][j] * D[i][j] + offset;
74 C[i][j] = C[j][i] = scale;
82 m->errorOut(e, "LinearAlgebra", "recenter");
87 /*********************************************************************************************************************************/
89 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
91 int LinearAlgebra::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
93 double scale, hh, h, g, f;
100 for(int i=n-1;i>0;i--){
104 for(int k=0;k<l+1;k++){
105 scale += fabs(a[i][k]);
111 for(int k=0;k<l+1;k++){
113 h += a[i][k] * a[i][k];
116 g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
121 for(int j=0;j<l+1;j++){
122 a[j][i] = a[i][j] / h;
124 for(int k=0;k<j+1;k++){
125 g += a[j][k] * a[i][k];
127 for(int k=j+1;k<l+1;k++){
128 g += a[k][j] * a[i][k];
134 for(int j=0;j<l+1;j++){
136 e[j] = g = e[j] - hh * f;
137 for(int k=0;k<j+1;k++){
138 a[j][k] -= (f * e[k] + g * a[i][k]);
153 for(int i=0;i<n;i++){
156 for(int j=0;j<l;j++){
158 for(int k=0;k<l;k++){
159 g += a[i][k] * a[k][j];
161 for(int k=0;k<l;k++){
162 a[k][j] -= g * a[k][i];
168 for(int j=0;j<l;j++){
169 a[j][i] = a[i][j] = 0.0;
175 catch(exception& e) {
176 m->errorOut(e, "LinearAlgebra", "tred2");
181 /*********************************************************************************************************************************/
183 double LinearAlgebra::pythag(double a, double b) { return(pow(a*a+b*b,0.5)); }
185 /*********************************************************************************************************************************/
187 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
189 int LinearAlgebra::qtli(vector<double>& d, vector<double>& e, vector<vector<double> >& z) {
192 double s, r, p, g, f, dd, c, b;
195 for(int i=1;i<=n;i++){
200 for(int l=0;l<n;l++){
203 for(myM=l;myM<n-1;myM++){
204 dd = fabs(d[myM]) + fabs(d[myM+1]);
205 if(fabs(e[myM])+dd == dd) break;
208 if(iter++ == 3000) cerr << "Too many iterations in tqli\n";
209 g = (d[l+1]-d[l]) / (2.0 * e[l]);
211 g = d[myM] - d[l] + e[l] / (g + SIGN(r,g));
214 for(i=myM-1;i>=l;i--){
217 e[i+1] = (r=pythag(f,g));
226 r = (d[i] - g) * s + 2.0 * c * b;
227 d[i+1] = g + ( p = s * r);
229 for(int k=0;k<n;k++){
231 z[k][i+1] = s * z[k][i] + c * f;
232 z[k][i] = c * z[k][i] - s * f;
235 if(r == 0.00 && i >= l) continue;
244 for(int i=0;i<n;i++){
246 for(int j=i;j<n;j++){
254 for(int j=0;j<n;j++){
264 catch(exception& e) {
265 m->errorOut(e, "LinearAlgebra", "qtli");
269 /*********************************************************************************************************************************/
270 //groups by dimension
271 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes, int dimensions){
274 vector< vector<double> > dists; dists.resize(axes.size());
275 for (int i = 0; i < dists.size(); i++) { dists[i].resize(axes.size(), 0.0); }
277 if (dimensions == 1) { //one dimension calc = abs(x-y)
279 for (int i = 0; i < dists.size(); i++) {
281 if (m->control_pressed) { return dists; }
283 for (int j = 0; j < i; j++) {
284 dists[i][j] = abs(axes[i][0] - axes[j][0]);
285 dists[j][i] = dists[i][j];
289 }else if (dimensions > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
291 for (int i = 0; i < dists.size(); i++) {
293 if (m->control_pressed) { return dists; }
295 for (int j = 0; j < i; j++) {
297 for (int k = 0; k < dimensions; k++) {
298 sum += ((axes[i][k] - axes[j][k]) * (axes[i][k] - axes[j][k]));
301 dists[i][j] = sqrt(sum);
302 dists[j][i] = dists[i][j];
310 catch(exception& e) {
311 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
315 /*********************************************************************************************************************************/
316 //returns groups by dimensions from dimensions by groups
317 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes){
320 vector< vector<double> > dists; dists.resize(axes[0].size());
321 for (int i = 0; i < dists.size(); i++) { dists[i].resize(axes[0].size(), 0.0); }
323 if (axes.size() == 1) { //one dimension calc = abs(x-y)
325 for (int i = 0; i < dists.size(); i++) {
327 if (m->control_pressed) { return dists; }
329 for (int j = 0; j < i; j++) {
330 dists[i][j] = abs(axes[0][i] - axes[0][j]);
331 dists[j][i] = dists[i][j];
335 }else if (axes.size() > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
337 for (int i = 0; i < dists[0].size(); i++) {
339 if (m->control_pressed) { return dists; }
341 for (int j = 0; j < i; j++) {
343 for (int k = 0; k < axes.size(); k++) {
344 sum += ((axes[k][i] - axes[k][j]) * (axes[k][i] - axes[k][j]));
347 dists[i][j] = sqrt(sum);
348 dists[j][i] = dists[i][j];
356 catch(exception& e) {
357 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
361 /*********************************************************************************************************************************/
362 //assumes both matrices are square and the same size
363 double LinearAlgebra::calcPearson(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
366 //find average for - X
368 float averageEuclid = 0.0;
369 for (int i = 0; i < euclidDists.size(); i++) {
370 for (int j = 0; j < i; j++) {
371 averageEuclid += euclidDists[i][j];
375 averageEuclid = averageEuclid / (float) count;
377 //find average for - Y
379 float averageUser = 0.0;
380 for (int i = 0; i < userDists.size(); i++) {
381 for (int j = 0; j < i; j++) {
382 averageUser += userDists[i][j];
386 averageUser = averageUser / (float) count;
388 double numerator = 0.0;
389 double denomTerm1 = 0.0;
390 double denomTerm2 = 0.0;
392 for (int i = 0; i < euclidDists.size(); i++) {
394 for (int k = 0; k < i; k++) { //just lt dists
396 float Yi = userDists[i][k];
397 float Xi = euclidDists[i][k];
399 numerator += ((Xi - averageEuclid) * (Yi - averageUser));
400 denomTerm1 += ((Xi - averageEuclid) * (Xi - averageEuclid));
401 denomTerm2 += ((Yi - averageUser) * (Yi - averageUser));
405 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
406 double r = numerator / denom;
408 //divide by zero error
409 if (isnan(r) || isinf(r)) { r = 0.0; }
414 catch(exception& e) {
415 m->errorOut(e, "LinearAlgebra", "calcPearson");
419 /*********************************************************************************************************************************/
420 //assumes both matrices are square and the same size
421 double LinearAlgebra::calcSpearman(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
426 map<float, int> tableX;
427 map<float, int>::iterator itTable;
428 vector<spearmanRank> scores;
430 for (int i = 0; i < euclidDists.size(); i++) {
431 for (int j = 0; j < i; j++) {
432 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
433 scores.push_back(member);
435 //count number of repeats
436 itTable = tableX.find(euclidDists[i][j]);
437 if (itTable == tableX.end()) {
438 tableX[euclidDists[i][j]] = 1;
440 tableX[euclidDists[i][j]]++;
446 sort(scores.begin(), scores.end(), compareSpearman);
450 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
451 double tx = (double) itTable->second;
452 Lx += ((pow(tx, 3.0) - tx) / 12.0);
456 map<string, float> rankEuclid;
457 vector<spearmanRank> ties;
459 for (int j = 0; j < scores.size(); j++) {
461 ties.push_back(scores[j]);
463 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
464 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
466 for (int k = 0; k < ties.size(); k++) {
467 float thisrank = rankTotal / (float) ties.size();
468 rankEuclid[ties[k].name] = thisrank;
473 }else { // you are the last one
475 for (int k = 0; k < ties.size(); k++) {
476 float thisrank = rankTotal / (float) ties.size();
477 rankEuclid[ties[k].name] = thisrank;
484 map<float, int> tableY;
487 for (int i = 0; i < userDists.size(); i++) {
488 for (int j = 0; j < i; j++) {
489 spearmanRank member(toString(scores.size()), userDists[i][j]);
490 scores.push_back(member);
492 //count number of repeats
493 itTable = tableY.find(userDists[i][j]);
494 if (itTable == tableY.end()) {
495 tableY[userDists[i][j]] = 1;
497 tableY[userDists[i][j]]++;
503 sort(scores.begin(), scores.end(), compareSpearman);
507 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
508 double ty = (double) itTable->second;
509 Ly += ((pow(ty, 3.0) - ty) / 12.0);
513 map<string, float> rankUser;
516 for (int j = 0; j < scores.size(); j++) {
518 ties.push_back(scores[j]);
520 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
521 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
523 for (int k = 0; k < ties.size(); k++) {
524 float thisrank = rankTotal / (float) ties.size();
525 rankUser[ties[k].name] = thisrank;
530 }else { // you are the last one
532 for (int k = 0; k < ties.size(); k++) {
533 float thisrank = rankTotal / (float) ties.size();
534 rankUser[ties[k].name] = thisrank;
542 for (int i = 0; i < userDists.size(); i++) {
543 for (int j = 0; j < i; j++) {
545 float xi = rankEuclid[toString(count)];
546 float yi = rankUser[toString(count)];
548 di += ((xi - yi) * (xi - yi));
554 double n = (double) count;
556 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
557 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
559 r = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
561 //divide by zero error
562 if (isnan(r) || isinf(r)) { r = 0.0; }
567 catch(exception& e) {
568 m->errorOut(e, "LinearAlgebra", "calcSpearman");
573 /*********************************************************************************************************************************/
574 //assumes both matrices are square and the same size
575 double LinearAlgebra::calcKendall(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
580 vector<spearmanRank> scores;
581 for (int i = 0; i < euclidDists.size(); i++) {
582 for (int j = 0; j < i; j++) {
583 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
584 scores.push_back(member);
589 sort(scores.begin(), scores.end(), compareSpearman);
592 map<string, float> rankEuclid;
593 vector<spearmanRank> ties;
595 for (int j = 0; j < scores.size(); j++) {
597 ties.push_back(scores[j]);
599 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
600 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
602 for (int k = 0; k < ties.size(); k++) {
603 float thisrank = rankTotal / (float) ties.size();
604 rankEuclid[ties[k].name] = thisrank;
609 }else { // you are the last one
611 for (int k = 0; k < ties.size(); k++) {
612 float thisrank = rankTotal / (float) ties.size();
613 rankEuclid[ties[k].name] = thisrank;
618 vector<spearmanRank> scoresUser;
619 for (int i = 0; i < userDists.size(); i++) {
620 for (int j = 0; j < i; j++) {
621 spearmanRank member(toString(scoresUser.size()), userDists[i][j]);
622 scoresUser.push_back(member);
627 sort(scoresUser.begin(), scoresUser.end(), compareSpearman);
630 map<string, float> rankUser;
633 for (int j = 0; j < scoresUser.size(); j++) {
635 ties.push_back(scoresUser[j]);
637 if (j != (scoresUser.size()-1)) { // you are not the last so you can look ahead
638 if (scoresUser[j].score != scoresUser[j+1].score) { // you are done with ties, rank them and continue
640 for (int k = 0; k < ties.size(); k++) {
641 float thisrank = rankTotal / (float) ties.size();
642 rankUser[ties[k].name] = thisrank;
647 }else { // you are the last one
649 for (int k = 0; k < ties.size(); k++) {
650 float thisrank = rankTotal / (float) ties.size();
651 rankUser[ties[k].name] = thisrank;
660 vector<spearmanRank> user;
661 for (int l = 0; l < scores.size(); l++) {
662 spearmanRank member(scores[l].name, rankUser[scores[l].name]);
663 user.push_back(member);
667 for (int l = 0; l < scores.size(); l++) {
669 int numWithHigherRank = 0;
670 int numWithLowerRank = 0;
671 float thisrank = user[l].score;
673 for (int u = l+1; u < scores.size(); u++) {
674 if (user[u].score > thisrank) { numWithHigherRank++; }
675 else if (user[u].score < thisrank) { numWithLowerRank++; }
679 numCoor += numWithHigherRank;
680 numDisCoor += numWithLowerRank;
683 r = (numCoor - numDisCoor) / (float) count;
685 //divide by zero error
686 if (isnan(r) || isinf(r)) { r = 0.0; }
691 catch(exception& e) {
692 m->errorOut(e, "LinearAlgebra", "calcKendall");
696 /*********************************************************************************************************************************/
697 double LinearAlgebra::calcKendall(vector<double>& x, vector<double>& y, double& sig){
699 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
702 vector<spearmanRank> xscores;
703 for (int i = 0; i < x.size(); i++) {
704 spearmanRank member(toString(i), x[i]);
705 xscores.push_back(member);
709 sort(xscores.begin(), xscores.end(), compareSpearman);
711 //convert scores to ranks of x
712 vector<spearmanRank*> ties;
714 for (int j = 0; j < xscores.size(); j++) {
716 ties.push_back(&(xscores[j]));
718 if (j != xscores.size()-1) { // you are not the last so you can look ahead
719 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
720 for (int k = 0; k < ties.size(); k++) {
721 float thisrank = rankTotal / (float) ties.size();
722 (*ties[k]).score = thisrank;
727 }else { // you are the last one
728 for (int k = 0; k < ties.size(); k++) {
729 float thisrank = rankTotal / (float) ties.size();
730 (*ties[k]).score = thisrank;
737 vector<spearmanRank> yscores;
738 for (int j = 0; j < y.size(); j++) {
739 spearmanRank member(toString(j), y[j]);
740 yscores.push_back(member);
744 sort(yscores.begin(), yscores.end(), compareSpearman);
747 map<string, float> rank;
748 vector<spearmanRank> yties;
750 for (int j = 0; j < yscores.size(); j++) {
752 yties.push_back(yscores[j]);
754 if (j != yscores.size()-1) { // you are not the last so you can look ahead
755 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
756 for (int k = 0; k < yties.size(); k++) {
757 float thisrank = rankTotal / (float) yties.size();
758 rank[yties[k].name] = thisrank;
763 }else { // you are the last one
764 for (int k = 0; k < yties.size(); k++) {
765 float thisrank = rankTotal / (float) yties.size();
766 rank[yties[k].name] = thisrank;
776 vector<spearmanRank> otus;
777 for (int l = 0; l < xscores.size(); l++) {
778 spearmanRank member(xscores[l].name, rank[xscores[l].name]);
779 otus.push_back(member);
783 for (int l = 0; l < xscores.size(); l++) {
785 int numWithHigherRank = 0;
786 int numWithLowerRank = 0;
787 float thisrank = otus[l].score;
789 for (int u = l+1; u < xscores.size(); u++) {
790 if (otus[u].score > thisrank) { numWithHigherRank++; }
791 else if (otus[u].score < thisrank) { numWithLowerRank++; }
795 numCoor += numWithHigherRank;
796 numDisCoor += numWithLowerRank;
799 double p = (numCoor - numDisCoor) / (float) count;
801 //calc signif - zA - http://en.wikipedia.org/wiki/Kendall_tau_rank_correlation_coefficient#Significance_tests
802 double numer = 3.0 * (numCoor - numDisCoor);
803 int n = xscores.size();
804 double denom = n * (n-1) * (2*n + 5) / (double) 2.0;
808 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
812 catch(exception& e) {
813 m->errorOut(e, "LinearAlgebra", "calcKendall");
817 /*********************************************************************************************************************************/
818 double LinearAlgebra::calcSpearman(vector<double>& x, vector<double>& y, double& sig){
820 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
823 map<float, int> tableX;
824 map<float, int>::iterator itTable;
825 vector<spearmanRank> xscores;
827 for (int i = 0; i < x.size(); i++) {
828 spearmanRank member(toString(i), x[i]);
829 xscores.push_back(member);
831 //count number of repeats
832 itTable = tableX.find(x[i]);
833 if (itTable == tableX.end()) {
843 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
844 double tx = (double) itTable->second;
845 Lx += ((pow(tx, 3.0) - tx) / 12.0);
850 sort(xscores.begin(), xscores.end(), compareSpearman);
852 //convert scores to ranks of x
854 map<string, float> rankx;
855 vector<spearmanRank> xties;
857 for (int j = 0; j < xscores.size(); j++) {
859 xties.push_back(xscores[j]);
861 if (j != xscores.size()-1) { // you are not the last so you can look ahead
862 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
863 for (int k = 0; k < xties.size(); k++) {
864 float thisrank = rankTotal / (float) xties.size();
865 rankx[xties[k].name] = thisrank;
870 }else { // you are the last one
871 for (int k = 0; k < xties.size(); k++) {
872 float thisrank = rankTotal / (float) xties.size();
873 rankx[xties[k].name] = thisrank;
879 vector<spearmanRank> yscores;
880 map<float, int> tableY;
881 for (int j = 0; j < y.size(); j++) {
882 spearmanRank member(toString(j), y[j]);
883 yscores.push_back(member);
885 itTable = tableY.find(member.score);
886 if (itTable == tableY.end()) {
887 tableY[member.score] = 1;
889 tableY[member.score]++;
896 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
897 double ty = (double) itTable->second;
898 Ly += ((pow(ty, 3.0) - ty) / 12.0);
901 sort(yscores.begin(), yscores.end(), compareSpearman);
904 map<string, float> rank;
905 vector<spearmanRank> yties;
907 for (int j = 0; j < yscores.size(); j++) {
909 yties.push_back(yscores[j]);
911 if (j != yscores.size()-1) { // you are not the last so you can look ahead
912 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
913 for (int k = 0; k < yties.size(); k++) {
914 float thisrank = rankTotal / (float) yties.size();
915 rank[yties[k].name] = thisrank;
920 }else { // you are the last one
921 for (int k = 0; k < yties.size(); k++) {
922 float thisrank = rankTotal / (float) yties.size();
923 rank[yties[k].name] = thisrank;
929 for (int k = 0; k < x.size(); k++) {
931 float xi = rankx[toString(k)];
932 float yi = rank[toString(k)];
934 di += ((xi - yi) * (xi - yi));
939 double n = (double) x.size();
940 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
941 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
943 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
945 //signifigance calc - http://en.wikipedia.org/wiki/Spearman%27s_rank_correlation_coefficient
946 double temp = (x.size()-2) / (double) (1- (p*p));
949 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
953 catch(exception& e) {
954 m->errorOut(e, "LinearAlgebra", "calcSpearman");
958 /*********************************************************************************************************************************/
959 double LinearAlgebra::calcPearson(vector<double>& x, vector<double>& y, double& sig){
961 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
964 float averageX = 0.0;
965 for (int i = 0; i < x.size(); i++) { averageX += x[i]; }
966 averageX = averageX / (float) x.size();
970 for (int j = 0; j < y.size(); j++) { sumY += y[j]; }
971 float Ybar = sumY / (float) y.size();
974 double numerator = 0.0;
975 double denomTerm1 = 0.0;
976 double denomTerm2 = 0.0;
978 for (int j = 0; j < x.size(); j++) {
982 numerator += ((Xi - averageX) * (Yi - Ybar));
983 denomTerm1 += ((Xi - averageX) * (Xi - averageX));
984 denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
987 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
989 r = numerator / denom;
991 //signifigance calc - http://faculty.vassar.edu/lowry/ch4apx.html
992 double temp = (1- (r*r)) / (double) (x.size()-2);
995 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
999 catch(exception& e) {
1000 m->errorOut(e, "LinearAlgebra", "calcPearson");
1004 /*********************************************************************************************************************************/
1006 vector<vector<double> > LinearAlgebra::getObservedEuclideanDistance(vector<vector<double> >& relAbundData){
1008 int numSamples = relAbundData.size();
1009 int numOTUs = relAbundData[0].size();
1011 vector<vector<double> > dMatrix(numSamples);
1012 for(int i=0;i<numSamples;i++){
1013 dMatrix[i].resize(numSamples);
1016 for(int i=0;i<numSamples;i++){
1017 for(int j=0;j<numSamples;j++){
1020 for(int k=0;k<numOTUs;k++){
1021 d += pow((relAbundData[i][k] - relAbundData[j][k]), 2.0000);
1023 dMatrix[i][j] = pow(d, 0.50000);
1024 dMatrix[j][i] = dMatrix[i][j];
1032 /*********************************************************************************************************************************/