linearalgebra.cpp

   1 /*
   2  *  linearalgebra.cpp
   3  *  mothur
   4  *
   5  *  Created by westcott on 1/7/11.
   6  *  Copyright 2011 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "linearalgebra.h"
  11
  12 // This class references functions used from "Numerical Recipes in C++" //
  13
  14 /*********************************************************************************************************************************/
  15 inline double SQR(const double a)
  16 {
  17     return a*a;
  18 }
  19 /*********************************************************************************************************************************/
  20
  21 inline double SIGN(const double a, const double b)
  22 {
  23     return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a);
  24 }
  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){
  28     try {
  29         double t,z,ans;
  30         z=fabs(x);
  31         t=1.0/(1.0+0.5*z);
  32
  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)))))))));
  36
  37         //cout << "in erfcc " << t << '\t' << ans<< endl;
  38         return (x >= 0.0 ? ans : 2.0 - ans);
  39     }
  40         catch(exception& e) {
  41                 m->errorOut(e, "LinearAlgebra", "betai");
  42                 exit(1);
  43         }
  44 }
  45 /*********************************************************************************************************************************/
  46 //Numerical Recipes pg. 232
  47 double LinearAlgebra::betai(const double a, const double b, const double x) {
  48     try {
  49         double bt;
  50         double result = 0.0;
  51
  52         if (x < 0.0 || x > 1.0) { m->mothurOut("[ERROR]: bad x in betai.\n"); m->control_pressed = true; return 0.0; }
  53
  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));  }
  56
  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; }
  59
  60         return result;
  61     }
  62         catch(exception& e) {
  63                 m->errorOut(e, "LinearAlgebra", "betai");
  64                 exit(1);
  65         }
  66 }
  67 /*********************************************************************************************************************************/
  68 //Numerical Recipes pg. 219
  69 double LinearAlgebra::gammln(const double xx) {
  70     try {
  71         int j;
  72         double x,y,tmp,ser;
  73         static const double cof[6]={76.18009172947146,-86.50532032941677,24.01409824083091,
  74             -1.231739572450155,0.120858003e-2,-0.536382e-5};
  75
  76         y=x=xx;
  77         tmp=x+5.5;
  78         tmp -= (x+0.5)*log(tmp);
  79         ser=1.0;
  80         for (j=0;j<6;j++) {
  81             ser += cof[j]/++y;
  82         }
  83         return -tmp+log(2.5066282746310005*ser/x);
  84     }
  85         catch(exception& e) {
  86                 m->errorOut(e, "LinearAlgebra", "gammln");
  87                 exit(1);
  88         }
  89 }
  90 /*********************************************************************************************************************************/
  91 //Numerical Recipes pg. 223
  92 double LinearAlgebra::gammp(const double a, const double x) {
  93     try {
  94         double gamser,gammcf,gln;
  95
  96         if (x < 0.0 || a <= 0.0) { m->mothurOut("[ERROR]: Invalid arguments in routine GAMMP\n"); m->control_pressed = true; return 0.0;}
  97         if (x < (a+1.0)) {
  98             gser(gamser,a,x,gln);
  99             return gamser;
 100         } else {
 101             gcf(gammcf,a,x,gln);
 102             return 1.0-gammcf;
 103         }
 104
 105         return 0;
 106     }
 107         catch(exception& e) {
 108                 m->errorOut(e, "LinearAlgebra", "gammp");
 109                 exit(1);
 110         }
 111 }
 112 /*********************************************************************************************************************************/
 113 //Numerical Recipes pg. 223
 114 double LinearAlgebra::gammq(const double a, const double x) {
 115     try {
 116         double gamser,gammcf,gln;
 117
 118         if (x < 0.0 || a <= 0.0) { m->mothurOut("[ERROR]: Invalid arguments in routine GAMMQ\n"); m->control_pressed = true; return 0.0; }
 119         if (x < (a+1.0)) {
 120             gser(gamser,a,x,gln);
 121             return 1.0-gamser;
 122         } else {
 123             gcf(gammcf,a,x,gln);
 124             return gammcf;
 125         }
 126
 127         return 0;
 128     }
 129         catch(exception& e) {
 130                 m->errorOut(e, "LinearAlgebra", "gammp");
 131                 exit(1);
 132         }
 133 }
 134 /*********************************************************************************************************************************/
 135 //Numerical Recipes pg. 224
 136 double LinearAlgebra::gcf(double& gammcf, const double a, const double x, double& gln){
 137     try {
 138         const int ITMAX=100;
 139         const double EPS=numeric_limits<double>::epsilon();
 140         const double FPMIN=numeric_limits<double>::min()/EPS;
 141         int i;
 142         double an,b,c,d,del,h;
 143
 144         gln=gammln(a);
 145         b=x+1.0-a;
 146         c=1.0/FPMIN;
 147         d=1.0/b;
 148         h=d;
 149         for (i=1;i<=ITMAX;i++) {
 150             an = -i*(i-a);
 151             b += 2.0;
 152             d=an*d+b;
 153             if (fabs(d) < FPMIN) { d=FPMIN; }
 154             c=b+an/c;
 155             if (fabs(c) < FPMIN) { c=FPMIN; }
 156             d=1.0/d;
 157             del=d*c;
 158             h *= del;
 159             if (fabs(del-1.0) <= EPS) break;
 160         }
 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;
 163
 164         return 0.0;
 165     }
 166         catch(exception& e) {
 167                 m->errorOut(e, "LinearAlgebra", "gcf");
 168                 exit(1);
 169         }
 170
 171 }
 172 /*********************************************************************************************************************************/
 173 //Numerical Recipes pg. 223
 174 double LinearAlgebra::gser(double& gamser, const double a, const double x, double& gln) {
 175     try {
 176         int n;
 177         double sum,del,ap;
 178         const double EPS = numeric_limits<double>::epsilon();
 179
 180         gln=gammln(a);
 181         if (x <= 0.0) {
 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;
 184         } else {
 185             ap=a;
 186             del=sum=1.0/a;
 187             for (n=0;n<100;n++) {
 188                 ++ap;
 189                 del *= x/ap;
 190                 sum += del;
 191                 if (fabs(del) < fabs(sum)*EPS) {
 192                     gamser=sum*exp(-x+a*log(x)-gln);
 193                     return 0.0;
 194                 }
 195             }
 196
 197             m->mothurOut("[ERROR]: a too large, ITMAX too small in routine GSER\n");
 198             return 0.0;
 199         }
 200         return 0;
 201     }
 202         catch(exception& e) {
 203                 m->errorOut(e, "LinearAlgebra", "gser");
 204                 exit(1);
 205         }
 206 }
 207 /*********************************************************************************************************************************/
 208 //Numerical Recipes pg. 233
 209 double LinearAlgebra::betacf(const double a, const double b, const double x) {
 210     try {
 211         const int MAXIT = 100;
 212         const double EPS = numeric_limits<double>::epsilon();
 213         const double FPMIN = numeric_limits<double>::min() / EPS;
 214         int m1, m2;
 215         double aa, c, d, del, h, qab, qam, qap;
 216
 217         qab=a+b;
 218         qap=a+1.0;
 219         qam=a-1.0;
 220         c=1.0;
 221         d=1.0-qab*x/qap;
 222         if (fabs(d) < FPMIN) d=FPMIN;
 223         d=1.0/d;
 224         h=d;
 225         for (m1=1;m1<=MAXIT;m1++) {
 226             m2=2*m1;
 227             aa=m1*(b-m1)*x/((qam+m2)*(a+m2));
 228             d=1.0+aa*d;
 229             if (fabs(d) < FPMIN) d=FPMIN;
 230             c=1.0+aa/c;
 231             if (fabs(c) < FPMIN) c=FPMIN;
 232             d=1.0/d;
 233             h *= d*c;
 234             aa = -(a+m1)*(qab+m1)*x/((a+m2)*(qap+m2));
 235             d=1.0+aa*d;
 236             if (fabs(d) < FPMIN) d=FPMIN;
 237             c=1.0+aa/c;
 238             if (fabs(c) < FPMIN) c=FPMIN;
 239             d=1.0/d;
 240             del=d*c;
 241             h *= del;
 242             if (fabs(del-1.0) < EPS) break;
 243         }
 244
 245         if (m1 > MAXIT) { m->mothurOut("[ERROR]: a or b too big or MAXIT too small in betacf."); m->mothurOutEndLine(); m->control_pressed = true; }
 246         return h;
 247
 248     }
 249         catch(exception& e) {
 250                 m->errorOut(e, "LinearAlgebra", "betacf");
 251                 exit(1);
 252         }
 253 }
 254 /*********************************************************************************************************************************/
 255
 256 vector<vector<double> > LinearAlgebra::matrix_mult(vector<vector<double> > first, vector<vector<double> > second){
 257         try {
 258                 vector<vector<double> > product;
 259
 260                 int first_rows = first.size();
 261                 int first_cols = first[0].size();
 262                 int second_cols = second[0].size();
 263
 264                 product.resize(first_rows);
 265                 for(int i=0;i<first_rows;i++){
 266                         product[i].resize(second_cols);
 267                 }
 268
 269                 for(int i=0;i<first_rows;i++){
 270                         for(int j=0;j<second_cols;j++){
 271
 272                                 if (m->control_pressed) { return product; }
 273
 274                                 product[i][j] = 0.0;
 275                                 for(int k=0;k<first_cols;k++){
 276                                         product[i][j] += first[i][k] * second[k][j];
 277                                 }
 278                         }
 279                 }
 280
 281                 return product;
 282         }
 283         catch(exception& e) {
 284                 m->errorOut(e, "LinearAlgebra", "matrix_mult");
 285                 exit(1);
 286         }
 287
 288 }
 289
 290 /*********************************************************************************************************************************/
 291
 292 void LinearAlgebra::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
 293         try {
 294                 int rank = D.size();
 295
 296                 vector<vector<double> > A(rank);
 297                 vector<vector<double> > C(rank);
 298                 for(int i=0;i<rank;i++){
 299                         A[i].resize(rank);
 300                         C[i].resize(rank);
 301                 }
 302
 303                 double scale = -1.0000 / (double) rank;
 304
 305                 for(int i=0;i<rank;i++){
 306                         A[i][i] = 0.0000;
 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;
 311                         }
 312                 }
 313
 314                 A = matrix_mult(C,A);
 315                 G = matrix_mult(A,C);
 316         }
 317         catch(exception& e) {
 318                 m->errorOut(e, "LinearAlgebra", "recenter");
 319                 exit(1);
 320         }
 321
 322 }
 323 /*********************************************************************************************************************************/
 324
 325 //  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
 326
 327 int LinearAlgebra::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
 328         try {
 329                 double scale, hh, h, g, f;
 330
 331                 int n = a.size();
 332
 333                 d.resize(n);
 334                 e.resize(n);
 335
 336                 for(int i=n-1;i>0;i--){
 337                         int l=i-1;
 338                         h = scale = 0.0000;
 339                         if(l>0){
 340                                 for(int k=0;k<l+1;k++){
 341                                         scale += fabs(a[i][k]);
 342                                 }
 343                                 if(scale == 0.0){
 344                                         e[i] = a[i][l];
 345                                 }
 346                                 else{
 347                                         for(int k=0;k<l+1;k++){
 348                                                 a[i][k] /= scale;
 349                                                 h += a[i][k] * a[i][k];
 350                                         }
 351                                         f = a[i][l];
 352                                         g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
 353                                         e[i] = scale * g;
 354                                         h -= f * g;
 355                                         a[i][l] = f - g;
 356                                         f = 0.0;
 357                                         for(int j=0;j<l+1;j++){
 358                                                 a[j][i] = a[i][j] / h;
 359                                                 g = 0.0;
 360                                                 for(int k=0;k<j+1;k++){
 361                                                         g += a[j][k] * a[i][k];
 362                                                 }
 363                                                 for(int k=j+1;k<l+1;k++){
 364                                                         g += a[k][j] * a[i][k];
 365                                                 }
 366                                                 e[j] = g / h;
 367                                                 f += e[j] * a[i][j];
 368                                         }
 369                                         hh = f / (h + h);
 370                                         for(int j=0;j<l+1;j++){
 371                                                 f = a[i][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]);
 375                                                 }
 376                                         }
 377                                 }
 378                         }
 379                         else{
 380                                 e[i] = a[i][l];
 381                         }
 382
 383                         d[i] = h;
 384                 }
 385
 386                 d[0] = 0.0000;
 387                 e[0] = 0.0000;
 388
 389                 for(int i=0;i<n;i++){
 390                         int l = i;
 391                         if(d[i] != 0.0){
 392                                 for(int j=0;j<l;j++){
 393                                         g = 0.0000;
 394                                         for(int k=0;k<l;k++){
 395                                                 g += a[i][k] * a[k][j];
 396                                         }
 397                                         for(int k=0;k<l;k++){
 398                                                 a[k][j] -= g * a[k][i];
 399                                         }
 400                                 }
 401                         }
 402                         d[i] = a[i][i];
 403                         a[i][i] = 1.0000;
 404                         for(int j=0;j<l;j++){
 405                                 a[j][i] = a[i][j] = 0.0;
 406                         }
 407                 }
 408
 409                 return 0;
 410         }
 411         catch(exception& e) {
 412                 m->errorOut(e, "LinearAlgebra", "tred2");
 413                 exit(1);
 414         }
 415
 416 }
 417 /*********************************************************************************************************************************/
 418
 419 double LinearAlgebra::pythag(double a, double b)        {       return(pow(a*a+b*b,0.5));       }
 420
 421 /*********************************************************************************************************************************/
 422
 423 //  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
 424
 425 int LinearAlgebra::qtli(vector<double>& d, vector<double>& e, vector<vector<double> >& z) {
 426         try {
 427                 int myM, i, iter;
 428                 double s, r, p, g, f, dd, c, b;
 429
 430                 int n = d.size();
 431                 for(int i=1;i<=n;i++){
 432                         e[i-1] = e[i];
 433                 }
 434                 e[n-1] = 0.0000;
 435
 436                 for(int l=0;l<n;l++){
 437                         iter = 0;
 438                         do {
 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;
 442                                 }
 443                                 if(myM != l){
 444                                         if(iter++ == 3000) cerr << "Too many iterations in tqli\n";
 445                                         g = (d[l+1]-d[l]) / (2.0 * e[l]);
 446                                         r = pythag(g, 1.0);
 447                                         g = d[myM] - d[l] + e[l] / (g + SIGN(r,g));
 448                                         s = c = 1.0;
 449                                         p = 0.0000;
 450                                         for(i=myM-1;i>=l;i--){
 451                                                 f = s * e[i];
 452                                                 b = c * e[i];
 453                                                 e[i+1] = (r=pythag(f,g));
 454                                                 if(r==0.0){
 455                                                         d[i+1] -= p;
 456                                                         e[myM] = 0.0000;
 457                                                         break;
 458                                                 }
 459                                                 s = f / r;
 460                                                 c = g / r;
 461                                                 g = d[i+1] - p;
 462                                                 r = (d[i] - g) * s + 2.0 * c * b;
 463                                                 d[i+1] = g + ( p = s * r);
 464                                                 g = c * r - b;
 465                                                 for(int k=0;k<n;k++){
 466                                                         f = z[k][i+1];
 467                                                         z[k][i+1] = s * z[k][i] + c * f;
 468                                                         z[k][i] = c * z[k][i] - s * f;
 469                                                 }
 470                                         }
 471                                         if(r == 0.00 && i >= l) continue;
 472                                         d[l] -= p;
 473                                         e[l] = g;
 474                                         e[myM] = 0.0;
 475                                 }
 476                         } while (myM != l);
 477                 }
 478
 479                 int k;
 480                 for(int i=0;i<n;i++){
 481                         p=d[k=i];
 482                         for(int j=i;j<n;j++){
 483                                 if(d[j] >= p){
 484                                         p=d[k=j];
 485                                 }
 486                         }
 487                         if(k!=i){
 488                                 d[k]=d[i];
 489                                 d[i]=p;
 490                                 for(int j=0;j<n;j++){
 491                                         p=z[j][i];
 492                                         z[j][i] = z[j][k];
 493                                         z[j][k] = p;
 494                                 }
 495                         }
 496                 }
 497
 498                 return 0;
 499         }
 500         catch(exception& e) {
 501                 m->errorOut(e, "LinearAlgebra", "qtli");
 502                 exit(1);
 503         }
 504 }
 505 /*********************************************************************************************************************************/
 506 //groups by dimension
 507 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes, int dimensions){
 508         try {
 509                 //make square matrix
 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); }
 512
 513                 if (dimensions == 1) { //one dimension calc = abs(x-y)
 514
 515                         for (int i = 0; i < dists.size(); i++) {
 516
 517                                 if (m->control_pressed) { return dists; }
 518
 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];
 522                                 }
 523                         }
 524
 525                 }else if (dimensions > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
 526
 527                         for (int i = 0; i < dists.size(); i++) {
 528
 529                                 if (m->control_pressed) { return dists; }
 530
 531                                 for (int j = 0; j < i; j++) {
 532                                         double sum = 0.0;
 533                                         for (int k = 0; k < dimensions; k++) {
 534                                                 sum += ((axes[i][k] - axes[j][k]) * (axes[i][k] - axes[j][k]));
 535                                         }
 536
 537                                         dists[i][j] = sqrt(sum);
 538                                         dists[j][i] = dists[i][j];
 539                                 }
 540                         }
 541
 542                 }
 543
 544                 return dists;
 545         }
 546         catch(exception& e) {
 547                 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
 548                 exit(1);
 549         }
 550 }
 551 /*********************************************************************************************************************************/
 552 //returns groups by dimensions from dimensions by groups
 553 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes){
 554         try {
 555                 //make square matrix
 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); }
 558
 559                 if (axes.size() == 1) { //one dimension calc = abs(x-y)
 560
 561                         for (int i = 0; i < dists.size(); i++) {
 562
 563                                 if (m->control_pressed) { return dists; }
 564
 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];
 568                                 }
 569                         }
 570
 571                 }else if (axes.size() > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
 572
 573                         for (int i = 0; i < dists[0].size(); i++) {
 574
 575                                 if (m->control_pressed) { return dists; }
 576
 577                                 for (int j = 0; j < i; j++) {
 578                                         double sum = 0.0;
 579                                         for (int k = 0; k < axes.size(); k++) {
 580                                                 sum += ((axes[k][i] - axes[k][j]) * (axes[k][i] - axes[k][j]));
 581                                         }
 582
 583                                         dists[i][j] = sqrt(sum);
 584                                         dists[j][i] = dists[i][j];
 585                                 }
 586                         }
 587
 588                 }
 589
 590                 return dists;
 591         }
 592         catch(exception& e) {
 593                 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
 594                 exit(1);
 595         }
 596 }
 597 /*********************************************************************************************************************************/
 598 //assumes both matrices are square and the same size
 599 double LinearAlgebra::calcPearson(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
 600         try {
 601
 602                 //find average for - X
 603                 int count = 0;
 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];
 608                                 count++;
 609                         }
 610                 }
 611                 averageEuclid = averageEuclid / (float) count;
 612
 613                 //find average for - Y
 614                 count = 0;
 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];
 619                                 count++;
 620                         }
 621                 }
 622                 averageUser = averageUser / (float) count;
 623
 624                 double numerator = 0.0;
 625                 double denomTerm1 = 0.0;
 626                 double denomTerm2 = 0.0;
 627
 628                 for (int i = 0; i < euclidDists.size(); i++) {
 629
 630                         for (int k = 0; k < i; k++) { //just lt dists
 631
 632                                 float Yi = userDists[i][k];
 633                                 float Xi = euclidDists[i][k];
 634
 635                                 numerator += ((Xi - averageEuclid) * (Yi - averageUser));
 636                                 denomTerm1 += ((Xi - averageEuclid) * (Xi - averageEuclid));
 637                                 denomTerm2 += ((Yi - averageUser) * (Yi - averageUser));
 638                         }
 639                 }
 640
 641                 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
 642                 double r = numerator / denom;
 643
 644                 //divide by zero error
 645                 if (isnan(r) || isinf(r)) { r = 0.0; }
 646
 647                 return r;
 648
 649         }
 650         catch(exception& e) {
 651                 m->errorOut(e, "LinearAlgebra", "calcPearson");
 652                 exit(1);
 653         }
 654 }
 655 /*********************************************************************************************************************************/
 656 //assumes both matrices are square and the same size
 657 double LinearAlgebra::calcSpearman(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
 658         try {
 659                 double r;
 660
 661                 //format data
 662                 map<float, int> tableX;
 663                 map<float, int>::iterator itTable;
 664                 vector<spearmanRank> scores;
 665
 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);
 670
 671                                 //count number of repeats
 672                                 itTable = tableX.find(euclidDists[i][j]);
 673                                 if (itTable == tableX.end()) {
 674                                         tableX[euclidDists[i][j]] = 1;
 675                                 }else {
 676                                         tableX[euclidDists[i][j]]++;
 677                                 }
 678                         }
 679                 }
 680
 681                 //sort scores
 682                 sort(scores.begin(), scores.end(), compareSpearman);
 683
 684                 //calc LX
 685                 double Lx = 0.0;
 686                 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
 687                         double tx = (double) itTable->second;
 688                         Lx += ((pow(tx, 3.0) - tx) / 12.0);
 689                 }
 690
 691                 //find ranks of xi
 692                 map<string, float> rankEuclid;
 693                 vector<spearmanRank> ties;
 694                 int rankTotal = 0;
 695                 for (int j = 0; j < scores.size(); j++) {
 696                         rankTotal += (j+1);
 697                         ties.push_back(scores[j]);
 698
 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
 701
 702                                         for (int k = 0; k < ties.size(); k++) {
 703                                                 float thisrank = rankTotal / (float) ties.size();
 704                                                 rankEuclid[ties[k].name] = thisrank;
 705                                         }
 706                                         ties.clear();
 707                                         rankTotal = 0;
 708                                 }
 709                         }else { // you are the last one
 710
 711                                 for (int k = 0; k < ties.size(); k++) {
 712                                         float thisrank = rankTotal / (float) ties.size();
 713                                         rankEuclid[ties[k].name] = thisrank;
 714                                 }
 715                         }
 716                 }
 717
 718
 719                 //format data
 720                 map<float, int> tableY;
 721                 scores.clear();
 722
 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);
 727
 728                                 //count number of repeats
 729                                 itTable = tableY.find(userDists[i][j]);
 730                                 if (itTable == tableY.end()) {
 731                                         tableY[userDists[i][j]] = 1;
 732                                 }else {
 733                                         tableY[userDists[i][j]]++;
 734                                 }
 735                         }
 736                 }
 737
 738                 //sort scores
 739                 sort(scores.begin(), scores.end(), compareSpearman);
 740
 741                 //calc LX
 742                 double Ly = 0.0;
 743                 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
 744                         double ty = (double) itTable->second;
 745                         Ly += ((pow(ty, 3.0) - ty) / 12.0);
 746                 }
 747
 748                 //find ranks of yi
 749                 map<string, float> rankUser;
 750                 ties.clear();
 751                 rankTotal = 0;
 752                 for (int j = 0; j < scores.size(); j++) {
 753                         rankTotal += (j+1);
 754                         ties.push_back(scores[j]);
 755
 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
 758
 759                                         for (int k = 0; k < ties.size(); k++) {
 760                                                 float thisrank = rankTotal / (float) ties.size();
 761                                                 rankUser[ties[k].name] = thisrank;
 762                                         }
 763                                         ties.clear();
 764                                         rankTotal = 0;
 765                                 }
 766                         }else { // you are the last one
 767
 768                                 for (int k = 0; k < ties.size(); k++) {
 769                                         float thisrank = rankTotal / (float) ties.size();
 770                                         rankUser[ties[k].name] = thisrank;
 771                                 }
 772                         }
 773                 }
 774
 775
 776                 double di = 0.0;
 777                 int count = 0;
 778                 for (int i = 0; i < userDists.size(); i++) {
 779                         for (int j = 0; j < i; j++) {
 780
 781                                 float xi = rankEuclid[toString(count)];
 782                                 float yi = rankUser[toString(count)];
 783
 784                                 di += ((xi - yi) * (xi - yi));
 785
 786                                 count++;
 787                         }
 788                 }
 789
 790                 double n = (double) count;
 791
 792                 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
 793                 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
 794
 795                 r = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
 796
 797                 //divide by zero error
 798                 if (isnan(r) || isinf(r)) { r = 0.0; }
 799
 800                 return r;
 801
 802         }
 803         catch(exception& e) {
 804                 m->errorOut(e, "LinearAlgebra", "calcSpearman");
 805                 exit(1);
 806         }
 807 }
 808
 809 /*********************************************************************************************************************************/
 810 //assumes both matrices are square and the same size
 811 double LinearAlgebra::calcKendall(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
 812         try {
 813                 double r;
 814
 815                 //format data
 816                 vector<spearmanRank> scores;
 817                 for (int i = 0; i < euclidDists.size(); i++) {
 818                         for (int j = 0; j < i; j++) {
 819                                 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
 820                                 scores.push_back(member);
 821                         }
 822                 }
 823
 824                 //sort scores
 825                 sort(scores.begin(), scores.end(), compareSpearman);
 826
 827                 //find ranks of xi
 828                 map<string, float> rankEuclid;
 829                 vector<spearmanRank> ties;
 830                 int rankTotal = 0;
 831                 for (int j = 0; j < scores.size(); j++) {
 832                         rankTotal += (j+1);
 833                         ties.push_back(scores[j]);
 834
 835                         if (j != (scores.size()-1)) { // you are not the last so you can look ahead
 836                                 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
 837
 838                                         for (int k = 0; k < ties.size(); k++) {
 839                                                 float thisrank = rankTotal / (float) ties.size();
 840                                                 rankEuclid[ties[k].name] = thisrank;
 841                                         }
 842                                         ties.clear();
 843                                         rankTotal = 0;
 844                                 }
 845                         }else { // you are the last one
 846
 847                                 for (int k = 0; k < ties.size(); k++) {
 848                                         float thisrank = rankTotal / (float) ties.size();
 849                                         rankEuclid[ties[k].name] = thisrank;
 850                                 }
 851                         }
 852                 }
 853
 854                 vector<spearmanRank> scoresUser;
 855                 for (int i = 0; i < userDists.size(); i++) {
 856                         for (int j = 0; j < i; j++) {
 857                                 spearmanRank member(toString(scoresUser.size()), userDists[i][j]);
 858                                 scoresUser.push_back(member);
 859                         }
 860                 }
 861
 862                 //sort scores
 863                 sort(scoresUser.begin(), scoresUser.end(), compareSpearman);
 864
 865                 //find ranks of yi
 866                 map<string, float> rankUser;
 867                 ties.clear();
 868                 rankTotal = 0;
 869                 for (int j = 0; j < scoresUser.size(); j++) {
 870                         rankTotal += (j+1);
 871                         ties.push_back(scoresUser[j]);
 872
 873                         if (j != (scoresUser.size()-1)) { // you are not the last so you can look ahead
 874                                 if (scoresUser[j].score != scoresUser[j+1].score) { // you are done with ties, rank them and continue
 875
 876                                         for (int k = 0; k < ties.size(); k++) {
 877                                                 float thisrank = rankTotal / (float) ties.size();
 878                                                 rankUser[ties[k].name] = thisrank;
 879                                         }
 880                                         ties.clear();
 881                                         rankTotal = 0;
 882                                 }
 883                         }else { // you are the last one
 884
 885                                 for (int k = 0; k < ties.size(); k++) {
 886                                         float thisrank = rankTotal / (float) ties.size();
 887                                         rankUser[ties[k].name] = thisrank;
 888                                 }
 889                         }
 890                 }
 891
 892                 int numCoor = 0;
 893                 int numDisCoor = 0;
 894
 895                 //order user ranks
 896                 vector<spearmanRank> user;
 897                 for (int l = 0; l < scores.size(); l++) {
 898                         spearmanRank member(scores[l].name, rankUser[scores[l].name]);
 899                         user.push_back(member);
 900                 }
 901
 902                 int count = 0;
 903                 for (int l = 0; l < scores.size(); l++) {
 904
 905                         int numWithHigherRank = 0;
 906                         int numWithLowerRank = 0;
 907                         float thisrank = user[l].score;
 908
 909                         for (int u = l+1; u < scores.size(); u++) {
 910                                 if (user[u].score > thisrank) { numWithHigherRank++; }
 911                                 else if (user[u].score < thisrank) { numWithLowerRank++; }
 912                                 count++;
 913                         }
 914
 915                         numCoor += numWithHigherRank;
 916                         numDisCoor += numWithLowerRank;
 917                 }
 918
 919                 r = (numCoor - numDisCoor) / (float) count;
 920
 921                 //divide by zero error
 922                 if (isnan(r) || isinf(r)) { r = 0.0; }
 923
 924                 return r;
 925
 926         }
 927         catch(exception& e) {
 928                 m->errorOut(e, "LinearAlgebra", "calcKendall");
 929                 exit(1);
 930         }
 931 }
 932 /*********************************************************************************************************************************/
 933 double LinearAlgebra::calcKendall(vector<double>& x, vector<double>& y, double& sig){
 934         try {
 935                 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
 936
 937                 //format data
 938                 vector<spearmanRank> xscores;
 939                 for (int i = 0; i < x.size(); i++) {
 940                         spearmanRank member(toString(i), x[i]);
 941                         xscores.push_back(member);
 942                 }
 943
 944                 //sort xscores
 945                 sort(xscores.begin(), xscores.end(), compareSpearman);
 946
 947                 //convert scores to ranks of x
 948                 vector<spearmanRank*> ties;
 949                 int rankTotal = 0;
 950                 for (int j = 0; j < xscores.size(); j++) {
 951                         rankTotal += (j+1);
 952                         ties.push_back(&(xscores[j]));
 953
 954                         if (j != xscores.size()-1) { // you are not the last so you can look ahead
 955                                 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
 956                                         for (int k = 0; k < ties.size(); k++) {
 957                                                 float thisrank = rankTotal / (float) ties.size();
 958                                                 (*ties[k]).score = thisrank;
 959                                         }
 960                                         ties.clear();
 961                                         rankTotal = 0;
 962                                 }
 963                         }else { // you are the last one
 964                                 for (int k = 0; k < ties.size(); k++) {
 965                                         float thisrank = rankTotal / (float) ties.size();
 966                                         (*ties[k]).score = thisrank;
 967                                 }
 968                         }
 969                 }
 970
 971
 972                 //format data
 973                 vector<spearmanRank> yscores;
 974                 for (int j = 0; j < y.size(); j++) {
 975                         spearmanRank member(toString(j), y[j]);
 976                         yscores.push_back(member);
 977                 }
 978
 979                 //sort yscores
 980                 sort(yscores.begin(), yscores.end(), compareSpearman);
 981
 982                 //convert to ranks
 983                 map<string, float> rank;
 984                 vector<spearmanRank> yties;
 985                 rankTotal = 0;
 986                 for (int j = 0; j < yscores.size(); j++) {
 987                         rankTotal += (j+1);
 988                         yties.push_back(yscores[j]);
 989
 990                         if (j != yscores.size()-1) { // you are not the last so you can look ahead
 991                                 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
 992                                         for (int k = 0; k < yties.size(); k++) {
 993                                                 float thisrank = rankTotal / (float) yties.size();
 994                                                 rank[yties[k].name] = thisrank;
 995                                         }
 996                                         yties.clear();
 997                                         rankTotal = 0;
 998                                 }
 999                         }else { // you are the last one
1000                                 for (int k = 0; k < yties.size(); k++) {
1001                                         float thisrank = rankTotal / (float) yties.size();
1002                                         rank[yties[k].name] = thisrank;
1003                                 }
1004                         }
1005                 }
1006
1007
1008                 int numCoor = 0;
1009                 int numDisCoor = 0;
1010
1011                 //associate x and y
1012                 vector<spearmanRank> otus;
1013                 for (int l = 0; l < xscores.size(); l++) {
1014                         spearmanRank member(xscores[l].name, rank[xscores[l].name]);
1015                         otus.push_back(member);
1016                 }
1017
1018                 int count = 0;
1019                 for (int l = 0; l < xscores.size(); l++) {
1020
1021                         int numWithHigherRank = 0;
1022                         int numWithLowerRank = 0;
1023                         float thisrank = otus[l].score;
1024
1025                         for (int u = l+1; u < xscores.size(); u++) {
1026                                 if (otus[u].score > thisrank) { numWithHigherRank++; }
1027                                 else if (otus[u].score < thisrank) { numWithLowerRank++; }
1028                                 count++;
1029                         }
1030
1031                         numCoor += numWithHigherRank;
1032                         numDisCoor += numWithLowerRank;
1033                 }
1034
1035                 double p = (numCoor - numDisCoor) / (float) count;
1036
1037                 sig = calcKendallSig(x.size(), p);
1038
1039                 return p;
1040         }
1041         catch(exception& e) {
1042                 m->errorOut(e, "LinearAlgebra", "calcKendall");
1043                 exit(1);
1044         }
1045 }
1046 double LinearAlgebra::ran0(int& idum)
1047 {
1048     const int IA=16807,IM=2147483647,IQ=127773;
1049     const int IR=2836,MASK=123459876;
1050     const double AM=1.0/double(IM);
1051     int k;
1052     double ans;
1053
1054     idum ^= MASK;
1055     k=idum/IQ;
1056     idum=IA*(idum-k*IQ)-IR*k;
1057     if (idum < 0) idum += IM;
1058     ans=AM*idum;
1059     idum ^= MASK;
1060     return ans;
1061 }
1062
1063 double LinearAlgebra::ran1(int &idum)
1064 {
1065         const int IA=16807,IM=2147483647,IQ=127773,IR=2836,NTAB=32;
1066         const int NDIV=(1+(IM-1)/NTAB);
1067         const double EPS=3.0e-16,AM=1.0/IM,RNMX=(1.0-EPS);
1068         static int iy=0;
1069         static vector<int> iv(NTAB);
1070         int j,k;
1071         double temp;
1072
1073         if (idum <= 0 || !iy) {
1074                 if (-idum < 1) idum=1;
1075                 else idum = -idum;
1076                 for (j=NTAB+7;j>=0;j--) {
1077                         k=idum/IQ;
1078                         idum=IA*(idum-k*IQ)-IR*k;
1079                         if (idum < 0) idum += IM;
1080                         if (j < NTAB) iv[j] = idum;
1081                 }
1082                 iy=iv[0];
1083         }
1084         k=idum/IQ;
1085         idum=IA*(idum-k*IQ)-IR*k;
1086         if (idum < 0) idum += IM;
1087         j=iy/NDIV;
1088         iy=iv[j];
1089         iv[j] = idum;
1090         if ((temp=AM*iy) > RNMX) return RNMX;
1091         else return temp;
1092 }
1093
1094 double LinearAlgebra::ran2(int &idum)
1095 {
1096         const int IM1=2147483563,IM2=2147483399;
1097         const int IA1=40014,IA2=40692,IQ1=53668,IQ2=52774;
1098         const int IR1=12211,IR2=3791,NTAB=32,IMM1=IM1-1;
1099         const int NDIV=1+IMM1/NTAB;
1100         const double EPS=3.0e-16,RNMX=1.0-EPS,AM=1.0/double(IM1);
1101         static int idum2=123456789,iy=0;
1102         static vector<int> iv(NTAB);
1103         int j,k;
1104         double temp;
1105
1106         if (idum <= 0) {
1107                 idum=(idum==0 ? 1 : -idum);
1108                 idum2=idum;
1109                 for (j=NTAB+7;j>=0;j--) {
1110                         k=idum/IQ1;
1111                         idum=IA1*(idum-k*IQ1)-k*IR1;
1112                         if (idum < 0) idum += IM1;
1113                         if (j < NTAB) iv[j] = idum;
1114                 }
1115                 iy=iv[0];
1116         }
1117         k=idum/IQ1;
1118         idum=IA1*(idum-k*IQ1)-k*IR1;
1119         if (idum < 0) idum += IM1;
1120         k=idum2/IQ2;
1121         idum2=IA2*(idum2-k*IQ2)-k*IR2;
1122         if (idum2 < 0) idum2 += IM2;
1123         j=iy/NDIV;
1124         iy=iv[j]-idum2;
1125         iv[j] = idum;
1126         if (iy < 1) iy += IMM1;
1127         if ((temp=AM*iy) > RNMX) return RNMX;
1128         else return temp;
1129 }
1130
1131 double LinearAlgebra::ran3(int &idum)
1132 {
1133         static int inext,inextp;
1134         static int iff=0;
1135         const int MBIG=1000000000,MSEED=161803398,MZ=0;
1136         const double FAC=(1.0/MBIG);
1137         static vector<int> ma(56);
1138         int i,ii,k,mj,mk;
1139
1140         if (idum < 0 || iff == 0) {
1141                 iff=1;
1142                 mj=labs(MSEED-labs(idum));
1143                 mj %= MBIG;
1144                 ma[55]=mj;
1145                 mk=1;
1146                 for (i=1;i<=54;i++) {
1147                         ii=(21*i) % 55;
1148                         ma[ii]=mk;
1149                         mk=mj-mk;
1150                         if (mk < int(MZ)) mk += MBIG;
1151                         mj=ma[ii];
1152                 }
1153                 for (k=0;k<4;k++)
1154                         for (i=1;i<=55;i++) {
1155                                 ma[i] -= ma[1+(i+30) % 55];
1156                                 if (ma[i] < int(MZ)) ma[i] += MBIG;
1157                         }
1158                 inext=0;
1159                 inextp=31;
1160                 idum=1;
1161         }
1162         if (++inext == 56) inext=1;
1163         if (++inextp == 56) inextp=1;
1164         mj=ma[inext]-ma[inextp];
1165         if (mj < int(MZ)) mj += MBIG;
1166         ma[inext]=mj;
1167         return mj*FAC;
1168 }
1169
1170 double LinearAlgebra::ran4(int &idum)
1171 {
1172 #if defined(vax) || defined(_vax_) || defined(__vax__) || defined(VAX)
1173         static const unsigned long jflone = 0x00004080;
1174         static const unsigned long jflmsk = 0xffff007f;
1175 #else
1176         static const unsigned long jflone = 0x3f800000;
1177         static const unsigned long jflmsk = 0x007fffff;
1178 #endif
1179         unsigned long irword,itemp,lword;
1180         static int idums = 0;
1181
1182         if (idum < 0) {
1183                 idums = -idum;
1184                 idum=1;
1185         }
1186         irword=idum;
1187         lword=idums;
1188         psdes(lword,irword);
1189         itemp=jflone | (jflmsk & irword);
1190         ++idum;
1191         return (*(float *)&itemp)-1.0;
1192 }
1193
1194 void LinearAlgebra::psdes(unsigned long &lword, unsigned long &irword)
1195 {
1196         const int NITER=4;
1197         static const unsigned long c1[NITER]={
1198                 0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L};
1199         static const unsigned long c2[NITER]={
1200                 0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L};
1201         unsigned long i,ia,ib,iswap,itmph=0,itmpl=0;
1202
1203         for (i=0;i<NITER;i++) {
1204                 ia=(iswap=irword) ^ c1[i];
1205                 itmpl = ia & 0xffff;
1206                 itmph = ia >> 16;
1207                 ib=itmpl*itmpl+ ~(itmph*itmph);
1208                 irword=lword ^ (((ia = (ib >> 16) |
1209                           ((ib & 0xffff) << 16)) ^ c2[i])+itmpl*itmph);
1210                 lword=iswap;
1211         }
1212 }
1213 /*********************************************************************************************************************************/
1214 double LinearAlgebra::calcKendallSig(double n, double r){
1215     try {
1216
1217         double sig = 0.0;
1218         double svar=(4.0*n+10.0)/(9.0*n*(n-1.0));
1219         double z= r/sqrt(svar);
1220         sig=erfcc(fabs(z)/1.4142136);
1221
1222                 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1223
1224         return sig;
1225     }
1226         catch(exception& e) {
1227                 m->errorOut(e, "LinearAlgebra", "calcKendallSig");
1228                 exit(1);
1229         }
1230 }
1231
1232 /*********************************************************************************************************************************/
1233 double LinearAlgebra::calcSpearman(vector<double>& x, vector<double>& y, double& sig){
1234         try {
1235                 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1236
1237                 //format data
1238         double sf = 0.0; //f^3 - f where f is the number of ties in x;
1239         double sg = 0.0; //f^3 - f where f is the number of ties in y;
1240                 map<float, int> tableX;
1241                 map<float, int>::iterator itTable;
1242                 vector<spearmanRank> xscores;
1243
1244                 for (int i = 0; i < x.size(); i++) {
1245                         spearmanRank member(toString(i), x[i]);
1246                         xscores.push_back(member);
1247
1248                         //count number of repeats
1249                         itTable = tableX.find(x[i]);
1250                         if (itTable == tableX.end()) {
1251                                 tableX[x[i]] = 1;
1252                         }else {
1253                                 tableX[x[i]]++;
1254                         }
1255                 }
1256
1257
1258                 //calc LX
1259                 double Lx = 0.0;
1260                 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
1261                         double tx = (double) itTable->second;
1262                         Lx += ((pow(tx, 3.0) - tx) / 12.0);
1263                 }
1264
1265
1266                 //sort x
1267                 sort(xscores.begin(), xscores.end(), compareSpearman);
1268
1269                 //convert scores to ranks of x
1270                 //convert to ranks
1271                 map<string, float> rankx;
1272                 vector<spearmanRank> xties;
1273                 int rankTotal = 0;
1274                 for (int j = 0; j < xscores.size(); j++) {
1275                         rankTotal += (j+1);
1276                         xties.push_back(xscores[j]);
1277
1278                         if (j != xscores.size()-1) { // you are not the last so you can look ahead
1279                                 if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
1280                                         for (int k = 0; k < xties.size(); k++) {
1281                                                 float thisrank = rankTotal / (float) xties.size();
1282                                                 rankx[xties[k].name] = thisrank;
1283                                         }
1284                     int t = xties.size();
1285                     sf += (t*t*t-t);
1286                                         xties.clear();
1287                                         rankTotal = 0;
1288                                 }
1289                         }else { // you are the last one
1290                                 for (int k = 0; k < xties.size(); k++) {
1291                                         float thisrank = rankTotal / (float) xties.size();
1292                                         rankx[xties[k].name] = thisrank;
1293                                 }
1294                         }
1295                 }
1296
1297                 //format x
1298                 vector<spearmanRank> yscores;
1299                 map<float, int> tableY;
1300                 for (int j = 0; j < y.size(); j++) {
1301                         spearmanRank member(toString(j), y[j]);
1302                         yscores.push_back(member);
1303
1304                         itTable = tableY.find(member.score);
1305                         if (itTable == tableY.end()) {
1306                                 tableY[member.score] = 1;
1307                         }else {
1308                                 tableY[member.score]++;
1309                         }
1310
1311                 }
1312
1313                 //calc Ly
1314                 double Ly = 0.0;
1315                 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
1316                         double ty = (double) itTable->second;
1317                         Ly += ((pow(ty, 3.0) - ty) / 12.0);
1318                 }
1319
1320                 sort(yscores.begin(), yscores.end(), compareSpearman);
1321
1322                 //convert to ranks
1323                 map<string, float> rank;
1324                 vector<spearmanRank> yties;
1325                 rankTotal = 0;
1326                 for (int j = 0; j < yscores.size(); j++) {
1327                         rankTotal += (j+1);
1328                         yties.push_back(yscores[j]);
1329
1330                         if (j != yscores.size()-1) { // you are not the last so you can look ahead
1331                                 if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
1332                                         for (int k = 0; k < yties.size(); k++) {
1333                                                 float thisrank = rankTotal / (float) yties.size();
1334                                                 rank[yties[k].name] = thisrank;
1335                                         }
1336                     int t = yties.size();
1337                     sg += (t*t*t-t);
1338                                         yties.clear();
1339                                         rankTotal = 0;
1340                                 }
1341                         }else { // you are the last one
1342                                 for (int k = 0; k < yties.size(); k++) {
1343                                         float thisrank = rankTotal / (float) yties.size();
1344                                         rank[yties[k].name] = thisrank;
1345                                 }
1346                         }
1347                 }
1348
1349                 double di = 0.0;
1350                 for (int k = 0; k < x.size(); k++) {
1351
1352                         float xi = rankx[toString(k)];
1353                         float yi = rank[toString(k)];
1354
1355                         di += ((xi - yi) * (xi - yi));
1356                 }
1357
1358                 double p = 0.0;
1359
1360                 double n = (double) x.size();
1361                 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
1362                 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
1363
1364                 p = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
1365
1366                 //Numerical Recipes 646
1367         sig = calcSpearmanSig(n, sf, sg, di);
1368
1369                 return p;
1370         }
1371         catch(exception& e) {
1372                 m->errorOut(e, "LinearAlgebra", "calcSpearman");
1373                 exit(1);
1374         }
1375 }
1376 /*********************************************************************************************************************************/
1377 double LinearAlgebra::calcSpearmanSig(double n, double sf, double sg, double d){
1378     try {
1379
1380         double sig = 0.0;
1381         double probrs = 0.0;
1382         double en=n;
1383         double en3n=en*en*en-en;
1384         double aved=en3n/6.0-(sf+sg)/12.0;
1385         double fac=(1.0-sf/en3n)*(1.0-sg/en3n);
1386         double vard=((en-1.0)*en*en*SQR(en+1.0)/36.0)*fac;
1387         double zd=(d-aved)/sqrt(vard);
1388         double probd=erfcc(fabs(zd)/1.4142136);
1389         double rs=(1.0-(6.0/en3n)*(d+(sf+sg)/12.0))/sqrt(fac);
1390         fac=(rs+1.0)*(1.0-rs);
1391         if (fac > 0.0) {
1392             double t=rs*sqrt((en-2.0)/fac);
1393             double df=en-2.0;
1394             probrs=betai(0.5*df,0.5,df/(df+t*t));
1395         }else {
1396             probrs = 0.0;
1397         }
1398
1399         //smaller of probd and probrs is sig
1400         sig = probrs;
1401         if (probd < probrs) { sig = probd; }
1402
1403                 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1404
1405         return sig;
1406     }
1407         catch(exception& e) {
1408                 m->errorOut(e, "LinearAlgebra", "calcSpearmanSig");
1409                 exit(1);
1410         }
1411 }
1412 /*********************************************************************************************************************************/
1413 double LinearAlgebra::calcPearson(vector<double>& x, vector<double>& y, double& sig){
1414         try {
1415                 if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
1416
1417                 //find average X
1418                 float averageX = 0.0;
1419                 for (int i = 0; i < x.size(); i++) { averageX += x[i];  }
1420                 averageX = averageX / (float) x.size();
1421
1422                 //find average Y
1423                 float sumY = 0.0;
1424                 for (int j = 0; j < y.size(); j++) { sumY += y[j]; }
1425                 float Ybar = sumY / (float) y.size();
1426
1427                 double r = 0.0;
1428                 double numerator = 0.0;
1429                 double denomTerm1 = 0.0;
1430                 double denomTerm2 = 0.0;
1431
1432                 for (int j = 0; j < x.size(); j++) {
1433                         float Yi = y[j];
1434                         float Xi = x[j];
1435
1436                         numerator += ((Xi - averageX) * (Yi - Ybar));
1437                         denomTerm1 += ((Xi - averageX) * (Xi - averageX));
1438                         denomTerm2 += ((Yi - Ybar) * (Yi - Ybar));
1439                 }
1440
1441                 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
1442
1443                 r = numerator / denom;
1444
1445                 //Numerical Recipes pg.644
1446         sig = calcPearsonSig(x.size(), r);
1447
1448                 return r;
1449         }
1450         catch(exception& e) {
1451                 m->errorOut(e, "LinearAlgebra", "calcPearson");
1452                 exit(1);
1453         }
1454 }
1455 /*********************************************************************************************************************************/
1456 double LinearAlgebra::calcPearsonSig(double n, double r){
1457     try {
1458
1459         double sig = 0.0;
1460         const double TINY = 1.0e-20;
1461         double z = 0.5*log((1.0+r+TINY)/(1.0-r+TINY)); //Fisher's z transformation
1462
1463         //code below was giving an error in betacf with sop files
1464         //int df = n-2;
1465         //double t = r*sqrt(df/((1.0-r+TINY)*(1.0+r+TINY)));
1466         //sig = betai(0.5+df, 0.5, df/(df+t*t));
1467
1468         //Numerical Recipes says code below gives approximately the same result
1469         sig = erfcc(fabs(z*sqrt(n-1.0))/1.4142136);
1470                 if (isnan(sig) || isinf(sig)) { sig = 0.0; }
1471
1472         return sig;
1473     }
1474         catch(exception& e) {
1475                 m->errorOut(e, "LinearAlgebra", "calcPearsonSig");
1476                 exit(1);
1477         }
1478 }
1479 /*********************************************************************************************************************************/
1480
1481 vector<vector<double> > LinearAlgebra::getObservedEuclideanDistance(vector<vector<double> >& relAbundData){
1482     try {
1483
1484         int numSamples = relAbundData.size();
1485         int numOTUs = relAbundData[0].size();
1486
1487         vector<vector<double> > dMatrix(numSamples);
1488         for(int i=0;i<numSamples;i++){
1489             dMatrix[i].resize(numSamples);
1490         }
1491
1492         for(int i=0;i<numSamples;i++){
1493             for(int j=0;j<numSamples;j++){
1494
1495                 if (m->control_pressed) { return dMatrix; }
1496
1497                 double d = 0;
1498                 for(int k=0;k<numOTUs;k++){
1499                     d += pow((relAbundData[i][k] - relAbundData[j][k]), 2.0000);
1500                 }
1501                 dMatrix[i][j] = pow(d, 0.50000);
1502                 dMatrix[j][i] = dMatrix[i][j];
1503
1504             }
1505         }
1506         return dMatrix;
1507         }
1508         catch(exception& e) {
1509                 m->errorOut(e, "LinearAlgebra", "getObservedEuclideanDistance");
1510                 exit(1);
1511         }
1512 }
1513
1514 /*********************************************************************************************************************************/
1515 vector<double> LinearAlgebra::solveEquations(vector<vector<double> > A, vector<double> b){
1516     try {
1517         int length = (int)b.size();
1518         vector<double> x(length, 0);
1519         vector<int> index(length);
1520         for(int i=0;i<length;i++){  index[i] = i;   }
1521         double d;
1522
1523         ludcmp(A, index, d);  if (m->control_pressed) { return b; }
1524         lubksb(A, index, b);
1525
1526         return b;
1527     }
1528         catch(exception& e) {
1529                 m->errorOut(e, "LinearAlgebra", "solveEquations");
1530                 exit(1);
1531         }
1532 }
1533 /*********************************************************************************************************************************/
1534 vector<float> LinearAlgebra::solveEquations(vector<vector<float> > A, vector<float> b){
1535     try {
1536         int length = (int)b.size();
1537         vector<double> x(length, 0);
1538         vector<int> index(length);
1539         for(int i=0;i<length;i++){  index[i] = i;   }
1540         float d;
1541
1542         ludcmp(A, index, d);  if (m->control_pressed) { return b; }
1543         lubksb(A, index, b);
1544
1545         return b;
1546     }
1547         catch(exception& e) {
1548                 m->errorOut(e, "LinearAlgebra", "solveEquations");
1549                 exit(1);
1550         }
1551 }
1552
1553 /*********************************************************************************************************************************/
1554
1555 void LinearAlgebra::ludcmp(vector<vector<double> >& A, vector<int>& index, double& d){
1556     try {
1557         double tiny = 1e-20;
1558
1559         int n = (int)A.size();
1560         vector<double> vv(n, 0.0);
1561         double temp;
1562         int imax;
1563
1564         d = 1.0;
1565
1566         for(int i=0;i<n;i++){
1567             double big = 0.0;
1568             for(int j=0;j<n;j++){   if((temp=fabs(A[i][j])) > big ) big=temp;  }
1569             if(big==0.0){   m->mothurOut("Singular matrix in routine ludcmp\n");    }
1570             vv[i] = 1.0/big;
1571         }
1572
1573         for(int j=0;j<n;j++){
1574             if (m->control_pressed) { break; }
1575             for(int i=0;i<j;i++){
1576                 double sum = A[i][j];
1577                 for(int k=0;k<i;k++){   sum -= A[i][k] * A[k][j];   }
1578                 A[i][j] = sum;
1579             }
1580
1581             double big = 0.0;
1582             for(int i=j;i<n;i++){
1583                 double sum = A[i][j];
1584                 for(int k=0;k<j;k++){   sum -= A[i][k] * A[k][j];   }
1585                 A[i][j] = sum;
1586                 double dum;
1587                 if((dum = vv[i] * fabs(sum)) >= big){
1588                     big = dum;
1589                     imax = i;
1590                 }
1591             }
1592             if(j != imax){
1593                 for(int k=0;k<n;k++){
1594                     double dum = A[imax][k];
1595                     A[imax][k] = A[j][k];
1596                     A[j][k] = dum;
1597                 }
1598                 d = -d;
1599                 vv[imax] = vv[j];
1600             }
1601             index[j] = imax;
1602
1603             if(A[j][j] == 0.0){ A[j][j] = tiny; }
1604
1605             if(j != n-1){
1606                 double dum = 1.0/A[j][j];
1607                 for(int i=j+1;i<n;i++){ A[i][j] *= dum; }
1608             }
1609         }
1610     }
1611         catch(exception& e) {
1612                 m->errorOut(e, "LinearAlgebra", "ludcmp");
1613                 exit(1);
1614         }
1615 }
1616
1617 /*********************************************************************************************************************************/
1618
1619 void LinearAlgebra::lubksb(vector<vector<double> >& A, vector<int>& index, vector<double>& b){
1620     try {
1621         double total;
1622         int n = (int)A.size();
1623         int ii = 0;
1624
1625         for(int i=0;i<n;i++){
1626             if (m->control_pressed) { break; }
1627             int ip = index[i];
1628             total = b[ip];
1629             b[ip] = b[i];
1630
1631             if (ii != 0) {
1632                 for(int j=ii-1;j<i;j++){
1633                     total -= A[i][j] * b[j];
1634                 }
1635             }
1636             else if(total != 0){  ii = i+1;   }
1637             b[i] = total;
1638         }
1639         for(int i=n-1;i>=0;i--){
1640             total = b[i];
1641             for(int j=i+1;j<n;j++){ total -= A[i][j] * b[j];  }
1642             b[i] = total / A[i][i];
1643         }
1644     }
1645         catch(exception& e) {
1646                 m->errorOut(e, "LinearAlgebra", "lubksb");
1647                 exit(1);
1648         }
1649 }
1650 /*********************************************************************************************************************************/
1651
1652 void LinearAlgebra::ludcmp(vector<vector<float> >& A, vector<int>& index, float& d){
1653     try {
1654         double tiny = 1e-20;
1655
1656         int n = (int)A.size();
1657         vector<float> vv(n, 0.0);
1658         double temp;
1659         int imax;
1660
1661         d = 1.0;
1662
1663         for(int i=0;i<n;i++){
1664             float big = 0.0;
1665             for(int j=0;j<n;j++){   if((temp=fabs(A[i][j])) > big ) big=temp;  }
1666             if(big==0.0){   m->mothurOut("Singular matrix in routine ludcmp\n");    }
1667             vv[i] = 1.0/big;
1668         }
1669
1670         for(int j=0;j<n;j++){
1671             if (m->control_pressed) { break; }
1672             for(int i=0;i<j;i++){
1673                 float sum = A[i][j];
1674                 for(int k=0;k<i;k++){   sum -= A[i][k] * A[k][j];   }
1675                 A[i][j] = sum;
1676             }
1677
1678             float big = 0.0;
1679             for(int i=j;i<n;i++){
1680                 float sum = A[i][j];
1681                 for(int k=0;k<j;k++){   sum -= A[i][k] * A[k][j];   }
1682                 A[i][j] = sum;
1683                 float dum;
1684                 if((dum = vv[i] * fabs(sum)) >= big){
1685                     big = dum;
1686                     imax = i;
1687                 }
1688             }
1689             if(j != imax){
1690                 for(int k=0;k<n;k++){
1691                     float dum = A[imax][k];
1692                     A[imax][k] = A[j][k];
1693                     A[j][k] = dum;
1694                 }
1695                 d = -d;
1696                 vv[imax] = vv[j];
1697             }
1698             index[j] = imax;
1699
1700             if(A[j][j] == 0.0){ A[j][j] = tiny; }
1701
1702             if(j != n-1){
1703                 float dum = 1.0/A[j][j];
1704                 for(int i=j+1;i<n;i++){ A[i][j] *= dum; }
1705             }
1706         }
1707     }
1708         catch(exception& e) {
1709                 m->errorOut(e, "LinearAlgebra", "ludcmp");
1710                 exit(1);
1711         }
1712 }
1713
1714 /*********************************************************************************************************************************/
1715
1716 void LinearAlgebra::lubksb(vector<vector<float> >& A, vector<int>& index, vector<float>& b){
1717     try {
1718         float total;
1719         int n = (int)A.size();
1720         int ii = 0;
1721
1722         for(int i=0;i<n;i++){
1723             if (m->control_pressed) { break; }
1724             int ip = index[i];
1725             total = b[ip];
1726             b[ip] = b[i];
1727
1728             if (ii != 0) {
1729                 for(int j=ii-1;j<i;j++){
1730                     total -= A[i][j] * b[j];
1731                 }
1732             }
1733             else if(total != 0){  ii = i+1;   }
1734             b[i] = total;
1735         }
1736         for(int i=n-1;i>=0;i--){
1737             total = b[i];
1738             for(int j=i+1;j<n;j++){ total -= A[i][j] * b[j];  }
1739             b[i] = total / A[i][i];
1740         }
1741     }
1742         catch(exception& e) {
1743                 m->errorOut(e, "LinearAlgebra", "lubksb");
1744                 exit(1);
1745         }
1746 }
1747
1748
1749 /*********************************************************************************************************************************/
1750
1751 vector<vector<double> > LinearAlgebra::getInverse(vector<vector<double> > matrix){
1752     try {
1753         int n = (int)matrix.size();
1754
1755         vector<vector<double> > inverse(n);
1756         for(int i=0;i<n;i++){   inverse[i].assign(n, 0.0000);   }
1757
1758         vector<double> column(n, 0.0000);
1759         vector<int> index(n, 0);
1760         double dummy;
1761
1762         ludcmp(matrix, index, dummy);
1763
1764         for(int j=0;j<n;j++){
1765             if (m->control_pressed) { break; }
1766
1767             column.assign(n, 0);
1768
1769             column[j] = 1.0000;
1770
1771             lubksb(matrix, index, column);
1772
1773             for(int i=0;i<n;i++){   inverse[i][j] = column[i];  }
1774         }
1775
1776         return inverse;
1777     }
1778         catch(exception& e) {
1779                 m->errorOut(e, "LinearAlgebra", "getInverse");
1780                 exit(1);
1781         }
1782 }
1783
1784 /*********************************************************************************************************************************/